The PSciCo Project: Parallel
Scientific Computing in ML
Robert Harper
Carnegie Mellon University
August, 2000
Acknowledgements
• This talk represents collaboration with
Guy Blelloch (CS), Gary Miller (CS), and
Noel Walkington (Math).
• Much of the work is done by our
students: Umut Acar, Hal Burch,
Franklin Chen, Perry Cheng, David
Garmire, Aleks Nanevski, Leaf Petersen,
Seth Porter, and Chris Stone.
10/3/2015
Robert Harper - PSciCo Project
2
Project Goal
Use high-level languages for scientific
computing to achieve ...
• Smooth integration of production and
prototyping.
• Richer data structures and more complex
algorithms.
• Free exchange of re-usable components.
• Expressiveness with acceptable efficiency.
10/3/2015
Robert Harper - PSciCo Project
3
Project Thesis
We may achieve these goals by exploiting
• Functional programming.
• Persistent data structures.
• Higher-order functions, especially staging
(currying) and parameterization.
• Implicit coarse-grained parallelism.
• Parallel operations on aggregates.
• Modular programming.
• Interfaces, parameterization, hierarchy.
10/3/2015
Robert Harper - PSciCo Project
4
How We Can Win
• Use methods that are “unthinkable” in
current languages.
• persistence, parallelism, higher-order
functions, modularity
• Change the terms of the debate.
• Efficiency matters, but it’s not everything.
• Expressiveness, concision, robustness,
flexibility matter too.
10/3/2015
Robert Harper - PSciCo Project
5
How We Can Win
• Eliminate the distinction between
• Prototyping languages, for experimenting
with new ideas.
• Encourage “quick and dirty” hacks.
• Poor software engineering properties.
• Production languages, for building systems
“for real”.
• Sacrifice ease of use and maintainability for the
sake of efficiency.
10/3/2015
Robert Harper - PSciCo Project
6
The State of the Art
• Prototyping languages
•
•
•
•
•
•
10/3/2015
eg, Mathematica, MatLab
Convenient, especially interactively
Support numeric and symbolic computation
High-level, language-oriented
Packages / toolkits for various domains
Easy to build an application
Robert Harper - PSciCo Project
7
The State of the Art
• Production languages
•
•
•
•
•
•
10/3/2015
eg, Fortran, C, C++
Batch-oriented, not interactive
Low-level, machine-oriented
Poor support for symbolic computation
Large collection of libraries, packages
More effort to build a system
Robert Harper - PSciCo Project
8
Starting Point: ML + NESL
• ML
• Excellent support for complex data
structures and symbolic computing
• Excellent support for modularity
• Poor support for numeric computation
• Poor support for parallelism
10/3/2015
Robert Harper - PSciCo Project
9
Starting Point: ML + NESL
• NESL
• First-order functional language.
• Excellent support for implicit parallelism.
• “scan vector” model of data parallelism
• No support for complex data stuctures or
symbolic computing.
• No support for modularity.
10/3/2015
Robert Harper - PSciCo Project
10
Current Work
• Shared libraries.
• Adaptive precision arithmetic.
• General purpose geometry and physics.
• Triangulated surfaces (simplicial complexes).
• Applications.
•
•
•
•
•
10/3/2015
Convex hull.
Delaunay triangulation.
n-Body problem.
Terrain modelling.
Finite element method.
Robert Harper - PSciCo Project
11
Current Work
• Language technology.
• TILT compiler for Standard ML.
• typed intermediate languages for efficient data
representation
• Parallel run-time system with space- and
time-efficient thread scheduling.
• exploited by data parallel sequence primitives.
• Parallel, concurrent, real-time copying GC.
• Multiple mutators and collectors on an SMP
10/3/2015
Robert Harper - PSciCo Project
12
Shared Libraries
• Geometry primitives
• Points, vectors, figures.
• Line-side test, in-circle test.
• Arithmetic.
• adaptive precision floating point
• interval arithmetic
• Simplicial complexes.
• Persistent representation of surfaces.
10/3/2015
Robert Harper - PSciCo Project
13
Representing Geometry
• Single dimension-independent interface:
signature GEOMETRY = ...
• Several dimension-dependent
implementations:
structure Geo2D :> GEOMETRY = ...
structure Geo3D :> GEOMETRY = ...
• Different types for different dimensions!
• avoids run-time checking for conformance
10/3/2015
Robert Harper - PSciCo Project
14
Representing Geometry
signature GEOMETRY = sig
structure Vec : VEC
structure Point : POINT
structure HalfSpace : HALF_SPACE
structure Box : BOX
sharing Point.Vec = Vec
and HalfSpace.Point = Point
end
10/3/2015
Robert Harper - PSciCo Project
15
Representing Geometry
signature VEC = sig
structure Scalar : NUMBER
type t
val + : t * t -> t
val scale : t * Scalar.t -> t
val dot : t * t -> Scalar.t
val norm : t -> Scalar.t
end
10/3/2015
Robert Harper - PSciCo Project
16
Representing Geometry
signature POINT = sig
structure Vec : VEC
structure Number : NUMBER
type t
val origin : t
val translate : t * Vec.t -> t
val - : t * t -> Vec.t
end
10/3/2015
Robert Harper - PSciCo Project
17
Representing Geometry
signature HALF_SPACE = sig
exception Degenerate
structure Pt : POINT
structure Vec : VEC
type t
val from_pts :
Pt.t Seq.t -> t
val side_test :
t -> Pt.t -> Side.t (* 3-valued *)
end
10/3/2015
Robert Harper - PSciCo Project
18
An Object-Oriented Approach
• CGAL uses OOP techniques.
• No common interface for any dimension.
• One implementation for all dimensions.
Point translate (Point p, Vector v)
• Requires run-time conformance checks.
• Uses templates for numeric precision (float,
double, etc)
10/3/2015
Robert Harper - PSciCo Project
19
ML vs OO
• In ML conformance is ensured by sharing
specifications.
• Vectors and points must agree on common types.
• Run-time checks are avoided.
• ML supports “coherent integration” of
abstractions better than do C++ or Java.
• Geometry library is a good example of this.
10/3/2015
Robert Harper - PSciCo Project
20
Arithmetic
• Key questions in geometric algorithms:
• Side test: is a point above or below a given
line or plane?
• Circle test: is a point within a circle or
sphere?
• These functions convert a floating point
number to a boolean.
• Check sign of a determinant.
10/3/2015
Robert Harper - PSciCo Project
21
Arithmetic
• It is essential that these questions be
answered consistently and accurately.
• Otherwise higher-level results are
nonsense (eg, “convex hull” not convex).
• It is essential that these be
implemented efficiently.
• Lie at the heart of just about any algorithm
in CG.
10/3/2015
Robert Harper - PSciCo Project
22
Arithmetic
• Adaptive precision via error analysis (Priest,
Shewchuk).
• Calculate using single precision.
• Calculate “fuzz factor” by forward error analysis.
• If reliable, answer, otherwise re-do in double
precision, then arbitrary precision.
• Key weakness: must derive and implement the
error calculation; not a general arithmetic
package.
10/3/2015
Robert Harper - PSciCo Project
23
Arithmetic
• “Lazy” interval arithmetic.
• Use a general interval arithmetic package
for all arithmetic.
• IA determines reliability of result (eg, for
positivity test, interval cannot span origin).
• Carry a symbolic representation of the
“history” of the computation of the point.
• Re-do using full precision if unreliable.
10/3/2015
Robert Harper - PSciCo Project
24
Arithmetic
• Intervals work well for polynomials.
• Enough to do CG tests, assuming “virgin data”.
• Polynomials are not enough, in general!
• eg, projections onto a sphere require square roots
• How much can we feasibly support?
• algebraic numbers?
• transcendentals?
10/3/2015
Robert Harper - PSciCo Project
25
Real Numbers
• Why not represent all (computable) real
numbers?
• eg, Edalat and Potts representation as infinite
sequences of convergent intervals
• Basic limitation: equality and ordering of reals
is not decidable.
• eg, cannot decide whether x>0 in general
• most algorithms assume a decidable equality test
10/3/2015
Robert Harper - PSciCo Project
26
Real Numbers
• CG algorithms assume that geometric objects
are “discrete”.
• eg, line-side and in-circle tests
• Computationally, this is nonsense!
• cannot decide in general if x lies on L or within C
• Edalat suggests to treat geometric objects
themselves as sequences of intervals.
• successive refinement to interior and exterior
• But we must re-work CG in this framework.
10/3/2015
Robert Harper - PSciCo Project
27
Representing Surfaces
• Triangulated surfaces arise frequently.
• Finite element method for solving
boundary-value problems.
• Computation of convex hull of a set of
points.
• Terrain modeling from altitude data.
• Putting clothes on the Toy Story
characters.
10/3/2015
Robert Harper - PSciCo Project
28
Representing Surfaces
• A triangulated surface is a simplicial complex.
• An n-complex is a set of n-simplexes that “fit
together” well.
• Intersection of two n-simplexes is an (n-1)simplex.
• A 1-simplex is a point, a 2-simplex a line, a 3simplex a triangle, etc.
• Vertices are affinely independent when embedded
in n-space.
10/3/2015
Robert Harper - PSciCo Project
29
Representing Surfaces
signature SIMPLEX = sig
structure Vertex : VERTEX
type simplex
val compare : simplex * simplex -> order
val dim : simplex -> int
val vertices : simplex -> Vertex.vertex seq
val simplex : Vertex.vertex seq -> simplex
val
val
val
val
end
10/3/2015
down : simplex -> Vertex.vertex * simplex
join : Vertex.vertex * simplex -> simplex
faces : simplex -> simplex seq
flip : simplex -> simplex
Robert Harper - PSciCo Project
30
Representing Surfaces
signature VERTEX = sig
type vertex
val compare : vertex * vertex -> order
type state
type point
val new : state * point -> state * vertex
val loc : vertex -> point
end
10/3/2015
Robert Harper - PSciCo Project
31
Representing Surfaces
signature SIMPCOMP = sig
structure Simplex : SIMPLEX
type ‘a complex
val dim : int
val empty : ‘a complex
val isempty : ‘a complex -> bool
val simplices :
‘a complex -> int -> Simplex.simplex seq
val data :
‘a complex -> Simplex.simplex -> ‘a option
...
end
10/3/2015
Robert Harper - PSciCo Project
32
Representing Surfaces
signature SIMPCOMP = sig
...
val grep :
‘a complex -> int * Simplex.simplex ->
Simplex.simplex seq
val add :
‘a complex -> Simplex.simplex * ‘a ->
‘a complex
...
end
10/3/2015
Robert Harper - PSciCo Project
33
Representing Surfaces
• Signature of complexes is dimensionindependent.
• Allows dimension-independent code to
manipulate surfaces.
• Implementations are dimensiondependent.
• Different dimensions have different types.
10/3/2015
Robert Harper - PSciCo Project
34
Implementing Surfaces
• Usual implementation is imperative.
• Doubly-linked edge list.
• Constant-time traversal operations.
• Ephemeral, not easily parallelizable.
• Our implementation is functional.
• Simplex is a sequence of vertices.
• Logarithmic-time traversal of complex.
• Persistent, easily parallelizable.
10/3/2015
Robert Harper - PSciCo Project
35
Cost of Persistence?
• Lower bound on 3D hull: (n lg n).
• Matching upper bounds are known.
• eg, Quickhull algorithm.
• But they rely on constant-time surface
operations!
• eg, find all three neighbors of a triangle
• Can we match this in the persistent
case?
10/3/2015
Robert Harper - PSciCo Project
36
Cost of Persistence?
• Naively, we’d expect O(n lg2 n) to pay for
persistence.
• Surface operations are O(lg n) instead of O(1).
• Burch and Miller discovered an optimal
randomized algorithm for the persistent case!
• Touches surface O(n) times, at a cost of O(n lg n).
• O(n lg n) floating point operations.
• Plus it is persistent and parallelizable!
10/3/2015
Robert Harper - PSciCo Project
37
Delaunay Triangulation
• Generate a triangulation of a region
with boundary (in 2-space).
• Must include the boundary segments.
• Configuration of triangles must be “good”
in several senses (we’ll discuss one).
• A triangulation is a 2-complex.
• Use simplicial complex package to
represent it.
10/3/2015
Robert Harper - PSciCo Project
38
Delaunay Triangulation
10/3/2015
Robert Harper - PSciCo Project
39
Delaunay Triangulation
• Goal: avoid “skinny” triangles.
• Lead to ill-conditioned numeric problems.
• Method: “improve” a triangulation by
adding new vertexes.
• Turn a “skinny” triangle into two “fat” ones.
• Where to add vertices?
• Circumcenter of triangle “usually” works.
• Other choices are sometimes forced.
10/3/2015
Robert Harper - PSciCo Project
40
Delaunay Triangulation
skinny triangle
new edges
remove base
new vertex
10/3/2015
Robert Harper - PSciCo Project
41
Rupert’s Algorithm
• Adding a vertex may violate the
boundary!
• New vertex may introduce a new “skinny”
triangle somewhere else in the region.
• If the new triangle is on the boundary,
further splitting may remove a boundary
segment!
10/3/2015
Robert Harper - PSciCo Project
42
Rupert’s Algorithm
cannot be removed
10/3/2015
new skinny triangle
Robert Harper - PSciCo Project
43
Rupert’s Algorithm
put new vertex here instead
10/3/2015
Robert Harper - PSciCo Project
44
Rupert’s Algorithm
• Typical solution: search whenever a
vertex is about to be added.
• Lots of ad hoc code for a rare case.
• Introduces overhead at each step.
• A simpler solution: back out when a
boundary violation occurs.
• Can be arbitrarily far in the future.
10/3/2015
Robert Harper - PSciCo Project
45
Rupert’s Algorithm in ML
• Use dynamic exceptions and persistent
surface representation.
• Associate an exception with each “new”
vertex.
• If the boundary is violated, raise the
exception of some “new” vertex, passing a
better location for it.
• Handler restarts triangulation from there.
10/3/2015
Robert Harper - PSciCo Project
46
Rupert’s Algorithm in ML
• For the 2D case this is nice.
• simplifies coding
• easier to explain
• For the 3D case this is a big win.
• conditions much harder to state and enforce
• persistence helps enormously
• first implementation of 3D Ruppert nearly
complete
10/3/2015
Robert Harper - PSciCo Project
47
N-Body Problem
• Problem: given n bodies, compute the
potentials between them.
• Gravitational potential between stars.
• Electrical potential between charged
bodies.
• Naive solution: O(n2).
• For each body, sum potential due to all
other bodies in the system
10/3/2015
Robert Harper - PSciCo Project
48
N-Body Problem
• Approximate potentials by clustering
distant bodies.
• eg, think of Andromeda as one body
relative to the Sun
• Several choices:
• How to form clusters?
• How far is “distant”?
• How to approximate potentials?
10/3/2015
Robert Harper - PSciCo Project
49
Clustering
• Build a tree by recursively dividing
space into regions.
• Each node is a cluster.
• Each leaf is a body.
• Two division methods:
• Equal-sized regions of space (quaternary).
• Split bounding box on longest dimension
(binary).
10/3/2015
Robert Harper - PSciCo Project
50
Clustering: Even Division
10/3/2015
Robert Harper - PSciCo Project
51
Clustering: Split Bounding Box
10/3/2015
Robert Harper - PSciCo Project
52
Approximation
• Treat each node as a “pseudo-body”.
• Potential based on approximation method.
• Valid approximation for “distant” objects.
• Several approximation methods.
• Monopole = first-order term (center of gravity).
• Multipole expansion = higher-order terms.
• Spherical distribution = k points evenly distributed
on the surface of a sphere.
10/3/2015
Robert Harper - PSciCo Project
53
Force Calculation
• For each body, traverse the tree.
• If far enough away from the cluster, use
the approximation.
• Otherwise, sum potentials due to each
child.
• Distance criterion: s/d > .
• s = size of cell, d = distance away
•  = threshold factor
10/3/2015
Robert Harper - PSciCo Project
54
Force Calculation
• There is a trade-off between the
approximation and the separation
criterion.
• Smaller  implies more interactions to be
calculated.
• More terms in the expansion makes
approximation harder to calculate, but
improves accuracy.
10/3/2015
Robert Harper - PSciCo Project
55
A Modular Implementation
• Use a functor parameterized on separation,
partition, and approximation method:
functor NBody
(structure Sep : SEPARATION
structure Part : PARTITION
structure Approx : APPROX
sharing Sep.Geom = Part.Geom
and Part.Geom = Approx.Geom
) = …
10/3/2015
Robert Harper - PSciCo Project
56
A Modular Implementation
• We can mix-and-match methods to
compare speed, accuracy, etc.
• Facilitates prototyping and comparison of
methods.
• A single framework for all tree-based
methods.
• Standard methods use one “trajectory”
through the space of options.
10/3/2015
Robert Harper - PSciCo Project
57
Opportunities For Parallelism
• Tree-based methods are very amenable
to parallelization.
• During construction, can construct cells in
parallel.
• During evaluation, can compute body-cell
interactions in parallel.
• Can calculate effect on all bodies
simultaneously.
10/3/2015
Robert Harper - PSciCo Project
58
Summary
• Scientific computing demands a
programming language with
• A flexible module system,
• Numeric and symbolic computation,
• Dynamically-allocated data structures
• ... while providing
• Competitive efficiency on stock hardware.
• Super reliability and maintainability.
10/3/2015
Robert Harper - PSciCo Project
59
Summary
• Merging ML and NESL is a good start ...
• superior expressiveness, modularity
• natural support for parallelism
• ... but there’s much more to be done!
• need a “killer app” (3D Ruppert?)
• performance comparisons
10/3/2015
Robert Harper - PSciCo Project
60
Further Information
Visit the PSciCo home page:
http://www.cs.cmu.edu/~pscico
10/3/2015
Robert Harper - PSciCo Project
61
Descargar

Parallel Scientific Computing: The PSciCo Project at CMU