We want to study a physical system which evolves towards equilibrium, from som given initial conditions. Recall the simple example of particles in a box. At an initial time \( t_0 \) all particles are in the left half of the box. Thereafter they are allowed to diffuse into the two halves of the box.
Use Markov chain Monte Carlo
The Markov process is used repeatedly in Monte Carlo simulations in order to generate new random states.
The reason for choosing a Markov process is that when it is run for a long enough time starting with a random state, we will eventually reach the most likely state of the system.
In thermodynamics, this means that after a certain number of Markov processes we reach an equilibrium distribution.
This mimicks the way a real system reaches its most likely state at a given temperature of the surroundings.
To reach this distribution, the Markov process needs to obey two important conditions, that of ergodicity and detailed balance. These conditions impose then constraints on our algorithms for accepting or rejecting new random states.
The Metropolis algorithm discussed here abides to both these constraints.
The Metropolis algorithm is widely used in Monte Carlo simulations and the understanding of it rests within the interpretation of random walks and Markov processes.
In a random walk one defines a mathematical entity called a walker, whose attributes completely define the state of the system in question.
The state of the system can refer to any physical quantities, from the vibrational state of a molecule specified by a set of quantum numbers, to the brands of coffee in your favourite supermarket.
The walker moves in an appropriate state space by a combination of deterministic and random displacements from its previous position.
This sequence of steps forms a chain.
$$
\frac{\partial V}{\partial t}+\frac{1}{2}\sigma^{2}S^{2}\frac{\partial^{2} V}{\partial S^{2}}+rS\frac{\partial V}{\partial S}-rV=0.
$$
The Black and Scholes equation is a partial differential equation, which describes the price
of the option over time. It is a diffusion equation with a random term.
The list of applications is endless
The obvious case is that of a random walker on a one-, or two- or three-dimensional lattice (dubbed coordinate space hereafter).
Consider a system whose energy is defined by the orientation of single spins. Consider the state \( i \), with given energy \( E_i \) represented by the following \( N \) spins
$$
\begin{array}{cccccccccc}
\uparrow&\uparrow&\uparrow&\dots&\uparrow&\downarrow&\uparrow&\dots&\uparrow&\downarrow\\
1&2&3&\dots& k-1&k&k+1&\dots&N-1&N\end{array}
$$
We may be interested in the transition with one single spinflip to a new state \( j \) with energy \( E_j \)
$$
\begin{array}{cccccccccc}
\uparrow&\uparrow&\uparrow&\dots&\uparrow&\uparrow&\uparrow&\dots&\uparrow&\downarrow\\
1&2&3&\dots& k-1&k&k+1&\dots&N-1&N\end{array}
$$
This change from one microstate \( i \) (or spin configuration) to another microstate \( j \) is the
configuration space analogue to a random walk on a lattice. Instead of jumping from
one place to another in space, we 'jump' from one microstate to another.
A Markov process allows in principle for a microscopic description of Brownian motion. As with the random walk studied in the previous section, we consider a particle which moves along the \( x \)-axis in the form of a series of jumps with step length \( \Delta x = l \). Time and space are discretized and the subsequent moves are statistically independent, i.e., the new move depends only on the previous step and not on the results from earlier trials. We start at a position \( x=jl=j\Delta x \) and move to a new position \( x =i\Delta x \) during a step \( \Delta t=\epsilon \), where \( i\ge 0 \) and \( j\ge 0 \) are integers. The original probability distribution function (PDF) of the particles is given by \( w_i(t=0) \) where \( i \) refers to a specific position on the grid in
For the Markov process we have a transition probability from a position \( x=jl \) to a position \( x=il \) given by
$$
\begin{equation*}
W_{ij}(\epsilon)=W(il-jl,\epsilon)=\left\{\begin{array}{cc}\frac{1}{2} & |i-j| = 1\\
0 & \mathrm{else} \end{array} \right. ,
\end{equation*}
$$
where \( W_{ij} \) is normally called
the transition probability and we can represent it, see below,
as a matrix.
Here we have specialized to a case where the transition probability is known.
Our new PDF \( w_i(t=\epsilon) \) is now related to the PDF at \( t=0 \) through the relation
$$
\begin{equation*}
w_i(t=\epsilon) =\sum_{j} W(j\rightarrow i)w_j(t=0).
\end{equation*}
$$
This equation represents the discretized time-development of an original PDF with equal probability of jumping left or right.
Since both \( W \) and \( w \) represent probabilities, they have to be normalized, i.e., we require that at each time step we have
$$
\begin{equation*}
\sum_i w_i(t) = 1,
\end{equation*}
$$
and
$$
\begin{equation*}
\sum_j W(j\rightarrow i) = 1,
\end{equation*}
$$
which applies for all \( j \)-values.
The further constraints are
\( 0 \le W_{ij} \le 1 \) and \( 0 \le w_{j} \le 1 \).
Note that the probability for remaining at the same place is in general
not necessarily equal zero.
The time development of our initial PDF can now be represented through the action of the transition probability matrix applied \( n \) times. At a time \( t_n=n\epsilon \) our initial distribution has developed into
$$
\begin{equation*}
w_i(t_n) = \sum_jW_{ij}(t_n)w_j(0),
\end{equation*}
$$
and defining
$$
\begin{equation*}
W(il-jl,n\epsilon)=(W^n(\epsilon))_{ij}
\end{equation*}
$$
we obtain
$$
\begin{equation*}
w_i(n\epsilon) = \sum_j(W^n(\epsilon))_{ij}w_j(0),
\end{equation*}
$$
or in matrix form
$$
\begin{equation} \tag{1}
\hat{w}(n\epsilon) = \hat{W}^n(\epsilon)\hat{w}(0).
\end{equation}
$$
The following simple example may help in understanding the meaning of the transition matrix \( \hat{W} \) and the vector \( \hat{w} \). Consider the \( 4\times 4 \) matrix \( \hat{W} \)
$$
\begin{equation*}
\hat{W} = \left(\begin{array}{cccc} 1/4 & 1/9 & 3/8 & 1/3 \\
2/4 & 2/9 & 0 & 1/3\\
0 & 1/9 & 3/8 & 0\\
1/4 & 5/9& 2/8 & 1/3 \end{array} \right),
\end{equation*}
$$
and we choose our initial state as
$$
\begin{equation*}
\hat{w}(t=0)= \left(\begin{array}{c} 1\\
0\\
0 \\
0 \end{array} \right).
\end{equation*}
$$
We note that both the vector and the matrix are properly normalized. Summing the vector elements gives one and summing over columns for the matrix results also in one. Furthermore, the largest eigenvalue is one. We act then on \( \hat{w} \) with \( \hat{W} \). The first iteration is
$$
\begin{equation*}
\hat{w}(t=\epsilon) = \hat{W}\hat{w}(t=0),
\end{equation*}
$$
resulting in
$$
\begin{equation*}
\hat{w}(t=\epsilon)= \left(\begin{array}{c} 1/4\\
1/2 \\
0 \\
1/4 \end{array} \right).
\end{equation*}
$$
The next iteration results in
$$
\begin{equation*}
\hat{w}(t=2\epsilon) = \hat{W}\hat{w}(t=\epsilon),
\end{equation*}
$$
resulting in
$$
\begin{equation*}
\hat{w}(t=2\epsilon)= \left(\begin{array}{c} 0.201389\\
0.319444 \\
0.055556 \\
0.423611 \end{array} \right).
\end{equation*}
$$
Note that the vector \( \hat{w} \) is always normalized to \( 1 \).
We find the steady state of the system by solving the set of equations
$$
\begin{equation*}
w(t=\infty) = Ww(t=\infty),
\end{equation*}
$$
which is an eigenvalue problem with eigenvalue equal to one!
This set of equations reads
$$
\begin{align}
W_{11}w_1(t=\infty) +W_{12}w_2(t=\infty) +W_{13}w_3(t=\infty)+ W_{14}w_4(t=\infty)=&w_1(t=\infty) \nonumber \\
W_{21}w_1(t=\infty) + W_{22}w_2(t=\infty) + W_{23}w_3(t=\infty)+ W_{24}w_4(t=\infty)=&w_2(t=\infty) \nonumber \\
W_{31}w_1(t=\infty) + W_{32}w_2(t=\infty) + W_{33}w_3(t=\infty)+ W_{34}w_4(t=\infty)=&w_3(t=\infty) \nonumber \\
W_{41}w_1(t=\infty) + W_{42}w_2(t=\infty) + W_{43}w_3(t=\infty)+ W_{44}w_4(t=\infty)=&w_4(t=\infty) \nonumber \\
\tag{2}
\end{align}
$$
with the constraint that
$$
\begin{equation*}
\sum_i w_i(t=\infty) = 1,
\end{equation*}
$$
yielding as solution
$$
\begin{equation*}
\hat{w}(t=\infty)= \left(\begin{array}{c}0.244318 \\
0.319602 \\ 0.056818 \\ 0.379261 \end{array} \right).
\end{equation*}
$$
The table here demonstrates the convergence as a function of the number of iterations or time steps. After twelve iterations we have reached the exact value with six leading digits.
| Iteration | \( w_1 \) | \( w_2 \) | \( w_3 \) | \( w_4 \) |
|---|---|---|---|---|
| 0 | 1.000000 | 0.000000 | 0.000000 | 0.000000 |
| 1 | 0.250000 | 0.500000 | 0.000000 | 0.250000 |
| 2 | 0.201389 | 0.319444 | 0.055556 | 0.423611 |
| 3 | 0.247878 | 0.312886 | 0.056327 | 0.382909 |
| 4 | 0.245494 | 0.321106 | 0.055888 | 0.377513 |
| 5 | 0.243847 | 0.319941 | 0.056636 | 0.379575 |
| 6 | 0.244274 | 0.319547 | 0.056788 | 0.379391 |
| 7 | 0.244333 | 0.319611 | 0.056801 | 0.379255 |
| 8 | 0.244314 | 0.319610 | 0.056813 | 0.379264 |
| 9 | 0.244317 | 0.319603 | 0.056817 | 0.379264 |
| 10 | 0.244318 | 0.319602 | 0.056818 | 0.379262 |
| 11 | 0.244318 | 0.319602 | 0.056818 | 0.379261 |
| 12 | 0.244318 | 0.319602 | 0.056818 | 0.379261 |
| \( \hat{w}(t=\infty) \) | 0.244318 | 0.319602 | 0.056818 | 0.379261 |
We have after \( t \)-steps
$$
\begin{equation*}
\hat{w}(t) = \hat{W}^t\hat{w}(0),
\end{equation*}
$$
with \( \hat{w}(0) \) the distribution at \( t=0 \) and \( \hat{W} \) representing the
transition probability matrix.
We can always expand \( \hat{w}(0) \) in terms of the right eigenvectors \( \hat{v} \) of \( \hat{W} \) as
$$
\begin{equation*}
\hat{w}(0) = \sum_i\alpha_i\hat{v}_i,
\end{equation*}
$$
resulting in
$$
\begin{equation*}
\hat{w}(t) = \hat{W}^t\hat{w}(0)=\hat{W}^t\sum_i\alpha_i\hat{v}_i=
\sum_i\lambda_i^t\alpha_i\hat{v}_i,
\end{equation*}
$$
with \( \lambda_i \) the \( i^{\mathrm{th}} \) eigenvalue corresponding to
the eigenvector \( \hat{v}_i \).
If we assume that \( \lambda_0 \) is the largest eigenvector we see that in the limit \( t\rightarrow \infty \), \( \hat{w}(t) \) becomes proportional to the corresponding eigenvector \( \hat{v}_0 \). This is our steady state or final distribution.
#include <iostream>
#include "armadillo"
using namespace arma;
using namespace std;
int main()
{
int dim = 5;
mat W = zeros<mat>(dim,dim);
vec wold = zeros<mat>(dim);
vec wnew = zeros<mat>(dim);
vec eigenvector = zeros<mat>(dim);
// Initializing the first vector
wold(0) = 1.0;
// Setting up the stochastic matrix W
W(0,0) = 0.; W(0,1) = 0.; W(0,2) = 0.25; W(0,3) = 0.0; W(0,4) = 0.;
W(1,0) = 0.; W(1,1) = 0.; W(1,2) = 0.25; W(1,3) = 0.; W(1,4) = 0.0;
W(2,0) = 0.5; W(2,1) = 1.0; W(2,2) = 0.; W(2,3) = 0.5; W(2,4) = 0.;
W(3,0) = 0.0; W(3,1) = 0.; W(3,2) = 0.25; W(3,3) = 0.; W(3,4) = 0.;
W(4,0) = 0.5; W(4,1) = 0.; W(4,2) = 0.25; W(4,3) = 0.5; W(4,4) = 1.0;
double eps = 1.0E-10;
W.print("W =");
double difference = norm(wold-wnew, 2);
int count = 0;
do{
// Multiplying the old vector with the transition probability
count += 1;
wnew = W*wold;
difference = norm(wold-wnew, 2);
wold = wnew;
cout << "Iteration number = " << count << endl;
wnew.print("New vector =");
} while(difference > eps);
// Getting the eigenvectors and eigenvalues of the stochastic matrix
cx_vec eigval;
eig_gen(eigval, W);
eigval.print("Eigenvalues=");
return 0;
}
Let us recapitulate some of our results about Markov chains and random walks.
one time-step from \( t=0 \) is given by
$$
\begin{equation*}
w_i(t=\epsilon) = W(j\rightarrow i)w_j(t=0).
\end{equation*}
$$
This equation represents the discretized time-development of an original PDF. We can rewrite this as a
$$
\begin{equation*}
w_i(t=\epsilon) = W_{ij}w_j(t=0).
\end{equation*}
$$
with the transition matrix \( W \) for a random walk given by
$$
\begin{equation*}
W_{ij}(\epsilon)=W(il-jl,\epsilon)=\left\{\begin{array}{cc}\frac{1}{2} & |i-j| = 1\\
0 & \mathrm{else} \end{array} \right.
\end{equation*}
$$
We call \( W_{ij} \) for the transition probability and we represent it as a matrix.
$$
\begin{equation*}
\sum_i w_i(t) = 1,
\end{equation*}
$$
and
$$
\begin{equation*}
\sum_j W(j\rightarrow i) = 1.
\end{equation*}
$$
Here we have written the previous matrix \( W_{ij}=W(j\rightarrow i) \).
The further constraints are \( 0 \le W_{ij} \le 1 \) and \( 0 \le w_{j} \le 1 \).
$$
\begin{equation*}
w_i(t+1) = \sum_jW_{ij}w_j(t),
\end{equation*}
$$
or as vector-matrix relation
$$
\begin{equation*}
\hat{w}(t+1) = \hat{W\hat{w}}(t),
\end{equation*}
$$
and if we have that \( ||\hat{w}(t+1)-\hat{w}(t)||\rightarrow 0 \), we say that
we have reached the most likely state of the system, the so-called steady state or equilibrium state.
Another way of phrasing this is
$$
\begin{equation}
w(t=\infty) = Ww(t=\infty).
\tag{3}
\end{equation}
$$
The question then is how can we model anything under such a severe lack of knowledge? The Metropolis algorithm comes to our rescue here. Since \( W(j\rightarrow i) \) is unknown, we model it as the product of two probabilities, a probability for accepting the proposed move from the state \( j \) to the state \( j \), and a probability for making the transition to the state \( i \) being in the state \( j \). We label these probabilities \( A(j\rightarrow i) \) and \( T(j\rightarrow i) \), respectively. Our total transition probability is then
$$
\begin{equation*}
W(j\rightarrow i)=T(j\rightarrow i)A(j\rightarrow i).
\end{equation*}
$$
The algorithm can then be expressed as
We wish to derive the required properties of the probabilities \( T \) and \( A \) such that \( w_i^{(t\rightarrow \infty)} \rightarrow w_i \), starting from any distribution, will lead us to the correct distribution.
We can now derive the dynamical process towards equilibrium. To obtain this equation we note that after \( t \) time steps the probability for being in a state \( i \) is related to the probability of being in a state \( j \) and performing a transition to the new state together with the probability of actually being in the state \( i \) and making a move to any of the possible states \( j \) from the previous time step.
We can express this as, assuming that \( T \) and \( A \) are time-independent,
$$
\begin{equation*}
w_i(t+1) = \sum_j \left [
w_j(t)T_{j\rightarrow i} A_{j\rightarrow i}
+w_i(t)T_{i\rightarrow j}\left ( 1- A_{i\rightarrow j} \right)
\right ] \,.
\end{equation*}
$$
All probabilities are normalized, meaning that \( \sum_j T_{i\rightarrow j} = 1 \). Using the latter, we can rewrite the previous equation as
$$
\begin{equation*}
w_i(t+1) = w_i(t) +
\sum_j \left [
w_j(t)T_{j\rightarrow i} A_{j\rightarrow i}
-w_i(t)T_{i\rightarrow j}A_{i\rightarrow j}\right ] \,,
\end{equation*}
$$
which can be rewritten as
$$
\begin{equation*}
w_i(t+1)-w_i(t) = \sum_j \left [w_j(t)T_{j\rightarrow i} A_{j\rightarrow i}
-w_i(t)T_{i\rightarrow j}A_{i\rightarrow j}\right ] .
\end{equation*}
$$
The last equation is very similar to the so-called Master equation, which relates the temporal dependence of a PDF \( w_i(t) \) to various transition rates. The equation can be derived from the so-called Chapman-Einstein-Enskog-Kolmogorov equation. The equation is given as
$$
\begin{equation}
\tag{4}
\frac{d w_i(t)}{dt} = \sum_j\left[ W(j\rightarrow i)w_j-W(i\rightarrow j)w_i\right],
\end{equation}
$$
which simply states that the rate at which the systems moves from a state \( j \)
to a final state \( i \) (the first term on the right-hand side of the last equation) is balanced by the rate at which the system undergoes transitions from the state \( i \) to a state \( j \) (the second term). If we have reached the so-called steady state, then the temporal development is zero. This means that in equilibrium we have
$$
\begin{equation*}
\frac{d w_i(t)}{dt} = 0.
\end{equation*}
$$
In the limit \( t\rightarrow \infty \) we require that the two distributions \( w_i(t+1)=w_i \) and \( w_i(t)=w_i \) and we have
$$
\begin{equation*}
\sum_j w_jT_{j\rightarrow i} A_{j\rightarrow i}= \sum_j w_iT_{i\rightarrow j}A_{i\rightarrow j},
\end{equation*}
$$
which is the condition for balance when the most likely state (or steady state) has been reached.
We see also that the right-hand side can be rewritten as
$$
\begin{equation*}
\sum_j w_iT_{i\rightarrow j}A_{i\rightarrow j}= \sum_j w_iW_{i\rightarrow j},
\end{equation*}
$$
and using the property that \( \sum_j W_{i\rightarrow j}=1 \), we can rewrite our equation
as
$$
\begin{equation*}
w_i= \sum_j w_jT_{j\rightarrow i} A_{j\rightarrow i}= \sum_j w_j W_{j\rightarrow i},
\end{equation*}
$$
which is nothing but the standard equation for a Markov chain when the steady state has been reached.
However, the condition that the rates should equal each other is in general not sufficient to guarantee that we, after many simulations, generate the correct distribution. We may risk to end up with so-called cyclic solutions. To avoid this we therefore introduce an additional condition, namely that of detailed balance
$$
\begin{equation*} W(j\rightarrow i)w_j= W(i\rightarrow j)w_i. \end{equation*}
$$
These equations were derived by Lars Onsager when studying irreversible processes.
At equilibrium detailed balance gives thus
$$
\begin{equation*} \frac{W(j\rightarrow i)}{W(i\rightarrow j)}=\frac{w_i}{w_j}. \end{equation*}
$$
Rewriting the last equation in terms of our transition probabilities \( T \) and
acceptance probobalities \( A \) we obtain
$$
\begin{equation*}
w_j(t)T_{j\rightarrow i}A_{j\rightarrow i}= w_i(t)T_{i\rightarrow j}A_{i\rightarrow j}.
\end{equation*}
$$
Since we normally have an expression for the probability distribution functions \( w_i \), we can rewrite the last equation as
$$
\begin{equation*}
\frac{T_{j\rightarrow i}A_{j\rightarrow i}}{T_{i\rightarrow j}A_{i\rightarrow j}}= \frac{w_i}{w_j}.
\end{equation*}
$$
In statistical physics this condition ensures that it is e.g., the Boltzmann distribution which is generated when equilibrium is reached.
We introduce now the Boltzmann distribution
$$
\begin{equation*}
w_i= \frac{\exp{(-\beta(E_i))}}{Z},
\end{equation*}
$$
which states that the probability of finding the system in a state \( i \) with energy \( E_i \)
at an inverse temperature \( \beta = 1/k_BT \) is \( w_i\propto \exp{(-\beta(E_i))} \).
The denominator \( Z \) is a normalization constant which ensures that the sum of all
probabilities is normalized to one. It is defined as the sum of probabilities over all microstates
\( j \) of the system
$$
\begin{equation*}
Z=\sum_j \exp{(-\beta(E_i))}.
\end{equation*}
$$
From the partition function we can in principle generate all interesting quantities for a given system in equilibrium with its surroundings at a temperature \( T \).
With the probability distribution given by the Boltzmann distribution we are now in a position where we can generate expectation values for a given variable \( A \) through the definition
$$
\begin{equation*}
\langle A \rangle = \sum_jA_jw_j=
\frac{\sum_jA_j\exp{(-\beta(E_j)}}{Z}.
\end{equation*}
$$
In general, most systems have an infinity of microstates making thereby the computation
of \( Z \) practically impossible and
a brute force Monte Carlo calculation over a given number of randomly selected microstates
may therefore not yield those microstates which are important
at equilibrium.
To select the most important contributions we need to
use the condition for detailed balance. Since this is just given by the ratios of probabilities,
we never need to evaluate the partition function \( Z \).
For the Boltzmann distribution, detailed balance results in
$$
\begin{equation*} \frac{w_i}{w_j}= \exp{(-\beta(E_i-E_j))}. \end{equation*}
$$
Let us now specialize to a system whose energy is defined by the orientation of single spins. Consider the state \( i \), with given energy \( E_i \) represented by the following \( N \) spins
$$
\begin{equation*}
\begin{array}{cccccccccc}
\uparrow&\uparrow&\uparrow&\dots&\uparrow&\downarrow&\uparrow&\dots&\uparrow&\downarrow\\
1&2&3&\dots& k-1&k&k+1&\dots&N-1&N\end{array}
\end{equation*}
$$
We are interested in the transition with one single spinflip to a new state \( j \) with energy \( E_j \)
$$
\begin{equation*}
\begin{array}{cccccccccc}
\uparrow&\uparrow&\uparrow&\dots&\uparrow&\uparrow&\uparrow&\dots&\uparrow&\downarrow\\
1&2&3&\dots& k-1&k&k+1&\dots&N-1&N\end{array}
\end{equation*}
$$
This change from one microstate \( i \) (or spin configuration) to another microstate \( j \) is the
configuration space analogue to a random walk on a lattice. Instead of jumping from
one place to another in space, we 'jump' from one microstate to another.
However, the selection of states has to generate a final distribution which is the Boltzmann distribution. This is again the same we saw for a random walker, for the discrete case we had always a binomial distribution, whereas for the continuous case we had a normal distribution. The way we sample configurations should result, when equilibrium is established, in the Boltzmann distribution. Else, our algorithm for selecting microstates is wrong.
As stated above, we do in general not know the closed-form expression of the transition rate and we are free to model it as \( W(i\rightarrow j)=T(i\rightarrow j)A(i\rightarrow j) \). Our ratio between probabilities gives us
$$
\begin{equation*}
\frac{A_{j\rightarrow i}}{A_{i\rightarrow j}}= \frac{w_iT_{i\rightarrow j}}{w_jT_{j\rightarrow i}}.
\end{equation*}
$$
The simplest form of the Metropolis algorithm (sometimes called for brute force Metropolis) assumes that
the transition probability \( T(i\rightarrow j) \) is symmetric, implying that \( T(i\rightarrow j)=T(j\rightarrow i) \).
We obtain then (using the Boltzmann distribution)
$$
\begin{equation*}
\frac{A(j\rightarrow i)}{A(i\rightarrow j)}= \exp{(-\beta(E_i-E_j))} .
\end{equation*}
$$
We are in this case interested in a new state \( E_j \) whose energy is lower than
\( E_i \), viz., \( \Delta E = E_j-E_i \le 0 \). A simple test would then be to accept only those
microstates which lower the energy.
Suppose we have ten microstates with energy \( E_0 \le E_1 \le E_2 \le E_3 \le \dots \le E_9 \).
Our desired energy is \( E_0 \).
At a given temperature \( T \) we start our simulation by randomly choosing state \( E_9 \). Flipping spins we may then find a path from \( E_9\rightarrow E_8 \rightarrow E_7 \dots \rightarrow E_1 \rightarrow E_0 \). This would however lead to biased statistical averages since it would violate the ergodic hypothesis discussed in the previous section. This principle states that it should be possible for any Markov process to reach every possible state of the system from any starting point if the simulations is carried out for a long enough time.
Any state in a Boltzmann distribution has a probability different from zero and if such a state cannot be reached from a given starting point, then the system is not ergodic. This means that another possible path to \( E_0 \) could be \( E_9\rightarrow E_7 \rightarrow E_8 \dots \rightarrow E_9 \rightarrow E_5 \rightarrow E_0 \) and so forth. Even though such a path could have a negligible probability it is still a possibility, and if we simulate long enough it should be included in our computation of an expectation value.
Thus, we require that our algorithm should satisfy the principle of detailed balance and be ergodic. The problem with our ratio
$$
\begin{equation*}
\frac{A(j\rightarrow i)}{A(i\rightarrow j)}= \exp{(-\beta(E_i-E_j))},
\end{equation*}
$$
is that we do not know the acceptance probability. This equation only specifies the ratio of pairs of probabilities. Normally we want an algorithm which is as efficient as possible and maximizes the number of accepted moves.
Moreover, we know that the acceptance probability has \( 0 \) as its smallest value and \( 1 \) as its largest.
If we assume that the largest possible acceptance probability is \( 1 \), we adjust thereafter the other acceptance probability
to this constraint.
To understand this better, assume that we have two energies, \( E_i \) and \( E_j \), with \( E_i < E_j \). This means that the largest acceptance value must be \( A(j\rightarrow i) \) since we move to a state with lower energy. It follows from also from the fact that the probability \( w_i \) is larger than \( w_j \). The trick then is to fix this value to \( A(j\rightarrow i)=1 \). It means that the other acceptance probability has to be
$$
\begin{equation*}
A(i\rightarrow j)= \exp{(-\beta(E_j-E_i))}.
\end{equation*}
$$
One possible way to encode this equation reads
$$
\begin{equation*}
A(j\rightarrow i)=\left\{\begin{array}{cc}
\exp{(-\beta(E_i-E_j))} & E_i-E_j > 0 \\ 1 & else \end{array} \right.,
\end{equation*}
$$
implying that if we move to a state with a lower energy, we always accept
this move with acceptance probability \( A(j\rightarrow i)=1 \). If the energy is higher, we need to check
this acceptance probability with the ratio between the probabilities from our PDF. From a practical point of view,
the above ratio is compared with a random number.
If the ratio is smaller than a given random number we accept the move to a higher energy, else we stay in the same state.
Nothing hinders us obviously in choosing another acceptance ratio, like a weighting of the two energies via
$$
\begin{equation*}
A(j\rightarrow i)=\exp{(-\frac{1}{2}\beta(E_i-E_j))}.
\end{equation*}
$$
However, it is easy to see that such an acceptance ratio would result in
fewer accepted moves.
The Monte Carlo approach, combined with the theory for Markov chains can be summarized as follows: A Markov chain Monte Carlo method for the simulation of a distribution \( w \) is any method producing an ergodic Markov chain of events \( x \) whose stationary distribution is \( w \). The Metropolis algorithm can be phrased as
$$
\begin{equation*}
x^{(i+1)}= \left\{\begin{array}{cc} y_t & \mathrm{with\hspace{0.1cm}probability} = A(x^{(i)}\rightarrow y_t) \\ x^{(i)} & \mathrm{with \hspace{0.1cm}probability} = 1-A(x^{(i)}\rightarrow y_t)\end{array}\right .
\end{equation*}
$$
$$
\begin{equation*} A(x\rightarrow y)= \mathrm{min}\left\{\frac{w(y)T(x|y)}{w(x)T(y|x)},1\right\}.
\end{equation*}
$$
# Program to test the Metropolis algorithm with one particle at given temp in one dimension
import numpy as np
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import random
from math import sqrt, exp, log
# initialize the rng with a seed
random.seed()
# Hard coding of input parameters
MCcycles = 100000
Temperature = 2.0
beta = 1./Temperature
InitialVelocity = -2.0
CurrentVelocity = InitialVelocity
Energy = 0.5*InitialVelocity*InitialVelocity
VelocityRange = 10*sqrt(Temperature)
VelocityStep = 2*VelocityRange/10.
AverageEnergy = Energy
AverageEnergy2 = Energy*Energy
VelocityValues = np.zeros(MCcycles)
# The Monte Carlo sampling with Metropolis starts here
for i in range (1, MCcycles, 1):
TrialVelocity = CurrentVelocity + (2.0*random.random() - 1.0)*VelocityStep
EnergyChange = 0.5*(TrialVelocity*TrialVelocity -CurrentVelocity*CurrentVelocity);
if random.random() <= exp(-beta*EnergyChange):
CurrentVelocity = TrialVelocity
Energy += EnergyChange
VelocityValues[i] = CurrentVelocity
AverageEnergy += Energy
AverageEnergy2 += Energy*Energy
#Final averages
AverageEnergy = AverageEnergy/MCcycles
AverageEnergy2 = AverageEnergy2/MCcycles
Variance = AverageEnergy2 - AverageEnergy*AverageEnergy
print(AverageEnergy, Variance)
n, bins, patches = plt.hist(VelocityValues, 400, facecolor='green')
plt.xlabel('$v$')
plt.ylabel('Velocity distribution P(v)')
plt.title(r'Velocity histogram at $k_BT=2$')
plt.axis([-5, 5, 0, 600])
plt.grid(True)
plt.show()
Diffusion and the diffusion equation are central topics in both Physics and Mathematics, and their ranges of applicability span from stellar dynamics to the diffusion of particles governed by Schroedinger's equation. The latter is, for a free particle, nothing but the diffusion equation in complex time!
Let us consider the one-dimensional diffusion equation. We study a large ensemble of particles performing Brownian motion along the \( x \)-axis. There is no interaction between the particles.
We define \( w(x,t)dx \) as the probability of finding a given number of particles in an interval of length \( dx \) in \( x\in [x, x+dx] \) at a time \( t \). This quantity is our probability distribution function (PDF).
From experiment there are strong indications that the flux of particles \( j(x,t) \), viz., the number of particles passing \( x \) at a time \( t \) is proportional to the gradient of \( w(x,t) \). This proportionality is expressed mathematically through
$$
\begin{equation*}
j(x,t) = -D\frac{\partial w(x,t)}{\partial x},
\end{equation*}
$$
where \( D \) is the so-called diffusion constant, with dimensionality length$^2$ per time.
If the number of particles is conserved, we have the continuity equation
$$
\begin{equation*}
\frac{\partial j(x,t)}{\partial x} = -\frac{\partial w(x,t)}{\partial t},
\end{equation*}
$$
which leads to
$$
\begin{equation} \tag{5}
\frac{\partial w(x,t)}{\partial t} =
D\frac{\partial^2w(x,t)}{\partial x^2},
\end{equation}
$$
which is the diffusion equation in one dimension.
With the probability distribution function \( w(x,t)dx \) we can evaluate expectation values such as the mean distance
$$
\langle x(t)\rangle = \int_{-\infty}^{\infty}xw(x,t)dx,
$$
or
$$
\langle x^2(t)\rangle = \int_{-\infty}^{\infty}x^2w(x,t)dx,
$$
which allows for the computation of the variance
\( \sigma^2=\langle x^2(t)\rangle-\langle x(t)\rangle^2 \). Note well that
these expectation values are time-dependent.
In a similar way we can also define expectation values of functions \( f(x,t) \) as
$$
\begin{equation*}
\langle f(x,t)\rangle = \int_{-\infty}^{\infty}f(x,t)w(x,t)dx.
\end{equation*}
$$
The normalization condition
$$
\begin{equation*}
\int_{-\infty}^{\infty}w(x,t)dx=1
\end{equation*}
$$
imposes significant constraints on \( w(x,t) \).
We have
$$
\begin{equation*}
w(x=\pm \infty,t)=0 \hspace{1cm}
\frac{\partial^{n}w(x,t)}{\partial x^n}|_{x=\pm\infty} = 0,
\end{equation*}
$$
implying that when we study the time-derivative
\( \partial\langle x(t)\rangle/\partial t \), we obtain after integration by parts and using
Eq. (5)
$$
\frac{\partial \langle x\rangle}{\partial t} =
\int_{-\infty}^{\infty}x\frac{\partial w(x,t)}{\partial t}dx=
D\int_{-\infty}^{\infty}x\frac{\partial^2w(x,t)}{\partial x^2}dx,
$$
leading to
$$
\frac{\partial \langle x\rangle}{\partial t} =
Dx\frac{\partial w(x,t)}{\partial x}|_{x=\pm\infty}-
D\int_{-\infty}^{\infty}\frac{\partial w(x,t)}{\partial x}dx.
$$
The result is
$$
\frac{\partial \langle x\rangle}{\partial t} = 0.
$$
This means in turn that \( \langle x\rangle \) is independent of time.
If we choose the initial position \( x(t=0)=0 \),
the average displacement \( \langle x\rangle= 0 \).
If we link this discussion to a random walk in one dimension with equal probability
of jumping to the left or right and with an initial position \( x=0 \), then our probability
distribution remains centered around \( \langle x\rangle= 0 \) as function of time.
The variance is not necessarily 0. Consider first
$$
\frac{\partial \langle x^2\rangle}{\partial t} =
Dx^2\frac{\partial w(x,t)}{\partial x}|_{x=\pm\infty}-
2D\int_{-\infty}^{\infty}x\frac{\partial w(x,t)}{\partial x}dx,
$$
where we have performed an integration by parts as we did
for \( \frac{\partial \langle x\rangle}{\partial t} \).
Integration by parts results in
$$
\frac{\partial \langle x^2\rangle}{\partial t} =
-Dxw(x,t)|_{x=\pm\infty}+
2D\int_{-\infty}^{\infty}w(x,t)dx=2D,
$$
leading to
$$
\langle x^2\rangle = 2Dt,
$$
and the variance as
$$
\begin{equation} \tag{6}
\langle x^2\rangle-\langle x\rangle^2 = 2Dt.
\end{equation}
$$
The root mean square displacement after a time \( t \) is then
$$
\begin{equation*}
\sqrt{\langle x^2\rangle-\langle x\rangle^2} = \sqrt{2Dt}.
\end{equation*}
$$
This should be contrasted to the displacement of a free particle with initial velocity \( v_0 \). In that case the distance from the initial position after a time \( t \) is \( x(t) = vt \) whereas for a diffusion process the root mean square value is \( \sqrt{\langle x^2\rangle-\langle x\rangle^2} \propto \sqrt{t} \). Since diffusion is strongly linked with random walks, we could say that a random walker escapes much more slowly from the starting point than would a free particle.
$$
\begin{equation*}
w(x,t)dx = \frac{1}{\sqrt{4\pi Dt}}\exp{(-\frac{x^2}{4Dt})}dx.
\end{equation*}
$$
At a time $t=2$s the new variance is $\sigma^2=4D$s, implying that the root mean square value
is \( \sqrt{\langle x^2\rangle-\langle x\rangle^2} = 2\sqrt{D} \).
At a further time \( t=8 \) we have \( \sqrt{\langle x^2\rangle-\langle x\rangle^2} = 4\sqrt{D} \).
While time has elapsed by a factor of \( 4 \), the root mean square has only changed by a factor
of 2.
Consider
$$
\begin{equation*}
\langle x(n)\rangle = \sum_{i}^{n} \Delta x_i = 0 \hspace{1cm} \Delta x_i=\pm l,
\end{equation*}
$$
since we have an equal probability of jumping either to the left or to right.
The value of \( \langle x(n)^2\rangle \) is
$$
\begin{equation*}
\langle x(n)^2\rangle = \left(\sum_{i}^{n} \Delta x_i\right)\left(\sum_{j}^{n} \Delta x_j\right)=\sum_{i}^{n} \Delta x_i^2+
\sum_{i\ne j}^{n} \Delta x_i\Delta x_j=l^2n.
\end{equation*}
$$
For many enough steps the non-diagonal contribution is
$$
\begin{equation*}
\sum_{i\ne j}^{N} \Delta x_i\Delta x_j=0,
\end{equation*}
$$
since \( \Delta x_{i,j} = \pm l \).
The variance is then
$$
\begin{equation}
\langle x(n)^2\rangle - \langle x(n)\rangle^2 = l^2n.
\tag{7}
\end{equation}
$$
It is also rather straightforward to compute the variance for \( L\ne R \). The result is
$$
\begin{equation*}
\langle x(n)^2\rangle - \langle x(n)\rangle^2 = 4LRl^2n.
\end{equation*}
$$
In Eq. (7) the variable \( n \) represents the number of time steps. If we define \( n=t/\Delta t \), we can then couple the variance result from a random walk in one dimension with the variance from the diffusion equation of Eq. (6) by defining the diffusion constant as
$$
\begin{equation*}
D = \frac{l^2}{\Delta t}.
\end{equation*}
$$
The main program reads the name of the output file from screen and sets up the arrays containing the walker's position after a given number of steps. The corresponding program for a two-dimensional random walk (not listed in the main text) is found under programs/chapter12/program2.cpp
/*
1-dim random walk program.
A walker makes several trials steps with
a given number of walks per trial
*/
#include <iostream>
#include <fstream>
#include <iomanip>
#include "lib.h"
using namespace std;
// Function to read in data from screen, note call by reference
void initialise(int&, int&, double&) ;
// The Mc sampling for random walks
void mc_sampling(int, int, double, int *, int *);
// prints to screen the results of the calculations
void output(int, int, int *, int *);
int main()
{
int max_trials, number_walks;
double move_probability;
// Read in data
initialise(max_trials, number_walks, move_probability) ;
int *walk_cumulative = new int [number_walks+1];
int *walk2_cumulative = new int [number_walks+1];
for (int walks = 1; walks <= number_walks; walks++){
walk_cumulative[walks] = walk2_cumulative[walks] = 0;
} // end initialization of vectors
// Do the mc sampling
mc_sampling(max_trials, number_walks, move_probability,
walk_cumulative, walk2_cumulative);
// Print out results
output(max_trials, number_walks, walk_cumulative,
walk2_cumulative);
delete [] walk_cumulative; // free memory
delete [] walk2_cumulative;
return 0;
} // end main function
\end{lstlisting}
The input and output functions are
\begin{lstlisting}
void initialise(int& max_trials, int& number_walks, double& move_probability)
{
cout << "Number of Monte Carlo trials =";
cin >> max_trials;
cout << "Number of attempted walks=";
cin >> number_walks;
cout << "Move probability=";
cin >> move_probability;
} // end of function initialise
The algorithm tests the probability of moving to the left or to the right by generating a random number.
void mc_sampling(int max_trials, int number_walks,
double move_probability, int *walk_cumulative,
int *walk2_cumulative)
{
long idum;
idum=-1; // initialise random number generator
for (int trial=1; trial <= max_trials; trial++){
int position = 0;
for (int walks = 1; walks <= number_walks; walks++){
if (ran0(&idum) <= move_probability) {
position += 1;
}
else {
position -= 1;
}
walk_cumulative[walks] += position;
walk2_cumulative[walks] += position*position;
} // end of loop over walks
} // end of loop over trials
} // end mc_sampling function


--- layout: default ---
I'm a conversion-driven marketer... In other words, I help turn more prospects into leads, more leads into customers, and more customers into diehard fans of your company.
Copy is the biggest tool I use to make that happen.
But it's not the only tool...
My goal is to be a "full stack" marketer... Someone who can write effective copy, but also handle the other necessary factors for that copy to be successful.
I've written for several markets, but I specialize in health & finance.
If you're looking for someone like me, use the form below to tell me about your project.
You can also see my testimonials by clicking here
Cheers,
Wil


Netflix's The Confession Killer could be most sinister true crime series yet
Netflix knows how to feed our true crime obsession but upcoming series The Confession Killer might be the most baffling and sinister story yet.
The streaming service has dropped the trailer for the new show, which documents the case of murderer Henry Lee Lucas.
Lucas was convicted of the murder of his mother in 1960, and then two others in 1983. But at this time, he started declaring he had killed hundreds of other people.
Texas Rangers worked closely with Lucas, who ‘solved’ dozens of seemingly unconnected cases as he claimed he would kill ‘at random’ over the course of eight years.
Advertisement Advertisement
However, people are soon quick to question both his confessions, which were often highly detailed, and whether or not they’re the truth.
If he is to be believed, his kill count goes beyond that of some of the world’s most known serial killers…but to this day, he’s only truly been convicted of 11 more people.
The new series examines just what went on in the investigation, and how he darkly embraced an almost celebrity persona as his case became more famous.
Befriending the police force who were trying to convict him, he would be given fast food, strawberry milkshakes and lived a much more comfortable life than the one he had as a drifter.
Speaking to relatives of the victims he claimed to have killed, it becomes clear that some, while they want to believe the case is solved, simply don’t believe that he did it.

Is he the world’s biggest serial killer? Or is it the biggest hoax in US history?
We’re already curious to see what the documentary reveals of the cases, and how they plan to conclude the hundreds of cases Lucas confessed to but may not have done.
Advertisement Advertisement
The Confession Killer is released 6 December on Netflix.
Got a showbiz story? If you've got a story, video or pictures get in touch with the Metro.co.uk Entertainment team by emailing us<EMAIL_ADDRESS>calling 020 3615 2145 or by visiting our Submit Stuff page - we'd love to hear from you.


SDS associate Scrap Metal is continually training its staff on environmental procedures. We follow the guidelines of Occupational Safety and Health Administration (OSHA) and Institute of Scrap Metal Recyclers (ISRI) in every scrap yard and recycling operation. Now certifies our scrap yards and recycling facility as compliant. Our other facilities are in the process of this certification.
Sustainability
We have an obligation to operate with sustainable practices, and are committed to improve the environment through recycling, and to create a clean and safe place for our colleagues to work. We value sustainable and renewable solutions to create a greener tomorrow.
Importance of Recycling Metal
Metals are recycled at a surprisingly low rate throughout the world. In fact, one report published by the United Nations Environment Programme (UNEP) claims less than a third of metals have an end-of-life recycling rate above 50 percent. We are working to grow scrap metal recycling – metal recycling conserves natural resources, saves energy, and reduces greenhouse gas emissions.
Recycling just one ton of steel saves 2,500 lbs. of iron ore, 120 lbs. of limestone, and 1,4000 lbs. of coal. Recycling one ton of aluminum has a similar effect – saving more than a dozen megawatt hours’ worth of electricity. It’s estimated if more steel was recycled, enough energy would be conserved to power 18 million homes in the U.S. for a year.
Our Commitments
- We will seek ways to continually improve our recycling process and evaluate our management procedures for setting and reviewing environmental objectives and targets.
- We will practice pollution prevention through consistent maintenance, education and operation upgrades.
- We actively promote dialog with our neighbours on our environmental policies through participation in a number of different community groups and community notices.
- We keep a close working relationship with professional environmental consultants in order to maintain compliance on all existing and future environmental laws and standards.
- We will involve each of our employees in environmental stewardship and promotion of pollution prevention.


Education
2015 University of OregonBS - Product Design
Minor - Computer Science
Work
2017SSI Shredding Systems
Wilsonville, OR
Web Developer (freelance)
Web design and development: www.unishear.com (public, but in progress).
2017Portland Pet Food Company
Portland, OR
Web Developer & Photographer (freelance)
Web design, development, and studio photography: www.portlandpet
Andrew Gower Design
Remote
CAD Drafter
Create manufacture-ready technical drawings, CAD models, and renders for furniture designs sold globally.
2016Superfab
Portland, OR
Web Designer & Fabricator
Website design and development: www.superfabyardsale.org
Build furniture for clients such as Airbnb, Uber, and Simple.
On-site installation of cabinetry, tables, and shelving.
2012-Tilt Scooters
Mokena, IL
Designer
Website redesign: www.tiltscooters.com
Direction, post-production, and packaging for full-length DVD.
Artwork for wheels, grips, headsets, and apparel.
Tech-packs and contracting for headwear and apparel.
2014Emberex
Eugene, OR
Design Intern
UI/UX design: www.bloxi.com
Graphic design and product development.
Skills
Design software
Illustrator, Photoshop, InDesign, Premiere, Final Cut, and After Effects.
CAD software
Solidworks, Rhino, SketchUp, and KeyShot.
Languages
HTML, CSS, JavaScript, and jQuery. Experience building for Wordpress, Shopify, Squarespace, and ExpressionEngine.


 Nah, me neither
3/30/2019, 11:09:03 AM
Favs: 5
Retweets: 0
link


Google Analytics tracker demo
Pssst. View source to see tracker being initialised
Track a tap event!


[gecode-users] define a domain for a SetVar using a IntVar
Steven Villegas
sagrath23 at gmail.com
Thu Mar 10 04:44:50 CET 2011
Hi,
Actually, i'm trying to implement a musical CSP using GECODE; in this
CSP i have to define a set based on the result based on the result of
the multiplication of two IntVar (each one define the root and
interval of a chord) and int (this define the number of notes of the
chord).
My Question is: how i define the domain of a SetVar using a IntVar??
that's possible???
Thanks.
--
Steven Augusto Villegas Castellanos
Tecnologo en Sistemas.
Colaborador Parquesoft Tuluá.
Estudiante Ingenieria de Sistemas.
Universidad del Valle.
More information about the users
mailing list


This Galaxy is maintained by the Bioinformatics group of the department of Pathology.
For any questions regarding this Galaxy server, please contact us.


Installation¶
Requirements¶
Grove is compatible with Python 2.7.
Install Grove¶
The preferred method for installing Grove is easy_install or pip.
easy_install grove
or
pip install grove
Alternatively, you can download or clone the repository and build from sources.
python setup.py install
Grove also optionally uses the ETE Toolkit for graphing parse trees generated during grammatical evolution. Installation through conda is required to use this feature.
conda install grove


--- layout: archive title: Zulip Chat Archive permalink: /stream/113489-new-members/topic/two.20functions.20with.20the.20same.20name.html ---
Is there a reason why lean doesn't let me make two functions with the same name? Even if they take different arguments?
It's confusing, and could even lead to people thinking they've proved something, when they've actually proved something about a function with the same name. You can always use namespaces.
You can have two functions in different namespaces with the same name, where both namespaces are open at once, and lean will use the type to disambiguate
The basic reason why is because lean needs a name for the function, one that is not ambiguous, and if you overload it what will that name be? How can you refer to it precisely?
I know I should be using Pi types for this sort of thing, but my two functions take in a different amount of arguments, so I think I might have to use if then else.
You can also use default values and optional arguments to make a single definition have varying number of arguments
Also note that +
is kind of like several functions with the same name -- and looking at the types of the inputs is exactly how Lean decides which function to actually use.
Last updated: Apr 24 2020 at 19:02 UTC


UN Security Council reaffirms stance on Chemical Weapons
United Nations: The UN Security Council has reaffirmed its stance on chemical weapons after it was briefed by the head of the Organization for the Prohibition of Chemical Weapons (OPCW) on the use of chemical weapons in Syria. In a presidential statement issued after Fernando Arias’ November 5 briefing, the Security Council on Friday reaffirmed that the proliferation of chemical weapons, as well as their means of delivery, constitutes a threat to international peace and security, Xinhua news agency reported. The council reaffirmed that the use of chemical weapons constitutes a violation of international law and condemned in the strongest possible terms the use of chemical weapons. (IANS)
Also Read: Vietnam Elected To Asian Seat On UN Security Council
Also Watch: Actress Amrita Gogoi promotes ‘Pratighat’ at Namti in Sivasagar


Face Manager¶
Overview¶
The architecture of the face manager is built around the concept of interfaces, which allows for a modular and extensible deployment.
Interfaces are used to implement in isolation various sources of information which help with the construction of faces (such as network interface and service discovery), and with handling the heterogeneity of host platforms.
Platform and supported interfaces¶
Currently, Android, Linux and MacOS are supported through the following interfaces:
hicn-light [Linux, Android, MacOS, iOS] An interface to the hicn-light forwarder, and more specifically to the Face Table and FIB data structures. This component is responsible to effectively create, update and delete faces in the forwarder, based on the information provided by third party interfaces, plus adding default routes for each of the newly created face. The communication with the forwarder is based on the hicn control library (
libhicnctrl
).netlink [Linux, Android] The default interface on Linux systems (including Android) to communicate with the kernel and receive information from various sources, including link and address information (both IPv4 and IPv6) about network interfaces.
android_utility [Android only] Information available through Netlink is limited with respect to cellular interfaces. This component allows querying the Android layer through SDK functions to get the type of a given network interface (Wired, WiFi or Cellular).
bonjour [Linux, Android] This component performs remote service discovery based on the bonjour protocol to discover a remote hICN forwarder that might be needed to establish overlay faces.
network_framework [MacOS, iOS]
This component uses the recommended Network framework on Apple devices, which provided all required information to query faces in a unified API: link and address information, interface types, and bonjour service discovery.
Interface API¶
TODO
Developing a new interface¶
Dummy template¶
The face manager source code includes a template that can be used as a skeleton
to develop new faces. It can be found in src/interface/dummy/dummy.{h,c}
. Both
include guard and specific interface functions are prefixed by a (short)
identifier which acts as a namespace for interface specific code (in our case
the string ‘dummy_’).
Registration and instantiation of the different interfaces is currently done at
compile time in the file src/api.c
, and the appropriate hooks to use the dummy
interface are available in the code between #if 0/#endif
tags.
Interface template header; configuration parameters¶
All interfaces have a standard interface defined in src/interface.{h,c}
, and
as such the header file is only used to specify the configuration parameters of
the interface, if any.
In the template, these configuration options are empty:
/*
* Configuration data
*/
typedef struct {
/* ... */
} dummy_cfg_t;
Overview of the interface template¶
The file starts with useful includes:
- the global include `<hicn/facemgr.h>` : this provides public facing elements
of the face manager, such the standard definition of faces (`face_t` from
`libhicnctrl`), helper classes (such as `ip_address_t` from `libhicn`), etc.
- common.h
- facelet.h : facelets are the basic unit of communication between the face
manager and the different interfaces. They are used to construct the faces
incrementally.
- interface.h : the parent class of interfaces, such as the current dummy
interface.
Each interface can hold a pointer to an internal data structure, which is declared as follows:
/*
* Internal data
*/
typedef struct {
/* The configuration data will likely be allocated on the stack (or should
* be freed) by the caller, we recommend to make a copy of this data.
* This copy can further be altered with default values.
*/
dummy_cfg_t cfg;
/* ... */
int fd; /* Sample internal data: file descriptor */
} dummy_data_t;
We find here a copy of the configuration settings (which allows the called to instantiate the structure on the stack), as well as a file descriptor (assuming most interfaces will react on events on a file descriptor).
The rest of the file consists in the implementation of the interface, in
particular the different function required by the registration of a new
interface to the system. They are grouped as part of the interface_ops_t
data
structure declared at the end of the file:
interface_ops_t dummy_ops = {
.type = "dummy",
.initialize = dummy_initialize,
.finalize = dummy_finalize,
.callback = dummy_callback,
.on_event = dummy_on_event,
};
The structure itself is declared and documented in src/interface.h
/**
* \brief Interface operations
*/
typedef struct {
/** The type given to the interfaces */
char * type;
/* Constructor */
int (*initialize)(struct interface_s * interface, void * cfg);
/* Destructor */
int (*finalize)(struct interface_s * interface);
/* Callback upon file descriptor event (iif previously registered) */
int (*callback)(struct interface_s * interface);
/* Callback upon facelet events coming from the face manager */
int (*on_event)(struct interface_s * interface, const struct facelet_s * facelet);
} interface_ops_t;
Such an interface has to be registered first, then one (or multiple) instance(s)
can be created (see src/interface.c
for the function prototypes, and
src/api.c
for their usage).
interface registration:
extern interface\_ops\_t dummy\_ops;
/* [...] */
rc = interface\_register(&dummy\_ops);
if (rc < 0)
goto ERR_REGISTER;
interface instantiation:
#include "interfaces/dummy/dummy.h"
/* [...] */
rc = facemgr_create_interface(facemgr, "dummy0", "dummy", &facemgr->dummy);
if (rc < 0) {
ERROR("Error creating 'Dummy' interface\n");
goto ERR_DUMMY_CREATE;
}
Implementation of the Interface API¶
We now quickly go other the different functions, but their usage will be better understood through the hands-on example treated in the following section.
In the template, the constructor is the most involved as it need to:
initialize the internal data structure:
dummy_data_t * data = malloc(sizeof(dummy_data_t));
if (!data)
goto ERR_MALLOC;
interface->data = data;
process configuration parameters, eventually setting some default values:
/* Use default values for unspecified configuration parameters */
if (cfg) {
data->cfg = *(dummy_cfg_t *)cfg;
} else {
memset(&data->cfg, 0, sizeof(data->cfg));
}
open an eventually required file descriptor
For the sake of simplicity, the current API only supports a single file descriptor per-interface, and it has to be created in the constructor, and set as the return value so as to be registered by the system, and added to the event loop for read events. A return value of 0 means the interface does not require any file descriptor. As usual, a negative return value indicates an error.
data->fd = 0;
/* ... */
/*
* We should return a negative value in case of error, and a positive value
* otherwise:
* - a file descriptor (>0) will be added to the event loop; or
* - 0 if we don't use any file descriptor
*/
return data->fd;
While support for multiple file descriptors might be added in the future, an
alternative short-term implementation might consider the instanciation of
multiple interface, as is done for Bonjour in the current codebase, in
src/api.c
.
Data reception on the file descriptor will get the callback function called, in
our case dummy_callback
. Finally, the destructor dummy_finalize
should close
an eventual open file descriptor.
In order to retrieve the internal data structure, that should in particular store such a file descriptor, all other function but the constructor can dereference it from the interface pointer they receive as parameter:
dummy_data_t * data = (dummy_data_t*)interface->data;
Raising and Receiving Events¶
An interface will receive events in the form of a facelet through the *_on_event
function. It can then use the facelet API we have described above to read
information about the face.
As this information is declared const, the interface can either create a new facelet (identified by the same netdevice and protocol family), or eventually clone it.
The facelet event can then be defined and raised to the face maanger for further processing through the following code:
facelet_set_event(facelet, EVENT_TYPE_CREATE);
interface_raise_event(interface, facelet);
Here the event is a facelet creation (EVENT_TYPE_CREATE
). The full facelet API
and the list of possible event types is available in src/facelet.h
Integration in the Build System¶
The build system is based on CMake. Each interface should declare its source
files, private and public header files, as well as link dependencies in the
local CMakeLists.txt
file.
TODO: detail the structure of the file
Hands-On¶
Architecture¶
In order to better illustrate the development of a new interface, we will
consider the integration of a sample server providing a signal instructing the
face manager to alternatively use either the WiFi or the LTE interface. The code
of this server is available in the folder examples/updownsrv/
, and the
corresponding client code in examples/updowncli
.
Communication between client and server is done through unix sockets over an abstract namespace (thereby not using the file system, which would cause issues on Android). The server listens for client connections, and periodically broadcast a binary information to all connected clients, in the form of one byte equal to either \0 (which we might interpret as enable LTE, disable WiFi), or \1 (enable WiFi, disable LTE).
Our objective is to develop a new face manager interface that would listen to such event in order to update the administrative status of the current faces. This would thus alternatively set the different interfaces admnistratively up and down (which takes precedence over the actual status of the interface when the forwarder establishes the set of available next hops for a given prefix). The actual realization of such queries will be ultimately performed by the hicn-light interface.
Sample Server and Client¶
In the folder containing the source code of hICN, the following commands allow to run the sample server:
cd ctrl/facemgr/examples/updownsrv
make
./updownsrv
The server should display “Waiting for clients…”
Similar commands allow to run the sample client:
cd ctrl/facemgr/examples/updowncli
make
./updowncli
The client should display “Waiting for server data…”, then every couple of seconds display either “WiFi” or “LTE”.
Face Manager Interface¶
An example illustrating how to connect to the dummy service from updownsrv
is
provided as the updown
interface in the facemgr source code.
This interface periodically swaps the status of the LTE interface up and down. It is instantiated as part of the facemgr codebase when the code is compiled with the ``-DWITH_EXAMPLE_UPDOWN` cmake option.


Details: Brewed with well water (the original well at the Old Brewery, sunk in 1758, is still in use, with the hard well water being drawn from 85 feet underground); fermented in ‘stone Yorkshire squares’ to create an almost opaque, wonderfully silky and smooth textured ale with a complex medium dry palate and bittersweet finish.
Style: Oatmeal Stout
BA Score: 94
ABV: 5.0%
Notes: Not too much head. Very creamy texture, sweet flavor. Nice dark color for a stout. Goes down very nicely.
Details: Aecht Schlenkerla Rauchbier Maerzen is unpasteurized and made entirely from smoked malt. Heller Trum is one of the very few breweries left in Germany that does its own malting. During germination of barley, the smoke of beech wood fire is "forced" through the germinating barley for a total of between 24 and 36 hours at varying temperatures; therefore, Aecht Schlenkerla Rauchbier gains an unusual smoky dryness from this process.
Style: Rauchbier
BA Score: 91
ABV: 5.4%
Notes: Very smokey flavor, but not too much. First taste kind of shocks you, but it gets better with each drink. Not too bad of a finish, either.
Details: Award-winning Kentucky Bourbon Barrel Ale® is a unique sipping beer with the distinctive nose of a well-crafted bourbon. Our Kentucky Ale® is aged for up to six weeks in freshly decanted bourbon barrels from some of Kentucky’s finest distilleries. Subtle yet familiar flavors of vanilla and oak are imparted to this special ale as it rests in the charred barrels. Pleasantly smooth and robust.
Style: English Strong Ale
BA Score: 82
ABV: 8.2%
Notes: My experience with bourbon barrel ales has not been good. Every one I've had has been too strong on the bourbon flavor. This one is nearly perfect. Very smooth, with a tiny hint of bourbon on the finish. Extremely pleased.


Payment History
6/1/2016 6 session credits remaining
Number of sessions: 10
Paid by: Brian Skinner
Payment ID: 305254
5 sessions at $110 with 10% off
5 sessions at $130 with 10% off
Total: $1080
Charlie - Tuesday, October 5th 2015
Sophie - Tuesday, October 5th 2015
Charlie - Tuesday, October 5th 2015
Sophie - Tuesday, October 5th 2015
5/16/2016
Number of sessions: 2
Paid by: Brian Skinner
Payment ID: 305204
1 session at $110
1 session at $130
Total: $240
Charlie - Tuesday, October 5th 2015
Sophie - Tuesday, October 5th 2015
2/23/2016
Number of sessions: 10
Paid by: Brian Skinner
Payment ID: 305111
5 sessions at $110 with 10% off
5 sessions at $130 with 10% off
Total: $1080
Charlie - Tuesday, October 5th 2015
Sophie - Tuesday, October 5th 2015
Charlie - Tuesday, October 5th 2015
Sophie - Tuesday, October 5th 2015
Charlie - Tuesday, October 5th 2015
Sophie - Tuesday, October 5th 2015
Charlie - Tuesday, October 5th 2015
Sophie - Tuesday, October 5th 2015
Charlie - Tuesday, October 5th 2015
Sophie - Tuesday, October 5th 2015


Active speaker
Click one of the videos below to make that user the "active speaker".
Self view 1
Self view 2
Self view 3


Rocko's dog food makes feeding your dog easy with food delivered straight to your door weekly so you will never run out of food again.
At Rocko's dog food, we only use the finest quality of dog food so you can rest assured knowing that your dog will be eating the best food around.
Whenever you receive a bag of food from Rocko, you will also be helping a dog in need with support for food, medical and housing.


The scenario:
A group of people are working on a sensitive data set that for practical reasons needs to be stored in a place that we’re not 100% happy with the security (e.g., Dropbox), or we’re concerned that files stored in plain text on users computers (e.g. laptops) may lead to the data being compromised.
If the data can be stored encrypted but everyone in the group can still read and write the data then we’ve improved the situation somewhat. But organising for everyone to get a copy of the key to decrypt the data files is non-trivial. The workflow described here aims to simplify this procedure using lower-level functions in the cyphr
package.
The general procedure is this:
A person will set up a set of personal keys and a key for the data. The data key will be encrypted with their personal key so they have access to the data but nobody else does. At this point the data can be encrypted.
Additional users set up personal keys and request access to the data. Anyone with access to the data can grant access to anyone else.
Before doing any of this, everyone needs to have ssh keys set up. By default the package will use your ssh keys found at “~/.ssh”; see the main package vignette for how to use this.
For clarity here we will generate two sets of key pairs for two actors Alice and Bob:
path_key_alice <- cyphr::ssh_keygen(password = FALSE)
path_key_bob <- cyphr::ssh_keygen(password = FALSE)
These would ordinarily be on different machines (nobody has access to anyone else’s private key) and they would be password protected. In the function calls below, all the path_user
arguments would be omitted.
We’ll store data in the directory data
; at present there is nothing there (this is in a temporary directory for compliance with CRAN policies but would ordinarily be somewhere persistent and under version control ideally).
## character(0)
First, create a personal set of keys. These will be shared across all projects and stored away from the data. Ideally one would do this with ssh-keygen
at the command line, following one of the many guides available. A utility function ssh_keygen
(which simply calls ssh-keygen
for you) is available in this package though. You will need to generate a key on each computer you want access from. Don’t copy the key around. If you lose your user key you will lose access to the data!
Second, create a key for the data and encrypt that key with your personal key. Note that the data key is never stored directly - it is always stored encrypted by a personal key.
## Generating data key
## Authorising ourselves
## Adding key 91:e2:bf:0b:1d:d0:0b:c9:81:af:c7:9a:8a:60:c6:0d:ff:ac:17:b6:63:0c:bc:56:e8:cc:cf:c4:bd:cf:70:25
## user: rich
## host: wpia-dide136
## date: 2022-06-20 10:55:37
## Verifying
The data key is very important. If it is deleted, then the data cannot be decrypted. So do not delete the directory data_dir/.cyphr
! Ideally add it to your version control system so that it cannot be lost. Of course, if you’re working in a group, there are multiple copies of the data key (each encrypted with a different person’s personal key) which reduces the chance of total loss.
This command can be run multiple times safely; if it detects it has been rerun and the data key will not be regenerated.
## Already set up at /tmp/RtmpI8wXy7/data
## Verifying
Third, you can add encrypted data to the directory (or to anywhere really). When run, cyphr::config_data
will verify that it can actually decrypt things.
This object can be used with all the cyphr
functions (see the “cyphr” vignette; vignette("cyphr")
)
filename <- file.path(data_dir, "iris.rds")
cyphr::encrypt(saveRDS(iris, filename), key)
dir(data_dir)
## [1] "iris.rds"
The file is encrypted and so cannot be read with readRDS
:
## Error in readRDS(filename): unknown input format
But we can decrypt and read it:
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
Fourth, have someone else join in. Recall that to simulate another person here, I’m going to pass an argument path_user = path_key_bob
though to the functions. This contains the path to “Bob”’s ssh keypair. If run on an actually different computer this would not be needed; this is just to simulate two users in a single session for this vignette (see minimal example below where this is simulated). Again, typically this user would also not use the cyphr::ssh_keygen
function but use the ssh-keygen
command from their shell.
We’re going to assume that the user can read and write to the data. This is the case for my use case where the data are stored on dropbox and will be the case with GitHub based distribution, though there would be a pull request step in here.
This user cannot read the data, though trying to will print a message explaining how you might request access:
But bob
is your collaborator and needs access! What they need to do is run:
## A request has been added
## Email someone with access to add you
##
## hash: 21:54:9b:4a:d9:d1:68:37:7c:cd:9e:31:80:be:78:fe:54:29:8c:07:7e:83:e2:1a:7b:6c:5a:d8:9a:11:80:98
##
## If you are using git, you will need to commit and push first:
##
## git add .cyphr
## git commit -m "Please add me to the dataset"
## git push
(again, ordinarily you would not need the bob
bit here)
The user should the send an email to someone with access and quote the hash in the message above.
Fifth, back on the first computer we can authorise the second user. First, see who has requested access:
## 1 key:
## 21:54:9b:4a:d9:d1:68:37:7c:cd:9e:31:80:be:78:fe:54:29:8c:07:7e:83:e2:1a:7b:6c:5a:d8:9a:11:80:98
## user: rich
## host: wpia-dide136
## date: 2022-06-20 10:55:37
We can see the same hash here as above (21549b4ad9d168377ccd9e3180be78fe54298c077e83e21a7b6c5ad89a118098
)
…and then grant access to them with the cyphr::data_admin_authorise
function.
## There is 1 request for access
## Adding key 21:54:9b:4a:d9:d1:68:37:7c:cd:9e:31:80:be:78:fe:54:29:8c:07:7e:83:e2:1a:7b:6c:5a:d8:9a:11:80:98
## user: rich
## host: wpia-dide136
## date: 2022-06-20 10:55:37
## Added 1 key
## If you are using git, you will need to commit and push:
##
## git add .cyphr
## git commit -m "Authorised rich"
## git push
If you do not specify yes = TRUE
will prompt for confirmation at each key added.
This has cleared the request queue:
## (empty)
and added it to our set of keys:
## 2 keys:
## 21:54:9b:4a:d9:d1:68:37:7c:cd:9e:31:80:be:78:fe:54:29:8c:07:7e:83:e2:1a:7b:6c:5a:d8:9a:11:80:98
## user: rich
## host: wpia-dide136
## date: 2022-06-20 10:55:37
## 91:e2:bf:0b:1d:d0:0b:c9:81:af:c7:9a:8a:60:c6:0d:ff:ac:17:b6:63:0c:bc:56:e8:cc:cf:c4:bd:cf:70:25
## user: rich
## host: wpia-dide136
## date: 2022-06-20 10:55:37
Finally, as soon as the authorisation has happened, the user can encrypt and decrypt files:
key_bob <- cyphr::data_key(data_dir, path_user = path_key_bob)
head(cyphr::decrypt(readRDS(filename), key_bob))
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
As above, but with less discussion:
Setup, on Alice’s computer:
## Generating data key
## Authorising ourselves
## Adding key 91:e2:bf:0b:1d:d0:0b:c9:81:af:c7:9a:8a:60:c6:0d:ff:ac:17:b6:63:0c:bc:56:e8:cc:cf:c4:bd:cf:70:25
## user: rich
## host: wpia-dide136
## date: 2022-06-20 10:55:37
## Verifying
Get the data key key:
Encrypt a file:
Request access, on Bob’s computer:
## A request has been added
## Email someone with access to add you
##
## hash: 21:54:9b:4a:d9:d1:68:37:7c:cd:9e:31:80:be:78:fe:54:29:8c:07:7e:83:e2:1a:7b:6c:5a:d8:9a:11:80:98
##
## If you are using git, you will need to commit and push first:
##
## git add .cyphr
## git commit -m "Please add me to the dataset"
## git push
Alice authorises this request::
## There is 1 request for access
## Adding key 21:54:9b:4a:d9:d1:68:37:7c:cd:9e:31:80:be:78:fe:54:29:8c:07:7e:83:e2:1a:7b:6c:5a:d8:9a:11:80:98
## user: rich
## host: wpia-dide136
## date: 2022-06-20 10:55:38
## Added 1 key
## If you are using git, you will need to commit and push:
##
## git add .cyphr
## git commit -m "Authorised rich"
## git push
Bob can get the data key:
Bob can read the secret data:
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
Encryption does not work through security through obscurity; it works because we can rely on the underlying maths enough to be open about how things are stored and where.
Most encryption libraries require some degree of security in the underlying software. Because of the way R works this is very difficult to guarantee; it is trivial to rewrite code in running packages to skip past verification checks. So this package is not designed to (or able to) avoid exploits in your running code; an attacker could intercept your private keys, the private key to the data, or skip the verification checks that are used to make sure that the keys you load are what they say they are. However, the data are safe; only people who have keys to the data will be able to read it.
cyphr
uses two different encryption algorithms; it uses RSA encryption via the openssl
package for user keys, because there is a common file format for these keys so it makes user configuration easier. It uses the modern sodium package (and through that the libsodium library) for data encryption because it is very fast and simple to work with. This does leave two possible points of weakness as a vulnerability in either of these libraries could lead to an exploit that could allow decryption of your data.
Each user has a public/private key pair. Typically this is in ~/.ssh/id_rsa.pub
and ~/.ssh/id_rsa
, and if found these will be used. Alternatively the location of the keypair can be stored elsewhere and pointed at with the USER_KEY
or USER_PUBKEY
environment variables. The key may be password protected (and this is recommended!) and the password will be requested without ever echoing it to the terminal.
The data directory has a hidden directory .cyphr
in it.
## [1] ".cyphr" "iris.rds"
This does not actually need to be stored with the data but it makes sense to (there are workflows where data is stored remotely where storing this directory might make sense). The “keys” directory contains a number of files; one for each person who has access to the data.
## [1] "21549b4ad9d168377ccd9e3180be78fe54298c077e83e21a7b6c5ad89a118098"
## [2] "91e2bf0b1dd00bc981afc79a8a60c60dffac17b6630cbc56e8cccfc4bdcf7025"
## [1] "21549b4ad9d168377ccd9e3180be78fe54298c077e83e21a7b6c5ad89a118098"
## [2] "91e2bf0b1dd00bc981afc79a8a60c60dffac17b6630cbc56e8cccfc4bdcf7025"
(the file test
is a small file encrypted with the data key used to verify everything is working OK).
Each file is stored in RDS format and is a list with elements:
h <- names(cyphr::data_admin_list_keys(data_dir))[[1]]
readRDS(file.path(data_dir, ".cyphr", "keys", h))
## $user
## [1] "rich"
##
## $host
## [1] "wpia-dide136"
##
## $date
## [1] "2022-06-20 10:55:38 BST"
##
## $pub
## [2048-bit rsa public key]
## md5: d85323cf26b481aff3975c96725807ad
##
## $key
## [1] ac 36 2f 00 66 02 3a 62 46 83 b4 e6 4f ae eb 79 86 48 81 10 fb c3 8a 03 e1
## [26] 12 c5 9e 9f 42 c0 4d 5b 69 f8 d2 bd 98 d8 7d 30 66 92 6a 60 0d 99 df e9 67
## [51] 0e d9 dd ee 55 94 b4 62 db d2 b1 47 58 2b e5 dd ae 5d 83 22 4c ca e4 2b 5f
## [76] c5 4c fd 46 c3 43 50 91 6b 92 bd 3c d2 11 e5 32 71 21 c9 24 24 64 82 6c cf
## [101] a0 56 c4 42 c5 a1 9c 73 42 d4 6a a2 4d cb 9d 82 b4 05 5b 59 24 ea 16 28 c3
## [126] 8a cc a1 16 4e 13 f3 9f fc ce 75 cb 13 8a 5f 8a 82 5c 03 84 33 db 46 e0 44
## [151] 22 ec 85 32 3a 5e b8 b8 05 6c cf 4b 0c af e1 e1 78 90 4c 8b 2b b0 a6 f5 d5
## [176] df 5f bc 61 48 3f a4 f7 03 f0 09 08 0d 0f cd 4c 1b 1e f4 b1 8f 3d 72 57 36
## [201] db 09 11 e0 46 da 8f e5 90 a9 b9 23 9b 0c ff a4 9a d9 21 b5 1c 70 7f 77 f8
## [226] 24 c8 54 4d b1 b3 7d 18 4e 28 ea 9c 1e 8b 4b a3 a7 bd 62 89 83 83 22 74 c7
## [251] 6d dd 14 2d 2e ad
You can see that the hash of the public key is the same as name of the stored file here (which is used to prevent collisions when multiple people request access at the same time).
## [1] "21549b4ad9d168377ccd9e3180be78fe54298c077e83e21a7b6c5ad89a118098"
When a request is posted it is an RDS file with all of the above except for the key
element, which is added during authorisation.
(Note that the verification relies on the package code not being attacked, and given R’s highly dynamic nature an attacker could easily swap out the definition for the verification function with something that always returns TRUE
.)
When an authorised user creates the data_key
object (which allows decryption of the data) secret
will:
~/.ssh/id_rsa
)$key
element from the list above).In the Dropbox scenario, non-password protected keys will afford only limited protection. This is because even though the keys and data are stored separately on Dropbox, they will be in the same place on a local computer; if that computer is lost then the only thing preventing an attacker recovering the data is security through obscurity (the data would appear to be random junk but they will be able to run your analysis scripts as easily as you can). Password protected keys will improve this situation considerably as without a password the data cannot be recovered.
The data is not encrypted during a running R session. R allows arbitrary modification of code at runtime so this package provides no security from the point where the data can be decrypted. If your computer was compromised then stealing the data while you are running R should be assumed to be straightforward.


Toggle navigation
pubmedBatch
{{ user }}
{{ folder }}
Delete folder
{% for d in files %}
x
{{ d }}
{% endfor %}
Column to search
Terms to search
(separated by space)
AND
OR
Your email
CSV (with a header) to upload
Pred cutoff
Verbose
Genes to highlight
Genes to mask
Submit to pubmed
Export to CSV
Rename record
Warning
Are you sure to delete this folder?
Message


--- layout: page title: Authoring css: "/css/aboutme.css" ---
(2022) Data Science Summit Conference
Talk: "Self-supervised learning in medicine: X-ray analysis with few labeled examples". [SLIDES]
(2023 - present) Intermediate Deep Learning with PyTorch
Interactive online course I teach at DataCamp. Learn about fundamental deep learning architectures such as CNNs, RNNs, LSTMs, and GRUs for modeling image and sequential data.
(2022) Deep Learning for Image Analysis
I've been a trainer during the European Molecular Biology Laboratory's (EMBL) deep learning course, guiding particpants as they worked on their project involving segmenting and classifying cells in microscope images.
(2020 - present) Bayesian Data Analysis in Python
Interactive online course I teach at DataCamp. Learn all about the advantages of the Bayesian approach to data analysis, and apply it to a variety of real-world use cases.
(2019 - present) Handling Missing Data with Imputations in R
Interactive online course I teach at DataCamp. Diagnose, visualize and treat missing data with a range of imputation techniques with tips to improve your results.
NannyML: Bad Machine Learning Models Can Still Be Well-Calibrated
NannyML: Detecting Covariate Shift: A Guide to the Multivariate Approach
NannyML: How to detect data drift with hypothesis testing
neptune.ai: Feature Selection Methods and How to Choose Them
neptune.ai: Reducing Pipeline Debt With Great Expectations
neptune.ai: Organizing ML Monorepo With Pants
MLOps Community: Monitoring Regression Models Without Ground-Truth
Nota AI: AI Model Optimization Has Never Been Easier
DataCamp: Regularization in R Tutorial: Ridge, Lasso and Elastic Net
DataCamp: Visualize Missing Data with VIM Package


Time is: {{mytime | date:'shortTime' }}
A lightweight & configurable timepicker directive.
All settings can be provided as attributes in the <timepicker>
or globally configured through the timepickerConfig
.
ng-model
: The Date object that provides the time state.
hour-step
(Defaults: 1) : Number of hours to increase or decrease when using a button.
minute-step
(Defaults: 1) : Number of minutes to increase or decrease when using a button.
show-meridian
(Defaults: true) : Whether to display 12H or 24H mode.
meridians
(Defaults: null) : Meridian labels based on locale. To override you must supply an array like ['AM', 'PM'].
readonly-input
(Defaults: false) : Whether user can type inside the hours & minutes input.
mousewheel
(Defaults: true) : Whether user can scroll inside the hours & minutes input to increase or decrease it's values.
arrowkeys
(Defaults: true) : Whether user can use up/down arrowkeys inside the hours & minutes input to increase or decrease it's values.


{% extends "bootstrap/layouts/base_navbar_responsive.html" %} {%- from "collection/macros.html" import show_collected_things with context %} {%- from "frontend/macros.html" import follow_toggle %} {%- from "talk/macros.html" import show_threads_list with context %} {% block title %} {{collection.title}} {% endblock %} {% block content %}
{{collection.description|markdown}}
{{subcollection.short_description}}
{{subcollection.description}}


|
Ising_OPV
v4.0.2
Generates and analyzes model bulk heterojunction morphologies in a parallel computing environment
|
This class contains the parameters used by Ising_OPV and functions to import parameter files. More...
#include <Parameters.h>
Public Member Functions | |
| Parameters () | |
| Default constructor that creates an empty Parameters object. More... | |
| bool | checkParameters () const |
| Checks the current values of all parameters to ensure that they are valid. More... | |
| bool | importParameters (std::ifstream ¶meterfile) |
| Imports the parameters from the input parameter filestream. More... | |
Public Attributes | |
| int | Length = 0 |
| x-direction size of the lattice More... | |
| int | Width = 0 |
| y-direction size of the lattice More... | |
| int | Height = 0 |
| z-direction size of the lattice More... | |
| bool | Enable_periodic_z = true |
| Choose whether or not to enable z-direction periodic boundary conditions. More... | |
| double | Mix_fraction = 0.0 |
| volume fraction of type1 sites More... | |
| double | Interaction_energy1 = 0.0 |
| energetic favorability for type1-type1 interactions over type1-type2 interactions More... | |
| double | Interaction_energy2 = 0.0 |
| energetic favorability for type2-type2 interactions over type1-type2 interactions More... | |
| int | MC_steps = 0 |
| number of MC steps to be executed (determines number of Ising swapping iterations), which sets the duration of the phase separation process More... | |
| bool | Enable_smoothing = false |
| choose whether or not to perform domain smoothing More... | |
| double | Smoothing_threshold = 0.0 |
| cutoff threshold for the smoothing algorithm More... | |
| bool | Enable_rescale = false |
| choose whether or nmot to perform lattice rescaling More... | |
| int | Rescale_factor = 0 |
| rescale factor to be used More... | |
| bool | Enable_shrink = false |
| chose whether or not to shrink the lattice by 1/rescale_factor instead of expand it More... | |
| bool | Enable_interfacial_mixing = false |
| choose whether or not to perform interfacial mixing More... | |
| double | Interface_width = 0.0 |
| interfacial width of the mixed interface More... | |
| double | Interface_conc = 0.0 |
| concentration of type1 sites in the interfacial mixed region More... | |
| bool | Enable_analysis_only = false |
| choose whether or not to only perform analysis on an imported morphology (no modification of the morphology) More... | |
| bool | Enable_correlation_calc = false |
| choose whether or not to perform the domain size calculation using the pair-pair autocorrelation function method More... | |
| int | N_sampling_max = 0 |
| maximum number of sites to be sampled for calculating the autocorrelation function More... | |
| bool | Enable_mix_frac_method = false |
| choose whether or not to use the mix fraction method for determining the domain size from the autocorrelation data More... | |
| bool | Enable_e_method = false |
| choose whether or not to use the 1/e method for determining the domain size from the autocorrelation data More... | |
| bool | Enable_extended_correlation_calc = false |
| choose whether or not to extend the autocorrelation function calculation to the specified cutoff distance More... | |
| int | Extended_correlation_cutoff_distance = 0 |
| cutoff distance for the extended autocorrelation function calculation More... | |
| bool | Enable_interfacial_distance_calc = false |
| choose whether or not to calculate the interfacial distance histograms More... | |
| bool | Enable_tortuosity_calc = false |
| choose whether or not to calculate the end-to-end tortuosity histograms and island volume fraction More... | |
| bool | Enable_reduced_memory_tortuosity_calc = false |
| choose whether or not to perform the tortuosity calculation using an algorithm that takes longer, but uses less memory More... | |
| bool | Enable_depth_dependent_calc = false |
| choose whether or not to calculate and output the film depth dependent morphology characteristics More... | |
| bool | Enable_areal_maps_calc = false |
| choose whether or not to calculate and output areal mappings of the morphology characteristics More... | |
| bool | Enable_checkerboard_start = false |
| choose whether or not to start from a alternating checkerboard-like configuration instead of a random blend (creates 0.5 mix fraction) More... | |
| bool | Enable_growth_pref = false |
| choose whether or not to implement a directional-dependent interaction energy that caues directional-dependent domain growth More... | |
| int | Growth_direction = 0 |
| direction that has the modified interaction energy with 1=x,2=y,3=z More... | |
| double | Additional_interaction = 0.0 |
| amount that the interaction energy is modified by in the specified direction More... | |
| bool | Enable_export_compressed_files = true |
| choose whether or not the output morphology data file is in compressed format More... | |
| bool | Enable_export_cross_section = false |
| choose whether or not to output uncompressed data for a cross-section (x=Length/2 plane) of the morphology More... | |
| bool | Enable_import_morphologies = false |
| choose whether or not the import a morphology set from the working directory More... | |
| bool | Enable_import_tomogram = false |
| choose whether or not to import a tomogram dataset from the working directory More... | |
| std::string | Tomogram_name = "" |
| name of the tomogram dataset corresponding to the Tomogram_name.xml and Tomogram_name.raw files More... | |
| double | Desired_unit_size = 0.0 |
| desired unit size (resolution) of the output morphology dataset extracted from the tomogram More... | |
| double | Mixed_frac = 0.0 |
| volume fraction of the mixed phase in the tomogram dataset More... | |
| double | Mixed_conc = 0.0 |
| volume fraction of type1 sites in the mixed phase of the tomogram dataset More... | |
| int | N_extracted_segments = 0 |
| number of equal size cuboids to extract from the tomogram dataset More... | |
| int | N_variants = 0 |
| number of random variants to create from each extracted cuboid segment More... | |
This class contains the parameters used by Ising_OPV and functions to import parameter files.
| Ising_OPV::Parameters::Parameters | ( | ) |
Default constructor that creates an empty Parameters object.
| bool Ising_OPV::Parameters::checkParameters | ( | ) | const |
Checks the current values of all parameters to ensure that they are valid.
| bool Ising_OPV::Parameters::importParameters | ( | std::ifstream & | parameterfile | ) |
Imports the parameters from the input parameter filestream.
| parameterfile | is an open input filestream that points to the beginning of a correctly formatted parameter file. |
| double Ising_OPV::Parameters::Additional_interaction = 0.0 |
amount that the interaction energy is modified by in the specified direction
| double Ising_OPV::Parameters::Desired_unit_size = 0.0 |
desired unit size (resolution) of the output morphology dataset extracted from the tomogram
| bool Ising_OPV::Parameters::Enable_analysis_only = false |
choose whether or not to only perform analysis on an imported morphology (no modification of the morphology)
| bool Ising_OPV::Parameters::Enable_areal_maps_calc = false |
choose whether or not to calculate and output areal mappings of the morphology characteristics
| bool Ising_OPV::Parameters::Enable_checkerboard_start = false |
choose whether or not to start from a alternating checkerboard-like configuration instead of a random blend (creates 0.5 mix fraction)
| bool Ising_OPV::Parameters::Enable_correlation_calc = false |
choose whether or not to perform the domain size calculation using the pair-pair autocorrelation function method
| bool Ising_OPV::Parameters::Enable_depth_dependent_calc = false |
choose whether or not to calculate and output the film depth dependent morphology characteristics
| bool Ising_OPV::Parameters::Enable_e_method = false |
choose whether or not to use the 1/e method for determining the domain size from the autocorrelation data
| bool Ising_OPV::Parameters::Enable_export_compressed_files = true |
choose whether or not the output morphology data file is in compressed format
| bool Ising_OPV::Parameters::Enable_export_cross_section = false |
choose whether or not to output uncompressed data for a cross-section (x=Length/2 plane) of the morphology
| bool Ising_OPV::Parameters::Enable_extended_correlation_calc = false |
choose whether or not to extend the autocorrelation function calculation to the specified cutoff distance
| bool Ising_OPV::Parameters::Enable_growth_pref = false |
choose whether or not to implement a directional-dependent interaction energy that caues directional-dependent domain growth
| bool Ising_OPV::Parameters::Enable_import_morphologies = false |
choose whether or not the import a morphology set from the working directory
| bool Ising_OPV::Parameters::Enable_import_tomogram = false |
choose whether or not to import a tomogram dataset from the working directory
| bool Ising_OPV::Parameters::Enable_interfacial_distance_calc = false |
choose whether or not to calculate the interfacial distance histograms
| bool Ising_OPV::Parameters::Enable_interfacial_mixing = false |
choose whether or not to perform interfacial mixing
| bool Ising_OPV::Parameters::Enable_mix_frac_method = false |
choose whether or not to use the mix fraction method for determining the domain size from the autocorrelation data
| bool Ising_OPV::Parameters::Enable_periodic_z = true |
Choose whether or not to enable z-direction periodic boundary conditions.
| bool Ising_OPV::Parameters::Enable_reduced_memory_tortuosity_calc = false |
choose whether or not to perform the tortuosity calculation using an algorithm that takes longer, but uses less memory
| bool Ising_OPV::Parameters::Enable_rescale = false |
choose whether or nmot to perform lattice rescaling
| bool Ising_OPV::Parameters::Enable_shrink = false |
chose whether or not to shrink the lattice by 1/rescale_factor instead of expand it
| bool Ising_OPV::Parameters::Enable_smoothing = false |
choose whether or not to perform domain smoothing
| bool Ising_OPV::Parameters::Enable_tortuosity_calc = false |
choose whether or not to calculate the end-to-end tortuosity histograms and island volume fraction
| int Ising_OPV::Parameters::Extended_correlation_cutoff_distance = 0 |
cutoff distance for the extended autocorrelation function calculation
| int Ising_OPV::Parameters::Growth_direction = 0 |
direction that has the modified interaction energy with 1=x,2=y,3=z
| int Ising_OPV::Parameters::Height = 0 |
z-direction size of the lattice
| double Ising_OPV::Parameters::Interaction_energy1 = 0.0 |
energetic favorability for type1-type1 interactions over type1-type2 interactions
| double Ising_OPV::Parameters::Interaction_energy2 = 0.0 |
energetic favorability for type2-type2 interactions over type1-type2 interactions
| double Ising_OPV::Parameters::Interface_conc = 0.0 |
concentration of type1 sites in the interfacial mixed region
| double Ising_OPV::Parameters::Interface_width = 0.0 |
interfacial width of the mixed interface
| int Ising_OPV::Parameters::Length = 0 |
x-direction size of the lattice
| int Ising_OPV::Parameters::MC_steps = 0 |
number of MC steps to be executed (determines number of Ising swapping iterations), which sets the duration of the phase separation process
| double Ising_OPV::Parameters::Mix_fraction = 0.0 |
volume fraction of type1 sites
| double Ising_OPV::Parameters::Mixed_conc = 0.0 |
volume fraction of type1 sites in the mixed phase of the tomogram dataset
| double Ising_OPV::Parameters::Mixed_frac = 0.0 |
volume fraction of the mixed phase in the tomogram dataset
| int Ising_OPV::Parameters::N_extracted_segments = 0 |
number of equal size cuboids to extract from the tomogram dataset
| int Ising_OPV::Parameters::N_sampling_max = 0 |
maximum number of sites to be sampled for calculating the autocorrelation function
| int Ising_OPV::Parameters::N_variants = 0 |
number of random variants to create from each extracted cuboid segment
| int Ising_OPV::Parameters::Rescale_factor = 0 |
rescale factor to be used
| double Ising_OPV::Parameters::Smoothing_threshold = 0.0 |
cutoff threshold for the smoothing algorithm
| std::string Ising_OPV::Parameters::Tomogram_name = "" |
name of the tomogram dataset corresponding to the Tomogram_name.xml and Tomogram_name.raw files
| int Ising_OPV::Parameters::Width = 0 |
y-direction size of the lattice


Donald Trump’s campaign was heavy with anti-immigration rhetoric. Analysis of the election data shows that low Immigrant Population correlated with higher GOP vote share in 2016 US Presidential Election. The purpose of this visualization is to show that Donald Trump was elected by populations that have less exposure to immigrants and other foreign-born people living in the United States. To demonstrate this relationship, I use a scatterplot and choropleth map in an interactive visualization.
After making preliminary sketches, I created a few quick visualizations in tableau to get a better idea of how to lay out the various parts of the visualization. I ended up deciding on a top-down approach because it can be easily navigated by scrolling up and down or automatically moving between places on the visualization.


I had chronic insomnia for years and nothing was working, infact it was ruining my life because I had to drop out of school and had to find a job that didn't require a set schedual. Well I started stripping to pay the bills because I can stay up all night and I am reasonably attractive and even used to model. I started taking Elavil at 50mg then 100mg and finally at 150mg. I was able to sleep everyday and wake up feeling great.
Well because of this drug I gained 25lbs and lost my job, but I had begun Real Estate so I didn't care and I was back in school.  I started drinking with my pill because if I could still sleep every night I could deal with the side effects. Alcohol did work for a short time, but eventually I needed more and more to go to sleep. The thing with Elavil though was I got a different type of drunk when I mixed this pill with booze.
This drunk consisted of irrational and unhealthy behavior which I had never had. This new kind of drunk gave me a feeling of euphoria and invincibility and because of this I got arrested. My designated driver was more drunk than I was and normally I would have called a cab, but instead I wanted to have fun. I got in my car to drive my friend home and thought screw the speed limit, I can drive fine, which I could. Well I got pulled doing 60 in a 35 and didn't care or even worry. When the cop came up to my car he asked how much I had to drink and I said too much. I passed all of the field tests with flying colors but the cop still wanted to breathalize me and I blew a .19. He was shocked because my physical reactions were perfect, but my attitude was so happy and care free that the cop felt bad for arresting me.
I was never put into handcuffs, the cop took a 30 min walk with me to try to sober me up and we even began to talk about our personal lives and I got him to give me very private information of his life because I was so laid back. He never put me in jail because I was so nice and he even drove me home after the arrest to ensure that I would be safe. Well I'm charged with a .17 and a speeding ticket but I'm waiting for my court date. The cop will testify in my behalf, but I still have to get a breathalizer on my car for a year, go to counciling, pay a fine, lose my license, do community service and hopefully not be givin any jail time (possible for me up to 3 years).
I cheeted on my boyfriend but don't remember a single thing the whole night, I tried to jog around my city at 1:30 in the morning in a bad neighborhood because I could and even lost $900 at work stripping because I drank champange and let the man walk out with out paying me because I didn't know what was going on. For some people the pill may work under a doctors supervision and with out alcohol.
Experience Reports are the writings and opinions of the individual authors who submit them.
Some of the activities described are dangerous and/or illegal and none are recommended by Erowid.


Build Your Website!
Home
Build Your Website
Get the Code
Developers
Website Input
Select Element
Choose an Element
one image
two image
three image
text
left image right article
left article right image
Placeholder
Provided Image
I have a URL
Placeholder
Provided Image
I have a URL
Placeholder
Provided Image
I have a URL
Add
preview themes:
plain
light
dark
Red
Colorful
Blue
Green
Preview


| Top |
Functions
| NcXcor * | nc_xcor_new () |
| NcXcor * | nc_xcor_ref () |
| void | nc_xcor_free () |
| void | nc_xcor_clear () |
| NcXcorKinetic * | nc_xcor_kinetic_copy () |
| void | nc_xcor_kinetic_free () |
| void | nc_xcor_prepare () |
| void | nc_xcor_limber () |
Properties
| NcDistance * | distance | Read / Write / Construct Only |
| NcXcorLimberMethod | meth | Read / Write / Construct |
| NcmPowspec * | power-spec | Read / Write / Construct Only |
Description
The angular power spectrum between observables $A$ and $B$ with kernels $W^A$ and $W^B$ is \begin{equation} C_{\ell}^{AB} &= \int dz W^A(z) \int dz^\prime W^B (z^\prime) \times \int dk \frac{2}{\pi} k^2 P(k, z, z^\prime) j_{\ell}(k\chi(z)) j_{\ell} (k\chi(z^\prime)). \end{equation} In the Limber approximation, it reduces to \begin{equation} C_{\ell}^{AB} = \int_0^{z_*} dz \frac{H(z)}{c \chi^2(z)} W^A(z) W^B (z) P\left(k = \frac{\ell +1/2}{\chi(z)} , z \right), \end{equation} where $P\left(k = \frac{\ell +1/2}{\chi(z)} , z \right)$ is the power spectrum (a NcmPowspec) at redshift $z$ and $chi(z)$ the comoving distance (a NcDistance).
Functions
nc_xcor_new ()
NcXcor * nc_xcor_new (NcDistance *dist
,NcmPowspec *ps
,NcXcorLimberMethod meth
);
Two methods are available to compute Limber-approximated integrals: independent GSL numerical integration or vector integration using Sundials's CVode algorithm.
nc_xcor_limber ()
void nc_xcor_limber (NcXcor *xc
,NcXcorLimberKernel *xclk1
,NcXcorLimberKernel *xclk2
,NcHICosmo *cosmo
,guint lmin
,guint lmax
,NcmVector *vp
);
Performs the computation of the power spectrum $C_{\ell}^{AB}$ in the Limber approximation. The kernels of observables A and B are xclk1
and xclk2
. If xclk2
is NULL, the auto power spectrum is computed. The result for multipoles lmin to lmax (included) is stored in the NcmVector vp
.
Property Details
The “distance”
property
“distance” NcDistance *
This property keeps the distance object.
Owner: NcXcor
Flags: Read / Write / Construct Only
The “meth”
property
“meth” NcXcorLimberMethod
This property keeps the method enumerator to compute the integrals.
Owner: NcXcor
Flags: Read / Write / Construct
Default value: NC_XCOR_LIMBER_METHOD_GSL
The “power-spec”
property
“power-spec” NcmPowspec *
This property keeps the matter power spectrum object.
Owner: NcXcor
Flags: Read / Write / Construct Only