text:
Similar radiation mechanism in gamma-ray bursts and blazars: evidence from two luminosity correlations

# Similar radiation mechanism in gamma-ray bursts and blazars: evidence from two luminosity correlations

## Abstract

Active galactic nuclei (AGNs) and gamma-ray bursts (GRBs) are powerful astrophysical events with relativistic jets. In this Letter the broadband spectral properties are compared between GRBs and the well-observed blazars. The distribution of GRBs are consistent with the well-known blazar sequence including the and correlations, where is defined as the broadband spectral slope in radio-to-X-ray bands, and is defined as the spectral peak frequency. Moreover, GRBs occupy the low radio luminosity end of these sequences. These two correlations suggest that GRBs could have a similar radiation process with blazars both in the prompt emission and afterglow phases, i.e., synchrotron radiation.

gamma-ray burst: general - BL Lacertae objects: general - methods: statistical
1

## 1 Introduction

Gamma-ray bursts (GRBs) and active galactic nuclei (AGNs) are both powered by relativistic jets from accreting black holes (Gehrels et al. 2009; Urry & Padovani 1995). The central engines of GRBs are argued to be stellar-mass black holes (Woosley 1993) and for AGNs the central engines are supermassive black holes. GRBs are most powerful explosions with isotropic-equivalent energy erg in the universe (Zhang 2011), and can be detected out to very high-redshift universe (Lamb & Reichart 2000; Wang et al. 2012). So GRBs can probe high-redshift universe, including dark energy (Dai, Liang & Xu 2004; Schaefer 2007; Wang & Dai 2011). Blazars include two subtype of AGNs, i.e., flat-spectrum radio quasars (FSRQs) and BL Lac objects (BL Lacs). A subclass of AGNs, e.g., super-Eddington accreting supermassive black holes, are also proposed to be a standard candle (Wang et al. 2013).

The radiation mechanism of balzars is well constrained. The spectral energy distribution (SED) of blazars is well understood, including the low-energy (infrared-soft X-ray) bump and the high-energy (MeV-GeV) bump. The synchrotron radiation can account for the low-energy peak, while the MeV-GeV peak is produced by inverse Compton radiation. But for GRBs, the radiation mechanism for the prompt emission is still highly debated. The spectrum of prompt emission can be modeled by the “Band function” (Band et al. 1993), whose origin is still unknown (but see Lucas Uhm & Zhang 2013). Some studies (Mészáros et al. 1994; Daigne & Mochkovitch 1998) proposed that synchrotron radiation is the leading mechanism. Other mechanisms are also proposed (Pe’er et al. 2006; Rees & Mészáros 2005; Beloborodov 2010). The radiation mechanism of afterglows is well understood (Sari et al. 1998). The observed afterglow radiation is well explained by synchrotron radiation (Sari et al. 1998; Panaitescu & Kumar 2001).

Some studies (Zhang 2007; Wang & Dai 2013; Wang et al. 2014) have proposed that the mechanisms in different scale outflow or jet systems may be the same. Some works have been done on comparison between GRBs and AGNs. Wang & Wei (2011) compared the spectral properties of blazars and optically bright GRB afterglows, and found that GRB afterglows have the same radiation mechanism as BL Lac objects. A similar correlation of the synchrotron luminosity and Doppler factor between GRBs and AGNs has been found (Wu et al. 2011). Nemmen et al. (2012) suggested that the relativistic jets in AGNs and GRBs have a similar energy dissipation efficiency. Ma et al. (2014) extended the analysis of Nemmen et al. (2012) by adding X-ray binaries and low-luminosity AGNS. Wang & Dai (2013) found that the GRB X-ray flares and solar X-ray flares have similar distributions, which indicate that the X-ray flares of GRBs are due to a magnetic reconnection process. These similar distributions also exist in X-ray flares from black hole systems with (Wang et al. 2014). Zhang et al. (2013) found that the prompt emission of GRBs may be produced by magnetic dominated jets.

In this Letter we compare the broadband spectral properties of GRBs and blazars, including the and correlations, where is the radio-to-X-ray spectral slope. For a GRB, is the peak frequency of spectrum of prompt emission, while is the low peak of spectrum for a blazar. The aim of this Letter is to explore a possible similarity in radiation mechanism between GRBs and blazars. this Letter is organized as follows. In section 2, we present the sample of blazars and GRBs. The fitting results are given in section 3. Section 4 gives conclusions and discussions.

## 2 Samples

Chandra & Frail (2012) compiled a sample of GRB radio afterglow observations from 1997 to 2011. This catalog consists of 304 GRBs with radio observations. We select 43 GRBs with redshift measurements. The sample is listed in Table 1. For each GRB, the name of GRB and redshift are presented in Columns 1 and 2, respectively. X-ray flux measured in 0.3-10 keV energy band is given in Column 3 observed by Swift at the time of Column 6. For bursts observed by BASTE, the energy range is between 2 keV and 10 keV. We use the typical GRB spectrum to convert the flux to 1 keV. The observed radio flux (Column 4) and frequency (Column 5) at the time of Column 6 are also provided. Column 7 is the jet opening angle. For GRBs without opening angle determination, we assume a typical value 5 degree. The derived collimation-corrected radio luminosity at 5 GHz is given in Column 8. Column 9 gives , which is the peak energy of the prompt spectrum. The parameters from Column 2 to Column 7 are taken from Chandra & Frail (2012). We use the value of from Wang, Qi & Dai (2011). In the calculation, we use the peak frequency in the GRB rest frame. The collimation-corrected radio luminosity is calculated as

 νLν(5GHz)=4πd2LνfνFbeam(1+z)−α−1, (1)

where is the luminosity distance, is radio flux at 5 GHz, is the beaming factor and is the spectral slope (Sari et al. 1998). We adopt in the slow cooling case. The observed radio flux is converted to flux at 5 GHz using spectral slope . In this work, we assume the cosmological parameters: km s Mpc, , and .

We use the spectral properties of balzars from Fossati et al. (1998). This sample consists of all the parameters that we require, including redshift, X-ray flux at 1 keV, radio flux at 5 GHz, and synchrotron peak frequency. Fossati et al. (1998) found a power spectral sequence for the blazars despite the difference in the continuum shapes among different sub-classes of blazars. This sequence indicates that the radio luminosity is anti-correlated with the synchrotron peak. A plausible interpretation is that relativistic jets radiate via synchrotron and inverse Compton processes if the physical parameters (i.e. magnetic field) vary with luminosity.

## 3 Results

### 3.1 νLν(5GHz)−νpeak correlation in blazars and GRBs

The radio luminosities at 5 GHz are anti-correlated with the synchrotron peaks of blazars, as found by Fossati et al. (1998). This correlation has not been studied in GRBs so far. We investigate this correlation in GRBs for the first time. Figure 1 shows the correlation of blazars and GRBs. The black and open dots represent blazars and GRBs, respectively. There is a tight correlation between and as expected from the blazar sequence (Fossati et al. 1998). The correlation coefficient is at a significance level from Spearman rank-order statistical test. From this figure, we can see that the GRBs occupy the low-luminosity region of this correlation. For both blazars and GRBs, the best fitting result is

 logνLν(5 GHz)=(−0.91±0.03)logνpeak+56.65±0.46. (2)

The correlation coefficient is improved to with probability . The correlation coefficient has an obvious enhancement after adding the GRB sample. GRB 060218 may deviate from this correlation. The possible reason is that this GRB usually called X-ray flash has a low peak energy (Soderberg et al. 2006; Pian et al. 2006).

### 3.2 νLν(5GHz)−αRX correlation in blazars and GRBs

The broad-band spectral slope is also correlated with luminosity at 5 GHz in blazars (Fossati et al. 1998). We also investigate this correlation in GRBs. The broad-band spectral slope is defined as

 αRX=−log(fνR)/log(fνx)log(νR/νx), (3)

where GHz, and keV. The X-ray flux at keV can be obtained as follows. From Column 3 of Table 1, X-ray flux measured in 0.3-10 keV energy range can be obtained. The X-ray flux usually evolves as (Sari et al. 1998), with is the power-law index of accelerated electrons distribution. The X-ray flux in 0.3-10 keV is . After obtaining the value of , the flux at keV can be derived.

Figure 2 shows a correlation between the spectral slope and radio luminosity at 5 GHz for blazars (black dots) and GRBs (open dots). The correlation coefficient is with probability using Spearman rank-order statistical test for blazars. After combining GRBs and blazars, the fitting result is

 αRX=(0.084±0.006)log(νLν(5GHz))+(−2.87±0.24). (4)

The Spearman’s rank correlation coefficient is with probability . GRBs also occupy the the low-luminosity end of this correlation.

## 4 Conclusions and Discussions

The physics of GRBs are poorly understood, i.e., the radiation mechanism of prompt emission, the value of Lorentz factor, jet composition, and central engine (Zhang 2011). Wang & Dai (2013) found similar frequency distributions between X-ray flares of GRBs and solar X-ray flares, which may indicate the magnetically dominated jets in GRBs. In this Letter we compile 43 GRBs with well X-ray and radio observations. Two new correlations between GRBs and blazars may provide a new clue as to the radiation mechanism of GRB prompt emission and afterglows. For example, our clear and correlations suggest that the radiation mechanism of GRBs in prompt and afterglow phases and blazars is similar, namely, synchrotron radiation. Moreover, GRBs are occupy the low-luminosity region of these correlations.

Fossati et al. (1998) found that is anti-correlated with the synchrotron peak luminosity for blazars. For GRBs, Liang et al. (2004) found that is positively correlated with the isotropic-equivalent luminosity (about total luminosity). But the luminosity in the correlation for blazars is the synchrotron peak luminosity, not the total luminosity. Because there are two peaks in a blazar spectral energy distribution and there is no correlation between and in GRBs (Chandra & Frail 2012), from our simple analysis above we cannot conclude that GRBs have a different correlation compared with that for blazars. In this paper, we find that GRBs occupy the low radio luminosity end of the blazar sequence, which is similar to that of Wang & Wei (2011).

We thank the referee for detailed and very constructive suggestions that have allowed us to improve our manuscript. We acknowledge helpful discussions with Y. C. Zou and X. F. Wu. This work is supported by the National Basic Research Program of China (973 Program, grant No. 2014CB845800) and the National Natural Science Foundation of China (grants 11373022, 11103007, and 11033002).

### Footnotes

1. affiliationtext: E-mail: fayinwang@nju.edu.cn (FYW); dzg@nju.edu.cn (ZGD)

### References

1. Amati, L., Frontera, F., Tavani, M, et al. 2002, A&A, 390, 81
2. Band, D., Matteson, J., Ford, L., et al. 1993, ApJ, 413, 281
3. Bannister, K. W., Murphy, T., Gaensler, B. M. & Reynolds, J. E., 2012, ApJ, 757, 38
4. Blandford, R. D., & Levinson, A. 1995, ApJ, 441, 79
5. Chandra, P. & Frail, D. A., 2012, ApJ, 746, 156
6. Dai, Z. G., Liang, E. W., & Xu, D. 2004, ApJ, 612, L101
7. Daigne, F., & Mochkovitch, R. 1998, MNRAS, 296, 275
8. Beloborodov, A. M., 2010, MNRAS, 407, 1033
9. Fardal, M. A., Katz, N., Weinberg, D. H. & Dave, R., 2007, MNRAS, 379, 985
10. Ferguson, H. C., Dickinson, M., & Papovich, C., 2002, ApJ, 569, L65
11. Fossati, G., Maraschi, L., Celotti, A., Comastri, A., & Ghisellini, G. 1998, MNRAS, 299, 433
12. Friedmann, A. S., & Bloom, J. S. 2005, ApJ, 627, 1
13. Gehrels, N., Ramirez-Ruiz, E. & Fox, D. B., 2009, ARA&A, 47, 567
14. Ghirlanda, G., Ghisellini, G. & Lazzati, D., 2004, ApJ, 616, 331
15. Ghirlanda, G., Nava, L. & Ghisellini, G., 2010, A&A, 511, A43
16. Harris, D. E., & Krawczynski, H., 2006, ARA&A, 44, 463
17. Jorstad, S. G., Marscher, A. P., Mattox, J. R., et al. 2001, ApJ, 556, 738
18. Kharb, P., Lister, M. L. & Cooper, N. J., 2010, ApJ, 710, 764
19. Lamb, D. Q., & Reichart, D. E. 2000, ApJ, 536, 1
20. Liang, E. W., Dai, Z. G. & Wu, X. F., 2004, ApJ, 606, L29
21. Lucas Uhm, Z., & Zhang, B., 2013, arXiv: 1303.2704
22. Lyubarsky, Y. 2008, ApJ, 682, 1443
23. Ma, R. Y., Xie, F. G., & Hou, S. J., 2014, ApJL, 780, L14
24. Macquart, J. P., 2007, ApJL, 658, L1
25. Marscher, A. P., Jorstad, S. G., Larionov, V. M., et al. 2010, ApJ, 710, L126
26. Mészáros, P., Rees, M. J., & Papathanassiou, H. 1994, ApJ, 432, 181
27. Nemmen, R. S., Georganopoulos, M., Guiriec, S., et al. 2012, Sci, 338, 1445
28. Panaitescu, A., & Kumar, P., 2001, ApJ, 560, L49
29. Pe’er, A., Mészáros, P., & Rees, M. J., 2006, ApJ, 642, 995
30. Pian, E., Mazzali, P. A., Masetti, N., et al., 2006, Nature, 442, 1011
31. Rees, M. J., & Mészáros, P., 2005, ApJ, 628, 847
32. Sagiv, A., & Waxman, E. 2002, ApJ, 574, 861
33. Sari, R., Piran, T., & Narayan, R. 1998, ApJ, 497, L17
34. Schaefer, B. E., 2007, ApJ, 660, 16
35. Soderberg, A. M., Kulkarni, S. R., Nakar, E., et al., 2006, Nature, 442, 1014
36. Urry, C. M., & Padovani, P. 1995, PASP, 107, 803
37. Usov, V. V., & Katz, J. I. 2000, A&A, 364, 655
38. Wang, F. Y., Bromm, V., Greif. T. H., Stacy, A., Dai, Z. G., Loeb, A. & Cheng, K. S. 2012, ApJ, 760, 27
39. Wang, F. Y., & Dai, Z. G., 2013, Nature Phys., 9, 465
40. Wang, F. Y., Dai, Z. G., & Yi, S. X., 2014, submitted
41. Wang, F. Y., & Dai, Z. G. 2011, A&A, 536, A96
42. Wang, F. Y., Qi, S., & Dai, Z. G., 2011, MNRAS, 415, 3423
43. Wang, J., & Wei, J. Y. 2011, ApJL, 726, L4
44. Wang, J. M., Du, P., Valls-Gabaud, D., et al. 2013, PRL, 110, 081301
45. Woosley, S. E., ApJ, 1993, 405, 273
46. Wu, Q., Zou, Y. C., Cao, X., Wang, D. X., & Chen, L. 2011, ApJL, 740, L21
47. Yonetoku, D., Murakami, T., Nakamura, T. et al., 2004, ApJ, 609, 935
48. Zhang, B. 2011, CRPhy, 12, 206
49. Zhang, J., Liang, E. W., Sun, X. N., et al., 2013, ApJL, 774, L5
50. Zhang, S. N. 2007, HiA, 14, 41
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters


text:
# Rubik's cube array rotation

I'm about to make a 3D Rubik's cube based game in Flash AS3 and Away3d.

I don't really know how to manage the 2D arrays of the Rubik's cube.

For example, how do I rotate the corresponding arrays if I rotate a side, or just rotate a middle part?

In this stage I also don't know how to rotate those smaller cube parts all together if a side is rotating. First I was thinking of "groups" ( like in sketchup or 3ds max, blender), but that would be tricky, because the group components would change every time.

So I was thinking of just rotating each individual piece along a global axis. However, I just know the Away3d functions to rotate the cube of his local X , Y or Z axis, but how to rotate in global axis?

Does anyone know of a algorithm for doing these types of rotations?

• Hey Ace. I edited your question to remove the additional question not related to the rotation. I suggest you ask that in a different question. – MichaelHouse Jun 27 '13 at 13:30

There may be very clever ways to do this, but there's also brute force simple ways. Create a function for each axis X,Y and Z. Each of these RotateAround functions will take an input for which layer to rotate (1, 2, 3).
So for example, the RotateAroundY(1) call will do exactly as the image above shows. It will take the first row (from the top) of the Xplus, Zminus, Xminus, Zplus faces, and shift all their values through each other.
Another example, RotateAroundX(2) will take the second row of the (from the left) of the Yplus, Zminus, Yminus, ZPlus faces and shift all their values through each other.

format: plain


text:
Saturday, December 14, 2013

There is no regular approximately invariant finitely additive probability measure on all subsets of a cube or ball

For a totally ordered field K, say a hyperreal one, write xy (and say that they are approximately equal) provided that xy is 0 or infinitesimal. A K-valued probability P defined for all subsets of Ω is said to be regular provided that P(A)>0 whenever A is non-empty. It is approximately rigid motion invariant provided that P(A)≈P(gA) for every rigid motion g and set A such that AgA⊆Ω. The following can be proved in Zermelo-Fraenkel (ZF) set theory without any Axiom of Choice:

Theorem 1. There is no totally ordered field K and a regular K-valued approximately rigid motion invariant finitely additive probability on all subsets of a ball or cube Ω.

If we delete "approximately", this follows from this.

The result follows from this post. Given such a regular measure we can define a preorder ≤ by letting AB if and only if P(A)≤P(B). By the Theorem from that post, it follows in ZF that Banach-Tarski is true. But Banach-Tarski implies that there is no approximately rigid motion invariant finitely additive probability on all subsets of a ball or cube.

(Why ball or cube? This saves me from having to worry about some edge effects given our definition of invariance.)

Another result, proved by similar methods:

Theorem 2. Let Ω be a subset of three-dimensional Euclidean space invariant under rotations about the origin 0. If K is a totally ordered field and P is a regular K-valued finitely additive probability on all subsets of Ω approximately invariant under rotations about the origin, then P({0})≈1.

Suppose now that we have a particle undergoing Brownian motion released at time t0 at the origin, and then observed at time t1. The probability of its being in some set at t1 should be at least approximately invariant under rotations, and of course it is unacceptable to say that the probability that it is at the origin is approximately one—on the contrary, with approximately unit probability it is going to be away from the origin.

Update: Similar things hold for full-conditional probabilities, where approximate invariance is replaced with
invariance conditionally on the whole space (but there is no requirement of invariance conditionally on subsets of the space).


text:
# Risk Management for Enterprises and Individuals, v. 1.0

by Etti Baranoff, Patrick L. Brockett, and Yehuda Kahane

## 3.2 Uncertainty, Expected Value, and Fair Games

### Learning Objectives

• In this section we discuss the notion of uncertainty. Mathematical preliminaries discussed in this section form the basis for analysis of individual decision making in uncertain situations.
• The student should pick up the tools of this section, as we will apply them later.

As we learned in the chapters Chapter 1 "The Nature of Risk: Losses and Opportunities" and Chapter 2 "Risk Measurement and Metrics", risk and uncertainty depend upon one another. The origins of the distinction go back to the Mr. Knight,See Jochen Runde, “Clarifying Frank Knight’s Discussion of the Meaning of Risk and Uncertainty,” Cambridge Journal of Economics 22, no. 5 (1998): 539–46. who distinguished between risk and uncertainty, arguing that measurable uncertainty is risk. In this section, since we focus only on measurable uncertainty, we will not distinguish between risk and uncertainty and use the two terms interchangeably.

As we described in Chapter 2 "Risk Measurement and Metrics", the study of uncertainty originated in games of chance. So when we play games of dice, we are dealing with outcomes that are inherently uncertain. The branch of science of uncertain outcomes is probability and statistics. Notice that the analysis of probability and statistics applies only if outcomes are uncertain. When a student registers for a class but does not attend any lectures nor does any assigned work or test, only one outcome is possible: a failing grade. On the other hand, if the student attends all classes and scores 100 percent on all tests and assignments, then too only one outcome is possible, an “A” grade. In these extreme situations, no uncertainty arises with the outcomes. But between these two extremes lies the world of uncertainty. Students often do research on the instructor and try to get a “feel” for the chance that they will make a particular grade if they register for an instructor’s course.

Even though we covered some of this discussion of probability and uncertainty in Chapter 2 "Risk Measurement and Metrics", we repeat it here for reinforcement. Figuring out the chance, in mathematical terms, is the same as calculating the probability of an event. To compute a probability empirically, we repeat an experiment with uncertain outcomes (called a random experiment) and count the number of times the event of interest happens, say n, in the N trials of the experiment. The empirical probability of the event then equals n/N. So, if one keeps a log of the number of times a computer crashes in a day and records it for 365 days, the probability of the computer crashing on a day will be the sum of all of computer crashes on a daily basis (including zeroes for days it does not crash at all) divided by 365.

For some problems, the probability can be calculated using mathematical deduction. In these cases, we can figure out the probability of getting a head on a coin toss, two aces when two cards are randomly chosen from a deck of 52 cards, and so on (see the example of the dice in Chapter 2 "Risk Measurement and Metrics"). We don’t have to conduct a random experiment to actually compute the mathematical probability, as is the case with empirical probability.

Finally, as strongly suggested before, subjective probability is based on a person’s beliefs and experiences, as opposed to empirical or mathematical probability. It may also depend upon a person’s state of mind. Since beliefs may not always be rational, studying behavior using subjective probabilities belongs to the realm of behavioral economics rather than traditional rationality-based economics.

So consider a lottery (a game of chance) wherein several outcomes are possible with defined probabilities. Typically, outcomes in a lottery consist of monetary prizes. Returning to our dice example of Chapter 2 "Risk Measurement and Metrics", let’s say that when a six-faced die is rolled, the payoffs associated with the outcomes are $1 if a 1 turns up,$2 for a 2, …, and $6 for a 6. Now if this game is played once, one and only one amount can be won—$1, $2, and so on. However, if the same game is played many times, what is the amount that one can expect to win? Mathematically, the answer to any such question is very straightforward and is given by the expected value of the game. In a game of chance, if $W 1 , W 2 ,…, W N$ are the N outcomes possible with probabilities $π 1 , π 2 ,…, π N$ , then the expected value of the game (G) is $E(U)= ∑ i=1 ∞ π i U( W i )= 1 2 ×ln( 2 )+ 1 4 ×ln( 4 )+…= ∑ i=1 ∞ 1 2 i ln( 2 i ).$ The computation can be extended to expected values of any uncertain situation, say losses, provided we know the outcome numbers and their associated probabilities. The probabilities sum to 1, that is, $∑ i=1 N π i = π 1 +…+ π N =1.$ While the computation of expected value is important, equally important is notion behind expected values. Note that we said that when it comes to the outcome of a single game, only one amount can be won, either$1, $2, …,$6. But if the game is played over and over again, then one can expect to win $E( G )= 1 6 1+ 1 6 2…+ 1 6 6=3.50$ per game. Often—like in this case—the expected value is not one of the possible outcomes of the distribution. In other words, the probability of getting $3.50 in the above lottery is zero. Therefore, the concept of expected value is a long-run concept, and the hidden assumption is that the lottery is played many times. Secondly, the expected valueThe sum of the products of two numbers, the outcomes and their associated probabilities. is a sum of the products of two numbers, the outcomes and their associated probabilities. If the probability of a large outcome is very high then the expected value will also be high, and vice versa. Expected value of the game is employed when one designs a fair gameGame in which the cost of playing equals the expected winnings of the game, so that net value of the game equals zero.. A fair game, actuarially speaking, is one in which the cost of playing the game equals the expected winnings of the game, so that net value of the game equals zero. We would expect that people are willing to play all fair value games. But in practice, this is not the case. I will not pay$500 for a lucky outcome based on a coin toss, even if the expected gains equal $500. No game illustrates this point better than the St. Petersburg paradox. The paradox lies in a proposed game wherein a coin is tossed until “head” comes up. That is when the game ends. The payoff from the game is the following: if head appears on the first toss, then$2 is paid to the player, if it appears on the second toss then $4 is paid, if it appears on the third toss, then$8, and so on, so that if head appears on the nth toss then the payout is $2n. The question is how much would an individual pay to play this game? Let us try and apply the fair value principle to this game, so that the cost an individual is willing to bear should equal the fair value of the game. The expected value of the game E(G) is calculated below. The game can go on indefinitely, since the head may never come up in the first million or billion trials. However, let us look at the expected payoff from the game. If head appears on the first try, the probability of that happening is $1 2 ,$ and the payout is$2. If it happens on the second try, it means the first toss yielded a tail (T) and the second a head (H). The probability of TH combination $= 1 2 × 1 2 = 1 4 ,$ and the payoff is $4. Then if H turns up on the third attempt, it implies the sequence of outcomes is TTH, and the probability of that occurring is $1 2 × 1 2 × 1 2 = 1 8$ with a payoff of$8. We can continue with this inductive analysis ad infinitum. Since expected is the sum of all products of outcomes and their corresponding probabilities, $E( G )= 1 2 ×2+ 1 4 ×4+ 1 8 ×8+…=∞.$

It is evident that while the expected value of the game is infinite, not even the Bill Gateses and Warren Buffets of the world will give even a thousand dollars to play this game, let alone billions.

Daniel Bernoulli was the first one to provide a solution to this paradox in the eighteenth century. His solution was that individuals do not look at the expected wealth when they bid a lottery price, but the expected utility of the lottery is the key. Thus, while the expected wealth from the lottery may be infinite, the expected utility it provides may be finite. Bernoulli termed this as the “moral value” of the game. Mathematically, Bernoulli’s idea can be expressed with a utility function, which provides a representation of the satisfaction level the lottery provides.

Bernoulli used $U( W )=ln(W)$ to represent the utility that this lottery provides to an individual where W is the payoff associated with each event H, TH, TTH, and so on, then the expected utility from the game is given by

$E(U)= ∑ i=1 ∞ π i U( W i )= 1 2 ×ln( 2 )+ 1 4 ×ln( 4 )+…= ∑ i=1 ∞ 1 2 i ln( 2 i , )$

which can be shown to equal 1.39 after some algebraic manipulation. Since the expected utility that this lottery provides is finite (even if the expected wealth is infinite), individuals will be willing to pay only a finite cost for playing this lottery.

The next logical question to ask is, What if the utility was not given as natural log of wealth by Bernoulli but something else? What is that about the natural log function that leads to a finite expected utility? This brings us to the issue of expected utility and its central place in decision making under uncertainty in economics.

### Key Takeaways

• Students should be able to explain probability as a measure of uncertainty in their own words.
• Moreover, the student should also be able to explain that any expected value is the sum of product of probabilities and outcomes and be able to compute expected values.

### Discussion Questions

1. Define probability. In how many ways can one come up with a probability estimate of an event? Describe.
2. Explain the need for utility functions using St. Petersburg paradox as an example.
3. Suppose a six-faced fair die with numbers 1–6 is rolled. What is the number you expect to obtain?
4. What is an actuarially fair game?


text:
# The statistical physics of a simple continuous system

Let $$\mathbb{R}_{\geq 0}$$ be the set of real numbers greater or equal to zero. Assume an average value $$\overline{R} \in \mathbb{R}$$, called the prior. Then, the probability distribution $$p(r), \forall r \in \mathbb{R}_{\geq 0}$$ which maximizes the entropy based on the prior is the Gibbs distribution. Since $$\mathbb{R}$$ is an uncountable set, the partition function is an integral:

$$Z=\int_0^\infty e^{-\beta r}dr=\frac{1}{\beta}$$

The average value $$\overline{R}$$ is

$$\overline{R}=-\frac{1}{Z}\frac{\partial Z}{\partial \beta}=-\beta (- \beta^{-2})=\frac{1}{\beta}$$

and the entropy is

$$S=k_B(\ln [Z]+\beta \overline{R})=k_B\left(\ln[\beta^{-1}]+\beta \frac{1}{\beta}\right)=k_B \left( 1-\ln [\beta] \right)$$

Graphing the entropy yields:

Why is the entropy negative when $$1<\ln[\beta]$$? Should the entropy not be greater than zero for all values of $$\beta$$? Error somewhere?

EDIT:

As requested in the comments, the Log Plot of S is

and $$\beta$$ is a Lagrange multiplier.

EDIT-2:

As clarified in the comments, here is the plot of $$S$$ with respect to $$\ln{\beta}$$.

• Two things: 1) please define $\beta$, and 2) plot the entropy on a log scale. It's weird that you've defined the average to be $\bar{R}$ but then introduce $\beta$ as a separate parameter. – DanielSank Mar 2 '19 at 19:51
• @DanielSank I'm assuming $\beta=1/kT$ , at least this is what I have always seen. – BioPhysicist Mar 2 '19 at 20:00
• Alexandre, are you looking at a system whose energy varies linearly with respect to some parameter $r$, and where there is no degeneracy for any energies? – BioPhysicist Mar 2 '19 at 20:02
• Please make your question one cohesive post. An edit history is available for those who are interested. – BioPhysicist Mar 2 '19 at 20:08
• @Aaron correct, it physically corresponds to a system where the energy varies linearly and where there is no degeneracy. One way to think of it is an entropic force $F$ over a distance $x$. In this case $E=Fx$. – Alexandre H. Tremblay Mar 2 '19 at 20:13


text:
2015
04-13

# Cubic Maze

Consider an N*N*N cubic maze. We create the coordinate system as the picture below.
Now a robot is in the maze. It can be represented as a vector from his foot to head, and we call this vector the roboctor. To discribe the direction the robot will move to, we use a direct vector, which is prependicular to the roboctor. Both the roboctor and the direct vector must parallel to the coordinate axes. The robot can only be in the integral point. The robot can’t go beyond the bound of the maze(i.e. The (1,1,1)->(N,N,N) space).

The robot can execute the following 6 instructions:

• F instruction. Move 1 unit along the direct vector.
• L instruction. Turn the direct vector left 90 degree around the roboctor with the roboctor unchanged.
• R instruction. Turn the direct vector right 90 degree around the roboctor with the roboctor unchanged.
• U instruction. Change the direct vector to the original roboctor and change the roboctor to the opposite of the original direct vector.
• D instruction. Change the direct vector to the opposite of the original roboctor and change the roboctor to the original direct vector.
• B instruction. Change the direct vector to the opposite of the original direct vector.

The instructions are NUMBERED FROM 1 TO 6 AS ABOVE.

There are some bombs in some integral points in the maze. The robot must avoid to move to the bomb point. For each integral point, we use 1 to represent the bomb and 0 for not. We have the map of the maze, but it is encoded. The map is encoded as follows:

Encode the N integral points with the same X coordinate and the same Z coordinate along the negative Y direction as a 0/1 string, and convert it to a 8-based number. The 8-based number is the code we get.

The entire maze can be represented by N*N 8-based numbers.

At first, the robot is at an integral point A(xa,ya,za), with direct vector DA and roboctor RA. The destination condition of the robot is to be on an integral point B(xb,yb,zb) with the direct vector DB and the roboctor RB. Your mission is to give a series of instructions to guide the robot to get to the destination condition safely(i.e. Not to hit any single bomb) with the least number of instructions.

DA, RA, DB, RB will be represented as 1~6. 1~3 represents the positive direction of X, Y, Z axis respectively. 4~6 represents the negative direction of X, Y, Z axis respectively.

An integer T comes first indicating the number of test cases. Each case contains N+3 lines:
The first line contains an integer number N (3<=N<=30) representing the cubic space.
The second line contains 5 space separated integer numbers xa,ya,za,DA,RA (1<=xa, ya, za<=N,1<=DA, RA<=6) representing the start condition.
The third line contains 5 space separated integer numbers xb,yb,zb,DB,RB (1<=xb, yb, zb<=N, 1<=DB, RB<=6) representing the destination condition.
The following N lines each contains N 8-based space separated numbers, representing the maze. The jth number of the ith line represents a line whose X coordinate is j and Z coordinate is i in the maze.

There may exist a bomb even in the start or the destination.

An integer T comes first indicating the number of test cases. Each case contains N+3 lines:
The first line contains an integer number N (3<=N<=30) representing the cubic space.
The second line contains 5 space separated integer numbers xa,ya,za,DA,RA (1<=xa, ya, za<=N,1<=DA, RA<=6) representing the start condition.
The third line contains 5 space separated integer numbers xb,yb,zb,DB,RB (1<=xb, yb, zb<=N, 1<=DB, RB<=6) representing the destination condition.
The following N lines each contains N 8-based space separated numbers, representing the maze. The jth number of the ith line represents a line whose X coordinate is j and Z coordinate is i in the maze.

There may exist a bomb even in the start or the destination.

2
5
1 1 1 6 4
5 5 5 6 4
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
3
1 1 1 6 5
3 3 3 4 6
6 0 0
0 0 4
0 6 2

17
RFFFFRFFFFUBFFFFU
Sorry, I can't get there.

1. bottes vernies blanches

I appreciate the efforts you men and women place in to share blogs on such sort of matters, it was certainly useful. Keep Posting!

2. #include <cstdio>

int main() {


text:
10 minutes
9.1

# Estimating Costs

Keywords: Budget, Cost, Estimating, Project, Schedule

### Learning Objectives

1. Describe methods of estimating costs.
2. Identify the effects of project phase and complexity on the choice of estimating method.
3. Describe the method of combining cost estimates with a schedule to create a budget.

## ESTIMATING COSTS TO COMPARE AND SELECT PROJECTS

During the conceptual phase when project selection occurs, economic factors are an important consideration when choosing between competing projects. To compare the simple paybacks or internal rates of return between projects, an estimate of the cost of each project is made. The estimates must be accurate enough so that the comparisons are meaningful, but the amount of time and resources used to make the estimates should be appropriate to the size and complexity of the project. The methods used to estimate the cost of the project during the selection phase are generally faster and consume fewer resources than those used to create detailed estimates in later phases. They rely more on the expert judgment of experienced managers who can make accurate estimates with less detailed information. Estimates in the earliest stages of project selection are usually based from previous projects that can be adjusted—scaled—to match the size and complexity of the current project or by applying standardized formulas.

### Analogous Estimate

An estimate that is based on other project estimates is an analogous estimate. If a similar project cost a certain amount, then it is reasonable to assume that the current project will cost about the same. Few projects are exactly the same size and complexity, so the estimate must be adjusted upward or downward to account for the difference. The selection of projects that are similar and the amount of adjustment needed is up to the judgment of the person who makes the estimate. Normally, this judgment is based on many years of experience estimating projects, including incorrect estimates that were learning experiences for the expert.

#### Analogous Estimate for John’s Move

In the John’s move example, John asked a friend for advice about the cost of his move. His friend replied, “I moved from an apartment a little smaller than yours last year and the distance was about the same. I did it with a fourteen-foot truck. It cost about $575 for the truck rental, pads, hand truck, rope, boxes, and gas.” Because of the similarity of the projects, John’s initial estimate of the cost of the move was less than$700 and he decided that the cost would be affordable and the project could go forward.

Less experienced managers who are required to make analogous estimates can look through the documentation that is available from previous projects. As explained in Chapter 2, if those projects were evaluated using the Darnall-Preston Complexity Index (DPCI), the manager can quickly identify projects that have similar profiles to the project under consideration even if those projects were managed by other people. Comparing the original estimates with the final project costs on several previous projects with the same DPCI ratings gives a less-experienced manager the perspective that it would take many years to acquire by trial and error. It also provides references which the manager can use to justify the estimate.

### Parametric Estimate

If the project consists of activities that are common to many other projects, average costs are available per unit. For example, if you ask a construction company how much it would cost to build a standard office building, they will ask for the size of the building in square feet and the city in which the building will be built. From these two factors—size and location—the company’s estimator can predict the cost of the building. Factors like size and location are parameters—measurable factors that can be used in an equation to calculate a result. The estimator knows the average cost per square foot of a typical office building and adjustments for local labor costs. Other parameters such as quality of finishes are used to further refine the estimate. Estimates that are calculated by multiplying measured parameters by cost-per-unit values are parametric estimates.

#### Parametric Estimate for John’s Move

To estimate the size of the truck needed for John’s move, the parameter used by a truck rental company is the number of bedrooms, as shown in Figure 9.1.

Figure 9.1

Number of Bedrooms Used for Parametric Cost Estimate

The moving company assumes that the number of bedrooms is the important parameter in determining how big a truck is needed for a move. For John’s move, he has a one-bedroom apartment, so he chooses the fourteen-foot truck. Once the size is determined, other parameters, such as distance and days, are used to estimate the cost of the truck rental.

## ESTIMATING COSTS TO INITIATE PROJECTS

Once the project is selected, more accurate estimates are often needed to raise funds and agree on contracts with vendors in the initiation phase.

#### Estimate During the Initiation Phase of John’s Move

John recalled that his friend also told him how tiring it was to do all the packing, loading, and driving himself, and some items were damaged when the load shifted inside the truck during the trip. John decides to call in favors from two friends, Dion and Carlita, to help him pack in Chicago and to hire some of the skilled labor like that needed to load the truck properly.

### Vendor Bid Analysis

If services or products will be provided by vendors, the cost of those services can be determined by issuing a request for proposal (RFP). The RFP describes the work, service, or product to be provided by the vendor and the quality level required. The RFP is sent to a list of vendors who are qualified—meet standards of reliability and capability—to perform this type of work. They respond with a proposal for completing the work described in the RFP, including an estimate of the cost. Some government organizations are required to use the qualified vendor with the lowest bid. Other organizations are not bound to take the lowest bid but are usually required to justify their reasons for not doing so.

#### Using RFPs to Make Estimates on John’s Move

John wants to find out how much it would cost to hire a skilled crew to load and secure the furniture in the truck and then have another crew from the same company meet him in Atlanta to unload the truck and help him unpack. He is not sure if any companies offer this option, so he decides to ask three moving companies for bids. He also decides to ask for bids on a standard move that includes all phases of packing, loading, transportation, and unloading as a comparison to see if his cost-saving plan is worth the extra effort.

The project management team can review the responses by several vendors to the RFP to determine if their estimate of the cost of that aspect of the project is close to the estimate made during the project selection stage. If the estimates by the vendors are much higher than expected, and if the project cannot be completed for the cost that was used to select the project, the selection decision might have to be reconsidered. Reconsidering the selection of the project should take into consideration the economic ratings of the competing projects that were not chosen and who the project champions are for the projects that would be affected.

Some vendors may suggest an alternative way to meet the objective of the RFP in a more cost-effective manner that does not match the specifics of the RFP. Such alternatives can reduce costs if they are acceptable.

### Bottom-Up Estimating

The most accurate and time-consuming estimating method is to identify the cost of each item in each activity of the schedule, including labor and materials. If you view the project schedule as a hierarchy where the general descriptions of tasks are at the top and the lower levels become more detailed, finding the price of each item at the lowest level and then summing them to determine the cost of higher levels is called bottom-up estimating.

#### Bottom-Up Estimate for John’s Move

After evaluating the bids by the moving companies, John decides the savings are worth his time if he can get the packing done with the help of his friends. He decides to prepare a detailed estimate of costs for packing materials and use of a rental truck. He looks up the prices for packing materials and truck rental costs on company Web sites and prepares a detailed list of items, quantities, and costs, as shown in Figure 9.2.

Figure 9.2

Detailed Cost Estimate

This type of estimate is typically more accurate than an analogous or parametric estimate. In this example, the sum of packing materials and truck expenses is estimated to be \$661.25.

The detail can be rolled up—subtotaled—to display less detail. This process is made easier using computer software. On projects with low complexity, the cost estimates can be done on spreadsheet software. On larger projects, software that manages schedules can also manage costs and display costs by activity and by category. For example, the subtotal feature could be used in Excel and collapsed to show the subtotals for the two categories of costs, as shown in Figure 9.3.

Figure 9.3

Sum of Detailed Costs by Type

### Activity-Based Estimates

An activity can have costs from multiple vendors in addition to internal costs for labor and materials. Detailed estimates from all sources can be reorganized so those costs that are associated with a particular activity can be grouped by adding the activity code to the detailed estimate, as shown in Figure 9.4.

Figure 9.4

Detailed Costs Associated with Activities

The detailed cost estimates can be sorted by activity and then subtotaled by activity to determine the cost for each activity.

### Establishing a Budget

Once you have broken your project down into activities, you will be able to calculate your overall project costs by estimating and totaling the individual activity costs.

This process of subtotaling costs by category or activity is called cost aggregation.

### Budget Timeline

Costs are associated with activities and since each activity has a start date and a duration period, it is possible to calculate how much money will be spent by any particular date during the project. The money needed to pay for a project is usually transferred to the project account shortly before it is needed. These transfers must be timed so that the money is there to pay for each activity without causing a delay in the start of the activity. If the money is transferred too far in advance, the organization will lose the opportunity to use the money somewhere else, or they will have to pay unnecessary interest charges if the money is borrowed. A schedule of money transfers is created that should match the need to pay for the activities. The process of matching the schedule of transfers with the schedule of activity payments is called reconciliation. Refer to Figure 9.5 which shows the costs of ten major activities in a project. Funds are transferred into the project account four times. Notice that during most of the project, there were more funds available than were spent except at activity 6 when all the available funds were spent.

Figure 9.5

Fund Transfers and Expenditures

In the project budget profile shown in Figure 9.5 there is no margin for error if the total of the first six activities exceeds the amount of funding at that point in the project.

Contractual agreements with vendors often require partial payment of their costs during the project. Those contracts can be managed more conveniently if the unit of measure for partial completion is the same as that used for cost budgeting. For example, if a graphic designer is putting together several pieces of artwork for a textbook, their contract may call for partial payment after 25% of their total number of drawings is complete.

### Key Takeaways

• Analogous estimating scales an estimate from a similar project to match the current project. Parametric estimating multiplies a standard cost-per-unit value by the number of units in the project. Bids from contractors can be compared to estimate costs. Bottom-up estimating determines the cost of each detail and aggregates them to determine activity cost estimates.
• During the project selection and approval stage, rough estimates are used that are usually obtained using analogous and parametric methods. Vendor bid analysis and detailed bottom-up estimates are used in the initiation phase to estimate project costs.
• Detailed estimates are associated with activities and aggregated during the planning phase to create an activity-based budget. Funding transfers are arranged to reconcile money spent to money from funding sources in a timely manner.

CC BY-NC-SA: This work is released under a CC BY-NC-SA license, which means that you are free to do with it as you please as long as you (1) properly attribute it, (2) do not use it for commercial gain, and (3) share any subsequent works under the same or a similar license.

### End-of-Chapter Survey

: How would you rate the overall quality of this chapter?
1. Very Low Quality
2. Low Quality
3. Moderate Quality
4. High Quality
5. Very High Quality
Comments will be automatically submitted when you navigate away from the page.


text:
# Homework Help: What is the derivative of this

1. Jul 24, 2009

### JakePearson

hey guys
whats is the derivative of you are asked to integrate;

x(t) = pi / t2

i dont have a rule for this, so i cant integrate it :( help

2. Jul 24, 2009

### Pengwuino

pi is just a constant so it's simply the integral of $$\frac{1}{{t^2 }}$$ that you're looking at. What function, when differentiated, results in $$\frac{\pi }{{t^2 }}$$ ?

3. Jul 24, 2009

### n!kofeyn

Sure you have a rule: the power rule for integration. If c is a constant and n is not equal to -1, then
$$\int cx^n \,dx = \frac{cx^{n+1}}{n+1} + C.$$

Pi is just a constant. Can you write the function x(t) as x(t)=pi*tn for some n, and then apply the above rule?

Last edited: Jul 24, 2009
4. Jul 24, 2009

### Staff: Mentor

Someone looking at your question would probably be confused about what you are asking. At first appearances it looks like you want the derivative of the function you show, but the "asked to integrate" part would throw most people off.

For your own sake and for the sakes of the people who are trying to help you, try to ask your questions in the clearest way you can.

5. Jul 24, 2009

### JakePearson

what would be the derivative if you are asked to differentiate x(t) = pi / t^2

6. Jul 24, 2009

### n!kofeyn

Are you being asked to differentiate the function you gave, or integrate (i.e. anti-differentiate)?

To integrate, follow my advice above. To differentiate, apply the same advice to the power rule for derivatives. If c is a constant, n is any number (positive or negative), and f(x)=cxn, then f'(x)=cnxn-1. To differentiate your function x(t), can you write it as x(t)=pi*tn, for some number n?

7. Jul 24, 2009

### Staff: Mentor

I rest my case.

8. Jul 24, 2009

### slider142

What rules of differentiation do you know? Or do you just know the definition of the derivative of a function?

9. Jul 24, 2009

### HallsofIvy

Do you know that pi/t^2= pi t^{-2}? Do you know the derivative of t^n?

Share this great discussion with others via Reddit, Google+, Twitter, or Facebook


text:
# Allow \listoffigures captions to break across pages

I'm using the command \listoffigures to list all figure captions after the contents page. Many of my figure captions are quite long and as a result LaTeX is placing 2 captions on a page with a large gap in between them. I would like to prevent these gaps by allowing figure captions to be split across pages just like normal text would be - how can I do this? I am aware of the short caption option with \caption[short]{long} for long captions but this is not appropriate for my purposes.

Minimum example:

\documentclass{article}
\usepackage{blindtext}
\usepackage{capt-of}

\begin{document}

\listoffigures
\newpage
\centering

\section{Section 1}
\rule{0.5\linewidth}{0.5\linewidth}
\captionof{figure}{\blindtext}

\section{Section 2}
\rule{0.5\linewidth}{0.5\linewidth}
\captionof{figure}{\blindtext}

\section{Section 3}
\rule{0.5\linewidth}{0.5\linewidth}
\captionof{figure}{\blindtext}

\section{Section 4}
\rule{0.5\linewidth}{0.5\linewidth}
\captionof{figure}{\blindtext}

\end{document}

This creates an analogous problem to the one I have: in the list of figures, the third figure legend starts on a new page leaving a large space at the end of the second figure legend (bottom of page 1). How can I force the third figure legend to start on page 1 and finish on page 2 like normal text?

• Can you provide a MWE demonstrating the problem? May 24 '14 at 14:38
• @darthbith I have edited the question to provide an example of the problem. May 24 '14 at 16:22
• Page breaking of caption is difficult, I suppose ;-)
– user31729
May 24 '14 at 17:51
• If they are so complicated, why bother with a LOF at all? Or if you want to list the page with each figure, set the short caption to empty maybe so you don't have to worry about confusing readers? May 25 '14 at 14:33
• Some of them are complicated, not all of them. The list of figures is a requirement for the thesis I'm writing. May 26 '14 at 9:35

## 1 Answer

LaTeX sets \interlinepenalty to the maximum value, 10000 while setting each caption in the LOF, if you locally disable \interlinepenalty by defining it to be a scratch register, so that the value set is ignored then page breaking happens as usual.

\documentclass{article}
\usepackage{blindtext}
\usepackage{capt-of}

\begin{document}

{\countdef\interlinepenalty255
\listoffigures

}
\newpage
\centering

\section{Section 1}
\rule{0.5\linewidth}{0.5\linewidth}
\captionof{figure}{\blindtext}

\section{Section 2}
\rule{0.5\linewidth}{0.5\linewidth}
\captionof{figure}{\blindtext}

\section{Section 3}
\rule{0.5\linewidth}{0.5\linewidth}
\captionof{figure}{\blindtext}

\section{Section 4}
\rule{0.5\linewidth}{0.5\linewidth}
\captionof{figure}{\blindtext}

\end{document}


text:
Sample Size, $n$

Power, $1-\beta$

Type I error rate, $\alpha$

X-axis

min

max

Documentation
Test 1 Mean
Compare 2 Means
Compare k Means
Test 1 Proportion
Compare 2 Proportions
Compare Paired Proportions
Compare k Proportions
Test Time-To-Event Data
Test Odds Ratio
Test Relative Incidence in Self Controlled Case Series Studies
Other

## Calculate Sample Size Needed to Test 1 Proportion: 1-Sample, 1-Sided

This calculator is useful for tests concerning whether a proportion, $p$, is equal to a reference value, $p_0$. The Null and Alternative hypotheses are

$H_0:p=p_0$
$H_1:p\lt p_0$

or
$H_0:p=p_0$
$H_1:p\gt p_0$

### Formulas

This calculator uses the following formulas to compute sample size and power, respectively: $$n=p_0(1-p_0)\left(\frac{z_{1-\alpha}+z_{1-\beta}\sqrt{\frac{p(1-p)}{p_0(1-p_0)}}}{p-p_0}\right)^2$$
$$1-\beta=\Phi\left(\sqrt{\frac{p_0(1-p_0)}{p(1-p)}}\left(\frac{|p-p_0|\sqrt{n}}{\sqrt{p_0(1-p_0)}}-z_{1-\alpha})\right)\right)$$
where

### R Code

R code to implement these functions:

p=0.05
p0=0.02
alpha=0.05
beta=0.20
(n=p0*(1-p0)*((qnorm(1-alpha)+qnorm(1-beta)*sqrt(p*(1-p)/p0/(1-p0)))/(p-p0))^2)
ceiling(n) # 191
z=(p-p0)/sqrt(p0*(1-p0)/n)
(Power=pnorm(sqrt(p0*(1-p0)/p/(1-p))*(abs(z)-qnorm(1-alpha))))

### References

Chow S, Shao J, Wang H. 2008. Sample Size Calculations in Clinical Research. 2nd Ed. Chapman & Hall/CRC Biostatistics Series. page 85.


text:
## Proof by Induction

The simplest proofs by induction have two steps. This is maybe best illustrated with an example.

Suppose we wanted to prove that(1) for allW e sayis the proposition thatand we can write this as

The first step, writtenoris the statement thatoris true.

We can takeas the statementiewhich is true.

Now we assume thatare all true. If we can prove thatis true then we will have proved thatis true for allsince we can take any value for includinghence provingorhence provingetc.

To proveis true, we first write down

We work with the left hand side:

Now use

Ifthen

We have provedsois proved and we have proved (1).

This is the general procedure:

Showor P(1) then assumeand useto prove

Example: Prove

is true since both sides equal 3

Assume

Prove

because of

Ifthenhence

Example Prove(2)

Ifthen

Assumethen

We must prove

which we obtain from (2) by replacingbythroughout.

Hence (2) is proved.


text:
Journal topic
Biogeosciences, 15, 3967–3973, 2018
https://doi.org/10.5194/bg-15-3967-2018
Biogeosciences, 15, 3967–3973, 2018
https://doi.org/10.5194/bg-15-3967-2018

Technical note 03 Jul 2018

Technical note | 03 Jul 2018

# Technical Note: An efficient method for accelerating the spin-up process for process-based biogeochemistry models

Technical Note: An efficient method for accelerating the spin-up process for process-based biogeochemistry models
Yang Qu1, Shamil Maksyutov2, and Qianlai Zhuang1,3 Yang Qu et al.
• 1Department of Earth, Atmospheric, and Planetary Sciences, Department of Agronomy, Purdue University, West Lafayette, IN 47907, USA
• 2National Institute for Environmental Studies, 16-2 Onogawa, Tsukuba, Ibaraki, 305-8506, Japan
• 3Department of Agronomy, Purdue University, West Lafayette, IN 47907, USA

Correspondence: Qianlai Zhuang (qzhuang@purdue.edu)

Abstract

To better understand the role of terrestrial ecosystems in the global carbon cycle and their feedbacks to the global climate system, process-based biogeochemistry models need to be improved with respect to model parameterization and model structure. To achieve these improvements, the spin-up time for those differential equation-based models needs to be shortened. Here, an algorithm for a fast spin-up was developed by finding the exact solution of a linearized system representing the cyclo-stationary state of a model and implemented in a biogeochemistry model, the Terrestrial Ecosystem Model (TEM). With the new spin-up algorithm, we showed that the model reached a steady state in less than 10 years of computing time, while the original method requires more than 200 years on average of model run. For the test sites with five different plant functional types, the new method saves over 90 % of the original spin-up time in site-level simulations. In North American simulations, average spin-up time savings for all grid cells is 85 % for either the daily or monthly version of TEM. The developed spin-up method shall be used for future quantification of carbon dynamics at fine spatial and temporal scales.

1 Introduction

Biogeochemistry models contain state variables representing various pools of carbon and nitrogen and a set of flux variables representing the element and material transfers among different state variables. Model spin-up is a step to get biogeochemistry models to a steady state for those state and flux variables (McGuire et al., 1992; King, 1995; Johns et al., 1997; Dickinson et al., 1998). Spin-up normally uses cyclic forcing data to force the model run and reach a steady state, which will be used as initial conditions for model transient simulations. The steady state is reached when modeled state variables show a cyclic pattern or a constant value and often requires a significant amount of computation time, which needs to be accelerated for regional and global simulations at fine spatial and temporal scales.

Spin-up is normally achieved by running the model repeatedly using one or several decades of meteorological or climatic data until a steady state is reached. The step could require that the model repeatedly run for more than 2000 annual cycles in some extreme cases. Specifically, the model will check the stability of the simulated carbon and nitrogen fluxes as well as state variables with specified threshold values. For instance, the model will check if the simulated annual net ecosystem production (NEP) is less than 1 g C m−2 yr−1 (McGuire et al., 1992). Another method to reach a steady state is to obtain the analytical solutions (King et al., 1995; Comins, 1997), which might also take a significantly long time.

For different biogeochemistry models, spin-up could take hundreds and thousands of years to reach a stability, normally longer than the model projection period (Thornton and Rosenbloom, 2005). Therefore, a more efficient method to reach the steady state will speed up the entire model simulation. Recently, a semi-analytical method (Xia et al., 2012) has been adapted to a carbon–nitrogen coupled model to speed up the spin-up process. The idea is to obtain an analytical solution very close to a steady condition, then start spin-up from the solution, which could significantly reduce spin-up time. This technique did not reach a cyclic pattern for state and flux variables and required an additional spin-up process to achieve the steady state. However, Lardy et al. (2011) and Martin et al. (2007) have implemented their spin-up methods for a linear problem of soil carbon dynamics including their seasonal cycles.

Here we developed a method to accelerate the spin-up process in a nonlinear model. We tested the method for representative PFTs and North America with both daily and monthly versions of the Terrestrial Ecosystem Model (TEM; Zhuang et al., 2003). In addition, we compared the performance of our algorithms with the semi-analytical version of Xia et al. (2012). The new algorithms will help us conduct very high spatial and temporal resolution simulations with process-based biogeochemistry models in the future.

2 Method

## 2.1 TEM description

We used a process-based biogeochemistry model, TEM (Zhuang et al., 2003), as a test bed to demonstrate the performance of the new algorithms of spin-up. TEM simulates the dynamics of ecosystem carbon and nitrogen fluxes and pools (McGuire et al., 1992; Zhuang et al., 2010, 2003). It contains five state variables: carbon in living vegetation (Cv), nitrogen in living vegetation (Nv), organic carbon in detritus and soils (Cs), organic nitrogen in detritus and soils (Ns), and available inorganic soil nitrogen (Nav). Carbon and nitrogen dynamics in TEM are governed by the following equations:

$\begin{array}{}\text{(1)}& & \frac{\mathrm{d}{C}_{\mathrm{v}}}{\mathrm{d}t}=\text{GPP}-{R}_{\mathrm{A}}-{L}_{\mathrm{C}},\text{(2)}& & \frac{\mathrm{d}{N}_{\mathrm{v}}}{\mathrm{d}t}=\text{NUPTAKE}-{L}_{\mathrm{N}},\text{(3)}& & \frac{\mathrm{d}{C}_{\mathrm{s}}}{\mathrm{d}t}={L}_{\mathrm{c}}-{R}_{\mathrm{H}},\text{(4)}& & \frac{\mathrm{d}{N}_{\mathrm{s}}}{\mathrm{d}t}={L}_{\mathrm{N}}-\text{NETNMIN},\text{(5)}& & \frac{\mathrm{d}{N}_{\mathrm{av}}}{\mathrm{d}t}=\text{NINPUT}+\text{NETNMIN}-\text{NLOST}-\text{NUPTAKE},\end{array}$

where GPP is gross primary production, RA is autotrophic respiration, LCis carbon in litterfall, NUPTAKE is nitrogen uptake by vegetation, LN is nitrogen in litterfall, RH is heterotrophic respiration, NETNMIN is net rate of mineralization of soil nitrogen, NINPUT is nitrogen input from the outside ecosystem, and NLOST is nitrogen loss from the ecosystem. Key carbon fluxes are defined as

$\begin{array}{ll}& \text{GPP}={C}_{max}f\left(\text{PAR}\right)f\left(\text{PHENOLOGY}\right)f\left(\text{FOLIAGE}\right)\\ \text{(6)}& & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}f\left(T\right)f\left({C}_{\mathrm{a}},{G}_{\mathrm{v}}\right)f\left(\text{NA}\right)f\left(\text{FT}\right),\text{(7)}& & \text{NPP}=\text{GPP}-{R}_{\mathrm{A}},\text{(8)}& & \text{NEP}=\text{GPP}-{R}_{\mathrm{A}}-{R}_{\mathrm{H}}.\end{array}$

For detailed GPP definition, see Zhuang et al. (2003). NEP will be near zero when the ecosystem reaches a steady state. Therefore, the spin-up goal is to keep running the model driven with repeated climate forcing data until NEP is close to zero with a certain tolerance value (e.g., 0.1 g C m−2 yr−1).

## 2.2 Spin-up acceleration method

TEM can be reformulated as

$\begin{array}{}\text{(9)}& \frac{\mathrm{d}\mathbit{x}}{\mathrm{d}t}=g\left(\mathbit{x},t\right)+\mathbit{h}\left(t\right),\end{array}$

where x is a vector of state variables (e.g., CV); h(t) is the vector of carbon–nitrogen input from the atmosphere (such as nitrogen input), independent of x; and g(x,t) is the process rate function of element pools (e.g., GPP).

By linearizing the model in terms of pools, we could obtain

$\begin{array}{}\text{(10)}& g\left(\mathbit{x},t\right)=g\left({\mathbit{x}}_{\mathrm{0}},t\right)+\mathbf{J}\left(\mathbit{x}-{\mathbit{x}}_{\mathrm{0}}\right),\end{array}$

where x0 is initial pool sizes and J is the Jacobian matrix of the process rate:

$\begin{array}{}\text{(11)}& \mathbf{J}=\frac{\mathrm{d}g}{\mathrm{d}x}=\left(\begin{array}{ccc}\frac{\partial {g}_{\mathrm{1}}}{\partial {x}_{\mathrm{1}}}& \mathrm{\cdots }& \frac{\partial {g}_{\mathrm{1}}}{\partial {x}_{n}}\\ \mathrm{⋮}& \mathrm{⋮}\\ \frac{\partial {g}_{n}}{\partial {x}_{\mathrm{1}}}& \mathrm{\cdots }& \frac{\partial {g}_{n}}{\partial {x}_{n}}\end{array}\right),\end{array}$

where g represents g(x,t). xn represents each state variable in the TEM (e.g., VC). The numerical discretization of Eq. (9) is

$\begin{array}{ll}& {x}_{i,\phantom{\rule{0.125em}{0ex}}k}-{x}_{i,\phantom{\rule{0.125em}{0ex}}k-\mathrm{1}}=\mathit{\tau }\cdot {\mathbf{J}}_{k-\frac{\mathrm{1}}{\mathrm{2}}}\cdot {x}_{i,\phantom{\rule{0.125em}{0ex}}k-\mathrm{1}}\\ \text{(12)}& & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}+\mathit{\tau }\left(g\left({x}_{\mathrm{0}},\phantom{\rule{0.125em}{0ex}}k-\mathrm{1}\right)-{\mathbf{J}}_{k-\frac{\mathrm{1}}{\mathrm{2}}}\cdot {x}_{\mathrm{0},\phantom{\rule{0.125em}{0ex}}k-\mathrm{1}}+{h}_{k-\mathrm{1}}\right),\end{array}$

where τ is time step (month), xi, k is the pool xi size at time k, and Jk−½ is a Jacobian matrix at time step $k-\frac{\mathrm{1}}{\mathrm{2}}$. Here $\frac{\mathrm{1}}{\mathrm{2}}$ refers to the half time step in the middle of a month, at which values of J are calculated as the mean value at time steps k and k−1. xi,0 refers to the initial pool xi size.

We introduce

$\begin{array}{}\text{(13)}& {f}_{k-\mathrm{1}}=g\left({x}_{\mathrm{0}},k-\mathrm{1}\right)-{\mathbf{J}}_{k-\frac{\mathrm{1}}{\mathrm{2}}}\cdot {x}_{\mathrm{0},\phantom{\rule{0.125em}{0ex}}k-\mathrm{1}}+{h}_{k-\mathrm{1}}.\end{array}$

The Eq. (12) can then be written as

$\begin{array}{}\text{(14)}& {x}_{i,\phantom{\rule{0.125em}{0ex}}k}-{x}_{i,\phantom{\rule{0.125em}{0ex}}k-\mathrm{1}}=\mathit{\tau }\cdot {\mathbf{J}}_{k-\frac{\mathrm{1}}{\mathrm{2}}}\cdot {x}_{i,\phantom{\rule{0.125em}{0ex}}k-\mathrm{1}}+\mathit{\tau }\cdot {f}_{k-\mathrm{1}},\end{array}$

where Jk−½ is a Jacobian matrix at the time step $k-\frac{\mathrm{1}}{\mathrm{2}}$. After running a large number of annual cycles, the model approaches a cyclo-stationary state, which can be expressed by condition ${x}_{T+i}={x}_{i}$, where T is the number of time steps in one cycle, for example, when spin-up is made at a monthly time step using monthly climatology of temperature, precipitation and other forcing data, T=12, and $\stackrel{\mathrm{\to }}{{x}_{\mathrm{1}}}$ is the size of carbon pools on 1 January, while J1.5 is the matrix of mean process rate constants for January.

By introducing

${A}_{k}=\mathit{\tau }\cdot {\mathbf{J}}_{k-\frac{\mathrm{1}}{\mathrm{2}}},\phantom{\rule{1em}{0ex}}{y}_{k}=\mathit{\tau }{f}_{k-\mathrm{1}},\phantom{\rule{1em}{0ex}}{B}_{k}=\mathbf{I},\phantom{\rule{1em}{0ex}}{C}_{k}=\mathbf{I}+A,$

where I is an identity matrix, Eq. (12) can be further written as

$\begin{array}{}\text{(15)}& -{C}_{k}\cdot {x}_{i,\phantom{\rule{0.125em}{0ex}}k-\mathrm{1}}+{B}_{k}\cdot {x}_{i,k}={y}_{k}.\end{array}$

The cyclic boundary condition is ${x}_{\mathrm{1}}={x}_{T+\mathrm{1}}$.

Then Eq. (13) will become

$\begin{array}{}\text{(15a)}& -{C}_{\mathrm{1}}\cdot {x}_{i,\phantom{\rule{0.125em}{0ex}}T}+{B}_{\mathrm{1}}\cdot {x}_{i,\phantom{\rule{0.125em}{0ex}}\mathrm{1}}={y}_{\mathrm{1}}.\end{array}$

Thus Eq. (15) and (15a) become a formulation of a linear problem with T unknown vectors $\stackrel{\mathrm{\to }}{{x}_{T}}$, which can be solved using LU (lower and upper) decomposition or Gaussian elimination (Martin et al., 2007). Xia et al. (2012; see Eq. 4) and Kwon and Primeau (2006) also had linear equations for a steady state, but conducted the model simulation at an annual time step. Going for annual average form reduces the size of the problem and prevents Xia et al. (2012) from obtaining the exact solution of the problem including seasonal cycle (see their Eqs. 15, 15a). While our new approach runs the model at a monthly time step with the cyclic boundary conditions for state variables x, it still targets a steady state for the ecosystem at an annual time step instead of a monthly time step.

## 2.3 Numerical implementation

Equation (15a) is explicitly expressed as

$\begin{array}{ll}& \left(\begin{array}{cccccccc}B& \mathrm{0}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}& \mathrm{0}& \mathrm{0}& -C\\ -C& B& \mathrm{0}& \mathrm{\dots }& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -C& B& \mathrm{\dots }& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{\dots }& -C& B& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}& -C& B& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}& \mathrm{0}& -C& B\end{array}\right)\\ \text{(16)}& & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}×\left(\begin{array}{c}{x}_{\mathrm{1}}\\ {x}_{\mathrm{2}}\\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ {x}_{T}\end{array}\right)=\left(\begin{array}{c}{y}_{\mathrm{1}}\\ {y}_{\mathrm{2}}\\ \cdot \\ \cdot \\ \cdot \\ \cdot \\ {y}_{T}\end{array}\right).\end{array}$

Equation (16) can be shown in the form Mx=Y.

We apply the Gaussian elimination to the upper block that reduces M to a lower triangular form (Fig. 1). The resulting matrix is lower diagonal:

$\begin{array}{}\text{(17)}& {M}^{\prime }=\left(\begin{array}{cccccccc}B{}^{\prime }& \mathrm{0}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -C& B& \mathrm{0}& \mathrm{\dots }& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -C& B& \mathrm{\dots }& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{\dots }& -C& B& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}& -C& B& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}& \mathrm{0}& -C& B\end{array}\right).\end{array}$

Equation (16) is thus reduced to the form ${M}^{\prime }x={Y}^{\prime }$, where M is the lower diagonal, and solution of Eqs. (15a) and (16) will be readily obtained for x.

Figure 1Algorithms and procedures of the new spin-up method.

## 2.4 Algorithm implementation for TEM

In the original TEM, carbon fluxes can be defined as

$\begin{array}{}\text{(18)}& & \text{NPP}=\text{GPP}-\text{MR}-\text{GR},\text{(19)}& & \text{MR}={V}_{\mathrm{C}}\cdot {K}_{T},\text{(20)}& & \text{GR}=\left\{\begin{array}{ll}\mathrm{0.25}\cdot \left(\text{GPP}-\text{MR}\right),& \text{if}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\text{GPP}>\text{MR}\\ \mathrm{0},& \text{otherwise}\end{array}\right\,\end{array}$

where NPP is defined as the difference of GPP and plant maintenance respiration (MR) and growth respiration (GR). MR is assumed as a function of CV and temperature (KT). Here we revised MR calculation:

$\begin{array}{ll}& \text{MR}\\ \text{(21)}& & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}=\left\{\begin{array}{ll}{V}_{\mathrm{C}}\cdot {K}_{T},& \mathrm{\text{if}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\text{GPP}>{V}_{\mathrm{C}}\cdot {K}_{T}\\ \mathrm{0.75}\cdot {V}_{\mathrm{C}}\cdot {K}_{T}+\mathrm{0.25}\cdot \text{GPP},& \text{otherwise}.\end{array}\right\\end{array}$

Table 1Test sites for new spin-up algorithms.

Figure 2The time for NEP (g C yr−1 m−2) to reach a steady state with the original spin-up method at the Harvard Forest site. x represents model simulation years.

The NEP is defined as the difference between NPP and heterotrophic respiration (RH).

The basic workflow to implement the method is (1) linearizing TEM first to obtain a sparse matrix with n variable (for TEM, n = 5) system, (2) performing Gaussian elimination for the linear system and (3) solving the sparse matrix to acquire the state variable values (Fig. 1). To adapt this method to a daily version of TEM, we changed the cyclic condition T from 12 to 365. The other steps are the same as for the monthly version. We tested the new method for the carbon-only version and carbon–nitrogen coupled version of TEM for different PFTs (Table 1). Specifically, for the carbon-only version, we only solved the differential equations that govern the carbon dynamics, while for the carbon–nitrogen coupled version, we solved the differential equations that govern both carbon and nitrogen dynamics in the system. For both versions, the spin-up process strives to reach a steady state for carbon pools and fluxes.

Figure 3Simulated NEP (g C m−2 yr−1) with the original spin-up method after different spin-up times of (a) 50, (b) 100, (c) 150 and (d) 200 years. After these spin-up times, 63, 89, 93 and 98 % of grids reach their steady states, respectively.

3 Results and discussion

At the Harvard Forest site, the traditional spin-up method took 564 years to reach the steady state for both the carbon-only and coupled carbon–nitrogen simulations with an annual NEP of less than 0.1 g C m−2 yr−1 (Fig. 2). In contrast, the improved method took 72 years for the carbon only and 122 for the coupled carbon–nitrogen simulations. For carbon and nitrogen pools, it took another 45 years (equivalent cyclic time) to reach a steady state with a NEP of less than 0.1 g C m−2 yr−1. In comparison with the traditional spin-up method (Zhuang et al., 2003), the new method saved 65 % of computational time to reach the steady state in the carbon-only simulations (Table 2). The differences in steady-state carbon pools between using the new method and traditional spin-up methods were small (less than 0.85 %). Similarly, for the coupled carbon–nitrogen simulations, the new method saves a similar amount of time to reach the steady state.

For all seven test sites, the original spin-up method in TEM takes 204–564 years (1.1–2.5 s of computing time) to reach the steady state at different sites. In contrast, our new method only takes 0.3–0.6 s, while the semi-analytical method (Xia et al., 2012) will need 0.5–0.9 s to reach the steady state at different sites (Table 2). Compared to the original spin-up method, the new method is not only faster but also computationally stable.

The time of spin-up to reach a steady state of NEP varied for different PFT grids using the original method (Fig. 2). In general, to allow 98 % of grid cells to reach their steady states of NEP, it takes 250 annual model runs, while the new method will only need on average of 0.6 s (equivalent to 60-year annual model runs with the original method) (Fig. 3). For regional tests in North America, we found that the average saving time with the new method with monthly TEM is 25, 32 and 22 % for Alaska, Canada and the conterminous US, respectively.

Table 2Spin-up time comparison for different methods for seven study sites; seconds represent real computation time, and years refer to the annual spin-up cycles.

To compare the performance of the new method with other existing methods, we adapted the semi-analytical method (Xia et al., 2012) to the TEM model. To do that, we first revised the TEM model structure to

$\begin{array}{}\text{(22)}& \frac{\mathrm{d}P\left(t\right)}{\mathrm{d}t}=\mathit{\epsilon }\mathbf{AC}P\left(t\right),\end{array}$

where P(t) is a vector of pools in TEM (e.g., CV and CS). ε is a scalar. A is a pool transfer matrix (in which Aij represents the fraction of carbon transfer from pool j to i). C is a diagonal matrix with pool components (where diagonal components quantify the fraction of carbon left from the state variables after each time step). With this method, we obtained an analytical solution for the intermediate state. We then kept running TEM with the traditional spin-up process. Specifically, we started TEM simulation to estimate the state variable values. Based on these values, the spin-up runs were conducted to reach the final steady state. We found that the semi-analytical solution is better than the original spin-up method but slower than the new method proposed in this study (Table 2).

The TEM model has a relatively small set of state variables for carbon and nitrogen. The version we used is TEM 5.0, which has only five state variables (Zhuang et al., 2003). Thus, the linearization process is relatively easy and the matrix size is relatively small; consequently, the computing is not a burden. To accelerate the spin-up for multiple soil carbon pool models with relatively simple and linear decomposition processes, implementing our method will still be relatively easy but will take a great amount of computing time to equilibrate. For models such as CLM, multiple methods have been tested to accelerate their spin-up process (e.g., Fang et al., 2015), the direct analytical solution method introduced in this study might be time consuming to achieve.

4 Summary

We developed a new method to speed up the spin-up process in process-based biogeochemistry models. We found that the new method shortened 90 % of the spin-up time using the traditional method. For regional simulations in North America, average spin-up time saving is 85 % for either daily or monthly version of TEM. We consider our method is a general approach to accelerate the spin-up process for process-based biogeochemistry models. As long as the governing equations of the models can be formulated as the form in Eq. (9), the algorithm could be adopted accordingly. Our method will significantly help future carbon dynamics quantification with biogeochemistry models at fine spatial and temporal scales.

Data availability
Data availability.

All data used in this study are available from the authors upon request.

Author contributions
Author contributions.

QZ and SM designed and supervised the research. YQ performed model simulations and data analysis. All authors contributed to the paper writing.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

This study is supported through projects funding to Qianlai Zhuang from the Department of Energy (DESC0008092 and DE-SC0007007) and the NSF Division of Information and Intelligent Systems (NSF-1028291). The supercomputing resource is provided by the Rosen Center for Advanced Computing at Purdue University.

Edited by: Christopher A. Williams
Reviewed by: two anonymous referees

References

Baldocchi, D. D., Xu, L., and Kiang, N.: How plant functional-type, weather, seasonal drought, and soil physical properties alter water and energy fluxes of an oak–grass savanna and an annual grassland, Agr. Forest Meteorol., 123, 13–39, 2004.

Comins, H. N.: Analysis of nutrient-cycling dynamics, for predicting sustainability and CO2-response of nutrient-limited forest ecosystems, Ecol. Model., 99, 51–69, 1997.

Dickinson, R. E., Shaikh, M., Bryant, R., and Graumlich, L.: Interactive canopies for a climate model, J. Climate, 11, 2823–2836, 1998.

Fang, Y., Liu, C., and Leung, L. R.: Accelerating the spin-up of the coupled carbon and nitrogen cycle model in CLM4, Geosci. Model Dev., 8, 781–789, https://doi.org/10.5194/gmd-8-781-2015, 2015.

Gilmanov, T. G., Demment, M. W., Wylie, B. K., Laca, E. A., Akshalov, K., Baldocchi, D. D., and Emmerich, W. E.: Quantification of the CO2 exchange in grassland ecosystems of the world using tower measurements, modeling and remote sensing, in: XX International Grassland Congress: offered papers, 26, p. 587, 2005

Goulden, M. L., Winston, G. C., McMillan, A. M. S., Litvak, M. E., Read, E. L., Rocha, A. V., and Rob Elliot, J.: An eddy covariance mesonet to measure the effect of forest age on land–atmosphere exchange, Glob. Change Biol., 12, 2146–2162, 2006.

Gu, L., Meyers, T., Pallardy, S. G., Hanson, P. J., Yang, B., Heuer, M., and Wullschleger, S. D.: Influences of biomass heat and biochemical energy storages on the land surface fluxes and radiative temperature, J. Geophys. Res.-Atmos., 112, https://doi.org/10.1029/2006JD007425, 2007.

Gu, L., Massman, W. J., Leuning, R., Pallardy, S. G., Meyers, T., Hanson, P. J., and Yang, B.: The fundamental equation of eddy covariance and its application in flux measurements, 152, 135-148, Agr. Forest Meteorol., 2012.

Johns, T. C., Carnell, R. E., Crossley, J. F., Gregory, J. M., Mitchell, J. F., Senior, C. A., Tett, S. F. B., and Wood, R. A.: The second Hadley Centre coupled ocean-atmosphere GCM: model description, spinup and validation, Clim. Dynam., 13, 103–134, 1997.

King, D. A.: Equilibrium analysis of a decomposition and yield model applied to Pinus radiata plantations on sites of contrasting fertility, Ecol. Model., 83, 349–358, 1995.

Kwon, E. Y. and Primeau, F.: Optimization and sensitivity study of a biogeochemistry ocean model using an implicit solver and in situ phosphate data, Global Biogeochem. Cy., 20, GB4009, https://doi.org/10.1029/2005GB002631, 2006.

Lardy, R., Bellocchi, G., and Soussana, J. F.: A new method to determine soil organic carbon equilibrium, Environ. Modell. Softw., 26, 1759–1763, 2011.

Martin, M. P., Cordier, S., Balesdent, J., and Arrouays, D.: Periodic solutions for soil carbon dynamics equilibriums with time-varying forcing variables, Ecol. Model., 204, 523–530, 2007.

McGuire, A. D., Melillo, J. M., Joyce, L. A., Kicklighter, D. W., Grace, A. L., Moore, B. I. I. I., and Vorosmarty, C. J.: Interactions between carbon and nitrogen dynamics in estimating net primary productivity for potential vegetation in North America, Global Biogeochem. Cy., 6, 101–124, 1992.

NVIDIA CUDA Compute Unified Device Architecture Programming Guide, NVIDIA Corporation, 2007.

Ollinger, S. V. and Smith, M. L.: Net primary production and canopy nitrogen in a temperate forest landscape: an analysis using imaging spectroscopy, modeling and field data, Ecosystems, 8, 760–778, 2005.

Thornton, P. E. and Rosenbloom, N. A.: Ecosystem model spin-up: Estimating steady state conditions in a coupled terrestrial carbon and nitrogen cycle model, Ecol. Model., 189, 25–48, 2005.

Turnipseed, A. A., Anderson, D. E., Blanken, P. D., Baugh, W. M., and Monson, R. K.: Airflows and turbulent flux measurements in mountainous terrain: Part 1. Canopy and local effects, Agr. Forest Meteorol., 119, 1–21, 2003.

Turnipseed, A. A., Anderson, D. E., Burns, S., Blanken, P. D., and Monson, R. K.: Airflows and turbulent flux measurements in mountainous terrain: Part 2: Mesoscale effects, Agr. Forest Meteorol., 125, 187–205, 2004.

Van Gorsel, E., Delpierre, N., Leuning, R., Black, A., Munger, J. W., Wofsy, S., and Chen, B.: Estimating nocturnal ecosystem respiration from the vertical turbulent flux and change in storage of CO2, Agr. Forest Meteorol., 149(11), 1919–1930, 2009.

Xia, J. Y., Luo, Y. Q., Wang, Y.-P., Weng, E. S., and Hararuk, O.: A semi-analytical solution to accelerate spin-up of a coupled carbon and nitrogen land model to steady state, Geosci. Model Dev., 5, 1259–1271, https://doi.org/10.5194/gmd-5-1259-2012, 2012.

Zhuang, Q., McGuire, A. D., Melillo, J. M., Clein, J. S., Dargaville, R. J., Kicklighter, D. W., Myneni, R. B., Dong, J., Romanovsky, V. E., Harden, J., and Hobbie, J. E.: Carbon cycling in extratropical terrestrial ecosystems of the Northern Hemisphere during the 20th century: a modeling analysis of the influences of soil thermal dynamics, Tellus B, 55, 751–776, 2003.

Zhuang, Q., He, J., Lu, Y., Ji, L., Xiao, J., and Luo, T.: Carbon dynamics of terrestrial ecosystems on the Tibetan Plateau during the 20th century: an analysis with a process-based biogeochemical model, Global Ecol. Biogeogr., 19, 649–662, 2010.


text:
# Math Help - Factor each polynomial

1. ## Factor each polynomial

I need help with these three questions,
can someone explain to me CLEARLY how to factor each polynomial

1. 125x^6 - 8
2. 64x^12 - 27
3. 8/125x^3 - 64y^6

2. Originally Posted by skeske1234
I need help with these three questions,
can someone explain to me CLEARLY how to factor each polynomial

1. 125x^6 - 8
2. 64x^12 - 27
3. 8/125x^3 - 64y^6
Hi skeske1234,

Do you know the model for factoring the difference of two cubes?

Review this:

$x^3-y^3=(x-y)(x^2+xy+y^2)$

Apply that to your problems. First let's make them look a bit more like the difference of two cubes.

[1] $(5x^2)^3-(2)^3$

[2] $(4x^4)^3-(3)^3$

[3] $\left(\frac{2}{5}x\right)^3-(4y^2)^3$

3. Originally Posted by masters
Hi skeske1234,

Do you know the model for factoring the difference of two cubes?

Review this:

$x^3-y^3=(x-y)(x^2+xy+y^2)$

Apply that to your problems. First let's make them look a bit more like the difference of two cubes.

[1] $(5x^2)^3-(2)^3$

[2] $(4x^4)^3-(3)^3$

[3] $\left(\frac{2}{5}x\right)^3-(4y^2)^3$
after you do that, how do you plug it into the formula?

4. Originally Posted by masters
Hi skeske1234,

Do you know the model for factoring the difference of two cubes?

Review this:

$x^3-y^3=(x-y)(x^2+xy+y^2)$

Apply that to your problems. First let's make them look a bit more like the difference of two cubes.

[1] $(5x^2)^3-(2)^3$

[2] $(4x^4)^3-(3)^3$

[3] $\left(\frac{2}{5}x\right)^3-(4y^2)^3$
Originally Posted by skeske1234
after you do that, how do you plug it into the formula?
I'll do the first one, and you see if that helps you do the other two.

$x^3-y^3=(x-y)(x^2+xy+y^2)$

In my model, let's let $x=5x^2 \ \ and \ \ y=2$

$(5x^2-2)\bigg((5x^2)^2+(2)(5x^2)+(2)^2\bigg)$

$(5x^2-2)(25x^4+10x^2+4$


text:
Here's the question you clicked on:

ChrisV 3 years ago The speed of light is 3.00 \times 10^8{\rm m/s}. How long does it take for light to travel from Earth to the Moon and back again?

• This Question is Closed
1. ChrisV

The distance from Earth to the Moon is approximately 240,000 mi

Ask your own question

Sign Up
Find more explanations on OpenStudy
Privacy Policy


text:
# django-jsonfield /

Filename Size Date modified Message
jsonfield
69 B
29 B
1.7 KB
562 B

# django-jsonfield

I had a serious need for a JSON field for django. There were a couple out there, but none packaged up nicely on bitbucket/github that were usable with pip install -e.

So I took the code from David Cramer's blog, and packaged it up.

To use, just install the package, and then use the field:

from django.db import models
import jsonfield

class MyModel(models.Model):
the_json = jsonfield.JSONField()

Now, it will validate the JSON on entry, and store it as a string in the database. When you instantiate/fetch the object, it will be turned back into a python list/dict/string.

There are also a couple of other bits and bobs:

• jsonify templatetag: This allows you to convert a python data structure into JSON within a template.

{% load jsonify %}

<script> var foo = {{ bar|jsonify }}; </script>

• JSONTableWidget I strongly recommend you don't use this: it's an experimental way of having a generic JSON structure represented as a table, enabling you to let non-programmer types edit JSON data. I have replaced a lot of the code that is here with some new KnockoutJS-based stuff, and that simplifies it a whole lot. As such, I am not using this widget any longer, and it may not work.

## History

• 0.7 : First time I tagged releases.

## Todo

Convert date/time objects nicely to/from ISO strings (YYYY-mm-dd HH:MM:SS TZNAME). This is actually a bit tricky, as we don't know if we are expecting a date/time object. We may parse objects as we go, but there could be some performance issues with this.

Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.


text:
# A modern Python interface for the Generic Mapping Tools - Demo¶

Thank you for giving GMT/Python and GMT6 a try!

This is a Jupyter notebook. It's an interactive computing environment where you can mix text (like this), code, and figures. The notebook is organized into cells. This is a Markdown cell (click on it to see the source) and it can contain text, hyperlinks, images, and even Latex equations.

To execute the code cells, select it and type shift+enter or click on the Run button above.

The following are examples of what you can currently do with GMT/Python. We'll also make the image for the poster background. There are some empty cells for you to experiment on the bottom of the notebook along with an example using the command line GMT6 in the new "modern mode".

The GMT modules are available as functions and classes in the gmt Python package. So we'll start by importing it:

In [ ]:
import gmt

## A first map¶

All figure generation in GMT/Python is handled by the gmt.Figure class. It has methods to add layers to your figure, like a basemap, coastlines, and data.

We start a new figure by creating an instance of gmt.Figure:

In [ ]:
fig = gmt.Figure()

We add elements to the figure using its methods. For example, lets add the coastlines of Central America to a 6 inch wide map using the Mercator projection (M). Our figure will also have a nice frame with automatic ticks.

In [ ]:
fig.coast(region=[-90, -70, 0, 20], projection='M6i', land='chocolate',
frame=True)

You can see a preview of the figure directly in the Jupyter notebook using fig.show().

In [ ]:
fig.show()

## A note for experienced GMT users¶

You'll probably have noticed several things that are different from classic command-line GMT. Many of these changes reflect the new GMT modern execution mode that will be part of the future 6.0 release. A few are GMT/Python exclusive (like the long argument names).

1. The name of method is coast instead of pscoast. As a general rule, all ps* modules had their ps removed. The exceptions are: psxy == plot, psxyz == plot3d, and psscale == colorbar.
2. The arguments don't use the GMT 1-letter syntax (R, J, B, etc). Those are still available as aliases and the methods will accept them (see below).
3. Arguments like region can take lists instead of strings like 1/2/3/4. You can still use the string form but the list form is easier in Python.
4. If a GMT argument has no options (like -B instead of -Baf), use a True value instead. An empty string would also be acceptable.
5. There is no output redirecting to a PostScript file. The figure is generated in the background and will only be shown or saved when you ask for it.

We could have generated the figure above using the classic GMT argument names (but not the module names):

In [ ]:
fig_alias = gmt.Figure()
fig_alias.coast(R='-90/-70/0/20', J='M6i', G='chocolate', B=True)
fig_alias.show()

## Saving the figure¶

Unlike the GMT command-line interface, no figure file was generated until you ask for one.

Use method fig.savefig (based on the matplotlib function) to save your figure to a file:

In [ ]:
fig.savefig('first-steps-central-america.png')

If you're running a Python script, you can tell fig.savefig to open the figure in an external viewer:

fig.savefig('first-steps-central-america.png', show=True)

## Plot point data¶

We can use gmt.Figure.plot to plot data on our map.

First, lets create some fake data to plot using numpy:

In [ ]:
import numpy as np

In [ ]:
# See the random number generator so we always
# get the same numbers
np.random.seed(42)

ndata = 30
region = [150, 240, -30, 60]
lon = np.random.uniform(region[0], region[1], ndata)
lat = np.random.uniform(region[2], region[3], ndata)
magnitude = np.random.uniform(1, 9, ndata)
depth = np.random.uniform(0, 1, ndata)

There are 3 ways to pass data into Figure.plot:

1. x and y coordinates as 1d numpy arrays using the x and y arguments.
2. A whole data table as a 2d numpy array using the data argument.
3. A file name using the data argument.

Let's explore all of these options.

### Using x and y¶

Now we can plot the data using Figure.plot and passing the x and y coordinate arrays:

In [ ]:
fig = gmt.Figure()
fig.coast(region='g', projection='G200/30/6i', frame='ag',
resolution='i', area_thresh=5000, land='white',
water='DarkTurquoise')
# Plot using circles (c) of 0.3 cm
fig.plot(x=lon, y=lat, style='c0.3c', color='black')
fig.show()

We can make the size of the markers follow a data value by passing in the argument sizes:

In [ ]:
fig = gmt.Figure()
fig.coast(region='g', projection='G200/30/6i', frame='ag',
resolution='i', area_thresh=5000, land='white',
water='DarkTurquoise')
fig.plot(x=lon, y=lat, sizes=0.1*magnitude, style='cc', color='black')
fig.show()

We can also map the colors of the markers to data by passing color as an array and providing a colormap name (cmap). We can even use the new matplotlib colormap "viridis":

In [ ]:
fig = gmt.Figure()
fig.coast(region='g', projection='G200/30/6i', frame='ag',
resolution='i', area_thresh=5000, land='white',
water='DarkTurquoise')
fig.plot(x=lon, y=lat, style='cc', sizes=0.1*magnitude,
color=depth, cmap='viridis')
fig.show()

### Using a matrix¶

Sometimes, you'll have data that you loaded from a file in the form of a single numpy 2d array (like a table). You can plot these data by passing in the data argument.

In [ ]:
table = np.transpose([lon, lat, magnitude, depth])

We can use the columns argument to specify the columns that correspond to x, y, color, and size, respectively. We'll use it to specify that we want to set the colors using the depths:

In [ ]:
fig = gmt.Figure()
fig.coast(region=region, projection='M6i', frame='ag',
resolution='i', area_thresh=5000, land='grey')
fig.plot(data=table, style='c0.8c', cmap='viridis', columns='0,1,3')
fig.show()

### Using a file name¶

If you have data in a file and you don't need to do any calculations on it, you can pass the file name to plot directly. The syntax for plotting will be the same as the example for a data matrix.

First, we must save our sample data to a file:

In [ ]:
np.savetxt('first-steps-data.txt', table)

In [ ]:
fig = gmt.Figure()
fig.coast(region='g', projection='N-165/6i', frame='ag',
resolution='i', area_thresh=5000, land='chocolate')
fig.plot(data='first-steps-data.txt', style='c0.4c',
cmap='viridis', columns='0,1,3')
fig.show()

### Using sample data from GMT¶

GMT ships some sample data that can be accessed by passing @data_file as the data argument. For example, we can plot the earthquake epicenters from the tut_quakes.ngdc sample dataset:

In [ ]:
fig_quakes = gmt.Figure()
fig_quakes.coast(
region=[130, 150, 35, 50], projection='M6i',
frame='afg', shorelines=True, land='gray',
water='lightblue')
fig_quakes.plot(data='@tut_quakes.ngdc', style='c0.5c',
color='blue', pen='faint', columns=[4, 3])
fig_quakes.show()

## Making the background plot for the poster¶

I'll use the USGS quake data that comes with GMT. I downloaded it and saved to a numpy-friendly format in file usgs_quakes.txt.

In [ ]:
lon, lat, magnitude = np.loadtxt('usgs_quakes.txt', unpack=True)

Make a global Mercator map without any borders. Plot the quakes mapping the marker size and color to the magnitude.

In [ ]:
fig = gmt.Figure()
fig.coast(region=[-270, 90, -70, 70], projection='M10i', land='#aaaaaa',
water='white', resolution='l')
fig.plot(lon, lat, style='cc', sizes=0.005*2**magnitude,
color=magnitude/magnitude.max(), cmap='ocean')
fig.show(width=900)

If it looks OK, save it to a high resolution PNG.

In [ ]:
fig.savefig('poster_background.png', dpi=1200)

# Experiment for yourself¶

Type anything you want in the cells below.

In [ ]:

In [ ]:

In [ ]:

# GMT6 demo¶

You can even try the GMT6 command line programs. The following are some of the plots from Paul's talk "The Generic Mapping Tools 6: Classic versus Modern Mode".

Cells that start with %%bash run their code through the bash shell (use them to execute GMT commands). We'll use the IPython.display.Image class to display in the notebook the PNG figures that you generate on the %%bash cells.

In [ ]:
from IPython.display import Image

## Example 1¶

In [ ]:
%%bash
gmt begin Chile png
gmt pscoast -RCL+r2 -JM15c+ -BWSne -B -Gbeige -Sblue -N1,1p
gmt end

In [ ]:
Image('Chile.png', width=120)

## Example 2¶

In [ ]:
%%bash
gmt begin map png
gmt grdimage @earth_relief_05m -RMG+r2 -Cgeo -I+
gmt coast -Wthin -BWSne -B
gmt colorbar -DJTC -B -C+Uk
gmt end

In [ ]:
Image('map.png', width=300)

## Example 3¶

In [ ]:
%%bash
gmt begin islands png
gmt set MAP_ANNOT_OBLIQUE 34 FORMAT_GEO_MAP ddd:mmF
gmt subplot begin 2x2 -M0.05i -Fs3i/3i -BWSne -A+gwhite+p0.5p
gmt grdimage @earth_relief_03s -R-30/30/-30/30+uk -JA159:32W/22:03N -Ba20mf5m -Csrtm -c1,1
gmt grdimage @earth_relief_03s -R-15/15/-15/15+uk -JA109:20W/27:07S -Ba20mf5m -c1,2
gmt grdimage @earth_relief_03s -R-30/30/-30/30+uk -JA149:22W/17:43S -Ba20mf5m -c2,1
gmt grdimage @earth_relief_03s -R-10/10/-10/10+uk -JA138:39W/10:29S -Ba20mf5m -c2,2
gmt subplot end
gmt colorbar -B -Dx0/0.4i+jBC+w5i+h -Xw/2 -Yh
gmt end

In [ ]:
Image('islands.png', width=500)

## Example 4¶

In [ ]:
%%bash
gmt begin
gmt figure hotspots png
gmt grdimage @earth_relief_10m -JG200/30/6i -Cgeo -I+
gmt coast -W -Dc -Bafg
gmt plot @hotspots.txt -Sc0.2c -Gred
gmt end

In [ ]:
Image('hotspots.png', width=500)

# Experiment for yourself¶

Try running different things in the cells below or edit and rerun the cells above.

In [ ]:
%%bash

In [ ]:
Image('XXX.png', width=500)

In [ ]:
%%bash

In [ ]:
Image('XXX.png', width=500)

In [ ]:
%%bash

In [ ]:
Image('XXX.png', width=500)


text:
# How to find slope using a line that has an x-intercept of 3 and a y-intercept of 4?

##### 1 Answer
Apr 7, 2015

By definition an x-intercept of $3$ is the point $\left(x , y\right) = \left(3 , 0\right)$
and
a y-intercept of $4$ is the point $\left(x , y\right) = \left(0 , 4\right)$

The slope is defined as
$\frac{\Delta y}{\Delta x} = \frac{{y}_{2} - {y}_{1}}{{x}_{2} - {x}_{1}}$

In this case
slope $= \frac{4 - 0}{0 - 3} = - \frac{4}{3}$