text:
# Total pressure in a chamber calculation

I have a chamber containing two different gases and a liquid which needs to be heated at a certain temperature to bring it to vapor state. I have worked out that we can calculate the total pressure of a mixture of gas using Dalton's law using partial pressure of each gas. To calculate total vapor pressure we have Raoult's law. If I have both gases and vapors in the medium can we calculate the total pressure by adding them both, as in a "collecting gases over water experiment"?

• Okay, let me put it this way. If I know the partial pressure of 2 gases and vapor pressure of liquid turned vapor, can I directly add them to get total pressure in the chamber? – NidhiS Sep 9 '15 at 10:53
• This not my homework question. I am an electronics engineer working on an interdisciplinary project wherein I need some information on mixture of gases. I have worked out that we can calculate the total pressure of a mixture of gas using Dalton's law using partial pressure of each gas. To calculate total vapor pressure we have Raoult's law. If I have both gases and vapors in the medium I think we can calculate the total pressure by adding the both like 'collecting gases over water experiment'. I read these things long ago so I wanted to confirm it. – NidhiS Sep 9 '15 at 11:15
• Thank you for the clarification, could you please add this detail to the question – user15489 Sep 9 '15 at 11:17

Dalton's Law is only true for ideal compounds, which theoretically have no intermolecular forces between their molecules. Those forces become more and more significant when increasing the concentration of any compound, because all molecules become packed up closer.

This will lower the pressure, because if any molecule is about to hit the container wall, it gets pulled backwards $-$ and thus slowed $-$ due to the attraction of the other molecules behind it. But for high temperature this effect is negligible, because they can't be slowed down so easily.

Now lets come to the calculation using Kay's rule:

Since you only have one liquid component, you don't need Roult's Law to calculate the vapor pressure for a mixture of liquids.

1.) You need to find the values for the critical temperatures and critical pressures for your components (gases and water vapor).

2.) You need to know the relative amounts of each substance: $$y_i = \frac{n_i}{n_m}$$

$n_i$... amount of a specific component [mol]

$n_m$... total amount of substance of the mixture [mol]

3.) Calculate the pseudo-critical temperature and pressure: $$T_c' = y_1 \times T_{c1} + y_2 \times T_{c2} + y_3 \times T_{c3}$$

$$P_c' = y_1\times P_{c1} +y_2\times P_{c2} +y_3\times P_{c3}$$

$T_c$... critical temperature of a specific component [K]

$P_c$... critical pressure of a specific component [Pa]

4.) calculate the pseudo-reduced values for T and the volume: $$T_r = T / T_c'$$ $$v_r' = V / (R \times T_c' / P_c')$$

$T$... temperature in the container

$V$... volume of the gas phase [m$^3$]

$R$... universal gas constant = 8.314 J/(K mol)

5.) Now you need to look up $v_r$ and $T_r$ in a compressibility chart to get the value for a specific $Z$ factor, where their lines intersect.

6.) And finally calculate the pressure: $$P = ZRT / V$$

This will be more accurate than if you would use the standard way for ideal gases, especially for low temperatures

References: real gases

calculation

compressibility chart examples


text:
# 22a: Driven Damped Pendulum: Period Doubling, Chaos, Strange Attractors

Michael Fowler, UVa

Note: although this lecture follows naturally from the previous one (on the driven damped nonlinear oscillator), it is self-contained$—$it does not depend on that earlier material.

## Introduction

We’ve previously discussed the driven damped simple harmonic oscillator, and in the last lecture we extended that work (following Landau) to an anharmonic oscillator, adding a quartic term to the potential.  Here we make a different extension of the simple oscillator: we go to a driven damped pendulum.  (A weight at one end of a light rigid rod, the other end of the rod being at a fixed position but the rod free to rotate about that fixed point in a vertical plane).  That is, we replace the potential term $-{\omega }_{0}^{2}x$ in the linear oscillator with $-{\omega }_{0}^{2}\mathrm{sin}x,$ or rather $-{\omega }_{0}^{2}\mathrm{sin}\varphi ,$ to make clear we have an angular system.  Going to a driven damped pendulum leads to many surprises!

For a sufficiently weak driving force, the behavior of the driven damped pendulum is of course close to that of the driven damped linear oscillator, but on gradually increasing the driving force, at a certain strength the period of the response doubles, then, with further increase, it doubles again and again, at geometrically decreasing intervals, going to a chaotic (nonperiodic) response at a definite driving strength.  But that is not the end of the story$—$the chaotic response regime has a lot of structure: many points within a generally chaotic region are in fact nonchaotic, well-defined cyclical patterns. And, as we’ll see later, the response pattern itself can be fractal in nature, see for example the strange attractor discussed at the end of this lecture.  You can use the accompanying applet to generate this attractor and its cousins easily.

Obviously, this is a very rich subject, we provide only a brief overview.  We closely follow the treatment in Taylor’s Classical Mechanics, but with the addition of our applets for illustrative purposes, and also to encourage further exploration: the applets accurately describe the motion and exhibit the strange attractors in the chaotic regime. With the applets, it is easy to check how these strange attractors change (or collapse) on varying the driving force or the damping. We end with some discussion of the fractal dimension of the attractors, and how it relates to the dynamics, in particular to the rate of divergence of initially close trajectories, here following Baker and Gollub, Chaotic Dynamics.

### Equation of Motion

For the driven pendulum, the natural measure of the driving force is its ratio to the weight $mg,$ Taylor calls this the drive strength, so for driving force

$F\left(t\right)={F}_{0}\mathrm{cos}\omega t,$

the drive strength $\gamma$ is defined by

$\text{ }\gamma ={F}_{0}/mg.$

The equation of motion (with resistive damping force $-bv$ and hence resistive torque $-b{L}^{2}\stackrel{˙}{\varphi }$ )  is:

$m{L}^{2}\stackrel{¨}{\varphi }=-b{L}^{2}\stackrel{˙}{\varphi }-mgL\mathrm{sin}\varphi +LF\left(t\right),$

Dividing by $m{L}^{2}$ and writing the damping term $b/m=2\beta$ (to coincide with Taylor’s notation, his equation 12.12) we get (with ${\omega }_{0}^{2}=g/L$ )

$\stackrel{¨}{\varphi }+2\beta \stackrel{˙}{\varphi }+{\omega }_{0}^{2}\mathrm{sin}\varphi =\gamma {\omega }_{0}^{2}\mathrm{cos}\omega t.$

### Behavior on Gradually Increasing the Driving Force: Period Doubling

The driving force $\gamma ={F}_{0}/mg$ is the dimensionless ratio of drive strength to weight, so if this is small the pendulum will not be driven to large amplitudes, and indeed we find that after initial transients it settles to motion at the driving frequency, close to the linearized case. We would expect things to get more complicated when the oscillations have amplitude of order a radian, meaning a driving force comparable to the weight.  And indeed they do.

Here we’ll show that our applet reproduces the sequence found by Taylor as the driving strength is increased.

In the equation of motion

$\stackrel{¨}{\varphi }+2\beta \stackrel{˙}{\varphi }+{\omega }_{0}^{2}\mathrm{sin}\varphi =\gamma {\omega }_{0}^{2}\mathrm{cos}\omega t,$

we therefore choose Taylor’s values ${\omega }_{0}=1.5,\text{ }\text{ }2\beta =0.75,\text{ }\text{ }\omega =1,$ and gradually increase $\gamma$ from 0.9 to 1.0829, where chaos begins.

For $\gamma =0.9$ (see picture) the oscillation (after brief initial transients) looks like a sine wave, although it’s slightly flatter, and notice the amplitude (turquoise in box) is greater than $\pi /2,$ the positive and negative swings are equal in magnitude to five figure accuracy after five or so oscillations.

You can see this for yourself by opening the applet!  Click here. The applet plots the graph, and simultaneously shows the swinging pendulum (click bar at top right).

For $\gamma =1.06,$ there is a larger initial transient, vaguely resembling the three-cycle we’ll meet shortly.  In fact, this can be suppressed by starting at $\phi =-\pi /2,$ (bottom slider on the right), but there are still transients: the peaks are not uniform to five-figure accuracy until about forty cycles.

For $\gamma =1.0662,$ there are very long transients: not evident on looking at the graph (on the applet), but revealed by monitoring the amplitude readout (turquoise figures in the box on the graph), the value at the most recent extremum.  To get past these transients, set the applet speed to 50 (this does not affect accuracy, only prints a lot less dots).  run 350 cycles then pause and go to speed = 2.  You will find the peaks now alternate in height to the five-figure accuracy displayed, looks like period doubling$—$but it isn’t, run a few thousand cycles, if you have the patience, and you find all peaks have the same height. That was a long lived transient precursor of the period doubling transition.

Going to 1.0664, you’ll see both peaks and dips now alternating in amplitude, for a well-defined period 2, which lasts until 1.0792. (Look at 1.075 after 70 or so cycles$—$and set the initial angle at -90.)

Use the “ Red Line” slider just below the graph to check the heights of successive peaks or dips.

For $\gamma =1.0793,$ the period doubles again: successive peaks are now of order 0.4 radians apart, but successive “high” peaks are suddenly about 0.01 radians apart. (To see these better on the applet, magnify the vertical scale and move the red line.)

For $\gamma =1.0821$ there is a further doubling to an 8 cycle, then at 1.0827 to a 16 cycle.

Look at the graph for 1.0826: in particular, the red line near the bottom. Every fourth dip goes below the line, but the next lowest dips alternate between touching the line and not quite reaching it.  This is the 8 cycle.

It is known that the intervals between successive doublings decrease by a factor $\delta =4.6692,$ found universally in period doubling cascades, and called the Feigenbaum number, after its discoverer. Our five-figure accuracy is too crude to follow this sequence further, but we can establish (or at least make very plausible!) that beyond the geometric series limit at ${\gamma }_{c}=1.0829$ the periodicity disappears (temporarily, as we’ll see), the system is chaotic. (Of course, the values of $\gamma$ depend on the chosen damping parameter, etc., only the ratio of doubling intervals is universal.)

In fact, the full picture is complex: there are further intervals of periodicity, for example a 6 cycle at $\gamma =1.0845,$ pictured here.

### Different Attractors

The periodic solutions described above are called “attractors”: configurations where the system settles down after initially wandering around.

Clearly the attractors change with the driving strength, what is less obvious is that they may be different for different initial conditions.  Taylor shows that for $\gamma =1.077$ taking $\phi \left(0\right)=0$ gives a 3-cycle after transients, but $\phi \left(0\right)=-\pi /2$ gives a 2-cycle. (Easily checked with the applet!)

Looking at the initial wanderings, which can be quite different for very small changes in the driving strength (compare 1.0730 to 1.0729 and 1.0731, use speed 5, it doesn’t affect the accuracy). But you can see these initial wanderings include elements from both attractors.

Exercise: use the next applet to plot at the same time 1.0729 and 1.0730, all other parameters the same, speed 5.

Exercise:  Use the applet to find at what angle the transition from one attractor to the other takes place.  And, explore what happens over a range of driving strengths.

These are the simplest attractors: there are far more complex entities, called strange attractors, we’ll discuss later.

Exercises:  Try different values of the damping constant and see how this affects the bifurcation sequence.

### Sensitivity to Initial Conditions

Recall that for the linear damped driven oscillator, which can be solved exactly, we found that changing the initial conditions changed the path, of course, but the difference between paths decayed exponentially: the curves converged.

This illustration is from an earlier applet: the red and green curves correspond to different initial conditions, the white curve is the difference between the two, obviously exponentially decreasing, as can be verified analytically.

For the damped driven pendulum, the picture is more complicated.  For  $\gamma <{\gamma }_{c}=1.0829$ curves corresponding to slightly different initial conditions will converge (except, for example, $\gamma =1.077,$ where, as mentioned above, varying the initial angle at a certain point switches from a final three-cycle to a two-cycle).

### The Liapunov Exponent

For $\gamma >{\gamma }_{c},$ curves even with very small initial differences (say, 10-4 radians) separate exponentially, as ${e}^{\lambda t},\text{ }\lambda$ is called the Lyapunov exponent.

Bear in mind, though, that this is chaotic motion, the divergence is not as smooth as the convergence pictured above for the linear system. This graph is from our two-track applet.

Nevertheless, the (average) exponential nature can be discerned by plotting the logarithm of the difference against time:

This is from our log difference applet.  It is very close to Taylor’s fig 12.13. The red line slope gives the Lyapunov exponent.

(The red line has adjustable position and slope.)

### Plotting Velocity against Time

As discussed in Taylor, further increase in the driving force beyond the chaotic threshold can lead to brief nonchaotic intervals, such as that containing the six cycle at 1.0845 illustrated above, but there are two long stretches of nonchaotic behavior in Taylor’s parameter range, from 1.1098 to 1.1482 and from 1.3 to 1.48.

In the stronger driving force range, the pendulum is driven completely around in each cycle, so plotting the position against time gives a “rounded staircase” kind of graph.  Check this with the applet.

.

The solution is to plot velocity against time, and thereby discover that there is a repetition of the period doubling route to chaos at the upper end of this interval. Click to plot $d\varphi /dt$ instead of $\varphi .$

### State-Space Trajectories

It can be illuminating to watch the motion develop in time in the two-dimensional state space $\left(\varphi ,\stackrel{˙}{\varphi }\right).$ (Equally called phase space.)  See the State Space applet!

Now for a particle in a time-independent potential, specifying position and velocity at a given instant determines the future path$—$but that is not the case here, the acceleration is determined by the phase of the driving force, which is time varying, so the system really needs three parameters to specify its subsequent motion.

That means the phase space is really three-dimensional, the third direction being the driving phase, or, equivalently, time, but periodic with the period of the driving force. In this three-dimensional space, paths cannot cross, at any point the future path is uniquely defined. Our two-dimensional plots are projections of these three-dimensional paths on to the $\left(\varphi ,\stackrel{˙}{\varphi }\right)$ plane.

Here’s the 4-cycle at γ = 1.081, minus the first 20 cycles, to eliminate transients.

For γ = 1.0826, there is an 8-cycle. Run the applet for 40 cycles or so to get rid of transients, then look at the far-left end of the curve generated by running, say, to 200 cycles. It doesn’t look like there are 8 lines, but the outermost line and the two innermost lines are doubles (notice the thickness).  You can check this with the applet, by pausing at the axis and noticing the position readout: there are 8 different values that repeat.  (You don’t have to stop right on the axis$—$you can stop close to it, then use the step buttons to get there.)

For γ = 1.0830, the path is chaotic$—$but doesn’t look too different on this scale!  Check it out with the applet. The chaos becomes more evident on further increasing γ.   For γ =1.087 the pattern is “fattened” as repeated cycles differ slightly.

For further increase in γ, the orbital motion covers more territory: at γ = 1.1, here are the first three hundred or so cycles.

Plotting many orbits at high speed we find:

However, here the story gets complicated: it turns out that this chaotic series of orbits is in fact a transient, and after about 300 cycles the system skips to a 3-cycle, where it stays.  In fact, we have reached a range of γ (approximately 1.1098 to 1.1482) where after initial transients (which can be quite long) the motion is not chaotic, but a 3-cycle.

Obviously, there is a lot to explore for this system!   To get a better grasp of these complexities, we try a different approach, developed by Poincaré.

### Poincaré Sections

Looking at the above pictures, as we go from a single orbit to successive period doublings and then chaos, the general shape of the orbit doesn’t change dramatically (until that final three-cycle).  The interesting things are the doubling sequence, chaos, and attractors$—$perhaps we’re plotting too much information.

To focus on what’s essential, Poincaré plotted a single point from each cycle, this is now called a Poincaré section. To construct this, we begin with the t = 0 position, label it  ${P}_{0}=\left({\phi }_{0},{\stackrel{˙}{\phi }}_{0}\right)$.  Then add the point ${P}_{1}=\left({\phi }_{1},{\stackrel{˙}{\phi }}_{1}\right)$ precisely one cycle later, and so on$—$points one cycle apart, ${P}_{2}=\left({\phi }_{2},{\stackrel{˙}{\phi }}_{2}\right)$, etc.  Now, knowing the position $P=\left(\phi ,\stackrel{˙}{\phi }\right)$ in state space is not enough information to plot the future orbit$—$we also need to know the phase of the driving force. But by plotting points one cycle in time apart, they will all see the same phase starting force, so the transformation that takes us from ${P}_{0}$ to ${P}_{1}$ just repeats in going from ${P}_{1}$ to ${P}_{2},$ etc.

To see these single points in the State Space applet, click “Show Red Dot”: on running, the applet will show one dot per cycle red.

Thinking momentarily of the full three-dimensional phase space, the Poincaré section is a cross-section at the initial phase of the driving cycle. We’ve added to the applet a red dot phase option, to find the Poincaré section at a different phase. Doing this a few times, and looking at the movement of the red dots, gives a fuller picture of the motion.

So the Poincaré section, on increasing γ through the doubling sequence (and always omitting initial transients) goes from a single point to two points, to four, etc.

To see all this a bit more clearly, the next applet, called Poincaré Section, shows only one dot per cycle, but has a phase option so you can look at any stage in the cycle.

Exercise: check this all out with the Poincaré applet!  To see it well, click Toggle Origin/Scale, use the sliders to center the pattern, more or less, then scale it up.  Run for a few hundred cycles, then click Hide Trace followed by Show Trace to get rid of transients.

Start with γ = 1.0826.  The Poincaré section is eight points on a curve, at 1.0827 you can discern 16. By 1.0829, we have chaos, and the section has finite stretches of the curve, not just points.  It looks very similar at 1.0831, but$—$surprise$—$at 1.0830 it’s a set of points (32?) some doubled.  This tells us we’re in a minefield. There’s nothing smooth about chaos.

Apart from interruptions, as γ increases, the Poincaré section fills a curve, which then evolves towards the strange attractor shown by Taylor and others.  Going from 1.090 in steps to 1.095, the curve evolves a second, closely parallel, branch. At 1.09650 the two branches are equal in length, then at 1.097 the lower branch suddenly extends a long curve, looks much the same at 1.100, then by 1.15 we recognize the emergence of the strange attractor.

But in fact this is a simplified narrative: there are interruptions.  For example, at γ = 1.0960 there’s a 5-cycle, no chaos (it disappears on changing by 0.0001).  And at 1.12 we’re in the 3-cycle interval mentioned earlier.

Anyway, the strange attractor keeps reappearing, and at 1.5 it looks like this:

This looks rather thin compared to Taylor’s picture for the same γ: the reason is we have kept the strong damping.

Changing to $2\beta =0.375,$ we recover Taylor’s picture.  Like him, we show here (from the applet) the attractor and successive magnifications of one part, to give some hint of the fractal structure:

We have to admit that at the highest magnification Taylor’s picture is certainly superior to that generated by the applet, but does not contradict it$—$actually, it reassures us that the applet is reliable up to the level of definition that is self-evident on looking at it.

Exercise:  Use the last two applets to explore other regions of parameter space. What happens on varying the damping?  Going from $2\beta =0.375$ to $2\beta =0.75$ the attractor has the same general form but is much narrower. Then go in 0.001 steps to 0.76.  For some values you see the attractor, but for others a cycle, different lengths. There is a two-cycle from 0.76 to 0.766, then a four cycle, then at 0.78 back to a two-cycle, at 0.814 a one-cycle.  If you up the driving strength, and the damping, you can find extremely narrow attractors.

Exercise:  Baker and Gollub, Chaotic Dynamics, page 52, give a sequence of Poincaré sections at $0.2\pi$ intervals for $\gamma =1.5$ and $2\beta =0.75.$ You can check them (and intermediate ones) with the applet, and try to imagine putting them in a stack to visualize the full three-dimensional attractor!

## Lyapunov Exponents and Dimensions of Strange Attractors

### Scaling the Attractor

Looking at the strange attractor pictured in the previous section, we found that on magnifying a small part of it we saw the same kind of structure the attractor has as a whole: if we look at a $\gamma =1.5$ attractor, with damping 0.75, there are long thin stretches, they end by looping over. We notice that halving the damping to 0.375 fattens the previously quasi-one-dimensional stretches and reveals complicated looping at several levels. (Remember that what we are looking at here is a Poincaré section of the attractor, the other dimension is periodic time (or driver phase, the same thing) so a curve here is a section of some sheet).  If more computing power is used, going to smaller and smaller scales, it turns out that the magnified tiny part of the attractor looks much the same as the attractor.  This kind of scale invariance is a characteristic of a fractal.  A mysterious aspect of fractals is their dimensionality. Look at the strange attractor. There are no places where it solidly fills a stretch of two-dimensional space, this is clearer on going to greater and greater magnification: we see more and more one-dimensional structures, with no end, so it surely has dimension less than two, but greater than one$—$how do we make sense of that?

### Fractals: the Cantor Set

To try to find a generalized concept of dimension of a set (i.e. not just an integer), we begin with perhaps the simplest example of a fractal, the Cantor set:  take the numbers between 0 and 1, and cut out the middle third. You now have two strips of numbers, from 0 to 1/3, and from 2/3 to 1. For each of those strips, cut out the middle third. You now have four strips$—$cut out the middle third of each of them (it might help to draw this). Do this forever. What’s left is the Cantor set. You can see this is scale invariant: after doing this many times, take one of the remaining tiny strips, what happens to it on continuing the process is identical (scaled down suitably) to what happened to the initial strip.

How big is this Cantor set? At each step, we cut the total length of line included by 2/3. Since (2/3)n goes to zero as n goes to infinity, it clearly has size zero, right?  But clearly there’s more to the Cantor set than there is to a single point, or for that to matter a finite number of points.  What about a countably infinite number of points$—$for example, the rational numbers between 0 and 1? Well, you can write them out in a list, ordered by increasing denominators, and for one denominator by increasing numerators. Then you can put them one by one into tiny intervals, $½$ goes into an interval of length $\epsilon ,$ 1/3 in an interval $\epsilon /2,$ 2/3 in one of length $\epsilon /{2}^{2},$ $¼$ in $\epsilon /{2}^{3},$ and so on, the total length of the infinite number of intervals being $2\epsilon ,$ so all the rationals can be covered by an arbitrarily small set of intervals. Can we count in order the numbers in the Cantor set in the same way? The answer is no, and to see why think first about all the numbers between 0 and 1, rationals and irrationals.  If you make an infinite list of them, I can show you missed some out: I just take your list and write down a decimal that differs from your nth number in the nth place. So we can’t put all the numbers in the interval in little boxes that add to zero, which is obvious anyway!

But now to the Cantor set: suppose we write all numbers between 0 and 1 using base 3, instead of the traditional base 10. That is, each number is a string of 0’s, 1’s and 2’s. Then the Cantor set is all those numbers that don’t have any 1’s, such as 0.2, 0.02202, etc.  (Check this yourself.)  But the number of these numbers is exactly the same as all the numbers between 0 and 1 in binary notation! So surely the Cantor set has dimension 1?  (These infinities are tricky.)

The bottom line from the above argument is that we can plausibly argue both that the Cantor set has dimension 0, and that it has dimension 1. To understand and categorize fractals better, we need a working definition of dimension for fractals. One approach is the capacity dimension.

### Dimensions: Capacity and Correlation

Suppose we cover the interval 0,1 with a set of small boxes, length $\epsilon ,$ there are clearly $N\left(\epsilon \right)=1/\epsilon$ such boxes (assume it’s an integer).  Now consider a subset of the numbers between 0 and 1, choose $\epsilon ,$ and find how many boxes $N\left(\epsilon \right)$ are necessary to cover this subset. The capacity dimension is defined as

${d}_{C}=\underset{\epsilon \to 0}{\mathrm{lim}}\frac{\mathrm{log}N\left(\epsilon \right)}{\mathrm{log}\left(1/\epsilon \right)}.$

For simplicity, we choose   $1/\epsilon =3,9,27,\dots$ so the necessary numbers of boxes to cover the Cantor set described in the previous section are 2, 4, 8… out of total numbers of boxes 3, 9, 27, … . Therefore

${d}_{C}\left(\text{Cantor}\right)=\frac{n\mathrm{log}2}{n\mathrm{log}3}=\frac{\mathrm{log}2}{\mathrm{log}3},$

of course between 0 and 1.  (There are many ways to define dimensionality of sets of numbers$—$this definition gives zero for a finite set of points, and one for all the numbers between 0 and 1, but also 1 for the set of rationals, which we’ve shown above can be covered by an infinite set of intervals of arbitrarily small total length.)

Another measure used is the correlation dimension, in which for a large number of points (such as our representation of the attractor) a correlation integral $C\left(r\right)$ is defined as the total number of pairs of two points less than $r$ apart. For small $r,$ this goes as a power ${r}^{\nu },$ and it turns out that in many cases  $\nu$ is close to the capacity definition of the fractal dimension.  (Grassberger and Procaccia.)

### Time Development of Systems in Phase Space

Recall first that the state space or phase space we have been plotting is really a projection of the full orbit space into two dimensions, the third dimension necessary to predict future motion being the phase of the sine-wave driving force, so this is just time (although of course cyclic).

Suppose now we populate this three-dimensional space with many points, like a gas, each representing a driven damped pendulum.  As time goes on the gas will flow, each gas atom’s path completely determined, and no two will ever be at the same point in this full space (except perhaps asymptotically at infinite time).

Take now a small volume, say a cube with sides parallel to the axes, containing many points. Consider first the undamped system: then Liouville’s Theorem  (link to my lecture) tells us that as time goes on the cube will generally distort, but it will not change in volume.  In other words, the gas of systems flows like an incompressible fluid.  (Details of the derivation are given in the linked lecture$—$briefly, the motions of the sides in time $\Delta t$ come from the equations of motion, etc.)

However, if the system has damping$—$as ours does$—$the same analysis leads to the conclusion that the volume the systems occupy in phase space (remember, this is now three-dimensional) shrinks at a rate determined by the damping.  As a trivial example, think of undriven damped pendula$—$they will all tend to the low point rest position. Lightly driven pendula will go to a one-dimensional cycle.

We can prove this shrinkage from the equation of motion:

$\stackrel{¨}{\varphi }+2\beta \stackrel{˙}{\varphi }+{\omega }_{0}^{2}\mathrm{sin}\varphi =\gamma {\omega }_{0}^{2}\mathrm{cos}\omega t.$

In the three-dimensional phase space, a pendulum’s position can be written in coordinates $\left(\varphi ,\stackrel{˙}{\varphi },\psi \right)$ where $\psi =\omega t,$ the driving phase, between 0 and $2\pi .$

The local phase space velocity $\stackrel{\to }{F}$ can be written in terms of the coordinates (this is just the above equation rewritten!):

$\begin{array}{l}\frac{\partial \varphi }{\partial t}=\stackrel{˙}{\varphi },\\ \frac{\partial \stackrel{˙}{\varphi }}{\partial t}=-2\beta \stackrel{˙}{\varphi }-{\omega }_{0}^{2}\mathrm{sin}\varphi +\gamma {\omega }_{0}^{2}\mathrm{cos}\psi ,\\ \frac{\partial \psi }{\partial t}=\omega ,\end{array}$

This is therefore the local velocity of the atoms of the gas (meaning the systems), and it is trivial to check that

$\stackrel{\to }{\nabla }\cdot \stackrel{\to }{F}=-2\beta .$

This means that if we have a small sphere containing many points corresponding to systems (a “gas” of systems) then the volume of the (now distorting) sphere enclosing those points is decreasing in volume at an exponential rate $V\left(t\right)={V}_{0}{e}^{-2\beta t}.$

### Relating this Picture to Lyapunov Exponents

Continuing to think about the development of a small sphere (containing many points corresponding to systems) in phase space, it will be moving along an orbit, but at the same time distorting, let’s say to an ellipsoid as an initial first approximation, and tumbling around.  In the chaotic regime, we know it must be growing in some direction, at least on average (the rates will vary along the orbit) because we know that points initially close together separate on average at an exponential rate given by the first Lyapunov exponent, $\propto {e}^{{\lambda }_{1}t}.$ We’ll make the simplifying assumption that the ellipsoid has its axes initially varying in time as $\left({e}^{{\lambda }_{1}t},{e}^{{\lambda }_{2}t},{e}^{{\lambda }_{3}t}\right),$ with ${\lambda }_{1}>{\lambda }_{2}>{\lambda }_{3}.$ From the result above $\stackrel{\to }{\nabla }\cdot \stackrel{\to }{F}=-2\beta ,$ we conclude that

${\lambda }_{1}+{\lambda }_{2}+{\lambda }_{3}=-2\beta .$

We need to say something more about ${\lambda }_{1}.$  We’re taking it as defined by the growth rate of distance between trajectories after any initial transients but before the distance is comparable to the size of the system (finding this interval plausibly has been termed a “dark art”). We envision our initially small sphere of gas elongating and tumbling around as it moves along. Hopefully its rate of elongation correlates well with what we actually measure, that is, the rate of growth of net $\varphi$ displacement, the coordinate separation of two initially close orbits, which we plot and approximately fit with an exponential, $\propto {e}^{{\lambda }_{1}t}.$

For the pendulum, the $\psi$ direction is just time, not scaled, so ${\lambda }_{3}={\lambda }_{\psi }=0.$ Then necessarily ${\lambda }_{2}<0$ to satisfy the damping equation.

So taking a local (in phase space) collection of systems, those inside a given closed surface, like a little sphere, and following their evolution in time in the chaotic regime, the sphere will expand in one direction; a direction, however, that varies with time, but contract or stay constant in the other directions.  As the surface grows this gets more complicated because it’s confined to a finite total phase space.  And it continues to expand at the same rate as time goes on, so the continual increase in surface must imply tighter and tighter foldings to stay in the phase space.  And this is what the strange attractor looks like.

### A Fractal Conjecture

In 1979, Kaplan and Yorke conjectured that the dimensionality of the strange attractor followed from the Lyapunov exponents taking part in its creation. In our case$—$the driven damped pendulum$—$there are only two relevant exponents, ${\lambda }_{1}>0,\text{ }{\lambda }_{2}<0$ and  ${\lambda }_{1}+{\lambda }_{2}=-2\beta .$

A plausibility argument is given in Baker and Gollub’s book, Chaotic Dynamics.  They define a Lyapunov dimension ${d}_{L}$ of the attractor by

${d}_{L}=\underset{\epsilon \to 0}{\mathrm{lim}}\left[\frac{d\left(\mathrm{log}N\left(\epsilon \right)\right)}{d\left(\mathrm{log}\left(1/\epsilon \right)\right)}\right],$

exactly analogous to the definition of capacity dimension in the previous section.

Now, as time passes a small square element will have its area multiplied by a factor ${e}^{\left({\lambda }_{1}+{\lambda }_{2}\right)t}.$ (No scaling takes place in the third (time) direction.)  At the same time, they argue that the length unit $\epsilon$ changes as ${e}^{{\lambda }_{2}t}.$ Then $N\left(\epsilon \right)$ is the area ${\epsilon }_{0}^{2}{e}^{\left({\lambda }_{1}+{\lambda }_{2}\right)t}$ divided by the shrinking basic area ${\epsilon }_{0}^{2}{e}^{2{\lambda }_{2}t}.$  The differential of $\mathrm{log}N\left(\epsilon \right)$ is ${\lambda }_{1}-{\lambda }_{2},$ that of $\mathrm{log}\left(1/\epsilon \right)$ is $-{\lambda }_{2},$ so their argument gives

${d}_{L}=1-\frac{{\lambda }_{1}}{{\lambda }_{2}}.$

The Lyapunov applet is designed to measure ${\lambda }_{1}$ by tracking separation of initially close trajectories.  Try it a few times: it becomes clear that there is considerable uncertainty in this approach. Given ${\lambda }_{1},{\lambda }_{2}$ follows from ${\lambda }_{1}+{\lambda }_{2}=-2\beta .$ There are various ways to assign a dimension to the attractor, such as the capacity and correlation dimensions mentioned above. Various attempts to verify this relationship have been made, but the uncertainties are considerable, and although results seem to be in the right ballpark, the results are off by ten or twenty percent typically.  It seems there is work still to be done on this fascinating problem.

Recommended reading:  chapter 5 of Baker and Gollub, Chaotic Dynamics.  The brief discussion above is based on their presentation.

format: plain


text:
[texhax] multi-page framed environments (was Re: )

Wed May 18 15:30:49 CEST 2005

On May 18, 2005, at 9:02 AM, Anna Medvedovsky wrote:
> I'm trying to define an environment that will put a frame around its
> contents, which may be more than one paragraph long -- something not
> unlike a graybox in a book.

> This is such a simple problem that there must be a solution! What
> could it be?! I'm inexperienced with boxes; perhaps there's another
> command I should use? Perhaps I need to draw the box manually? Perhaps
> I should use a table environment?

The memoir documentclass affords a framed'' environment which does
this sort of thing.

\documentclass{memoir}
\begin{document}
\begin{framed}
This is my text. This is more text.

Some more text in a box.
\end{framed}
\end{document}

ISTR there was discussion of the code being based on a framed''
package on CTAN, but I may not have the name precisely.

William

--


text:
## Transient heat conduction - partial differential equations

| categories: pde | tags: heat transfer | View Comments

We solved a steady state BVP modeling heat conduction. Today we examine the transient behavior of a rod at constant T put between two heat reservoirs at different temperatures, again T1 = 100, and T2 = 200. The rod will start at 150. Over time, we should expect a solution that approaches the steady state solution: a linear temperature profile from one side of the rod to the other.

$$\frac{\partial u}{\partial t} = k \frac{\partial^2 u}{\partial x^2}$$

at $$t=0$$, in this example we have $$u_0(x) = 150$$ as an initial condition. with boundary conditions $$u(0,t)=100$$ and $$u(L,t)=200$$.

In Matlab there is the pdepe command. There is not yet a PDE solver in scipy. Instead, we will utilze the method of lines to solve this problem. We discretize the rod into segments, and approximate the second derivative in the spatial dimension as $$\frac{\partial^2 u}{\partial x^2} = (u(x + h) - 2 u(x) + u(x-h))/ h^2$$ at each node. This leads to a set of coupled ordinary differential equations that is easy to solve.

Let us say the rod has a length of 1, $$k=0.02$$, and solve for the time-dependent temperature profiles.

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

N = 100  # number of points to discretize
L = 1.0
X = np.linspace(0, L, N) # position along the rod
h = L / (N - 1)

k = 0.02

def odefunc(u, t):
dudt = np.zeros(X.shape)

dudt[0] = 0 # constant at boundary condition
dudt[-1] = 0

# now for the internal nodes
for i in range(1, N-1):
dudt[i] = k * (u[i + 1] - 2*u[i] + u[i - 1]) / h**2

return dudt

init = 150.0 * np.ones(X.shape) # initial temperature
init[0] = 100.0  # one boundary condition
init[-1] = 200.0 # the other boundary condition

tspan = np.linspace(0.0, 5.0, 100)
sol = odeint(odefunc, init, tspan)

for i in range(0, len(tspan), 5):
plt.plot(X, sol[i], label='t={0:1.2f}'.format(tspan[i]))

# put legend outside the figure
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.xlabel('X position')
plt.ylabel('Temperature')

# adjust figure edges so the legend is in the figure
plt.savefig('images/pde-transient-heat-1.png')

# Make a 3d figure
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()

SX, ST = np.meshgrid(X, tspan)
ax.plot_surface(SX, ST, sol, cmap='jet')
ax.set_xlabel('X')
ax.set_ylabel('time')
ax.set_zlabel('T')
ax.view_init(elev=15, azim=-124) # adjust view so it is easy to see
plt.savefig('images/pde-transient-heat-3d.png')

# animated solution. We will use imagemagick for this

# we save each frame as an image, and use the imagemagick convert command to
# make an animated gif
for i in range(len(tspan)):
plt.clf()
plt.plot(X, sol[i])
plt.xlabel('X')
plt.ylabel('T(X)')
plt.title('t = {0}'.format(tspan[i]))
plt.savefig('___t{0:03d}.png'.format(i))

import commands
print commands.getoutput('convert -quality 100 ___t*.png images/transient_heat.gif')
print commands.getoutput('rm ___t*.png') #remove temp files

This version of the graphical solution is not that easy to read, although with some study you can see the solution evolves from the initial condition which is flat, to the steady state solution which is a linear temperature ramp.

The 3d version may be easier to interpret. The temperature profile starts out flat, and gradually changes to the linear ramp.

Finally, the animated solution.

org-mode source

## Boundary value problem in heat conduction

| categories: bvp | tags: heat transfer | View Comments

For steady state heat conduction the temperature distribution in one-dimension is governed by the Laplace equation:

$$\nabla^2 T = 0$$

with boundary conditions that at $$T(x=a) = T_A$$ and $$T(x=L) = T_B$$.

The analytical solution is not difficult here: $$T = T_A-\frac{T_A-T_B}{L}x$$, but we will solve this by finite differences.

For this problem, lets consider a slab that is defined by x=0 to x=L, with $$T(x=0) = 100$$, and $$T(x=L) = 200$$. We want to find the function T(x) inside the slab.

We approximate the second derivative by finite differences as

$$f''(x) \approx \frac{f(x-h) - 2 f(x) + f(x+h)}{h^2}$$

Since the second derivative in this case is equal to zero, we have at each discretized node $$0 = T_{i-1} - 2 T_i + T_{i+1}$$. We know the values of $$T_{x=0} = \alpha$$ and $$T_{x=L} = \beta$$.

$A = \left [ \begin{array}{ccccc} % -2 & 1 & 0 & 0 & 0 \\ 1 & -2& 1 & 0 & 0 \\ 0 & \ddots & \ddots & \ddots & 0 \\ 0 & 0 & 1 & -2 & 1 \\ 0 & 0 & 0 & 1 & -2 \end{array} \right ]$

$x = \left [ \begin{array}{c} T_1 \\ \vdots \\ T_N \end{array} \right ]$

$b = \left [ \begin{array}{c} -T(x=0) \\ 0 \\ \vdots \\ 0 \\ -T(x=L) \end{array} \right]$

These are linear equations in the unknowns $$x$$ that we can easily solve. Here, we evaluate the solution.

import numpy as np

#we use the notation T(x1) = alpha and T(x2) = beta
x1 = 0; alpha = 100
x2 = 5; beta = 200

npoints = 100

# preallocate and shape the b vector and A-matrix
b = np.zeros((npoints, 1));
b[0] = -alpha
b[-1] = -beta

A = np.zeros((npoints, npoints));

#now we populate the A-matrix and b vector elements
for i in range(npoints ):
for j in range(npoints):
if j == i: # the diagonal
A[i,j] = -2
elif j == i - 1: # left of the diagonal
A[i,j] = 1
elif j == i + 1: # right of the diagonal
A[i,j] = 1

# solve the equations A*y = b for Y
Y = np.linalg.solve(A,b)

x = np.linspace(x1, x2, npoints + 2)
y = np.hstack([alpha, Y[:,0], beta])

import matplotlib.pyplot as plt

plt.plot(x, y)

plt.plot(x, alpha + (beta - alpha)/(x2 - x1) * x, 'r--')

plt.xlabel('X')
plt.ylabel('T(X)')
plt.legend(('finite difference', 'analytical soln'), loc='best')
plt.savefig('images/bvp-heat-conduction-1d.png')

org-mode source

## Finding the nth root of a periodic function

| categories: nonlinear algebra | tags: heat transfer | View Comments

There is a heat transfer problem where one needs to find the n^th root of the following equation: $$x J_1(x) - Bi J_0(x)=0$$ where $$J_0$$ and $$J_1$$ are the Bessel functions of zero and first order, and $$Bi$$ is the Biot number. We examine an approach to finding these roots.

First, we plot the function.

from scipy.special import jn, jn_zeros
import matplotlib.pyplot as plt
import numpy as np

Bi = 1

def f(x):
return x * jn(1, x) - Bi * jn(0, x)

X = np.linspace(0, 30, 200)
plt.plot(X, f(X))
plt.savefig('images/heat-transfer-roots-1.png')

You can see there are many roots to this equation, and we want to be sure we get the n^{th} root. This function is pretty well behaved, so if you make a good guess about the solution you will get an answer, but if you make a bad guess, you may get the wrong root. We examine next a way to do it without guessing the solution. What we want is the solution to $$f(x) = 0$$, but we want all the solutions in a given interval. We derive a new equation, $$f'(x) = 0$$, with initial condition $$f(0) = f0$$, and integrate the ODE with an event function that identifies all zeros of $$f$$ for us. The derivative of our function is $$df/dx = d/dx(x J_1(x)) - Bi J'_0(x)$$. It is known (http://www.markrobrien.com/besselfunct.pdf) that $$d/dx(x J_1(x)) = x J_0(x)$$, and $$J'_0(x) = -J_1(x)$$. All we have to do now is set up the problem and run it.

from pycse import *  # contains the ode integrator with events

from scipy.special import jn, jn_zeros
import matplotlib.pyplot as plt
import numpy as np

Bi = 1

def f(x):
"function we want roots for"
return x * jn(1, x) - Bi * jn(0, x)

def fprime(f, x):
"df/dx"
return x * jn(0, x) - Bi * (-jn(1, x))

def e1(f, x):
"event function to find zeros of f"
isterminal = False
value = f
direction = 0
return value, isterminal, direction

f0 = f(0)
xspan = np.linspace(0, 30, 200)

x, fsol, XE, FE, IE = odelay(fprime, f0, xspan, events=[e1])

plt.plot(x, fsol, '.-', label='Numerical solution')
plt.plot(xspan, f(xspan), '--', label='Analytical function')
plt.plot(XE, FE, 'ro', label='roots')
plt.legend(loc='best')
plt.savefig('images/heat-transfer-roots-2.png')

for i, root in enumerate(XE):
print 'root {0} is at {1}'.format(i, root)

plt.show()

root 0 is at 1.25578377377
root 1 is at 4.07947743741
root 2 is at 7.15579904465
root 3 is at 10.2709851256
root 4 is at 13.3983973869
root 5 is at 16.5311587137
root 6 is at 19.6667276775
root 7 is at 22.8039503455
root 8 is at 25.9422288192
root 9 is at 29.081221492

You can work this out once, and then you have all the roots in the interval and you can select the one you want.


text:
Indico search will be reestablished in the next version upgrade of the software: https://getindico.io/roadmap/

#### This search is only for public events. Restricted events are not available.

IMPORTANT! Indico has been upgraded. Please let us know as soon as possible if you find any issues and email indico-support@fnal.gov

# New Perspectives 2018

18-19 June 2018
Fermilab, Wilson Hall
US/Central timezone

## Calorimeter Clustering Studies for the Mu2e Experiment at Fermilab

Jun 19, 2018, 9:45 AM
15m
One West (Fermilab, Wilson Hall)

### One West

#### Fermilab, Wilson Hall

Oral Presentation

Emma Castiglia

### Description

Mu2e (muon-to-electron-conversion), an Intensity Frontier experiment at Fermilab set to begin data taking in 2022, will search for flavor violation in charged leptons with a muon decaying to only an electron without any neutrinos. Protons will be accelerated by the 8GeV proton beam before reaching Mu2e’s production target. From there, pions enter an S-shaped Transport Solenoid and travel through a series of absorbers and collimators to select low-energy negatively charged muons from pion decay that will hit the stopping target. Downstream from the stopping target are a curved straw tracker and two annular disk crystal calorimeters. The signal will be a mono-energetic $e^{-}$ with an energy close to the rest mass of a muon. The Mu2e expected single event sensitivity is $2.87e^{-17}$. The complimentary trigger strategy will take advantage of both online track searches from the tracker and calorimeter reconstruction. The trigger requires a rejection factor greater than 100 and must make its decision within a 3ms time budget. In order to provide high efficiency on signal events, our calorimeter clustering algorithms must be accurate and meet the timing performance. I will discuss the expected performance of two clustering algorithms: a fast algorithm to be used at the trigger level and a more accurate clustering algorithm to be used in the offline reconstruction. The fast algorithm defines clusters by only selecting crystals adjacent to a high energy crystal seed. The full-blown algorithm can find smaller clusters separated from the main cluster, within a distance and timing window, and combine them into one cluster. We have analyzed the topology of the reconstructed clusters by the two algorithms on events with only a conversion electron and combined signal and background events. The goal of these studies is to improve the performance of the fast algorithm to better reconstruct the cluster energy and work within the strict timing budget.

### Co-authors

Gianantonio Pezzullo (Yale University) Sarah Demers (Yale)

 Slides


text:
# 3.12. Hard Multiple Choice Questions¶

These problems are harder than most of those that you will usually see on the AP CS A exam.

3-11-1: Which of the following would be the correct result from the following expression: 150 (in decimal) - 21 (in octal) + 101 (in binary) - E (in hexadecimal)?
• 186
• This would be true if it was 150 + 17 (21 in octal) + 5 (101 in binary) + 14 (E in hexadecimal) which is 186.
• 123
• This would be true if it was 150 - 17 (21 in octal) + 5 (101 in binary) - 15 (F in hexadecimal) which is 123.
• 125
• This would be true if it was 150 - 17 (21 in octal) + 5 (101 in binary) - 13 (D in hexadecimal) which is 123.
• 168
• How did you get this? Maybe look up how to convert from octal and hexadecimal to decimal.
• 124
• This is 150 - 17 (21 in octal) + 5 (101 in binary) - 14 (E in hexadecimal) which is 124.

The following question assumes you know about arrays and methods. You can skip it if you haven’t covered these yet and come back when you have.

3-11-2: Given the following code are the contents of itemArray and val after a call of mod(itemArray,val)?

int[] itemArray = {9, 8, 7, 6};
int val = 5;

public static void mod(int[] a, int value)
{
for (int i=0; i < a.length; i++)
{
a[i] = i;
}
value = a[a.length-1];
}

• itemArray = {0, 1, 2, 3} and val = 3;
• This would be true if Java used pass by reference rather than pass by value (it creates copies of the values that are passed).
• itemArray = {0, 1, 2, 3} and val = 5;
• Java passes parameters by copying the value. With an array it creates a copy of the object reference. So, mod will change the itemArray, but val won't change since mod only changes the copy of the primitive value.
• itemArray = {0, 0, 0, 0} and val = 0;
• How could this have happened?
• itemArray = {9, 8, 7, 6} and val = 3;
• Java passes parameters by passing the values this means that the contents of itemArray will be changed by the mod method, but val won't change.
• itemArray = {9, 8, 7, 6} and val = 5;
• Java passes parameters by passing the values this means that the contents of itemArray will be changed by the mod method.
Next Section - 3.13. Variables - Summary


text:
## ArrayFire v3.5 Official Release

Today we are pleased to announce the release of ArrayFire v3.5, our open source library of parallel computing functions supporting CUDA, OpenCL, and CPU devices. This new version of ArrayFire improves features and performance for applications in machine learning, computer vision, signal processing, statistics, finance, and more. This release focuses on thread-safety, support for simple sparse-dense arithmetic operations, canny edge detector function, and a genetic algorithm example. A complete list of ArrayFire v3.5 updates and new features are found in the product Release Notes. Thread Safety ArrayFire now supports threading programming models. This is not intended to improve the performance since most of the parallelism is happening on the device, but it does allow you to use multiple devices in …

## Benchmarking parallel vector libraries

There are many open source libraries that implement parallel versions of the algorithms in the C++ standard template libraries. Inevitably we get asked questions about how ArrayFire compares to the other libraries out in the open. In this post we are going to compare the performance of ArrayFire to that of BoostCompute, HSA-Bolt, Intel TBB and Thrust. The benchmarks include the following commonly used vector algorithms across 3 different architectures. Reductions Scan Transform The following setup has been used for the benchmarking purposes. The code to reproduce the benchmarks is linked at the bottom of the post. The hardware used for the benchmarks is listed below: NVIDIA Tesla K20 AMD FirePro S10000 Intel Xeon E5-2560v2 Background ArrayFire ArrayFire provides high …

## ArrayFire v3.0 is here!

Today we are pleased to announce the release of ArrayFire v3.0. This new version features major changes to ArrayFire’s visualization library, a new CPU backend, and dense linear algebra for OpenCL devices. It also includes improvements across the board for ArrayFire’s OpenCL backend. A complete list ArrayFire v3.0 updates and new features can be found in the product Release Notes. With over 8 years of continuous development, the open source ArrayFire library is the top CUDA and OpenCL software library. ArrayFire supports CUDA-capable GPUs, OpenCL devices, and other accelerators. With its easy-to-use API, this hardware-neutral software library is designed for maximum speed without the hassle of writing time-consuming CUDA and OpenCL device code. With ArrayFire’s library functions, developers can maximize …

## Intel OpenCL performance: 3rd generation hardware

Brian Kloppenborg 1 Comment

Introduction With Intel CPUs making up nearly 80% of the CPU market and 66% of computers using integrated graphics one can easily argue that integrated graphics devices represent one of the greatest markets for GPU-accelerated computing. Here at ArrayFire, we have long recognized the potential of these devices and offer built-in support for Intel CPUs, GPUs, and AMD APUs in the OpenCL backend of our ArrayFire GPU computing library. Yet one common theme for debate in the office has been how the hardware performs on different operating systems with different drivers across hardware revisions. To answer these questions (and, perhaps, to win some intra-office geek cred) I decided to write a series of blog posts about Intel’s GPU OpenCL performance. In this first installment I will compare the performance …

## Using zero-copy buffers on integrated GPUs

Brian Kloppenborg 1 Comment

One of the most powerful aspects of parallel program on integrated GPUs is taking advantage of shared memory and caches. The best example of this is sharing common data between the CPU and GPU via. zero-copy buffers. This technique permits your program to avoid the O(N) cost of copying data to/from the GPU. This feature is particularly useful for applications that deal with real-time data streams, like video processing.

## Machine Learning with ArrayFire: Linear Classifiers

Linear classifiers perform classification based on the linear combinition of the component features. Some examples of Linear Classifiers include: Naive Bayes Classifier, Linear Discriminant Analysis, Logistic Regression and Perceptrons. ArrayFire’s easy to use API enables users to write such classifiers from scratch fairly easily. In this post, we show how you can map mathematical equations to ArrayFire code and implement them from scratch. Naive Bayes Classifier Perceptron Naive Bayes Classifier Naive bayes classifier is a probabilistic classifier that assumes all the features in a feature vector are independent of each other. This assumption simplifies the bayes rule to a simple multiplication of probabilities as show below. First we start with the simple Baye’s rule.  p(C_k | x) = \frac{p(C_k)}{p(x)} …


text:
Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

A risk stratification tool for hospitalisation in Australia using primary care data

Abstract

Predictive risk models using general practice (GP) data to predict the risk of hospitalisation have the potential to identify patients for targeted care. Effective use can help deliver significant reductions in the incidence of hospitalisation, particularly for patients with chronic conditions, the highest consumers of hospital resources. There are currently no published validated risk models for the Australian context using GP data to predict hospitalisation. In addition, published models for other contexts typically rely on a patient’s history of prior hospitalisations, a field not commonly available in GP information systems, as a predictor. We present a predictive risk model developed for use by GPs to assist in targeting coordinated healthcare to patients most in need. The algorithm was developed and validated using a retrospective primary care cohort, linked to records of hospitalisation in Victoria, Australia, to predict the risk of hospitalisation within one year. Predictors employed include demographics, prescription history, pathology results and disease diagnoses. Prior hospitalisation information was not employed as a predictor. Our model shows good performance and has been implemented within primary care practices participating in Health Care Homes, an Australian Government initiative being trialled for providing ongoing comprehensive care for patients with chronic and complex conditions.

Introduction

The growing burden of chronic conditions (also known as non-communicable diseases), is now responsible for 70% of deaths globally1. Chronic conditions affect one in two people in developed countries and the World Health Organisation attributes three quarters of global chronic condition-related deaths to developing countries1. There is a strong imperative to address this disproportionate burden of chronic conditions on health care globally.

Worldwide efforts to move from episodic to integrated and coordinated care have delivered significant improvements in the management of chronic conditions in the primary care sector, demonstrating benefits to patients and the healthcare system2,3,4,5,6,7. To ensure optimal use of limited healthcare resources, it is crucial that such programs and resources target patients who would otherwise be admitted for an unplanned or avoidable hospitalisation.

The best performing risk stratification algorithms for primary care settings have generally been developed in locations where large linked datasets are practical and available. These datasets can offer good coverage of data variables and include information from outside the primary care setting, such as the patient’s history of previous hospitalisations, a highly significant predictor in such models8,9.

However, such linked data sets are not available in most health systems throughout the world. For example, in Australia there are no simple mechanisms to link patient records across primary and acute care settings. Furthermore, patient information about previous hospitalisations is not available in most primary care information systems. In such cases, there is a need for new algorithms based on data from the setting in which they will be used8.

Health Care Homes is a new Australian Government initiative with ongoing comprehensive care provided for up to 65,000 patients with chronic and complex conditions in up to 200 primary care clinics and Aboriginal Community Controlled Health Services (their ‘Health Care Home’). Patients eligible for the 22-month trial are identified in the general practice (GP) using an algorithm that predicts the risk of a patient being hospitalised over the next 12 months. By targeting services to patients with an imminent risk of hospitalisation, a significant reduction in the burden on the acute care health system is expected while ensuring valuable resources are allocated to programs and patients in greatest need.

This study describes the development and validation of a prediction model using Australian hospital and GP data to identify patients at risk of hospitalisation for the Health Care Homes trial. The developed prediction model employs routinely collected data from primary care clinics. We assume patient history of hospitalisation is not available as a predictor as this is typically not available in GP clinic databases. The algorithm developed and validated in this study is designed for use in Australian primary care. However, it also illustrates how risk stratification algorithms can be developed and validated to perform adequately in settings where hospitalisation history linked to acute care data may not be readily available.

Methods

Study Design and Data Sources

The work described here uses a retrospective cohort comprised of primary care patients whose data is linked with records of hospitalisations. Primary care data for the calendar years 2007 to 2016 was obtained from 29 general practice medical centres in Victoria, Australia. Hospitalisation data for patients admitted to public hospitals during the calendar years 2012 to 2016 in Victoria was provided by the Department of Health and Human Services (DHHS) Victoria.

Ethics approval for this study was provided by the CSIRO Health and Medical Human Research Ethics Committee (CHM HREC) Low Risk Review Panel – (Proposal LR 7/2017). All methods were performed in accordance with the guidelines and regulations relevant to the ethics approval. The ethics committee approved a waiver of consent for several factors including the nature of study, the size of the dataset, the fact that the data had been previously collected, and because the research posed minimum intrusion on people’s privacy.

Both data extracts were de-identified by their source organisation by attaching a unique identifier created using a SHA-256 hash. This helped ensure privacy as no primary identifiers left each data custodian’s organisation. The hashing algorithm used name, date of birth and Medicare number (a unique identifier for patients under the publicly funded universal health care system in Australia). Consistent pre-processing of patient names was undertaken by the custodians of both the primary care and hospitalisation data to ensure a high match rate (e.g. removing punctuation, numbers and spaces from name fields). These datasets were then provided for the study and the SHA field was used to link the primary care and hospitalisation records.

Participants

The Victorian primary care patient database contained data for over 1.8 million patients. This was filtered to consider only patients who attended their primary care clinics at least once between 2011 and 2016, resulting in approximately 1.3 million patients. To ensure every patient in the resulting cohort had at least one year of history (“look back”) available, and that every patient visited their GP after the one-year prediction period, and was alive for the entire prediction period, two key dates were calculated for each patient. “One year ahead” is one year from the date of a patient’s first visit (minimum value 1 Jan 2012). “One year back” is one year before the date of a patient’s last visit (maximum value 31 December 2015). The patient cohort was subsequently filtered to include only patients where the “One year back” date occurred on or after the “One year ahead” date.

The cohort was further filtered to exclude patients below the age of 18 and above the age of 106 years, patients without Victorian postcodes, patients with gender other than “Male” or “Female”, and physiological observations and pathology records with inconsistent or incorrect units and values. Figure 1 presents the cohort population at each stage of the cohort selection process. The resulting cohort, hereafter referred to as the Primary Cohort, comprised 393,229 patients, of whom 7.2% had hospitalisations of interest. A chronic-only subset of the Primary Cohort was also created for sub-analysis and is described in following sections.

To avoid a bias arising from consistently choosing the earliest or latest possible date as a patient’s “Prediction Date”, a random date was chosen between the “One Year Ahead” and “One Year Back” dates. While ensuring that sufficient history and prediction windows were available, this suited the perceived consumer use of the model, where GP practices may want to predict on their cohort at an “ad hoc” time period, not linked to a particular prediction date. Figure 2 presents the methodology chosen for selecting a random prediction date for each patient.

The public hospitalisation extract revealed several patients that appeared to be different people in the primary care dataset but were determined to be the same person. Where this was identified, their primary care records were merged to improve data quality.

Emergency or Potentially Preventable Hospitalisations

The outcome variable of interest for study was whether a patient has an Emergency or Potentially Preventable Hospitalisation within one year of the prediction date. Potentially Preventable Hospitalisations were included as they are a defined indicator for ambulatory care sensitive conditions10 and the case for targeting these through general practice integrated care interventions is well justified11,12. An Admission Type flag was employed to identify Emergency Hospitalisation. International Classification of Diseases (ICD) and Procedure Code information was employed to identify Potentially Preventable Hospitalisations using the guideline provided by the Australian Institute for Health and Welfare13. A binary (1/0) “One Year Emergency or Potentially Preventable Hospitalisation” variable was created after calculating, for each patient, the days from their prediction date to their next Emergency or Potentially Preventable Hospitalisation, whichever came first. The data comprised an extract from the inpatient (VAED) and emergency department (VEMD) datasets linked to the death registry (Registry of Births, Deaths and Marriages, Victoria). The only records of interest were related to inpatient hospitalisation.

Predictors

The list of all predictors used to develop the model is shown in Table 1 with further description below.

Demographics

Age at cohort entry was expected to be a key predictor of risk and was implemented as a continuous variable instead of categorical to prevent loss of information within categories. To capture a non-linear relationship of risk on age, the square root of age (√age), age2 (age × age) and age3 (age × age × age) predictors were also created. This was preferred over Cubic Splines, given that these were more practical to implement within the production environment of software systems used in primary care practices participating in the Health Care Homes trial.

For the Australian context, ethnicity was defined as Indigenous Australians, non-Indigenous Australians, or Unknown (where no ethnicity data was available).

Postcode was used to match with the Australia Bureau of Statistics’ SEIFA Index of Relative Socio-economic Advantage and Disadvantage (IRSAD)14, and with the Australian Department of Health’s Modified Monash Model (MMM)15. The SEIFA index summarises information about the economic and social conditions of people and households within an area, including both relative advantage and disadvantage measures, and was employed as a categorical decile value (i.e. grouped into 10 bands with lower values signifying higher levels of disadvantage, plus “Unknown” for missing postcodes). MMM is a recently developed geographical classification system designed to better address the maldistribution of medical services across Australia, taking values ranging from 1 to 7, with increasing values representing higher levels of rurality, plus “Unknown” for missing postcodes.

Patient Observations and History

Recorded physiological observations were included following data preparation to remove non-numeric characters, conversion of units to standardised forms, and removal of values and units identified as being “incorrectly recorded”. Since all physiological observations are not repeated at each GP visit, the system looked back through all available physiological observations (going back to 2007) to retrieve the most recent available reading.

Smoking status was included as a four-category predictor. Four variations of alcohol consumption were included. In the Australian context, the main software systems collect the number of alcoholic drinks per day consumed by the patient on average, the number of days per week the patient consumed alcohol on average, both of these, or a related measure. To increase the usability of the model in the Australian context, we defined alcohol consumption as a binary condition indicated by any consumption of alcohol per day and alcohol days per week data. We note that “volume” information is lost when converting to a binary category which is commented on in the Discussion. An additional category was included to indicate no alcohol consumption data was available.

Height, weight and body mass index (BMI) were also employed for prediction. BMI was calculated where it was not specifically recorded, but where height and weight were available. Height and weight were used as continuous predictors while BMI was employed as a categorical predictor, either with four categories or with six categories. In the 4-category predictor, the categories corresponded to: underweight/normal, overweight, obese, and “not recorded”. In the 6-category predictor, the obese category was subdivided into moderately obese, severely obese, and very severely obese.

Medication Variables

The usage of certain medications as indicated in the GP prescription data was included in the modelling using 0/1 binary indicator variables for each of six medication groups (statins, anticoagulants, antidepressants, antipsychotics, anti-inflammatories, and steroids). Groupings of medications into each medication group was completed in consultation with clinical and pharmacy experts from the Australian Government Department of Health. Data was searched back to 2007 (see Supplementary Table S1). The number of prescribed medication groups was also included in modelling either as the total number of prescribed groups (range 0–6) or as a categorical variable with five categories.

Diagnosis Variables

Grouping of diagnosis conditions/diseases of interest into “diagnosis families” was completed in consultation with the Australian Government Department of Health. In total, 35 diagnosis families were identified as relevant for the model (e.g., type 1 and 2 diabetes, osteoporosis, rheumatoid arthritis, cancer, learning difficulties, falls etc.). Data was searched back to 2007. Due to sparsity in the hospitalisation data, these “diagnosis families” were then further combined into 16 “disease groups”. The grouping was logical and clinically relevant (e.g., hyperlipidaemia, hypercholesterolaemia, and hypertriglyceridaemia were combined into a “bloodfats” category), and was informed by the incidence of each of the diagnosis conditions to ensure counts of the employed categorical variable were sufficient for analysis. A 0/1 binary indicator predictor for each group indicated the presence of at least one of the constituent diagnoses of interest in the GP records. The mapping of diagnosed conditions into diagnosis families and diagnosis groups is presented in Supplementary Table S2. These diagnosis group variables were the primary means to include the occurrence of relevant diagnosed conditions/diseases in the models. Note that these predictors correspond to diagnosis, not occurrence; they provide no way to capture additional risk from undiagnosed conditions/diseases of interest. However, the morbidity risk groups from pathology test results do provide a limited way to do this. Comorbidity associated with each patient was also included in modelling as the total number of distinct diagnosis families, either as a continuous predictor (range 0–35) or as a categorical variable with nine categories.

A “Chronic-only” group of diagnosis families was also defined, being one of 33 chronic conditions (i.e. any of the 35, excluding learning difficulties and falls). A “Chronic-only subset” of the patient cohort was then defined as including only those patients who had at least one diagnosis belonging to one of the “Chronic-only” group of diagnosis families.

Pathology Variables

Twelve pathology results were included as predictors. Following removal of values and units identified as being “incorrectly recorded”, the most current value of the available pathology results was used. Normal ranges for pathology results were defined in consultation with the Australian Government Department of Health and used to calculate morbidity risk flags representing three categories of risk – Low, Medium and High (see Supplementary Table S3). While not strictly a pathology test result, blood pressure was handled in a similar way. Due to data sparsity, medium and high levels were each combined for Bilirubin, Creatinine and Triglycerides.

For a number of patients, no test results were available. If these were treated as missing values it would be very difficult to include those patients in models using test results as predictors. An alternative approach we employed was to add a category for “no test history” for each predictor which enabled inclusion of all patients. Note that the absence of a test result is not the same as a test result indicating low risk.

Methodology

A predictive model for the binary event of Emergency or Potentially Preventable hospitalisation within one year was developed. Several machine learning approaches were chosen for model development and validation. Logistic regression16 was chosen because it is an established method for prediction in binary problems and implementation of models in a range of production environments is straightforward using basic mathematical functions. For comparison, four alternative types of models were also considered (Naïve Bayes17, Artificial Neural Networks18, Random Forests19 and Generalised Boosting20). Random Forests, Artificial Neural Networks and Generalised Boosting algorithms were considered because of their established superiority in pattern recognition from large complex data. Naive Bayes was employed as it offers a different approach to model building. Multiple variants of the full model were explored. For example, BMI with four or six categories, alcohol consumption as drinks per day or non-drinker/drinker, and the number of diseases as a numeric or a categorical variable were each considered. For each model variant, the performance of the logistic regression model and the four machine learning approaches were assessed.

Models were compared and validated using the area under the receiver operating characteristic curve (AUC) with 10-fold cross-validation. Optimism was estimated using bootstrap validatation16 with 100 bootstrap samples. For additional results (e.g., receiver operating characteristic (ROC) curves, calibration curves), patients were randomly split into 70%/30% training and test portions.

The subgroup of patients with at least one chronic condition (Chronic-only subset) is an important subgroup of primary care patients, especially since many programs like Health Care Homes tend to focus on patients with established chronic conditions. Although models were developed using the entire cohort of patients, additional results for this subset of patients were obtained and are reported in the Supplementary material (see Supplementary Results for the Chronic-Only Subset). This approach was used instead of restricting all model development and validation to the chronic-only subset to take advantage of the information and large patient numbers amongst those without a diagnosed chronic condition.

Software

All analyses were performed using the R statistical computing environment. Logistic regression models used the “glm” command. Random forest models employed the “randomforest” package. Generalised boosting models used the “GBM” package. Naïve Bayes models employed the “klaR” package. Artificial Neural Networks models employed the “nnet” package.

Code Availabilty

The code used for model development and validation is available from the corresponding author on reasonable request for non-commercial purposes.

Results

Participants

Table 2 summarises the patient characteristics in the Primary Cohort all together, as split by the outcome variable, and as split randomly into 70% training (development) and 30% testing (validation) datasets. For categorical quantities, values used as the reference category in modelling are listed first.

Outcomes

Table 3 profiles hospitalisations for the Primary Cohort. The rate of emergency or potentially preventable hospitalisations within one year is 7.2%.

Model Performance

Table 4 presents a comparison of AUC results using 10-fold cross-validation for the best performing model (the “final” model), a simpler canonical model (age and number of diagnosis families) and a variant on the final model capturing the size of alcohol consumption using the 6-category variable for alcohol per day. (Additional results for the Chronic-only subset are provided in the Supplementary material). Two modelling approaches, Logistic Regression and Generalised boosting are presented as no other method had higher AUC values. Results from 10-fold AUC for these other models were never better than for the final model. The bootstrap estimate of optimism of the final logistic regression model is 0.0017.

As part of the model development process, deviance residuals for the logistic regression models were calculated. For seven variables (age, gender, ethnicity, 6-category BMI, smoking status, any alcohol consumption and SEIFA IRSAD), ANOVA was used to check whether there was a statistically significant relationship between the predictor and the deviance residuals, the presence of which would suggest the model was missing a key relationship. The checks were performed for the residuals of the entire training dataset and for just the Chronic-only subset of the training dataset. A highly significant relationship for gender was found in the Chronic-only subset (P = 0.002) but not in the entire training dataset (P = 0.168). Interaction terms between diagnosis group and gender were investigated by creating an interim model with all diagnosis group-gender interaction terms. Interaction terms in the model were retained for the final model if the P-value for the interaction term in the model output was below a generous threshold of 0.5. As a result, 9 interaction terms were included in the final model. Repeating the ANOVA checks for the final model, inclusion of these 9 terms largely removed the effect for gender (training dataset P = 0.043, Chronic-only P = 0.052). P-values for other variables were always larger than 0.03 for both datasets before and after adding the interaction terms.

Figure 3 shows the receiver operating characteristic (ROC) curve for the final model. Labelled points indicate the corresponding location of risk for indicated quantiles. The diagonal reference line is the line of no discrimination and corresponds to the performance of random guesses. Figure 4 shows a calibration curve for risk groups defined by deciles of predicted risk. For each group, the horizontal coordinate is the mean of predicted risk of hospitalisation and the vertical coordinate is the mean (with 95% confidence interval) for the observed proportion of hospitalisations. The dashed diagonal line indicates ideal performance where observed and predicted risk are equal.

Using a 70%/30% train/test data split, we checked the amount of bias in our unknown/not recorded categories by stratifying the test dataset into known and unknown parts for variables of interest (BMI, ethnicity, smoking, alcohol consumption, IRSAD advantage/disadvantage). Predictive performance on the 30% test dataset (AUC 0.665) is similar to that from 10-fold cross validation (AUC 0.663), so the results should be indicative. Calibration for the known parts are universally excellent. For the unknown parts, if there is bias it is small. Comparing mean predictions to observed proportions of hospitalisations, all results are within confidence intervals (BMI mean prediction 0.184, observed 0.190 (95% CI: 0.181–0.199; ethnicity mean prediction 0.212, observed 0.218 (95% CI: 0.202–0.234), smoking mean prediction 0.214, observed 0.233 (95% CI: 0.211–0.255); IRSAD unknown or AnyAlcohol not recorded mean prediction 0.232, observed 0.264 (95% CI: 0.228–0.300). For all of these results, the difference in mean with the neighbouring risk group decile greatly exceeds the difference between known/recorded and unknown/not recorded within the decile.

Model Specification and Computation

Table 5 shows the complete model specification for the final logistic model. The model coefficient for each predictor is shown in column 2. For categorical predictors, coefficients are for 0/1 binary indicator variables (dummy coding). With dummy coding, each category (except the reference category) gets its own binary variable, where 1 indicates the categorical predictor takes that value. Columns 3–4 shows the standard error and P-value for each coefficient. Column 5 shows the odds ratio with 95% confidence interval for each category compared to its reference level. Odds ratios are not provided for continuous predictors as their interpretation is made more complicated by quadratic and higher order terms.

As with other logistic regression models, predictions are calculated via the linear predictor. The linear predictor, y, is defined using the equation

$$y={\rm{Intercept}}+\,\sum _{i=1}^{N}({C}_{i}\times {V}_{i})$$

where for each i, Ci is the model coefficient and Vi is the corresponding variable value (e.g., 0 or 1 for binary indicator and dummy coding variables). The value of the intercept is also shown in the table. Note that for this model, $$N=95$$. The predicted probability, p, is calculated from the linear predictor according to the equation

$$p=\exp (y)/(1+\exp (y))$$

which necessarily satisfies $$0\le p\le 1$$.

Five test cases showing values of the predictor variables, the linear predictor and the predicted probability are presented in the Supplementary material to support computational reproducibility.

Discussion

This study presents a prediction model designed for Australian primary care practices to identify patients with chronic conditions in their patient population that are at high risk of hospitalisation over the next 12 months. Despite not having access to a rich linked data repository and high impact predictors such as previous hospitalisation history, the prediction model performs at similar levels to other state-of-the-art prediction models and demonstrates the efficacy of using this approach in similar settings.

Our model shows good performance by the area under the Receiver Operating Characteristic curve (AUC~0.66) and is well calibrated. As might be expected, obesity and smoking increase the risk of hospitalisation, while higher economic advantage lowers the risk of hospitalisation. Interestingly, alcohol consumption is also associated with decreased risk of hospitalisation. However, in our cohort many more drinkers have trivial to moderate consumption, so we would expect the category to be dominated by these kinds of drinkers. QAdmissions8 reported decreased risk of hospitalisations for this group, compared with increased risk for heavy drinkers, so our model appears consistent with other published results. We further refined our focus to the subset of patients with at least one diagnosed chronic condition, for which targeted healthcare is more important. On this subset, our model showed even better performance. Smoothed plots of predicted vs observed probabilities suggest that our model is doing well to capture increasing risk by age and number of diseases, two variables known to be important for this kind of risk stratification.

In current literature, variables related to prior health care utilisation have been shown to have a strong influence on the performance of such models. A review by Wallace et al.9 explored 18 risk prediction models that used routinely collected data for prediction. Seven of the 8 models in this group that were comparable to this study employed prior emergency admission in the final model, while the eighth created customised models that used a combination of data sources that included hospitalisation and community data. For example, the QAdmissions8 algorithm predicts emergency admission within 2 years (AUC 0.77–0.78) but, like most other algorithms, employs prior emergency hospitalisation information (sourced from both GP data and linked hospital episode data) as the most significant predictor.

Heterogeneity in the definitions of study populations and outcome measures limits the ability to directly compare our model’s performance with other work. Models such as the adapted PEONY21 and HARP22, and the models presented by Johnson et al.23 focus on older patients only and are therefore not directly comparable. In addition, while Johnson et al. present models based on “GP-like-data”, data for their study is drawn from a self-reported, English-only, voluntary survey and simulates conditions where data is accessed through primary care practice management systems. It is therefore likely to be more complete than real-world primary care data. Also, they use survey variables like ‘health condition count’ and ‘self-rated health’ to approximate the patient health history information that is commonly captured in GP patient management software. Our experience of data extracted directly from GP patient management systems suggests that this is far from representative of the true picture. Haas et al.24 evaluated the performance of six algorithms in predicting outcome measures including inpatient hospitalisation (AUC 0.67–073) and Emergency Department visits (AUC 0.58 to 0.67) within 12 months. While somewhat similar, both of these outcome measures are different enough from the one employed for this study to make direct comparison difficult. On the other hand, efforts such as the Gold Coast Integrated Care in Australia25 employ a purposely designed Risk of Hospitalisation scoring mechanism for which no validation or performance information is published, making comparison impossible. In general, when comparing similar outcome measures, our predictive risk model performed at a level that was similar to, if not better than, other models that did not include previous hospitalisation as a predictor variable.

As with other studies of this type, there are limitations to this work. We suspect that model performance was constrained by primary care data quality issues, especially data correctness, inconsistencies, and completeness. Significant effort was devoted to data preparation to ensure quality issues (e.g., data entry issues, unit inconsistencies) were appropriately handled, derived predictors were calculated in a consistent manner, and steps were taken to reduce missingness where possible (e.g., calculating BMI where height and weight could be found in the data). Nevertheless, missing data was a feature in several key variables (BMI, ethnicity, alcohol consumption, smoking history, index of advantage/disadvantage). The expected use for our model is in GP clinics using computerised systems similar to those from which the sample data was obtained. Thus, we anticipate similar patterns of missingness in the production environment.

Statistical approaches to missing data include complete case analysis and multiple imputation. We avoided complete case analysis due to concerns about introducing bias. Multiple imputation is a common approach but requires an assumption that the data is “missing at random” and development of additional models to impute missing data, possibly using additional covariates. We deliberately used a different approach, creating an unknown or “not recorded” category for each of these categorical predictors. We believe that in many cases the missingness reflects the nature of the missing value and possibly the reason for the GP visit. In short, we believe the data is missing not at random and has predictive power of its own. An advantage of our approach is that it easily extends to a production environment in which it is necessary to make predictions with data that similarly exhibits missingness. Unfortunately, the approach can also introduce bias. In our estimates that bias is small, but it means care should be exercised concerning the interpretation of model coefficients and standard errors for those predictors in particular. What remains clear is that complete data is preferable to missing data for these kinds of models, and we hope that data completeness in these key predictors increases.

Due to limitations in recording processes and available data fields, clinically relevant semantics of medication and diagnosis variables were not considered. For example, mode of delivery and dosage were not incorporated into medication predictors, diagnosis variables did not distinguish between treated and untreated conditions/diseases, and medications were not ascribed to particular conditions/diseases. We leave the study of such predictors for future work. Also, our choice to include potentially preventable hospitalisations in the response variable was guided by studies11,12 recommending the targeting of potentially preventable hospitalisations through integrated care and other reform measures. We recognise that not all such hospitalisations are actually preventable in the short-term through management in primary care10. In practice, the addition of clinical assessment of patients with high probability of hospitalisation would allow further tailoring of primary care appropriate to patient needs.

Our prediction model provides a validated tool for risk stratification in GP patient populations and is deployed within a software suite used in the Health Care Homes trial. It is implemented in a number of general practice clinics (up to 200) within a web application which sources data from five different general practice information systems and uses the model to identify patients at high risk of hospitalisation over the following 12 months. Combined with other elements of the Health Care Homes trial, which include a clinical assessment to assign identified patients to tiers and a patient-centric integrated care program, it is hoped that this model will help deliver significant improvements to patient outcomes while reducing the burden of chronic conditions on the Australian health system.

Data Availability

The datasets analysed in the current study are not publicly available. Due to reasonable privacy and security concerns and requirements imposed by the ethics approval process, they are not redistributable to researchers other than those engaged in the ethics committee approved research protocol.

References

1. World Health Organisation. Non communicable diseases. World Health Organization Available at: http://www.who.int/news-room/fact-sheets/detail/noncommunicable-diseases. (Accessed: 13th May 2018).

2. Loewenson, R. & Simpson, S. Strengthening Integrated Care Through Population-Focused Primary Care Services: International Experiences Outside the United States. Annu. Rev. Public Health 38, 413–429 (2017).

3. Crocker, R. L. et al. Integrative medicine primary care: assessing the practice model through patients’ experiences. BMC Complement. Altern. Med. 17, (2017).

4. Desmedt, M. et al. Economic Impact of Integrated Care Models for Patients with Chronic Diseases: A Systematic Review. Value Health 19, 892–902 (2016).

5. Bauer, U. E., Briss, P. A., Goodman, R. A. & Bowman, B. A. Prevention of chronic disease in the 21st century: elimination of the leading preventable causes of premature death and disability in the USA. The Lancet 384, 45–52 (2014).

6. Mahomed, O. H. & Asmall, S. Development and implementation of an integrated chronic disease model in South Africa: lessons in the management of change through improving the quality of clinical practice. Int. J. Integr. Care 15, e038 (2015).

7. de Bruin, S. R., Heijink, R., Lemmens, L. C., Struijs, J. N. & Baan, C. A. Impact of disease management programs on healthcare expenditures for patients with diabetes, depression, heart failure or chronic obstructive pulmonary disease: A systematic review of the literature. Health Policy 101, 105–121 (2011).

8. Predicting risk of emergency admission to hospital using primary care data: derivation and validation of QAdmissions score | BMJ Open. Available at: http://bmjopen.bmj.com/content/3/8/e003482. (Accessed: 14th May 2018).

9. Wallace, E. et al. Risk prediction models to predict emergency hospital admission in community-dwelling adults: a systematic review. Med. Care 52, 751–765 (2014).

10. Falster, M. O. & Jorm, L. R. A guide to the potentially preventable hospitalisations indicator in Australia. (Centre for Big Data Research in Health, University of New South Wales in consultation with Australian Commission on Safety and Quality in Health Care and Australian Institute of Health and Welfare, 2017).

11. Australian Commission on Safety and Quality in Health Care. The second Australian atlas of healthcare variation. (Sydney, NSW: Australian Commission on Safety and Quality in Health Care, 2017).

12. Department of Health. Primary Health Care Advisory Group Final Report: Better Outcomes for People with Chronic and Complex Health Conditions. (2016).

13. Australian Institute of Health and Welfare. National Healthcare Agreement: PI 18-Selected potentially preventable hospitalisations, 2015. Available at: http://meteor.aihw.gov.au/content/index.phtml/itemId/559032. (Accessed: 14th May 2018).

14. Australian Bureau of Statistics. Socio-Economic Indexes for Areas (SEIFA) - Technical Paper. (2011).

15. Mason, J. Review of Australian Government Health Workforce Programs. (2013).

16. Regression Modeling Strategies - With Applications to Linear Models, Logistic Regression, and Survival Analysis | Frank Harrell | Springer. Available at: https://www.springer.com/gp/book/9781441929181. (Accessed: 16th May 2018).

17. Rish, I. An Empirical Study of the Naive Bayes Classifier. in Proceedings of the IJCAI Workshop on Empirical Methods in AI (2001).

18. Hinton, G. E., Osindero, S. & Teh, Y.-W. A fast learning algorithm for deep belief nets. Neural Comput. 18, 1527–1554 (2006).

19. Breiman, L. Random Forests. Mach. Learn. 45, 5–32 (2001).

20. Friedman, J. H. Stochastic gradient boosting. Comput. Stat. Data Anal. 38, 367–378 (2002).

21. Tomlin, A. M., Lloyd, H. S. & Tilyard, M. W. Risk stratification of New Zealand general practice patients for emergency admissions in the next year: adapting the PEONY model for use in New Zealand. J. Prim. Health Care 8, 227–237 (2016).

22. Canadian Institute for Health Information (CIHI). Early Identification of People at-risk of Hospitalization: Hospital Admission Risk Prediction (HARP). (Health Quality, Ontario, 2013).

23. Johnson, T. L. et al. Predictive risk modelling under different data access scenarios: who is identified as high risk and for how long? BMJ Open 8, (2018).

24. Haas, L. R. et al. Risk-stratification methods for identifying patients for care coordination. Am. J. Manag. Care 19, 725–732 (2013).

25. Scuffham, P. A., Mihala, G., Ward, L., McMurray, A. & Connor, M. Evaluation of the Gold Coast Integrated Care for patients with chronic disease or high risk of hospitalisation through a non-randomised controlled clinical trial: a pilot study protocol. BMJ Open 7, e016776 (2017).

Acknowledgements

The authors would like to thank Roland Balodis, Imaina Widagdo and Dr. Tony Hobbs of the Australian Government Department of Health for useful suggestions on translating GP data into predictor variables. The development and validation of the algorithm was funded through the award of a competitive tender by the Australian Government Department of Health. While clinical context and support was provided, the funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. The manuscript was reviewed by the funding organisation to ensure confidential information was not released.

Author information

Authors

Contributions

Study concept and design: S.K., J.B., R.J., M.H., M.G. Data preparation and analysis: D.A.R., S.K., Y.X., M.H. Initial manuscript preparation: D.A.R., S.K., J.B., R.J. Manuscript review and editing: S.K., D.A.R., J.B., M.H., R.J., M.G.

Corresponding author

Correspondence to Sankalp Khanna.

Ethics declarations

Competing Interests

The authors declare no competing interests.

Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Reprints and Permissions

Khanna, S., Rolls, D.A., Boyle, J. et al. A risk stratification tool for hospitalisation in Australia using primary care data. Sci Rep 9, 5011 (2019). https://doi.org/10.1038/s41598-019-41383-y

• Accepted:

• Published:

• DOI: https://doi.org/10.1038/s41598-019-41383-y

• Exploring the ability of self-report measures to identify risk of high treatment burden in chronic disease patients: a cross-sectional study

• Ruth Hardman
• Stephen Begg
• Evelien Spelten

BMC Public Health (2022)

• Risk Stratification in Primary Care: Value-Based Contributions of Provider Adjudication

• Brian C. Ricci
• Jonathan Sachs
• David A. Dorr

Journal of General Internal Medicine (2022)


text:
Products
Rewards
from HOLOOLY

We are determined to provide the latest solutions related to all subjects FREE of charge!

Enjoy Limited offers, deals & Discounts by signing up to Holooly Rewards Program

HOLOOLY

HOLOOLY
TABLES

All the data tables that you may search for.

HOLOOLY
ARABIA

For Arabic Users, find a teacher/tutor in your City or country in the Middle East.

HOLOOLY
TEXTBOOKS

Find the Source, Textbook, Solution Manual that you are looking for in 1 click.

HOLOOLY
HELP DESK

Need Help? We got you covered.

## Q. 16.9

The inverse of the Hubble parameter has the units of time. Use the value of the Hubble constant 71 (km/s)/Mpc to find the upper limit t for the age of the universe.

Strategy We will need to understand the relationship between the various astronomical parameters to do this. We begin with the units given by $H_{0}^{-1}$ and convert the result to years.

## Verified Solution

We begin with

\begin{aligned}\tau &=\frac{1}{H_{0}}=\frac{1}{71 \frac{ km }{ s \cdot Mpc }} \\&=\frac{1 s \cdot Mpc }{71 km } \frac{1 y }{3.16 \times 10^{7} s } \frac{10^{6} pc }{ Mpc } \frac{3.26 ly }{1 pc } \frac{9.46 \times 10^{12} km }{ ly } \\&=13.7 \times 10^{9} y\end{aligned}

The upper limit of the age of the universe is 13.7 Gly.