This package is an implementation of Ulam's method ^{1 2} (see also Galerkin projection ³) for the discretization of a stochastic operator using pure Julia. Given a set of two-dimensional, one-step trajectories

defined in a domain, the essential goal of Ulam's method is to partition the domain into a series of non-intersecting regions and construct a transition probability matrix $P$ on these regions. In UlamMethod.jl, this is accomplished in two main steps

1. The user provides a polygon containing the data, and covering of the the domain is generated by polygons according to one of several different binning algorithms (rectangles, triangles, hexagons and Voronoi cells.)

2. The number of trajectories beginning in polygon $i$ and ending in polygon $j$ is used to create the entry $P_{i, j}$ of $P$ such that the edges of the domain are handled by a stochasticization algorithm which re-injects the data either according to interior-exterior trajectories or at a user-defined location.

The polygons which form the covering and the transition probability matrix are the main outputs.

In the Julia REPL, run the following code and follow the prompts:

using Pkg
Pkg.add("UlamMethod")

Make the functions in this package available to use in your code by including the following line:

using UlamMethod

The core functionality is provided by

ulam_method(traj::UlamTrajectories, domain::UlamDomain)

where traj contains information about trajectories and domain contains information about the domain and covering. Here are 10000 random trajectories in the domain $[0, 10]^2$

using UlamMethod

n_data = 10000
x0, y0, xT, yT = 10*rand(n_data), 10*rand(n_data), 10*rand(n_data), 10*rand(n_data)
traj = UlamTrajectories(x0 = x0, y0 = y0, xT = xT, yT = yT)

We will take our domain to be the rectangular subset $[3, 5] \times [4, 8]$ and generate a covering with 40 rectangles.

xmin, xmax, ymin, ymax = 3, 5, 4, 8
domain = UlamDomain(xmin, xmax, ymin, ymax, poly_type = "rec", poly_number = 40)

Run Ulam's method.

ulam = ulam_method(traj, domain)    # the main calculation

ulam.P_closed                       # the transition matrix
pt = PolyTable(ulam.polys)          # PolyTable makes a simple list of nodes and edges
pt.nodes                            # |x|y| table of polygon nodes
pt.edges[:,3]                       # the index of the polygon that the i'th node belongs to