text:
# Newton–Cartan theory

Newton–Cartan theory (or geometrized Newtonian gravitation) is a geometrical re-formulation, as well as a generalization, of Newtonian gravity first introduced by Élie Cartan[1][2] and Kurt Friedrichs[3] and later developed by Dautcourt,[4] Dixon,[5] Dombrowski and Horneffer, Ehlers, Havas,[6] Künzle,[7] Lottermoser, Trautman,[8] and others. In this re-formulation, the structural similarities between Newton's theory and Albert Einstein's general theory of relativity are readily seen, and it has been used by Cartan and Friedrichs to give a rigorous formulation of the way in which Newtonian gravity can be seen as a specific limit of general relativity, and by Jürgen Ehlers to extend this correspondence to specific solutions of general relativity.

## Classical spacetimes

In Newton–Cartan theory, one starts with a smooth four-dimensional manifold ${\displaystyle M}$ and defines two (degenerate) metrics. A temporal metric ${\displaystyle t_{ab}}$ with signature ${\displaystyle (1,0,0,0)}$, used to assign temporal lengths to vectors on ${\displaystyle M}$ and a spatial metric ${\displaystyle h^{ab}}$ with signature ${\displaystyle (0,1,1,1)}$. One also requires that these two metrics satisfy a transversality (or "orthogonality") condition, ${\displaystyle h^{ab}t_{bc}=0}$. Thus, one defines a classical spacetime as an ordered quadruple ${\displaystyle (M,t_{ab},h^{ab},\nabla )}$, where ${\displaystyle t_{ab}}$ and ${\displaystyle h^{ab}}$ are as described, ${\displaystyle \nabla }$ is a metrics-compatible covariant derivative operator; and the metrics satisfy the orthogonality condition. One might say that a classical spacetime is the analog of a relativistic spacetime ${\displaystyle (M,g_{ab})}$, where ${\displaystyle g_{ab}}$ is a smooth Lorentzian metric on the manifold ${\displaystyle M}$.

## Geometric formulation of Poisson's equation

In Newton's theory of gravitation, Poisson's equation reads

${\displaystyle \Delta U=4\pi G\rho \,}$

where ${\displaystyle U}$ is the gravitational potential, ${\displaystyle G}$ is the gravitational constant and ${\displaystyle \rho }$ is the mass density. The weak equivalence principle motivates a geometric version of the equation of motion for a point particle in the potential ${\displaystyle U}$

${\displaystyle m_{t}\,{\ddot {\vec {x}}}=-m_{g}{\vec {\nabla }}U}$

where ${\displaystyle m_{t}}$ is the inertial mass and ${\displaystyle m_{g}}$ the gravitational mass. Since, according to the weak equivalence principle ${\displaystyle m_{t}=m_{g}}$, the according equation of motion

${\displaystyle {\ddot {\vec {x}}}=-{\vec {\nabla }}U}$

does not contain anymore a reference to the mass of the particle. Following the idea that the solution of the equation then is a property of the curvature of space, a connection is constructed so that the geodesic equation

${\displaystyle {\frac {d^{2}x^{\lambda }}{ds^{2}}}+\Gamma _{\mu \nu }^{\lambda }{\frac {dx^{\mu }}{ds}}{\frac {dx^{\nu }}{ds}}=0}$

represents the equation of motion of a point particle in the potential ${\displaystyle U}$. The resulting connection is

${\displaystyle \Gamma _{\mu \nu }^{\lambda }=\gamma ^{\lambda \rho }U_{,\rho }\Psi _{\mu }\Psi _{\nu }}$

with ${\displaystyle \Psi _{\mu }=\delta _{\mu }^{0}}$ and ${\displaystyle \gamma ^{\mu \nu }=\delta _{A}^{\mu }\delta _{B}^{\nu }\delta ^{AB}}$ (${\displaystyle A,B=1,2,3}$). The connection has been constructed in one inertial system but can be shown to be valid in any inertial system by showing the invariance of ${\displaystyle \Psi _{\mu }}$ and ${\displaystyle \gamma ^{\mu \nu }}$ under Galilei-transformations. The Riemann curvature tensor in inertial system coordinates of this connection is then given by

${\displaystyle R_{\kappa \mu \nu }^{\lambda }=2\gamma ^{\lambda \sigma }U_{,\sigma [\mu }\Psi _{\nu ]}\Psi _{\kappa }}$

where the brackets ${\displaystyle A_{[\mu \nu ]}={\frac {1}{2!}}[A_{\mu \nu }-A_{\nu \mu }]}$ mean the antisymmetric combination of the tensor ${\displaystyle A_{\mu \nu }}$. The Ricci tensor is given by

${\displaystyle R_{\kappa \nu }=\Delta U\Psi _{\kappa }\Psi _{\nu }\,}$

which leads to following geometric formulation of Poisson's equation

${\displaystyle R_{\mu \nu }=4\pi G\rho \Psi _{\mu }\Psi _{\nu }}$

More explicitly, if the roman indices i and j range over the spatial coordinates 1, 2, 3, then the connection is given by

${\displaystyle \Gamma _{00}^{i}=U_{,i}}$

the Riemann curvature tensor by

${\displaystyle R_{0j0}^{i}=-R_{00j}^{i}=U_{,ij}}$

and the Ricci tensor and Ricci scalar by

${\displaystyle R=R_{00}=\Delta U}$

where all components not listed equal zero.

Note that this formulation does not require introducing the concept of a metric: the connection alone gives all the physical information.

## Bargmann lift

It was shown that four-dimensional Newton–Cartan theory of gravitation can be reformulated as Kaluza–Klein reduction of five-dimensional Einstein gravity along a null-like direction.[9] This lifting is considered to be useful for non-relativistic holographic models.[10]

## References

1. ^ Cartan, Élie (1923), "Sur les variétés à connexion affine et la théorie de la relativité généralisée (première partie)" (PDF), Ann. Ecole Norm., 40: 325, doi:10.24033/asens.751
2. ^ Cartan, Élie (1924), "Sur les variétés à connexion affine, et la théorie de la relativité généralisée (première partie)" (PDF), Ann. Ecole Norm., 41: 1, doi:10.24033/asens.753
3. ^ Friedrichs, K. O. (1927), "Eine Invariante Formulierung des Newtonschen Gravitationsgesetzes und der Grenzüberganges vom Einsteinschen zum Newtonschen Gesetz", Mathematische Annalen, 98: 566–575, doi:10.1007/bf01451608
4. ^ Dautcourt, G. (1964), "Die Newtonische Gravitationstheorie als strenger Grenzfall der allgemeinen Relativitätstheorie", Acta Physica Polonica, 65: 637–646
5. ^ Dixon, W. G. (1975), "On the uniqueness of the Newtonian theory as a geometric theory of gravitation", Communications in Mathematical Physics, 45: 167–182, Bibcode:1975CMaPh..45..167D, doi:10.1007/bf01629247
6. ^ Havas, P. (1964), "Four-dimensional formulations of Newtonian mechanics and their relation to the special and general theory of relativity", Reviews of Modern Physics, 36: 938–965, Bibcode:1964RvMP...36..938H, doi:10.1103/revmodphys.36.938
7. ^ Künzle, H. (1976), "Covariant Newtonian limts of Lorentz space-times", General Relativity and Gravitation, 7: 445–457, Bibcode:1976GReGr...7..445K, doi:10.1007/bf00766139
8. ^ Trautman, A. (1965), Deser, Jürgen; Ford, K. W. (eds.), Foundations and current problems of general relativity, 98, Englewood Cliffs, New Jersey: Prentice-Hall, pp. 1–248
9. ^ C. Duval, G. Burdet, H. P. Künzle, and M. Perrin, Bargmann structures and Newton-Cartan theory, Phys. Rev. D 31, 1841–1853 (1985)
10. ^ Walter D. Goldberger, AdS/CFT duality for non-relativistic field theory, JHEP03(2009)069 [1]


text:
## Exercise 1

### a)

Show mathematically that the marginal distribution on return values for these three programs is the same by directly computing the probability using the rules of probability (hint: write down each possible history of random choices for each program).

flip() ? flip(.7) : flip(.1)

flip(flip() ? .7 : .1)

flip(.4)

### b)

Check your answers by sampling from the programs, 1000 times each, and plotting the results.

### c)

Write another new program with the same marginal distribution on return values that looks different from all of the above programs.

## Exercise 2

### a)

Explain why (in terms of the evaluation process) these two programs give different answers (i.e. have different distributions on return values).

var foo = flip()
display([foo, foo, foo])

var foo = function() { return flip() }
display([foo(), foo(), foo()])

### b)

Modify the second program using mem so that it has the same distribution as the first program.

### c)

Change the program in Part B so that the first two elements in the list are always the same as each other, but the third element can be different. Optional challenge: try to do this by adding only these 4 characters: x, 0, 0, and 1.

## Exercise 3

### a)

Which of these programs would be more likely to generate the following proportions for 100 values of C? Justify your response.

// Program "A"
var A = flip()
var B = flip(0.9)
var C = flip() ? A && B : A || B
display([A, B, C])

// Program "B"
var A = flip(0.9);
var B = A && flip(0.9)
var C = B && flip(0.9)
display([A, B, C])

### b)

Could the program you did not choose in Part A have also generated those return values? Explain.

## Exercise 4

In the simple medical diagnosis example, we imagined a generative process for diseases and symptoms of a single patient. In this exercise, we’ll write a version of that model that represents the diseases and symptoms of many patients.

### a)

Let’s look at just two common conditions (a cold and allergies) and just two symptoms (sneeze and fever), and let’s assume that symptoms are deterministic.

var allergies = flip(0.3)
var cold = flip(0.2)

var sneeze = cold || allergies
var fever = cold

display([sneeze, fever])

Under this model, what is the probability that the patient is sneezing? What is the probability that the patient is sneezing and has a fever?

### b)

Inspect the joint probability distributions of sneeze and fever using Infer.

Infer({method: "forward", samples: 1000}, function() {
...
return { ... }
})

### c)

If we wanted to represent the diseases of many patients we might have tried to make each disease and symptom into a function from a person to whether they have that disease, like this:

var allergies = function(person) { return flip(.3) }
var cold = function(person) { return flip(.2) }

var sneeze = function(person) { return cold(person) || allergies(person) }

display([sneeze('bob'), sneeze('alice')])

Add fever to the program above, and use Infer to inspect the probability distribution over Bob’s symptoms. Is this the same probability distribution that you computed for the single patient in Part A? If not, what can you do to fix this program to capture our intuitions correctly?

## Exercise 5

Work through the evaluation process for the bend higher-order function in this example:

var makeCoin = function(weight) {
return function() {
return flip(weight) ? 'h' : 't'
}
}
var bend = function(coin) {
return function() {
return coin() == 'h' ? makeCoin(.7)() : makeCoin(.1)()
}
}

var fairCoin = makeCoin(.5)
var bentCoin = bend(fairCoin)

### a)

Directly compute the probability distribution of the bent coin in the example.

### b)

Check your answer by using Infer.

## Exercise 6

### a)

Directly compute the probability that the geometric distribution defined by the following stochastic recursion returns the number 5. Hint: What is the default parameter for flip()?

var geometric = function() {
return flip() ? 0 : 1 + geometric()
}

### b)

Check your answer by using Infer.

## Exercise 7

### a)

Convert the following probability table to a compact WebPPL program:

A B P(A,B)
F F 0.14
F T 0.06
T F 0.4
T T 0.4

Requirement: fix the probability of A first and then define the probability of B to depend on whether A is true or not.

var a = ...
var b = ...
display([a, b])

### b)

Run your WebPPL program and use Infer to check that you get the correct distribution.

## Exercise 8

Below we’ve defined a higher-order function flipSequence that takes a coin flipping function (e.g. trickCoin, below) and flips that coin until it gets a sequence of two heads in a row (in which case it returns heads 'h') or two tails in a row (in which case it returns tails 't'). Try out different weights for the trickCoin.

var makeCoin = function(weight) {
return function() {
return flip(weight) ? 'h' : 't'
}
}
var flipSequence = function(coin) {
return function() {
var flip1 = coin()
var flip2 = coin()
if (flip1 == flip2) {
return flip1
} else {
return flipSequence(coin)()
}
}
}

var trickCoin = makeCoin(.6)

var n_samples = 10000;
viz(Infer({method: "forward", samples: n_samples}, trickCoin))
viz(Infer({method: "forward", samples: n_samples}, flipSequence(trickCoin)))

### a)

How does flipSequence change the distribution over return values (qualitatively)? Explain why requiring two flips in a row to be the same has this effect.

### b)

What would happen if a fair coin (with weight 0.5) were input to flipSequence? Explain.

## Exercise 9

Box2D is a two dimensional simulation engine for simulating rigid bodies (those with constant shape). It allows for the construction of arbitray worlds and models important physical concepts including collisions, friction, gravity, momentum, and more.

We have provided a wrapper around Box2D that allows for the easy construction of worlds. A world consists of list of shapes. Shapes are created by JavaScript objects with the following properties:

 shape “circle” or “rect” dims [width, height] for rect or [radius] for circle x x_position_as_number (distance from left) y y_position_as_number (distance from top) static boolean (does the object move or stay still?) velocity [x_velocity, y_velocity]

The variables worldWidth and worldHeight are constants representing the visible size of the simulation window.

Here’s an example with a ground and a single rectangle. Add another object to bowlingWorld and give it an initial velocity so that it knocks the original rectangle down.

var ground = {
shape: 'rect',
static: true,
dims: [worldWidth, 10],
x: worldWidth/2,
y: worldHeight
}

var rect = {
shape: 'rect',
static: false,
dims: [10, 100],
x: worldWidth/2,
y: 390
}

var bowlingWorld = [ground, rect]
physics.animate(1000, bowlingWorld)

## Exercise 10

In Example: Intuitive physics we modeled instability of a tower as the probability that the tower falls when perturbed, and we modeled “falling” as getting shorter. It would be reasonable to instead measure how much shorter the tower gets.

### a)

Below, modify the stability/instability model by writing a continuous measure, towerFallDegree. Let this measure take different values between 0 and 1. That way, your continuous measure will be numerically comparable to the discrete measure, doesTowerFall (defined here as either 0 or 1). Explain what your continuous measure represents and why it might be a good continuous measure of instability.

///fold:
var listMin = function(xs) {
if (xs.length == 1) {
return xs[0]
} else {
return Math.min(xs[0], listMin(rest(xs)))
}
}

var highestY = function (w) { listMin(map(function(obj) { return obj.y }, w)) }
var ground = {
shape: 'rect',
static: true,
dims: [worldWidth, 10],
x: worldWidth / 2,
y: worldHeight + 6
}

var almostUnstableWorld = [
ground,
{shape: 'rect', static: false, dims: [24, 22], x: 175, y: 473},
{shape: 'rect', static: false, dims: [15, 38], x: 159.97995044874122, y: 413},
{shape: 'rect', static: false, dims: [11, 35], x: 166.91912737427202, y: 340},
{shape: 'rect', static: false, dims: [11, 29], x: 177.26195677111082, y: 276},
{shape: 'rect', static: false, dims: [11, 17], x: 168.51354470809122, y: 230}
]

var noisify = function (world) {
var perturbX = function (obj) {
var noiseWidth = 10
if (obj.static) {
return obj
} else {
return _.extend({}, obj, {x: uniform(obj.x - noiseWidth, obj.x + noiseWidth)})
}
}
map(perturbX, world)
}

///

// Returns height of tower
var getTowerHeight = function(world) {
return worldHeight - highestY(world)
}

var doesTowerFall = function (initialW, finalW) {
var approxEqual = function (a, b) { Math.abs(a - b) < 1.0 }
return 1 - approxEqual(highestY(initialW), highestY(finalW))
}

var towerFallDegree = function(initialW, finalW) {
...
}

var visualizeInstabilityMeasure = function(measureFunction) {
var initialWorld = noisify(almostUnstableWorld)
var finalWorld = physics.run(1000, initialWorld)
var measureValue = measureFunction(initialWorld, finalWorld)
print("Instability measure: " + measureValue)
print("Initial height: " + getTowerHeight(initialWorld))
print("Final height: " + getTowerHeight(finalWorld))
physics.animate(1000, initialWorld)
}

// Test binary doesTowerFall measure
// visualizeInstabilityMeasure(doesTowerFall)

// Test custom towerFallDegree measure
visualizeInstabilityMeasure(towerFallDegree)

### b)

Describe a tower with a very different doesTowerFall and towerFallDegree measures look like. Which measure captures the meaning of “unstable” better?

format: plain


text:
# List of experiments contradicting the expected utility model

I proposed this question as an example question for the Area51 proposal "Mathematical modeling". User Artem Kaznatcheev suggested that it be asked here too, which I thought was a good idea.

So I am looking for a list of experiments which cannot be accounted for by the expected utility model. By the expected utility model, I mean the model of individual preferences over vectors of uncertain events (e.g. $\Big(P(rain) = 0.4, P(sunshine) = 0.6\Big)$ and $\Big(P(rain) = 0.6, P(sunshine) = 0.4\Big)$) which satisfies a list of axioms proposed by Von Neuman and Morgernstern, namely

• Completeness
• Transitivity
• Continuity
• Independence

A rigorous formulation of these axioms can be found on page 8 of Axiomatic Foundations of Expected Utility and Subjective Probability, by Edi Karni, from the Handbook of Economics of risk and uncertainty.. The violations of these axioms I am most interested in are the ones related to the Independence axiom (violations of completeness, transitivity and continuity would deserve a separate question. See this question for an example of intransitivity.).

Alternatively, by Von-Neuman and Morgenstern's representation theorem (page 9 of the same reference), these axioms are know to be equivalent to the fact that the preferences of the agent can be represented by a utility function of the form (in the discrete case):

$U(L) = \sum_{all~possible~events "e"} P(e)u(e)$

where $P(e)$ is again the probability that $e$ occurs and $u(e)$ is the utility of getting event $e$ for sure.

I am looking for situations which cannot be accounted for by the expected utility model. Some well-known examples are the the Allais and Ellsberg paradoxes (although there is still a debate regarding Ellsberg paradox). On the other hand, I do not see the Saint-Peterborough paradox as contradicting expected utility theory, because it can be accounted for by the theory if one assumes an appropriate degree of risk aversion.

I hope this question can serve as a repository of famous experiments contradicting expected utility theory, so feel free to add many.

• Interesting question. When you say 'contradicting', I presume you mean cases where people's actual behaviour systematically departs from what UI theory would predict (like the Allais paradox), rather than actual logical paradoxes in UI (like the St. Petersburg paradox)? I'm not sure if 'paradox' is the appropriate term for these descriptive violations of normative theory. You might want to look here as a starting point (I'll try and provide a proper answer later). – Eoin May 14 '14 at 14:07
• Thanks for your comment. I am not sure what you mean by "UI". "Cases where people's actual behaviour systematically departs from what expected utility theory would predict" is precisely what I am looking for though. Equivalently, I could ask for situations which violate the Von-Neuman and Morgenstern independence axiom (or continuity, but that is maybe less interesting). I guess I used paradoxes because the most famous cases are called "paradoxes". Maybe I should reframe the question and the title... – Martin Van der Linden May 14 '14 at 14:55
• Also, it would be better if you expanded your question to be more self-contained for those that might not be familiar, by for instance giving a one sentence sketch of expected utility theory (because there are different levels of generality people use). I would also avoid "Most famous paradoxes" since that is a judgement call and experimental results are only 'paradoxical' if you make the mistake of assuming expected utility theory describes reality; just 'experiments contradicting EUT' seems more neutral. – Artem Kaznatcheev May 14 '14 at 17:14
• Sorry, I was being lazy, and also amazingly sloppy; I meant EU (for Expected Utility). I see there might be a clash of terminology here between those of us from psychology and economics backgrounds: your 'positive theory' seems to be our 'descriptive model'. Also +1 on @ArtemKaznatcheev's point about how the question is phrased. – Eoin May 14 '14 at 17:54
• I completely disagree with @NickStauner, there is absolutely no way in which EUT fits with Cross Validated, and it is perfectly on topic here because psychology was central to forming behavioral economics, and early behavioral economics formed largely as a way to explain failures of EUT. – Artem Kaznatcheev May 14 '14 at 21:13


text:
vantage6
Search…
2 | Harukas
Release Notes
Algorithm Development
How to Contribute
References
Creating a New Algorithm
In this section the basic steps for creating an algorithm for horizontal partitioned data are explained.
The final code of this tutorial is published on Github. The algorithm is also published in our Docker registry: harbor2.vantage6.ai/demo/average
It is assumed that it is mathematically possible to create a federated version of the algorithm you want to use. In the following sections we create a federated algorithm to compute the average of a distributed dataset. An overview of the steps that we are going though:
1.
Mathematically decompose the model
2.
Implement and test locally
3.
Vantage6 algorithm wrapper
4.
Dockerize and push to a registry

# 📝 The mathematical problem

We want to now the average of column X from a dataset Q which contains n samples. Dataset Q is horizontally partitioned in dataset
$A = [a_1, a_2 ... a_j] = [q_1, q_2 ... q_j]$
and
$B = [b_{1}, b_{2} ... b_k] = [q_{j+1}, q_{j+2}...q_{n}]$
. The average of dataset Q is computed as:
$Q_{mean} = \frac{1}{n} \sum \limits_{i=1}^{n} {q_i} = \frac{q_1 + q_2 + ... + q_n}{n}$
Now we would like to compute
$Q_{mean}$
from dataset A and B. This could be computed as:
$Q_{mean} = \frac{(a_1+a_2+...+a_j) + (b_1+b_2+...+b_k)}{j+k} = \frac{\sum A + \sum B }{j+k}$
We both need to count the number of samples in each dataset and we need the total sum of each dataset. Then we can compute the global average of dataset A and B.
We cannot simply compute the average on each node and combine them, as this would be mathematically incorrect. This would only work in the case if dataset A and B contain the same number of samples.

# 👨💻 Implementation

In this examples we use python, however you are free to choose any language. The only requirements are: 1) It needs to be able to create HTTP-requests, and 2) and needs to be able to read and write to files.
However is you use a different language you are not able to use our wrapper. Reach us on Discord to discuss how this works.
Now that we have figured out the maths, we can translate it to an implementation. A federated algorithm consist of two parts:
1.
Central part of the algorithm which is responsible for combining the partial results from the data station. In our case that would be dividing the sum of the totals with the sum of observations.
2.
Federated part of the algorithm which is responsible for the creating the partial results. In our case this would be computing the total (=sum) and number of observations.
The central part of the algorithm can either be run on the machine of the researcher himself or in a master container which runs on a node, the latter is the preferred method.
In case the researcher runs this part himself he needs to have a proper setup to do so (i.e. python 3.5+ and the necessary dependencies). This is useful when developing new algorithms.

## 💕 Federated Part

The node that runs this part contains a CSV-file with one column numbers which we want to use to compute the global mean. We assume that this column has no NaN values.
1
import pandas
2
3
def federated_part(path, column_name="numbers"):
4
"""Compute the sum and number of observations of a column"""
5
6
# extract the column numbers from the CSV
7
8
9
# compute the sum, and count number of rows
10
local_sum = numbers.sum()
11
local_count = len(numbers)
12
13
# return the values as a dict
14
return {
15
"sum": local_sum,
16
"count": local_count
17
}
Copied!

## ❤ The central algorithm

The central algorithm receives the sums and counts from all sites and combines these to a global mean. This could be from one or more sites.
1
def central_part(node_outputs):
2
"""Combine the partial results to a global average"""
3
global_sum = 0
4
global_count = 0
5
for output in node_outputs:
6
global_sum += output["sum"]
7
global_count += output["count"]
8
9
return {"average": global_sum / global_count}
Copied!

## 🧪 Local testing

To test simple create two datasets A and B, both having a numerical column numbers. Then run the following:
1
outputs = [
2
federated_part("path/to/dataset/A"),
3
federated_part("path/to/dataset/B")
4
]
5
Q_average = central_part(outputs)["average"]
6
print(f"global average = {Q_average}.")
Copied!

# 🌯 Algorithm wrapper

A good starting point would be to use the boilerplate from our Github. This section gives background on all the steps needed to get to this boilerplate but also provides some background information.
Now that we have a federated implementation of our algorithm we need to incorporate it in the vantage6 infrastructure. The infrastructure handles the communication with the server and provides data access to the algorithm.
The algorithm consumes a file containing the input. This contains both the method name to be triggered as well as the arguments provided to the method. The algorithm has also access to a CSV file (in the future this could also be a database) on which the algorithm can run. Finally when the algorithm is finished it writes back the output to a different file.
The central part of the algorithm needs to be able to create (sub)tasks. These subtasks are responsible to execute the federated part of the algorithm. The central part of the algorithm can either be executed on the machine of the researcher or also on one of the nodes in the vantage6 network. In this example we only show the case in which one of the nodes executes the central part of the algorithm. The node provides the algorithm with a JWT token so that the central part of the algorithm has access to the server to post these subtasks.
In this example the node uses a CSV-file as database 📔. There are implementations that use traditional databases and triple stores. We expect to support these use cases better in the future.

## 📂Package Structure

The algorithm need to be structured as a package. This way the algorithm can be installed within the Docker image. The minimal file-structure would be:
1
project_folder
2
├── Dockerfile
3
├── setup.py
4
└── algorithm_pkg
5
└── __init__.py
Copied!
We also recommend adding a README.md, LICENSE and requirements.txt to the project_folder.

### setup.py

Contains the setup method to create a package from your algorithm code. Here you specify some details about your package and the dependencies it requires.
1
from os import path
2
from codecs import open
3
from setuptools import setup, find_packages
4
5
# we're using a README.md, if you do not have this in your folder, simply
6
# replace this with a string.
7
here = path.abspath(path.dirname(__file__))
8
with open(path.join(here, 'README.md'), encoding='utf-8') as f:
9
10
11
# Here you specify the meta-data of your package. The name argument is
12
# needed in some other steps.
13
setup(
14
name='v6-average-py',
15
version="1.0.0",
16
description='vantage6 average',
17
long_description=long_description,
18
long_description_content_type='text/markdown',
19
url='https://github.com/IKNL/v6-average-py',
20
packages=find_packages(),
21
python_requires='>=3.6',
22
install_requires=[
23
'vantage6-client',
24
25
# pandas, ...
26
]
27
)
Copied!
The setup.py above is sufficient in most cases. However if you want to do more advanced stuff (like adding static data, or a CLI) you can use the extra arguments from setup.

### 🐳 Dockerfile

Contains the recipe for building the Docker image. Typically you only have to change the argument PKG_NAME to the name of you package. This name should be the same as as the name you specified in the setup.py. In our case that would be v6-average-py.
1
# This specifies our base image. This base image contains some commonly used
2
# dependancies and an install from all vantage6 packages. You can specify a
3
# different image here (e.g. python:3). In that case it is important that
4
# vantage6-client is a dependancy of you project as this contains the wrapper
5
# we are using in this example.
6
FROM harbor.vantage6.ai/algorithms/algorithm-base
7
8
# Change this to the package name of your project. This needs to be the same
9
# as what you specified for the name in the setup.py.
10
ARG PKG_NAME="v6-average-py"
11
12
# This will install your algorithm into this image.
13
COPY . /app
14
RUN pip install /app
15
16
# This will run your algorithm when the Docker container is started. The
17
# wrapper takes care of the IO handling (communication between node and
18
# algorithm). You dont need to change anything here.
19
ENV PKG_NAME=${PKG_NAME} 20 CMD python -c "from vantage6.tools.docker_wrapper import docker_wrapper; docker_wrapper('${PKG_NAME}')"
Copied!

### __init__.py

This contains the code for your algorithm. It is possible to split this into multiple files, however the methods that should be available to the researcher should be in this file. You can do that by simply importing them into this file (e.g. from .average import my_nested_method)
We can distinguish two types of methods that a user can trigger:
name
description
prefix
arguments
master
Central part of the algorithm. Recieves a client as argument which provides an interface to the central server. This way the master can create tasks and collect their results.
(client, data, *args, **kwargs)
Remote procedure call
Consumes the data at the node to compute the partial.
RPC_
(data, *args, **kwargs)
The client the master method receives is a ContainerClient which is different than the client you use as a user.
Everything that is behind a return statement is send back to the central server. This should never contain any privacy sensitive information
For our average algorithm the implementation will look as follows:
1
import time
2
3
from vantage6.tools.util import info
4
5
def master(client, data, column_name):
6
"""Combine partials to global model
7
8
First we collect the parties that participate in the collaboration.
9
Then we send a task to all the parties to compute their partial (the
10
row count and the column sum). Then we wait for the results to be
11
ready. Finally when the results are ready, we combine them to a
12
global average.
13
14
Note that the master method also receives the (local) data of the
15
node. In most usecases this data argument is not used.
16
17
The client, provided in the first argument, gives an interface to
18
the central server. This is needed to create tasks (for the partial
19
results) and collect their results later on. Note that this client
20
is a different client than the client you use as a user.
21
"""
22
23
# Info messages can help you when an algorithm crashes. These info
24
# messages are stored in a log file which is send to the server when
25
# either a task finished or crashes.
26
info('Collecting participating organizations')
27
28
# Collect all organization that participate in this collaboration.
29
30
organizations = client.get_organizations_in_my_collaboration()
31
ids = [organization.get("id") for organization in organizations]
32
33
# Request all participating parties to compute their partial. This
34
# will create a new task at the central server for them to pick up.
35
# We've used a kwarg but is is also possible to use args. Although
36
# we prefer kwargs as it is clearer.
37
info('Requesting partial computation')
38
39
input_={
40
'method': 'average_partial',
41
'kwargs': {
42
'column_name': column_name
43
}
44
},
45
organization_ids=ids
46
)
47
48
# Now we need to wait untill all organizations(/nodes) finished
49
# their partial. We do this by polling the server for results. It is
50
# also possible to subscribe to a websocket channel to get status
51
52
info("Waiting for resuls")
53
54
55
56
57
info("Waiting for results")
58
time.sleep(1)
59
60
# Once we now the partials are complete, we can collect them.
61
info("Obtaining results")
62
63
64
# Now we can combine the partials to a global average.
65
global_sum = 0
66
global_count = 0
67
for result in results:
68
global_sum += result["sum"]
69
global_count += result["count"]
70
71
return {"average": global_sum / global_count}
72
73
def RPC_average_partial(data, column_name):
74
"""Compute the average partial
75
76
The data argument contains a pandas-dataframe containing the local
77
data from the node.
78
"""
79
80
# extract the column_name from the dataframe.
81
info(f'Extracting column {column_name}')
82
numbers = data[column_name]
83
84
# compute the sum, and count number of rows
85
info('Computing partials')
86
local_sum = numbers.sum()
87
local_count = len(numbers)
88
89
# return the values as a dict
90
return {
91
"sum": local_sum,
92
"count": local_count
93
}
Copied!

## 🏡 Local testing

Now that we have a vantage6 implementation of the algorithm it is time to test it. Before we run it into a vantage6 setup we can test it locally by using the ClientMockProtocol which simulates the communication with the central server.
Before we can locally test it we need to (editable) install the algorithm package so that the Mock client can use it. Simply move to the root directory of your algorithm package (with the setup.py file) and run the following:
1
pip install -e .
Copied!
Then create a script to test the algorithm:
1
from vantage6.tools.mock_client import ClientMockProtocol
2
3
# Initialize the mock server. The datasets simulate the local datasets from
4
# the node. In this case we have two parties having two different datasets:
5
# a.csv and b.csv. The module name needs to be the name of your algorithm
6
# package. This is the name you specified in setup.py, in our case that
7
# would be v6-average-py.
8
client = ClientMockProtocol(
9
datasets=["local/a.csv", "local/b.csv"],
10
module="v6-average-py"
11
)
12
13
# to inspect which organization are in your mock client, you can run the
14
# following
15
organizations = client.get_organizations_in_my_collaboration()
16
org_ids = ids = [organization["id"] for organization in organizations]
17
18
# we can either test a RPC method or the master method (which will trigger the
19
# RPC methods also). Lets start by triggering an RPC method and see if that
20
# works. Note that we do *not* specify the RPC_ prefix for the method! In this
21
# example we assume that both a.csv and b.csv contain a numerical column age.
22
23
input_={
24
'method':'average_partial',
25
'kwargs': {
26
'column_name': 'age'
27
}
28
},
29
organization_ids=org_ids
30
)
31
32
# You can directly obtain the result (we dont have to wait for nodes to
33
34
35
print(results)
36
37
# To trigger the master method you also need to supply the master-flag
38
# to the input. Also note that we only supply the task to a single organization
39
# as we only want to execute the central part of the algorithm once. The master
40
# task takes care of the distribution to the other parties.
41
42
input_={
43
'master': 1,
44
'method':'master',
45
'kwargs': {
46
'column_name': 'age'
47
}
48
},
49
organization_ids=[org_ids[0]]
50
)
51
52
print(results)
Copied!

## 🏗 Building and 🚛 Distributing

Now that we have a full tested algorithm for the vantage6 infrastructure. We need to package it so that it can be distributed to the data-stations/nodes. Algorithms are delivered in Docker images. So that's where we need the Dockerfile for. To build an image from our algorithm (make sure you have docker installed and it's running) you can run the following command from the root directory of your algorithm project.
1
docker build -t harbor2.vantage6.ai/demo/average .
Copied!
The option -t specifies the (unique) identifier used by the researcher to use this algorithm. Usually this includes the registry address (harbor2.vantage6.ai) and the project name (demo).
In case you are using docker hub as registry, you do not have to specify the registry or project as these are set by default to the Docker hub and your docker hub username.
1
docker push harbor2.vantage6.ai/demo/average
Copied!
Reach out to us on Discord if you want to use our registries (harbor.vantage6.ai and harbor2.vantage6.ai).

# 🤞 Cross-language serialization

It is possible that a vantage6 algorithm is developed in one programming language, but you would like to run the task from another language. For these kinds of usecases the python algorithm wrapper and client support cross-language serialization. By default input to the algorithms and output back to the client are serialized using pickle. However, it is possible to define a different serialization format.
Input and output serialization can be specified as follows:
1
2
3
image='harbor2.vantage6.ai/testing/v6-test-py',
4
collaboration_id=COLLABORATION_ID,
5
organization_ids=ORGANIZATION_IDS,
6
data_format='json', # Specify input format to the algorithm
7
input_={
8
'method': 'column_names',
9
'kwargs': {'data_format': 'json'}, # Specify output format
10
}
11
)
Copied!


text:
# Gram-Schmidt Process in Three Dimensions

This application illustrates the Gram-Schmidt process in three dimensions with several steps.

Identifier:
http://demonstrations.wolfram.com/GramSchmidtProcessInThreeDimensions/
Rating:
Creator(s):
Abby Brown
Cataloger:
Publisher:
Wolfram Demonstrations Project
Rights:
Abby Brown and Wolfram Research
Format Other:
This application steps users through the Gram-Schmidt process applied to three vectors in $$\mathbf{R}^3$$. The user has a large number of choices for the initial vectors. Users can elect to see the formulas as they move through each step and the vectors are updated on the graph.


text:
Infoscience

Journal article

# A Complete Family of Scaling Functions: The (α, τ)-Fractional Splines

We describe a new family of scaling functions, the (α, τ)-fractional splines, which generate valid multiresolution analyses. These functions are characterized by two real parameters: α, which controls the width of the scaling functions; and τ, which specifies their position with respect to the grid (shift parameter). This new family is complete in the sense that it is closed under convolutions and correlations. We give the explicit time and Fourier domain expressions of these fractional splines. We prove that the family is closed under generalized fractional differentiations, and, in particular, under the Hilbert transformation. We also show that the associated wavelets are able to whiten $1⁄ƒ ^{ \lambda }$ -type noise, by an adequate tuning of the spline parameters. A fast (and exact) FFT-based implementation of the fractional spline wavelet transform is already available. We show that fractional integration operators can be expressed as the composition of an analysis and a synthesis iterated filterbank.


text:
## Scaling

If the sides of a cube are doubled fromtoas shown,

the surface area of the cube increases fromto– an increase by a factor ofand

the volume of the cube increases fromto– an increase by a factor of

These factors of 2 for the length,for the surface area andfor the volume are dirct consequences of the uniform scaling factor and the dimensionsfor length,for surface area andfor volume.

Any quantity directly related to the area or volume will scale in the same way.

The rate of thermal energy loss is directly related to the area so will increase as

The mass is directly related to the volume so will increase as

The two quantities above – rate of thermal energy loss and mass – can be directly related to living objects, both across and within species. We would expect a person 1.3 times as tall to lose thermal energytimes as fast and to betimes as heavy.

In addition, physical strength is directly related to surface area (physical strength is proportional to muscle area), so increases asbut weight increases aswhich implies there is an upper mass limit at which the body cannot support it's own weight. This applies across all species.


text:
### 3098

IVIM values in healthy brain
Steren Chabert1, Jorge Verdu1,2, Gamaliel Huerta1, Cristian Montalba3, Pablo Cox4, Rodrigo Riveros4,5, Sergio Uribe3,6, Rodrigo Salas1, and Alejandro Veloz1,7

1Biomedical Engineering Department, Universidad de Valparaiso, Valparaiso, Chile, 2Universidad Politécnica de Valencia, Valencia, Spain, 3Center for Biomedical Imaging, Pontificia Universidad Católica de Chile, Santiago, Chile, 4Servicio de Imagenología, Hospital Carlos van Buren, Valparaiso, Chile, 5Facultad de Medicina, Universidad de Valparaiso, Valparaiso, Chile, 6Radiology Department, Pontificia Universidad Católica de Chile, Santiago, Chile, 7Informatics Department, Universidad Técnica Federico Santa María, Valparaíso, Chile

### Synopsis

Even though there is much interest in brain IVIM imaging, it is difficult to get a clear view from literature on which values to expect. Our purpose is to obtain healthy brain D, D* and f, to add findings and get closer to reference values. Two distributions of 16 b-values were used to acquire data on 10 volunteers, at 1.5T: one commonly found in literature and the other considered as optimal. Values obtained from the “optimal distribution” were significantly different in all cases but D in white matter. This study emphasizes the dependence of IVIM results on the acquisition scheme applied.

### Introduction

It has become easier to obtain good quality diffusion images due to the recent improvements in MR equipment. As a consequence, there has been a renewed interest towards brain IVIM images and their relation with perfusion-related information. Various studies emphasize the relevance of these images in clinical applications, such as tumor differentiation1, grading gliomas2, stroke3 among other applications. Nevertheless, a revision of current literature does not give a clear view of what values to expect in healthy tissues, as summarized in table 1. The goal of this study is to collect a set of cerebral IVIM values in healthy subjects to provide additional findings to consolidate reference values, to be obtained in a context as close as commonly found in clinical set-ups.

### Materials and Methods

This study was approved by the institutional Ethical Committee and included ten healthy subjects (7 males, 24.7 ± 6.8 y.o.). Images were acquired on a 1.5T Philips scanner. Conventional PGSE-EPI sequence was used, with TR/TE of 4000/110ms, acceleration factor of 2, FOV 230 mm, matrix size 1282. Two b-value distributions were used, one commonly found in literature {0; 10; 20; 30; 40; 50; 60; 70; 80; 90; 100; 150; 200; 400; 800; 1000} s.mm-2 and the other one, “optimal” according to Lemke et al.4 {0; 40; 50; 60; 150; 160; 170; 190; 200; 260; 440; 560; 600; 700; 980; 1000} s.mm-2. Images were processed in MATLAB. We adjusted the perfusion fraction f, pseudo-diffusion coefficient D* and diffusion coefficient D according to equation 1, where S stands for the signal magnitude of diffusion-weighted image and S0 stands for the signal magnitude without diffusion weighting.

$$S = S_0 \left [fe^{-bD^*} + (1-f)e^{-bD} \right ]$$

IVIM parameters were fitted using trust-region-reflective algorithm in two steps. Four ROIs were positioned in GM and four ROI in WM in each volunteer’s images. Analysis was then undertaken over signal average over each ROI. We considered as outlier the ROI whose values were superior to the population mean plus three standard deviations. Heteroscedastic, two-tails t-test was applied to check difference between populations.

### Results

Quality of fit was confirmed in all cases. Even though a slightly higher number of outliers were detected in literature distribution than in optimal distribution (4 vs 1), no significant difference in quality of fit was found between optimal and literature distribution (p=0.058 in case of GM and p = 0.950 in case of WM).

Averaged values for perfusion fraction, pseudo-diffusion and diffusion coefficients are summarized in table 2, and visualized as box-plot in figure 1. All values obtained with the optimal distribution were significantly different from the ones obtained with literature distribution, except for D in white matter. f values are higher with the optimal distribution and D* are lower with the optimal distribution. f and D are different (p<0.01) in WM compared to GM, in both cases of optimal and literature distributions. D* variations were such that no difference is found between GM and WM in both literature and optimal distribution schemes.

Standard deviations of IVIM parameters adjusted using data obtained with the optimal b-value distribution are lower than the ones obtained with the literature distribution. There is more variance observed in analysis in GM. Of all parameters (f, D* and D), D* is the least stable.

### Discussion and conclusion

Discussion has been active about the impact of the methodology used to estimate IVIM parameters: different fitting methods have been used, using one or two-steps for parameters fitting, constraining fit or not, about a possible dependence over b-value threshold to be used when adjusting parameters by parts, showing among other point the influence of SNR over robustness of evaluation5-7.

This work offers new elements to answer two of the current questions in IVIM imaging: which values are to be expected and how to set up the acquisition in a clinical setting. This study shows the dependence of IVIM parameters estimation to b-value distribution scheme. Our recommendation would be to choose an optimal b-value distribution such as the one proposed by Lemke et al.4, as the parameters variance is lower. Much is to be explored still in IVIM acquisition, to develop its full potential interest as biomarkers in different pathologies.

### Acknowledgements

Grant support: CONICYT - PIA - Anillo ACT1416; PMI UVA 1402

### References

1. Federau C, Meuli R, O’Brien K, Maeder P, Hagmann P. Perfusion measurement in brain gliomas with intravoxel incoherent motion MRI. Am J Neuroradiol 2014;35:256–62
2. Lin Y, Li J, Zhang Z, et al. Comparison of Intravoxel Incoherent Motion Diffusion-Weighted MR Imaging and Arterial Spin Labeling MR Imaging in Gliomas. Biomed Res Int 2015;2015:234245.
3. Suo S, Cao M, Zhu W et al. Stroke assessment with intravoxel incoherent motion diffusion-weighted MRI. NMR Biomed 2016;29:320–8.
4. Lemke A, Stieltjes B, Schad LR, Laun FB. Toward an optimal distribution of b values for intravoxel incoherent motion imaging. Magn Reson Imaging 2011;29:766–76.
5. Stieb S, Boss A, Wurnig MC et al. Non-parametric intravoxel incoherent motion analysis in patients with intracranial lesions: Test-retest reliability and correlation with arterial spin labeling. NeuroImage Clin 2016;11:780–8.
6. Meeus EM, Novak J, Withey SB, Zarinabad N, Dehghani H, Peet AC. Evaluation of intravoxel incoherent motion fitting methods in low-perfused tissue. J Magn Reson Imaging 2017;45:1325–34.
7. Wurnig MC, Donati OF, Ulbrich E, et al. Systematic analysis of the intravoxel incoherent motion threshold separating perfusion and diffusion effects: Proposal of a standardized algorithm. Magn Reson Med 2015;74:1414–22.
8. Wirestam R, Borg M, Brockstedt S, Lindgren A, Holtås S, Ståhlberg F. Perfusion-related parameters in intravoxel incoherent motion MR imaging compared with CBV and CBF measured by dynamic susceptibility-contrast MR technique. Acta Radiol 2001;42:123–8.
9. Rydhög AS, van Osch MJP, Lindgren E, et al. Intravoxel incoherent motion (IVIM) imaging at different magnetic field strengths: What is feasible? Magn Reson Imaging 2014;32:1247–58.
10. Wu WC, Chen YF, Tseng HM, Yang SC, My PC. Caveat of measuring perfusion indexes using intravoxel incoherent motion magnetic resonance imaging in the human brain. Eur Radiol 2015;25:2485–92.
11. Filli L, Wurnig MC, Luechinger R, Eberhardt C, Guggenberger R, Boss A. Whole-body intravoxel incoherent motion imaging. Eur Radiol 2015;25:2049–58.
12. Federau C, O’Brien K, Birbaumer A, Meuli R, Hagmann P, Maeder P. Functional mapping of the human visual cortex with intravoxel incoherent motion MRI. PLoS One 2015;10:1–11.
13. Ahlgren A, Knutsson L, Wirestam R, et al. Quantification of microcirculatory parameters by joint analysis of flow-compensated and non-flow-compensated intravoxel incoherent motion (IVIM) data. NMR Biomed 2016;29:640–9.
14. Wang C, Ren D, Guo Y et al. Distribution of intravoxel incoherent motion MRI-related parameters in the brain: evidence of interhemispheric asymmetry. Clinical Radiology 2017;72:94.e1-94.e6.
15. Wong SM, Zhang CE, van Bussel FCG et al. Simultaneous investigation of microvasculature and parenchyma in cerebral small vessel disease using intravoxel incoherent motion imaging. NeuroImage Clin 2017;14:216–21.
16. Finkenstaedt T, Klarhoefer M, Eberhardt C et al. The IVIM signal in the healthy cerebral gray matter: A play of spherical and non-spherical components. NeuroImage 2017;152:340–347.

### Figures

Table 1: IVIM values from literature in healthy subjects

Table 2: IVIM values adjusted in White and Gray Matter, using the two b-values distributions. Values are given in average ± standard deviation (with standard deviation expressed in % of the mean within parenthesis). Last row in each case indicates the p-value obtained from t-test between results obtained from data from the two distributions. ** indicates significant difference between results obtained using data from literature distribution or from optimal distribution, alfa 1%. N value corresponds to the number of ROI included, after elimination of outliers.

Figure 1. Box plots over ROI in all volunteers, “opt” indicates results obtained with optimal b-value distribution and “lit” indicates results obtained with literature b-value distribution. Top row: f values (%), middle row: D* values (x 10-3 mm2/s), bottom row: D values (x 10-3 mm2/s). Left column: Grey Matter; Right column: White Matter.

Proc. Intl. Soc. Mag. Reson. Med. 26 (2018)
3098


text:
# Rank of a matrix when points are on the unit sphere

Consider the following matrix:

$\begin{bmatrix} 2x & -2y & 0 \\\ 0 & z & y \\ z & 0 & x \\ y & x & 0 \\ 2x & 2y & 2z\end{bmatrix}$

How can we prove that this matrix has always rank $3$ whenever $(x,y,z)$ satisfy $x^{2}+y^{2}+z^{2}=1$? i.e they lie on the unit sphere. (without using row echelon form as it seems messy)

-
Has it? Try $x=y=0$. –  Did Oct 7 '11 at 16:32
Hmmm... the second version seems better. –  Did Oct 7 '11 at 16:35

Assuming that $z$, $y$ and $z$ are real: The determinant of rows 1, 3, and 4 is $$2\begin{vmatrix}x&-y&0\\z&0&x\\y&x&0\end{vmatrix}=2(-xy^2-x^3)=-2x(x^2+y^2)$$ which is nonzero as long as $x$ is.

Similarly, rows 1, 2, and 4 have nonzero determinant whenever $y\ne 0$.

The only case left is $x=y=0$, $z\ne 0$ ...

-


text:
# How do solve 2x^3+x^2>6x and write the answer as a inequality and interval notation?

Feb 23, 2017

The solutions are $- 2 < x < 0$ and $x > \frac{3}{2}$
Or in interval notation, x in ]-2,0[uu]3/2,+oo[

#### Explanation:

Let's rewrite and factorise the inequality

$2 {x}^{3} + {x}^{2} > 6 x$

$2 {x}^{3} + {x}^{2} - 6 x > 0$

$x \left(2 {x}^{2} + x - 6\right) > 0$

$x \left(2 x - 3\right) \left(x + 2\right) > 0$

Let $f \left(x\right) = x \left(2 x - 3\right) \left(x + 2\right)$

We can build the sign chart

$\textcolor{w h i t e}{a a a a}$$x$$\textcolor{w h i t e}{a a a a}$$- \infty$$\textcolor{w h i t e}{a a a a}$$- 2$$\textcolor{w h i t e}{a a a a a}$$0$$\textcolor{w h i t e}{a a a a a}$$\frac{3}{2}$$\textcolor{w h i t e}{a a a a a a}$$+ \infty$

$\textcolor{w h i t e}{a a a a}$$x + 2$$\textcolor{w h i t e}{a a a a a a}$$-$$\textcolor{w h i t e}{a a a a}$$+$$\textcolor{w h i t e}{a a a a}$$+$$\textcolor{w h i t e}{a a a a a}$$+$

$\textcolor{w h i t e}{a a a a}$$x$$\textcolor{w h i t e}{a a a a a a a a a}$$-$$\textcolor{w h i t e}{a a a a}$$-$$\textcolor{w h i t e}{a a a a}$$+$$\textcolor{w h i t e}{a a a a a}$$+$

$\textcolor{w h i t e}{a a a a}$$2 x - 3$$\textcolor{w h i t e}{a a a a a}$$-$$\textcolor{w h i t e}{a a a a}$$-$$\textcolor{w h i t e}{a a a a}$$-$$\textcolor{w h i t e}{a a a a a}$$+$

$\textcolor{w h i t e}{a a a a}$$f \left(x\right)$$\textcolor{w h i t e}{a a a a a a a}$$-$$\textcolor{w h i t e}{a a a a}$$+$$\textcolor{w h i t e}{a a a a}$$-$$\textcolor{w h i t e}{a a a a a}$$+$

Therefore,

$f \left(x\right) > 0$ when $- 2 < x < 0$ and $x > \frac{3}{2}$

Or in interval notation, x in ]-2,0[uu]3/2,+oo[


text:
# Talk:Dynamical billiards

## Author's Second Response to Reviews (07/21/07)

All additional suggestions by the referees (including the new referee C), are met besides the following ones. I believe that "Dynamical Billiards" is a natural name for a WEB real life encyclopedia (otherwise real billiard will pop up when one clicks "billiards"). The paper Bunimovich, DelMagno (2006) deals with semi-focusing components rather than with focusing ones. It was supposed that some links to the existing on the WEB animations of dynamical billiards will be added. This will take care of the Referee A comment that the illustrations are poor.

## Author's First Response to Reviews

First, I'd like to thank both referees for their comments and suggestions. I have followed almost all of them with a very few exceptions.

The Reviewer A has two general remarks and I do disagree with both of them. I believe that both these remarks are based on quite different view of what Scholarpedia is (or is supposed to be). First of all it is Encyclopedia. It means that it is allowed (and should be encouraged) to use "technical terms" which are introduced in details in other articles of Scholarpedia. It is not a popular review article. Moreover, I had (and still have) a feeling that the articles should be kept within four pages or at least close enough to this bound. Secondly, it is a scientific Encyclopedia for scientists. Observe also that it is also very focused: it is not an Encyclopedia of Sciences, not a Mathematical Encyclopedia. It is focused on neuroscience and general ideas, notions and tools of the theory of dynamical systems. The emphasis here is on "general," rather than on very specific, although nice and difficult, mathematical questions.

I did include hard balls Boltzmann gas as one of the major motivations to study billiards as well as references to Boltzmann-Sinai ergodic hypothesis. It is a standard way to start a review paper on billiards. In fact, I thought that Scholarpedia will include a separate article on foundations of statistical mechanics which motivated the creation and development of Ergodic Theory. I still think that it'd be a good idea to include such article.

Observe that "Dynamical Billiards" appear in Scholarpedia in "Ergodic Theory" and "Hyperbolic Dynamics." Therefore it is clear and natural that the emphasis should be on typical billiards and their generic behavior. Typical billiards have chaotic (hyperbolic) components of positive measure. Thus the mechanisms of chaos should be the main points of discussion. So it necessarily should be discussed both mechanisms of chaos and natural obstacles like astigmatism to their existence. Therefore I completely disagree with the Reviewer A comment about "rather special multidimensional constructions." Actually I am very puzzled by this comment because in fact there is (and there was) no word in my article about these constructions besides that focusing components cannot be large in d > 2 because of astigmatism.

I admire the results, techniques and the area of polygonal billiards. In fact, I feel that mathematics there is more beautiful than in hyperbolic billiards. However, polygonal billiards, especially in rational polygons are extremely nontypical billiards. Their studies explore very sophisticated mathematics and they are concerned with a very fine and very specific structure of orbits. There is no room here to properly introduce and define the needed notions. Observe that interval exchange transformations, which are more general are not even discussed in Scholarpedia. Nevertheless, I did include a section on billiards in polygons and polyhedra. I didn't include outer (or external) billiards because they are not billiards. Those are very nice dynamical systems which have no relations to billiards in Euclidean spaces which are the only one considered in this article.

I consider the first book on billiards by S. Tabachnikov to be the best existing introduction to the theory. In fact it was in my list of references but was omitted (together with Lazutkin, 1973) in typing. I included all the references suggested by both referees besides two. The paper by Sinai is already outdated and covered better by the other references. I didn't know about Katok's lecture. However after reading I didn't include it because of some factual and historical inaccuracies. I didn't include numerous references to the studies of rational polygonal billiards but put a reference to a recent review. Also, I didn't include a mirror formula (although it is my favorite tool in billiards) because it is a technical tool and it'd need to include explanations why and where it is useful. Otherwise, all the suggestions and comments are taken into account. Finally, I'd like again to thank the reviewers for their thoughtful input.

Leonid Bunimovich

## Reviewer A

As I understand this article should be written for general public, not for experts. In my opinion it should introduce the topic in common (simple) terms, provide a few most basic definitions and facts, relate the topic in question to other areas in sciences, and supply references for further reading.

In this respect, I have two general remarks:

First, the author often uses special technical terms without defining or describing them, for example "configuration space", "phase space", "orbit", "singularity", etc. I suggest the author attempts to use simpler language to improve readability. Special terms should be defined or described.

Editor-in-chief has included some links to articles describing these notions.

Second, the presentation is heavily biased and focused on the author's own contribution. His focusing billiards are discussed in detail, including his rather special multidimensional constructions. On the other hand, mechanical models of hard colliding balls (hard ball gases), which originally motivated Sinai to study chaotic billiards and which still supply the most challenging billiard-related problems, are not even mentioned. The author could have mentioned intensive studies of periodic orbits (going back to Birkhoff), dual billiards, and profound recent results by Burago et al on counting reflections near corner points (that is related to hard balls again). I also suggest to write a separate (larger) section on polygonal billiards, including those with rational angles, rather than casually mentioning them in the middle of a section on focusing billiards.

Including outstanding open problem besides Birkhoff's conjecture (say, ergodicity of hard ball gases, existence of periodic orbits in any polygon, etc.) might be a good idea, too.

A few other remarks:

1. References should include one of the best books around on the subject: "Billiards" by S. Tabachnikov (Panorama and Syntheses, 1995).

2. The paper (Lasutkin, 1971) is mentioned in the text but missing from the reference list.

3. "American Mathematical Society Press" is an inaccurate description of the publishing house. For example: "Tabachnikov... AMS, Student Mathematical Library, Vol. 30". Same goes to the books by Kozlov & Treshchev and by Chernov & Markarian. (The book by Birkhoff is referred to correctly.)

I approve the article.

## Reviewer B

My general impression of the article is positive. What follows are a several suggestions to the author.

1. I would mention the "mirror equation" of geometric optics; this equation plays a major role in the study of billiards.

2. I would mention Mather's theorem on non-existence of caustics (a very natural place for that is when the existence of caustics is discussed).

3. I believe, billiards in polygons deserve a separate new section: there is a very substantial activity on rational polygons (Veech, Masur, Smillie, Eskin, Zorich, just to mention a few names).

4. I would mention Birkhoff's periodic trajectories and, possibly, recent multi-dimensional generalizations (Farber-Tabachnikov).

5. I would mention recent progress toward proof of the Boltzmann-Sinai ergodic hypothesis.

6. I would mention recent results on the number of collisions in systems of elastic balls (Burago and his coworkers).

7. I'd mention various mechanical systems with elastic collisions that motivate the study of billiards.

8. Editing.

(i) l. 7 of "General Properties": "... is equal TO zero..."

(ii) Section "Integrable Billiards": Berger's theorem holds in all dimensions not less than three; one should refer to a different, longer, paper of his on the subject (Bull. Soc. Math. France 123 (1995), 107--116.)

(iii) Section "Chaotic Billiards" with one line content makes a strange impression; delete it.

(iv) In the references, I'd add Ya. Sinai's ICM 1990 talk "Hyperbolic billiards"; S. Tabachnikov's book "Billiards" (SMF "Panoramas et Syntheses", No. 1, 1995); references to the recent "Handbook of Dynamical Systems"; and a survey by A. Katok "Billiard table as a playground for a mathematician" in Surveys in modern mathematics, 216--242, Cambridge Univ. Press, 2005.

I am satisfied with the changes made and endorse the paper. Two items for the author to correct:

1). 1st sentence of "Billiards in polygons..." should read: "Let ... HAVE a flat ..."

2). Birkhoff's theorem asserts that for every $n$ and every $r\leq n/2$, coprime with $n$, there exist at least two $n$-periodic billiard trajectories with rotation number $r$. The current formulation is too weak.

## Reviewer C

I completely agree with the remarks made by the Reviewers A and B. I add the following obvservations:

0) I do not like the title of the article. All billiards are dynamical.

1) Dispersing billiards. It is not true that "all these properties (ergodicity, ..., exponential decay of correlations) are ensured by the mechanism of dispersing." Dispersing billiards with cusps have polynomial decay of correlations (Chernov- Markarian).

2) Focusing billiards. I do not understand the exact meaning of the following sentence: "The mechanism of defocusing... average (over the invariant measure)time of convergence."

3) Chaos and astigmatism. The mechanism of defocusing works as in demension 2, in semi-focusing billiards with cylindrical components (Bunimovich-DelMagno).

4) There are a lot of unjustified rigid assertions: "most visual", "All these properties", "there is no other mechanism of chaos", "The most general class", "The only visual and rigorously".

I suggest that the author rewrite the article, following the recommendations of all three Reviewers.

First of all (this is not related to the article), the refereeing system is confusing and I'm not sure that I'm typing my comments in the right place. I hope the author will see them anyway.

The author made many improvements, and the article is much better now. In my view, it is still written for a narrow circle of experts, rather than for a broader science community. But the author argues (in his response to the reviewers) that his article's style is in line with others in Scholarpedia, so I drop my criticism on that.

I also feel that the author could have provided better illustrations: he states that "billiards form the most visual class of Hamiltonian systems" (by the way, such categorical claims in the article appear to be too rigid), but the visual images he provides are rather poorly drawn and simplistic. One could use colors, moving visual elements, or make an interactive picture by Java code.

Lastly I have a few minor suggestions:

"reflections from" --> "reftections off (or at)" "free passes" --> "free paths" "Let ... has .." --> "Let ... have ..." "a typical billiard in a polygon" --> "billiard in a typical polygon"

elliptical and hyperbolic caustics are NOT separated, they intersect each other all over in the ellipse; it is the corresponding trajectories in the phase space that are separated.

Sinai billiards have exponential decay of correlations only if the corner points are not cusps (this condition is missing in the article).


text:
## Optimal levels for modular mod 2 representations over totally real fields

### Summary

In this paper, we study the level lowering problem for mod 2 representations of the absolute Galois group of a totally real field ${F}$. In the case ${F}=\Bbb Q$, this was done by Buzzard; here, we generalise some of Buzzard's results to higher weight and arbitrary totally real fields, using Rajaei's generalisation of Ribet's theorem and previous work of Fujiwara and the author.

### Mathematics Subject Classification

11F33, 11F80, 11F41, 11G18, 14G35

### Keywords/Phrases

Hilbert modular forms, totally real fields, Galois representations


text:
# Algorithmic trading with Julia

A detailed version of this article appeared in the Automated Trader magazine.

What makes algorithmic trading particularly challenging is that it
needs to be a polymath to do it well. It requires a unique blend of
mathematics, finance, databases, operating systems, and street
smarts. Julia makes it easier. Prototype in Julia, use the most
sophisticated algorithms to trade once in a while, or the simplest
ones to trade with high frequency, backtest at scale, and deploy in
production extracting the most out of the best hardware. The
JuliaFin suite of
packages makes this entire workflow smooth and
easy. JuliaDB makes it easy to store and query
historical data and write in-database real-time
algorithms. Miletus is a domain
specific language to value complex financial instruments without
programming the mathematics from scratch. While the
Bloomberg API
wrapper package makes it easy to query real-time data from the
Bloomberg terminal, it is equally easy to consume data from other
feeds and sockets in a variety of formats. Finally, with
JuliaRun, all of this
can be put into production and scaled with a single click. This blog
will introduce each of these packages from a trading perspective and
tie it all up with an arbitrage example.

Now, you may say that is all great, but does Julia integrate with
Excel, because otherwise there is no way my users will use it? That’s
why we will start off with a quick overview of
JuliaInXL.

# Excel integration – Deployment to front office teams

Regardless of where one sits within the financial services industry,
one product that is ubiquitous is Microsoft® Excel®. Analysts and
traders often want library functionalities to be accessible from
Excel. Julia Computing’s JuliaInXL package provides an extension to
Excel that brings the power of the Julia language and its ecosystem
into the familiar spreadsheet work
environment. JuliaInXL
allows users to launch Julia processes, load modules and execute
functions directly from the Excel UI. In Figure 01 below, we show an
example of how a user can construct new functionality in the IDE,
launch a JuliaInXL server and execute Julia functions against that
JuliaInXL server session.

Figure 01: Julia and Excel working together

# Connectors for historical data sources for tick data and aggregated data

To access historical financial data from within Julia there are
options available to either use community-developed packages or
commercially-supported options from Julia Computing. In the set of
community-developed packages, the
Quandl.jl package allows
users to access Quandl’s online financial data feeds using your own
unique authentication token for the Quandl service. Quandl.jl provides
functionality for both searching and retrieving information from
Quandl’s database which covers a broad set of categories for financial
and economic datasets. Similarly, the
YStockData.jl package
allows for retrieving historical securities pricing information from
Yahoo! Finance. (Although the Yahoo! Finance API does not give access
to current data anymore, it is still possible to download historical
data.) As part of the JuliaFin suite of offerings from Julia
Computing, the
Blpapi.jl package
implements a full wrapper interface of Bloomberg Professional’s BLPAPI
for C/C++, as well as high level implementations of the popular bdp
and bdh functions for retrieving historical daily data and intraday
tick and bar functions for retrieving tick level data (Listing 1).

# JuliaDB – Store and query historical data. Write in-database real-time algorithms.

Developing an effective trading strategy normally involves testing the
strategy using historical data before going into production using live
data. Building a robust platform for backtesting requires integration
of a number of distinct components for accessing, storing and querying
data libraries for statistical modelling, numerical optimisation and
financial analytics. Additionally, incorporating visualisation
packages can help lead to clear insights. The Julia ecosystem includes
a variety of packages that address each step in developing and
deploying a backtesting workflow.

## Efficient warehousing and querying of time series and other structured data

An effective backtesting system not only requires access to historical
financial data, but also needs means of storing, persisting and
querying that data. Financial time series datasets are often highly
structured data, but can also be composed of messy datasets that
include missing values or multiple columns of temporally mismatched
data. Handling financial time series requires data structures that can
be efficiently read, stored, queried and extracted for this type of
naturally ordered data. Algorithms need to be able to handle missing
or non-overlapping data in an effective manner.

In this section, we focus on a package that has recently been
developed and open sourced by Julia Computing, aiming to provide a
user-friendly and high-performing table interface for interacting with
column oriented data
sets. JuliaDB.jl is a
distributed columnar data table implementation that builds on a
framework for parallel, in-memory and out- of-core execution and
provides a column store data table implementation along with fast data
input/output from csv-files.

JuliaDB enables large datasets spread through numerous files across a
cluster of Julia worker processes to be ingested into a single,
distributed table data structure. A JuliaDB table separates columns
into two distinct sets wherein one set forms a sorted index and the
remaining columns are the associated data columns. The table data
structure is equivalent to an N-dimensional sparse array with an
interface that follows Julia’s array API as closely as possible, but
also allows for index values to have non-integer data types. A JuliaDB
table provides a mapping of index tuples to individual elements in one
or more data columns, essentially the individual row elements in the
data columns.

Listing 1

# Create a Bloomberg session over a particular IP
address and port number
IP = “localhost”
Port = 8194
session = createSession(IP, Port)
### Reference data request
# Required parameters of reference data request:
# ticker names
tickers =  [“IBM US Equity”, “AAPL US Equity”]
#  elds requested
elds = [“PX_Open”, “PX_High”, “PX_Last”]
# Call the bdp function by providing session,
tickers and  elds variables
Response = bdp(session, tickers,  elds)
# The response from bdp function is a Julia type
ReferenceDataResponse object
# Initializing tickers and  elds arrays to be passed
to the bdp function call
tickers = [“IBM US Equity”, “AAPL US Equity”]
elds = [“PX_Last”, “PX_Open”]
# Getting the response in variable ‘Response’
Response = bdp(Session, tickers,  elds)
# extracting response data by providing ticker and
eld name
ibmLastPrice = Response[“IBM US Equity”, “PX_Last”]
ibmOpenPrice = Response[“IBM US Equity”, “PX_Open”]
appleLastPrice = Response[“AAPL US Equity”,
“PX_Last”]
appleOpenPrice = Response[“AAPL US Equity”,
“PX_Open”]

Listing 2

               _
_       _ _(_)_     |  A fresh approach to technical computing
(_)     | (_) (_)    |  Documentation: https://docs.julialang.org
_ _   _| |_  __ _   |  Type "?help" for help.
| | | | | | |/ _  |  |
| | |_| | | | (_| |  |  Version 0.6.0-rc1.0 (2017-05-07 00:00 UTC)
_/ |\__'_|_|_|\__'_|  |  Official http://julialang.org/ release
|__/                   |  x86_64-pc-linux-gnu

julia> addprocs(20);

julia> @everywhere using JuliaDB
julia> dir = "/home/juser/TrueFX/data";
julia> files = glob("*.csv",dir)
1380-element Array{String,1}:
"/home/juser/TrueFX/data/AUDJPY-2009-05.csv"
"/home/juser/TrueFX/data/AUDJPY-2009-06.csv"
"/home/juser/TrueFX/data/AUDJPY-2009-07.csv"
"/home/juser/TrueFX/data/AUDJPY-2009-08.csv"
"/home/juser/TrueFX/data/AUDJPY-2009-09.csv"
"/home/juser/TrueFX/data/AUDJPY-2009-10.csv"
"/home/juser/TrueFX/data/AUDJPY-2009-11.csv"
"/home/juser/TrueFX/data/AUDJPY-2009-12.csv"
"/home/juser/TrueFX/data/AUDJPY-2010-01.csv"
⋮
"/home/juser/TrueFX/data/USDJPY-2016-05.csv"
"/home/juser/TrueFX/data/USDJPY-2016-06.csv"
"/home/juser/TrueFX/data/USDJPY-2016-07.csv"
"/home/juser/TrueFX/data/USDJPY-2016-08.csv"
"/home/juser/TrueFX/data/USDJPY-2016-09.csv"
"/home/juser/TrueFX/data/USDJPY-2016-10.csv"
"/home/juser/TrueFX/data/USDJPY-2016-11.csv"
"/home/juser/TrueFX/data/USDJPY-2016-12.csv"

julia> a = ingest(files, dir, colnames=["pair","timestamp","bid","ask"], indexcols=[1,2],
colparsers=Dict("pair" => TextParse.StrRange, "timestamp" => Dates.DateTime,
"bid" => Float32, "ask" => Float32))

Reading 1380 csv files totalling 237.217 GiB...
DTable with 5539994782 rows in 1380 chunks:

pair       timestamp               │ bid     ask
───────────────────────────────────┼───────────────
"AUD/JPY"  2009-05-01T00:00:00.31  │ 72.061  72.092
"AUD/JPY"  2009-05-01T00:00:00.318 │ 72.062  72.09
"AUD/JPY"  2009-05-01T00:00:00.361 │ 72.066  72.092
"AUD/JPY"  2009-05-01T00:00:00.528 │ 72.061  72.087
"AUD/JPY"  2009-05-01T00:00:00.547 │ 72.061  72.087
...

julia> a["GBP/USD", DateTime(2016,6,15):Dates.Minute(5):DateTime(2016,7,1)]
DTable with 1 chunks:

pair       timestamp           │ bid      ask
───────────────────────────────┼─────────────────
"GBP/USD"  2016-06-15T07:40:00 │ 1.41645  1.41657
"GBP/USD"  2016-06-15T10:35:00 │ 1.41995  1.4201
"GBP/USD"  2016-06-17T11:45:00 │ 1.42951  1.42968
"GBP/USD"  2016-06-20T07:40:00 │ 1.45941  1.45954
"GBP/USD"  2016-06-20T18:10:00 │ 1.46716  1.46724

The index columns in JuliaDB tables are automatically sorted
lexicographically from left to right, which allows for extremely fast
data extraction on data sets that have a natural order to their index
values, a common scenario for financial time- series data. JuliaDB
provides a set of functions for efficient selection, aggregation,
permutation and conversion of data along one or more of the index
dimensions of a given table. Whenever possible, operations on the
columns of a JuliaDB table happen in the individual Julia process
where a subset of data is already located. While JuliaDB is fully
capable of automatically resorting data between processes when
necessary, providing the system with advanced knowledge of how on-disk
data sets should be partitioned can improve performance.

Listing 2 is an example of loading JuliaDB across a cluster of 20
Julia worker processes, ingesting more than seven years worth of
foreign exchange rate data that was obtained from
TrueFX. A simple indexing operation is used to
extract a subset of the table.

In an example presented in the final section, we will walk through how
one can go about integrating JuliaDB with other packages from the
Julia ecosystem to build an efficient and high-performing backtesting
system.

# Miletus – Model, price and analyze individual as well as baskets of securities

The ability to price and hedge a variety of securities requires access
to an extensive set of financial modelling and simulation tools. These
must reflect the diverse range of possible contract terms that can be
used in defining individual securities and baskets of
securities. Mathematical models need to reflect the conditions under
which the market can operate. When speed of development and execution
are both critical, traders need access to both pre-built libraries of
commonly used contracts and models, as well as functionality for
quickly building models for new contracts under changing market
conditions.

The Julia package ecosystem provides all of the foundational
mathematical and statistical functionality necessary to build
financial models of any desired complexity. In addition it includes a
set of packages specifically focused on the construction, modelling
and time-series analytics of financial securities. The JuliaStats
organisation includes packages covering basic probability and
statistics, model fitting, Monte Carlo analysis and foundational
machine learning tools. With specific regards to time series
analytics, the
MarketTechnicals.jl,
Indicators.jl and
TimeModels.jl packages
implement a variety of algorithms commonly utilised in technical
analysis, including moving averages, momentum and volatility
indicators, rolling statistical calculations and GARCH models.

In the realm of financial contract definition, modelling and pricing,
Julia Computing has developed
Miletus.jl
as part of its
JuliaFin
offering. Miletus is a package that consists of a domain specific
language for financial contract definition inspired by the research
work of Peyton-Jones and Eber. Miletus allows for complex financial
contracts to be constructed from a combination of a few simple
primitive components and operations. When viewed through the lens of
functional programming, this basic set of primitive objects and
operations form a set of ‘combinators’ that can be used to create more
complex financial instruments, including equities, options, currencies
and bonds of various kinds.

In addition, Miletus includes a decoupled set of valuation model
routines that can be applied to combinations of contract
primitives. These combinators and valuation routines are implemented
through the use of Julia’s user-defined types, generic programming and
multiple dispatch capabilities.

Some existing implementations of financial contract modelling
environments (created in languages such as Haskell or OCaml) rely
heavily on pure functional programming for the contract definition
language, but may then switch to a second language (such as C++, Java,
APL) for implementation of valuation processes. Miletus differs in
that it leverages Julia’s strong type system and multiple dispatch
capabilities to both express these contract primitive constructs and
provide for generation of efficient valuation code. As seen elsewhere
in the Julia ecosystem, Miletus solves the two-language problem with
regards to defining and modelling financial contracts.

In Listing 3 we show an example of how to construct and value a basic
European call option in terms of the primitive constructs available in
Miletus, as well as convenient, high- level constructors.

Now we can create a few basic operations using the type constructors
provided by Miletus, shown in Listing 4. These basic primitives can
be combined into higher level operations as displayed in Listing 5.

Listing 3

 # Import the library
julia> using Miletus
import Miletus: When, Give, Receive, Buy,
Both, At, Either, Zero

Listing 4

# Receive an amount of 100 USD
julia> x=Receive(100USD)
Amount
└─100USD
# Pay is the opposite of Receive
julia> x=Pay(100USD)
Give
└─Amount
└─100USD
# Models are constructed and valued on a generic
SingleStock type
julia> s=SingleStock()
SingleStock

Listing 5

# Acquisition of a stock by paying 100USD
julia> x=Both(s, Pay(100USD))
Both
├─SingleStock
└─Give
└─Amount
└─100USD
# Which is equivalent to buying the stock
julia> x=Buy(s, 100USD)
Both
├─SingleStock
└─Give
└─Amount
└─100USD
# The notion of optionality is expressed by a choice
of two outcomes
julia> x=Either(s, Zero())
Either
├─SingleStock
└─Zero

Listing 6

julia> x=When(At(Date("2017-12-25")),
Receive(100USD))
When
├─{==}
│ ├─DateObs
│ └─2017-12-25
└─Amount
└─100USD

Listing 7

julia> x=When(At(Date("2017-12-25")),
Either(Buy(s, 100USD), Zero()))

When
├─{==}
│  ├─DateObs
│  └─2017-12-25
└─Either
├─Both
│  ├─SingleStock
│  └─Give
│      └─Amount
│          └─100USD
└─Zero

julia> eucall = EuropeanCall(Date("2017-12-25"),
SingleStock(), 100USD)
When
├─{==}
│  ├─DateObs
│  └─2017-12-25
└─Either
├─Both
│  ├─SingleStock
│  └─Give
│      └─Amount
│          └─100USD
└─Zero

Another important aspect of any contract is the notion of time. In
Listing 6 we define a temporal condition on which to receive a
payment of 100 USD. Combining that temporal condition with optionality
defines a basic European call option, which is demonstrated in Listing 7.

Listing 8

julia> gbmm = GeomBMModel(today(), 100.0USD, 0.1,
0.05, .15)
Geometric Brownian Motion Model
-------------------------------
S_0 = 100.0USD
T = 2017-03-14
Yield Constant Continuous Curve with r = 0.1,
T = 2017-03-14
Carry Constant Continuous Curve with r = 0.05,
σ = 0.15
julia> value(gbmm, eucall)
7.054679704161716USD

Listing 9

 julia> using Miletus, Gadfly, Colors
import Miletus: Both, Give, Contract, WhenAt,
value

Listing 10

julia> expirydate = Date("2017-12-25")
startdate  = Date("2017-12-1")
interestrate = 0.05
carryrate    = 0.1
Volatility = 0.15 K1 = 98.0USD
K2 = 100.0USD
K3 = 102.0USD
L = 11 # Layers in the binomial lattice
price = K1-1USD:0.1USD:K3+1USD

Listing 11

julia> function payoff_curve(c, d::Date, prices)
payoff  = [value(GeomBMModel(d, x, 0.0, 0.0,
0.0), c) for x in prices]
p = [x.val for x in payoff]
r = [x.val for x in prices]
return r, p
end

To price an option requires a valuation model. In Listing 8 we define a simple Geometric Brownian Motion
model (commonly referred to as the Black-Scholes equation) and value our European call option using that
model.

Having defined a basic call option, more complex payoffs can be
created through combinations of multiple options. Option spreads
consist of various combinations of call and put options having
different strike prices and/or expiry dates. We give examples of a few
different common option spread strategies. First, as shown in Listing 9,
let’s load the packages we will need to construct, visualise and
value our spreads.

Next, let’s define a few parameters that we will use when constructing
the different component options of our spread payoffs and valuation
model, including a range of strike prices we will use to visualise the
payoff curves (Listing 10).Then, in Listing 11, we define a function
that allows for calculating the payoff at expiry.

Listing 12

julia> function butterfly_call(expiry::Date, K1, K2, K3)
@assert K1 < K2 < K3
c1 = EuropeanCall(expiry, SingleStock(), K1)
c2 = EuropeanCall(expiry, SingleStock(), K2)
c3 = EuropeanCall(expiry, SingleStock(), K3)
Both(Both(c1,c3), Give(Both(c2,c2)))
end

julia> bfly1 = butterfly_call(expirydate, K1, K2, K3)

julia> s, p_bfly1 = payoff_curve(bfly1, expirydate, price)
call₁ = EuropeanCall(expirydate, SingleStock(), K₁)
call₂ = EuropeanCall(expirydate, SingleStock(), K₂)
call₃ = EuropeanCall(expirydate, SingleStock(), K₃)
s₁,cp₁ = payoff_curve(call₁, expirydate, price)
s₂,cp₂ = payoff_curve(call₂, expirydate, price)
s₃,cp₃ = payoff_curve(call₃, expirydate, price)
blk = colorant"black"
red = colorant"red"
grn = colorant"green"
blu = colorant"blue"
plot(layer(x=s , y=p_bfly1, Geom.line, Theme(default_color=blk, line_width=1.5mm)),
layer(x=s1,     y=cp₁, Geom.line, Theme(default_color=red, line_width=1.0mm)),
layer(x=s₃,     y=cp₃, Geom.line, Theme(default_color=grn, line_width=1.0mm)),
layer(x=s₂,   y=-2cp₂, Geom.line, Theme(default_color=blu, line_width=1.0mm)), Guide.manual_color_key("", ["Butterfly Call", "call1", "call3", "-2call2"],
["black", "red", "green", "blue"]),
Guide.title("Butterfly Call Payoff Curve at Expiry"),
Guide.xlabel("Stock Price"), Guide.ylabel("Payoff"))

Listing 13

julia> volatility = 0.2
julia> gbmm = GeomBMModel(startdate, K2,
interestrate, carryrate, volatility)
Geometric Brownian Motion Model
-------------------------------
S0 = 100.0USD
T = 2017-12-1
Yield Constant Continuous Curve with r = 0.1,
T = 2017-12-1
Carry Constant Continuous Curve with r = 0.05,
T = 2017-12-1
σ = 0.15

julia> value(gbmm, b y1)
0.40245573232657295USD

Note: Julia allows for using the full unicode character set in its
source code. Special characters can be entered in many ways depending
on the input editor, but the primary method supported on the REPL is
to use TAB-completion with Tex expressions. So, for example, typing
\pi produces π, while typing K\_1 produces K1.

A call butterfly consists of four options. It can be built by buying
two call options at the high and low strike price and selling two call
options at the central strike price.

Figure 02 shows the payoffs for both the call butterfly (in black) and
its constituent call options, produced by the code in Listing 12. The
options can be valued using a variety of models. For example, Listing
13 shows code for valuing the contract using the simple Geometric
Brownian Motion model.

## Optimal trading strategies

When developing any trading strategy, optimising with regards to
timing, price, volume, risk and other metrics is key to ensuring
profitable execution. Julia and its package ecosystem have unique
strengths in the area of mathematical optimisation that are not found
in other technical computing language. Optimisation development in
Julia is coordinated through the JuliaOpt,
JuliaNLSolvers and
JuliaDiff communities.

The JuliaOpt and JuliaNLSolvers organisations contain packages, such
as Optim.jl,
LineSearches.jl,
LsqFit.jl and
NLsolve.jl, which are
implemented as liberally-licensed, pure Julia code, as well as
integrations with other best-of- breed commercial and open-source
optimisation libraries implemented in C, C++ and Fortran. The full
list of open- source and commercial optimisation solvers integrated
with Julia can be found on the JuliaOpt home page. Most of these
interface packages have implemented integrations into JuliaOpt’s
MathProgBase.jl
abstraction layer package. This provides users with both high-level,
one-shot functions for linear and mixed-integer programming, as well
as a solver- independent, low-level interface for implementing
advanced techniques requiring the efficient solution of sequential
linear programming problems.

Figure 02: Call butterfly payoff at expiry

MathProgBase provides the solver-independent abstractions needed for
implementation of the JuMP.jl domain specific language for numerical
optimisation. JuMP is an award-
winning, domain-specific language for optimisation in Julia
implemented in the tradition of existing open source and commercial
modelling languages such as AMPL, GAMS, AIMMS, OPL and MPL. Unlike
these existing modelling languages, JuMP’s implementation as a
standard Julia package allows for easy integration of high-level
optimisation modelling within larger analytic workflows. Further down,
we will give an example of a basic FX arbitrage strategy that walks
through the use of JuMP and uses data pulled from JuliaDB.

In addition to traditional linear and nonlinear programming approaches
to mathematical optimisation, the current movement in machine learning
and deep learning is developing along a similar trajectory within
Julia. The JuliaML organisation encompasses low-level interfaces to
existing learning frameworks, such as
MXnet.jl,
TensorFlow.jl and
Knet.jl. Model abstraction
packages allow users to easily switch between frameworks, and high
level modelling tools provide concise representations of machine
learning models.

While early in its development, the
Flux.jl package is a Keras like
package that includes its own DSL for the definition of machine
learning models. It allows for switching between backends like MXNet
or TensorFlow and easily fits into a larger Julia workflow.

# Tying it all together

In this section we will show how to tie together a number of distinct
Julia packages into a workflow for determining basic arbitrage
opportunities in the FX markets. As a simple trading strategy for FX
arbitrage, we will be implementing a Julia version of an example from
Cornuejols and Tütüncü,
2006
. In
this strategy, imbalances within exchange rates between multiple
currency pairs can be exploited by setting up a linear programming
optimisation problem defined with JuMP. In this problem structure, for
a given set of exchange rates at a particular point in time, a set of
linear constraint equations is constructed that define all of the
possible relationships for exchanging US dollars, Euros, British
pounds and Japanese yen for a trader starting with one US dollar in
net assets. For this example, additional constraints are added to
limit all currency amounts to be positive and to limit the number of
dollars returned to be 10,000 USD. Our optimisation problem will be
set to maximise the number of dollars returned.

This example builds upon the TrueFX data set that we loaded into a
distributed JuliaDB table previously.To construct our optimisation
problem in Julia, we first need to load a set of packages (Listing
14). In this case we are loading JuMP along with the GLPK.jl and
GLPKMathProgInterface.jl
packages. GLPK.jl is a wrapper
for the open-source Gnu Linear Programming Kit optimisation library
and the
GLPKMathProgInterface.jl
package implements the interface to MathProgBase.jl that enables use
of this solver from JuMP.

Listing 14

julia> using JuMP, GLPK, GLPKMathProgInterface;

Listing 15

julia> function fx_arbitrage(eurusd, eurgbp, eurjpy,
gbpusd, gbpjpy, usdjpy)
usdeur = 1.0/eurusd
usdgbp = 1.0/gbpusd
gbpeur = 1.0/eurgbp
jpyusd = 1.0/usdjpy
jpyeur = 1.0/eurjpy
jpygbp = 1.0/gbpjpy

m = JuMP.Model(solver = GLPKSolverLP(
msg_lev = GLPK.MSG_ERR))

@variables m begin
de; dp; dy; ed; ep; ey; pd; pe; py; yd;
ye; yp; d
end

@objective(m, Max, d)

@constraints(m, begin
d + de + dp + dy - eurusd*ed - gbpusd*pd
- jpyusd*yd == 1.0
ed + ep + ey - usdeur*de - gbpeur*pe
- jpyeur*ye == 0.0
pd + pe + py - usdgbp*dp - eurgbp*ep
- jpygbp*yp == 0.0
yd + ye + yp - usdjpy*dy - eurjpy*ey
- gbpjpy*py == 0.0
d  <= 10000.0
de >= 0.0
dp >= 0.0
dy >= 0.0
ed >= 0.0
ep >= 0.0
ey >= 0.0
pd >= 0.0
pe >= 0.0
py >= 0.0
yd >= 0.0
ye >= 0.0
yp >= 0.0
end)

solve(m)

DE,DP,DY,D = getvalue(de), getvalue(dp),
getvalue(dy), getvalue(d)
ED,EP,EY   = getvalue(ed), getvalue(ep),
getvalue(ey)
PD,PE,PY   = getvalue(pd), getvalue(pe),
getvalue(py)
YD,YE,YP   = getvalue(yd), getvalue(ye),
getvalue(yp)

return D, DE, DP, DY, ED, EP, EY, PD, PE,
PY, YD, YE, YP
end
fx_arbitrage (generic function with 1 method)

Next, we will define our optimisation problem within a function using
JuMP. Our function accepts a set of exchange rates as input arguments
and then constructs a JuMP model using GLPK as the linear programming
solver.To fit with the assumptions made in the Cornuejols and
Tütüncü example, we will work with only the bid prices for each
currency pair and also ignore transaction costs.

Listing 16

julia> bids = compute(map(IndexedTables.pick(:bid),
a),allowoverlap=false)
DTable with 5539994778 rows in 1380 chunks:

pair       timestamp               │
───────────────────────────────────┼───────
"AUD/JPY"  2009-05-01T00:00:00.31  │ 72.061
"AUD/JPY"  2009-05-01T00:00:00.318 │ 72.062
"AUD/JPY"  2009-05-01T00:00:00.361 │ 72.066
"AUD/JPY"  2009-05-01T00:00:00.528 │ 72.061
"AUD/JPY"  2009-05-01T00:00:00.547 │ 72.061
...

julia> @everywhere function  veminutes(t)
w = trunc(t, Dates.Minute)
m = Dates.Minute(w).value % 5
w + Dates.Minute(5-m)
end

julia> bids_5m = compute(convertdim(bids, 2,
veminutes, agg = max))
DTable with 8274322 rows in 1380 chunks:

───────────────────────────────┬───────
"AUD/JPY" 2009-05-01T00:05:00 │ 72.092
"AUD/JPY" 2009-05-01T00:10:00 │ 72.26
"AUD/JPY" 2009-05-01T00:15:00 │ 72.264
"AUD/JPY" 2009-05-01T00:20:00 │ 72.151
"AUD/JPY" 2009-05-01T00:25:00 │ 72.21

Figure 03: Currency volumes traded over time

To change our program for use with a different open-source or
commercial optimisation library, one only needs to change the solver
argument used when constructing the JuMP model. Creation of an
optimisation problem with JuMP is as simple as applying a few simple
macros to the model that create new symbolic variables, define an
objective function and apply a set constraints to those symbolic
variables. Once the problem is fully defined, executing the solve
function on the model variable transforms the defined fields of the
model into inputs suitable for use by GLPK. It then executes GLPK’s
linear programming solver and returns an optimal set of currency
values that meet the requirements of our objective function and
constraints. See the JuMP documentation for more information on its
capabilities for the construction and solution of optimisation
problems.

With our optimisation problem defined in Listing 15, we will now
transform and extract a subset of the data on which to execute our
arbitrage algorithm in Listing 16. First, we extract the bid data into
its own table, then we perform a windowing operation over a
five-minute interval for each currency pair to synchronise their
timestamps.

Listing 17

julia> t_5m  = DateTime(2016,6,15):Dates.Minute(5):DateTime(2016,7,1)
2016-06-15T00:00:00:5 minutes:2016-07-01T00:00:00

julia> TS = gather(bids_5m["EUR/USD",t_5m]).index.columns[2];

julia> EURUSD = gather(bids_5m["EUR/USD",t_5m]).data;

julia> EURGBP = gather(bids_5m["EUR/GBP",t_5m]).data;

julia> EURJPY = gather(bids_5m["EUR/JPY",t_5m]).data;

julia> GBPUSD = gather(bids_5m["GBP/USD",t_5m]).data;

julia> GBPJPY = gather(bids_5m["GBP/JPY",t_5m]).data;

julia> USDJPY = gather(bids_5m["USD/JPY",t_5m]).data;

julia> D,DE,DP,DY,ED,EP,EY,PD,PE,PY,YD,YE,YP = fx_arbitrage(EURUSD[1], EURGBP[1], EURJPY[1], GBPUSD[1],
GBPJPY[1], USDJPY[1])
(10000.0, 1.5005627439278532e7, 0.0, 0.0, 0.0, 1.3384616979325693e7, 0.0, 1.063233813089466e7, 0.0, 0.0, 0.0,
0.0, 0.0)

julia> FM  = map(fx_arbitrage, EURUSD, EURGBP, EURJPY, GBPUSD, GBPJPY, USDJPY);

julia> function fx_arbitrage_arrays(A,TS)

M = length(A)
D  = Array{Float64}(M)
DE = Array{Float64}(M)
DP = Array{Float64}(M)
DY = Array{Float64}(M)
ED = Array{Float64}(M)
EP = Array{Float64}(M)
EY = Array{Float64}(M)
PD = Array{Float64}(M)
PE = Array{Float64}(M)
PY = Array{Float64}(M)
YD = Array{Float64}(M)
YE = Array{Float64}(M)
YP = Array{Float64}(M)
for m = 1:M
D[m],DE[m],DP[m],DY[m],ED[m],EP[m],EY[m],PD[m],PE[m],PY[m],YD[m],YE[m], YP[m]  = (A[m]...)
end
IndexedTable(Columns(timestamp = TS),
end
Columns(DE=DE,DP=DP,DY=DY,ED=ED,EP=EP,PD=PD,PE=PE,PY=PY,YD=YD,YE=YE,YP=YP))
fx_arbitrage_arrays (generic function with 1 method)

julia> FM_arbitrage  = fx_arbitrage_arrays(FM, TS)

Next, we index into the table for each currency pair over a two-week
subset of five-minute intervals and extract the resulting data columns
into their own arrays.The rate values in each of these arrays can be
passed individually into a single execution of our arbitrage
optimisation routine, or the routine can be mapped over every element
of these arrays. Finally, as shown in Listing 17, we execute a
separate routine that constructs a new IndexedTable for display of the
currencies exchanged as part of each optimisation at every timestamp
and plot the currency volumes traded over time for each trade in
Figure 03.

# Summing up – Why Julia for algorithmic trading

• Julia provides a high productivity and high performance language for both financial analytics and infrastructure.
• Julia solves the two-language problem for both development and deployment of automated trading systems.
• Julia allows easy development of libraries both within the Julia package ecosystem and by using other language functionality from the outside.
• Julia’s package ecosystem covers the entire analytics and backtesting workflow.

We have shown how to create financial models in Julia and the
performance of their implementation. Julia’s approach to technical
computing has excited hundreds of thousands of quantitative
programmers worldwide and we hope you will join their ranks.


text:
# Bibliography

Export 6 results:
Author Title Type [ Year]
Filters: Author is Somesh Jha  [Clear All Filters]
2009
. AS-awareness in Tor path selection. In: Proceedings of the 2009 ACM Conference on Computer and Communications Security, CCS 2009, Chicago, Illinois, USA, November 9-13, 2009. Proceedings of the 2009 ACM Conference on Computer and Communications Security, CCS 2009, Chicago, Illinois, USA, November 9-13, 2009. ACM; 2009. pp. 380–389. Available from: http://portal.acm.org/citation.cfm?id=1653662.1653708
. The bayesian traffic analysis of mix networks. In: Proceedings of the 2009 ACM Conference on Computer and Communications Security, CCS 2009, Chicago, Illinois, USA, November 9-13, 2009. Proceedings of the 2009 ACM Conference on Computer and Communications Security, CCS 2009, Chicago, Illinois, USA, November 9-13, 2009. ACM; 2009. pp. 369–379. Available from: http://portal.acm.org/citation.cfm?id=1653662.1653707
. Membership-concealing overlay networks. In: Proceedings of the 2009 ACM Conference on Computer and Communications Security, CCS 2009, Chicago, Illinois, USA, November 9-13, 2009. Proceedings of the 2009 ACM Conference on Computer and Communications Security, CCS 2009, Chicago, Illinois, USA, November 9-13, 2009. ACM; 2009. pp. 390–399. Available from: http://portal.acm.org/citation.cfm?id=1653662.1653709
2008
. Dependent Link Padding Algorithms for Low Latency Anonymity Systems. In: {Proceedings of the 15th ACM Conference on Computer and Communications Security (CCS 2008)}. {Proceedings of the 15th ACM Conference on Computer and Communications Security (CCS 2008)}. Alexandria, Virginia, USA: ACM Press; 2008. pp. 323–332. Available from: http://portal.acm.org/citation.cfm?id=1455812
. Information Leaks in Structured Peer-to-peer Anonymous Communication Systems. In: Proceedings of the 15th ACM Conference on Computer and Communications Security (CCS 2008). Proceedings of the 15th ACM Conference on Computer and Communications Security (CCS 2008). Alexandria, Virginia, USA: ACM Press; 2008. pp. 267–278. Available from: http://portal.acm.org/citation.cfm?id=1455805
. PEREA: Towards Practical TTP-Free Revocation in Anonymous Authentication. In: {Proceedings of the 15th ACM Conference on Computer and Communications Security (CCS 2008)}. {Proceedings of the 15th ACM Conference on Computer and Communications Security (CCS 2008)}. Alexandria, Virginia, USA: ACM Press; 2008. pp. 333–345. Available from: http://portal.acm.org/citation.cfm?id=1455813