Julia implementation of GFI

First part - the polyhedra sampler

Stéphane Laurent

Goal

Implement Ciszewski & Hannig's sampler of the fiducial distribution for normal linear mixed models

Why Julia ?

The algorithm is computationnaly intensive
It requires a high numerical precision; we hope the BigFloat type in Julia will achieve this precision
Currently, the available Matlab implementation is not sufficient for large datasets

Why these slides ?

I'm not a Julia specialist; these slides should firstly help me to request some help
Ciszewski & Hannig's paper is not easy to read for non-mathematicians, and the algorithm may appear complicated to them
Because I like the slidify package

Problem addressed in these slides

These slides only address one part of the algorithm: the sampling of random polyhedra in the Euclidean space
This is the point requiring high numerical precision, because the polyhedra are sequentially sampled and become smaller and smaller and smaller...

Polyhedra construction: overview

Data: Some pairs of points are given on the $y$-axis
Sampling: Some ribbons issued from these points are sampled at random
Computation: The polyhedron at the intersection of the ribbons

Polyhedra construction: algorithm

Step 1: the slopes of the first two pairs of lines are sampled without restriction
Step 2: the polyhedron is computed
Step 3: the slope of the next pair is sampled in a constrained range, assuring there's an intersection
Step 4: the polyhedron is updated
Repeat steps 3 and 4

Line and polyhedron representation

A line has a ''type'': upper or lower
A line has $0$, $1$ or $2$ intersections with the active polyhedron

Only lines having two intersections are kept
A polyhedron is represented by its set of vertices, but it is more convienently handled by ordering its vertices

Line and polyhedron in Julia

I will use these types to handle lines and ribbons:

type Line
        a::Float64   # intercept
        b::BigFloat  # slope
        typ::Bool    # type of the line (true:upper, false:lower)
end

type Ribbon
        aLow::Float64
        aUpp::Float64
        b::BigFloat
end

The intercepts are provided as data by the user, hence they are Float64. The slopes are generated by the algorithm and are treated as BigFloat.

A polyhedron (or particle) will be treated as a BigFloat array, and the following function is used to order the vertices:

# order poly
function orderPart(poly::Array{BigFloat,2})
        # compute an interior point
        O = [ mean(poly[1,1:3]) , mean(poly[2,1:3]) ]
        # center the polyhedron around O
        cpoly = poly .- O
        # compute the angular parts of polar coordinates
        angles = atan2(cpoly[2,:], cpoly[1,:])
        # find the order
        ord  = sortperm(vec(angles))
        return poly[:,ord]
end

The first particle

In reality the first particle is sampled at random, but we define a fixed particle for illustration.

# extract lower/upper line from ribbon
function Dlow(ribbon::Ribbon)
        Line(ribbon.aLow, ribbon.b, false)
end
function Dupp(ribbon::Ribbon)
        Line(ribbon.aUpp, ribbon.b, true)
end

# returns the intersection of two lines given by (intercept, slope)
function intersect(D1, D2)
        x = (D1[1]-D2[1])/(D2[2]-D1[2])
        return [x, D1[1] + D1[2]*x]
end

# returns the polyhedron intersection of two ribbons
function ipart(R1::Ribbon, R2::Ribbon)
    D11 = Dlow(R1); D12 = Dupp(R1); D21 = Dlow(R2); D22 = Dupp(R2)
    A = intersect((D11.a, D11.b), (D21.a, D21.b));
    B = intersect((D11.a, D11.b), (D22.a, D22.b));
    C = intersect((D12.a, D12.b), (D21.a, D21.b));
    D = intersect((D12.a, D12.b), (D22.a, D22.b));
    return orderPart(hcat(A,B,C,D))
end

Example:

# first ribbon: 
R1 = Ribbon(0.4, 1.5, BigFloat("1.5"));
# second ribbon:
R2 = Ribbon(4.5, 5.9, BigFloat("-2"));
# find the intersection polyhedron:
poly0 = poly = ipart(R1, R2)

julia> float64(poly)
2x4 Array{Float64,2}:
 1.17143  1.57143  1.25714  0.857143
 2.15714  2.75714  3.38571  2.78571

Plotting a Javascript particle with Gadfly

using Gadfly
function plotPart(poly::Array{BigFloat,2})
        poly = orderPart(poly)
        x = float64(poly[1,[1:size(poly,2);1]])
        y = float64(poly[2,[1:size(poly,2);1]])
        p = plot(x = x, y = y, Geom.point, Geom.line(preserve_order=true))
        return p
end

p = plotPart(poly)
draw(D3("part01.js", 630px, 340px), p)

Don't forget to play with this graphic ! (zoom in/out and move it by maintaining the click of the mouse)

Computing the range and intersection

Recall the two steps, when a particle and a new pair of starting points on the $y$-axis is given:

Calculate the range of the possible slope of the new ribbon
Generate the new ribbon and compute the new particle

plot of chunk xxx

Computing the range

Two different situations are distinguished for the range calculation:

plot of chunk twosituations

The first situation is easier to handle. We will restrict to this situation in these slides.

The range in the simple situation

Denote by $P$ the current particle and by $\color{red}{\lbrace a^-, a^+\rbrace}$ the new pair of points on the $y$-axis. Then the possible range for the slope of the new ribbon is the interval $(m,M)$ where $$ m = \underset{(x,y) \in P}{\min} \left\lbrace\frac{y - a^-}{x}, \frac{y - a^+}{x}\right\rbrace \quad \text{and} \quad M = \max_{(x,y) \in P} \left\lbrace\frac{y - a^-}{x}, \frac{y - a^+}{x}\right\rbrace $$

function findRange(poly::Array{BigFloat,2}, lower::Float64, 
            upper::Float64)
        if ( (minimum(poly[1,:]) > 0) || (maximum(poly[1,:])<0) )
            slopes = [(poly[2,:]-lower)./poly[1,:] (poly[2,:]-upper)./poly[1,:]]
            return [ minimum(slopes) ; maximum(slopes) ]
        else # the y-axis cuts the paricle
            .... not shown in theses slides

It remains to write a function calculating the new particle once the new ribbon is sampled.

plot of chunk xxx

The intersection (1/6)

# new ribbon
a_low = 2.;
a_upp = 3.;
b = BigFloat("0.5");
R3 = Ribbon(a_low, a_upp, b)

plot of chunk plot_newribbon_01

We firstly calculate the intersection with the lower line, and then with the upper line

For each edge ($1$, $2$, $3$, $4$) of the current particle, here are the number of vertices above the new lower line:

D = Dlow(R3) # the lower line of R3
test1 = vec(poly[2,:]) .> D.a .+ D.b .* vec(poly[1,:])
test2 = test1[[2:size(opoly)[2]; 1]]
test = test1 + test2

julia> test
4-element Array{Int64,1}:
 0
 1
 2
 1

The intersection (2/6)

Given these results:

julia> test
4-element Array{Int64,1}:
 0
 1
 2
 1

plot of chunk plot_newribbon_01_copy

$\implies$ we know that:

The first edge has to be removed
For the third edge, there's nothing to do.
For the second and fourth edges, we calculate the intersection.

julia> Dinters = find(test.== 1) # should be 0 or 2 elements
2-element Array{Int64,1}:
 2
 4

The intersection (3/6)

There are four situations for the intersection:

updatePoly1() handles the case of two points (one edge) outside
updatePoly2() handles the case of one point outside
updatePoly() handles the general case:
- it determines the situation
- it runs updatePoly1() or updatePoly2() in case 1 or 2
- deletes some vertices and runs updatePoly1() in case of $>2$ points outside

The intersection (4/6)

We use the following function to get the $i$-th edge of the ordered particle:

# converts an edge to (intercept, slope)
function getLine(poly::Array{BigFloat,2}, index::Int)
        A = poly[:,index]
        B = poly[:,rem1(index+1,size(poly,2))]
        slope = (B[2]-A[2])/(B[1]-A[1])
        intercept = A[2] - slope*A[1]
        return (intercept, slope)
end

Intersection in case of two points (one edge) outside:

function updatePoly1(opoly::Array{BigFloat,2}, D::Line, toRemove::Int)
        opoly = deepcopy(opoly) # see https://groups.google.com/d/topic/julia-users/PfTZhZu6OMo/discussion
        # first edge
        index = if toRemove==1 size(opoly)[2] else toRemove-1 end
        M = intersect((D.a,D.b), getLine(opoly,index))
        # second edge
        index = if toRemove==size(opoly)[2] 1 else toRemove+1 end
        N = intersect((D.a,D.b), getLine(opoly,index))
        #
        opoly[:,toRemove] = M
        opoly[:,index] = N
        #
        return opoly
end

The intersection (5/6)

Intersection in case of one point outside:

function updatePoly2(opoly::Array{BigFloat,2}, D::Line, Dinters::Array{Int64,1}, test1::BitArray{1})
    # shift to put the two edges at first positions
    ncol = size(opoly)[2]
    if Dinters[2]-Dinters[1] != 1
            arrange = [ncol, [1:ncol-1]]
    else
            arrange = [rem1(i+Dinters[1]-1,ncol)::Int for i=1:ncol]  
    end
    opoly = opoly[:, arrange]
    # intersections
    M = intersect((D.a,D.b), getLine(opoly,1))
    N = intersect((D.a,D.b), getLine(opoly,2))
    #
    test = test1[arrange][2]
    if( (!D.typ && !test) || (D.typ && test) )
        return hcat(opoly[:,1], M, N, opoly[:, [3:ncol]])
    else
        return hcat(M, opoly[:,2], N)
    end
end

The intersection (6/6) - main function

function updatePoly(opoly::Array{BigFloat,2}, D::Line)
            test1 = vec(opoly[2,:]) .> D.a .+ D.b .* vec(opoly[1,:])
            test2 = test1[[2:size(opoly)[2]; 1]]
            test = test1 + test2
            if(D.typ==false)
                Remove = test .== 0
            else
                Remove = test .== 2
            end
            toRemove = find(Remove)
            if length(toRemove) == 1
                    return updatePoly1(opoly, D, toRemove[1])
            elseif length(toRemove) == 0
                    Dinters=find(test.== 1)
                    if length(Dinters) == 2
                        return updatePoly2(opoly, D, Dinters, test1)
                    else
                        return opoly
                    end
            else
                    if Remove[1] && last(Remove)
                        indices = find(!Remove)
                        torem =  size(indices)[1]+1
                        indices = [indices, last(indices)+1]
                    else
                        indices = [1:size(opoly)[2]]
                        indices = deleteat!(indices, toRemove[2]:last(toRemove))
                        torem = toRemove[1]
                    end
                    return updatePoly1(opoly[:,indices], D, torem)
            end
end

updatePoly() takes an ordered particle and returns an ordered particle - thus we only have to apply orderPart() to the initial particle