Hybrid Synchronous Languages

1
Hybrid Synchronous Languages

• Vijay Saraswat
• IBM TJ Watson Research Center
• March 2005

2
Outline

• CCP as a framework for concurrency theory
• Constraints
• CCP
• Defaults
• Discrete Time
• Hybrid Time
• Examples
• Themes
• Synchronous languages are widely applicable
• Space, as well as time, needs to be treated.

3
Constraint systems
A,B,D atomic formulae AB XA G
multiset of formulae (Id) A - A (Id) (Cut) G -
B G,B - D ? G,G - D (Weak) G - A ? G,B -
A (Dup) G, A, A - B ? G,A - B (Xchg) G,A,B,G
- D ? G,B,A,G - D (-R) G,A,B - D ? G, AB -
D (-L) G - A G- B ? G - AB (-R) G -
At/X ? G - XA (-L,) G,A - D ? G,XA - D

• Any (intuitionistic, classical) system of partial
  information
• For Ai read as logical formulae, the basic
  relationship is
• A1,, An - A
• Read as If each of the A1,, An hold, then A
  holds
• Require conjunction, existential quantification

4
Constraint system Examples

• Gentzen
• Herbrand
• Lists
• Finite domain
• Propositional logic (SAT)
• Arithmetic constraints
• Naïve,
• Linear
• Nonlinear
• Interval arithmetic
• Orders
• Temporal Intervals

• Hash-tables
• Arrays
• Graphs
• Constraint systems are ubiquitous in computer
  science
• Type systems
• Compiler analysis
• Symbolic computation
• Concurrent system analysis

5
Concurrent Constraint Programming
(Agents) A c c ? A
A,B
XA (Config) G A,, A G, AB ? G,A,B G, XA ?
G,A (X not free in G) G, c ? A ? G,A (s(G) -
c) A set of fixed points of a
clop Completeness for constraint entailment

• Use constraints for communication and control
  between concurrent agents operating on a shared
  store.
• Two basic operations
• tell c Add c to the store
• ask c then A If the store is strong enough to
  entail c, reduce to A.

6
Concurrent Constraint Programming

• A constraint expresses information about the
  possible values of variables. E.g. x 0, 7x
  3y 21.
• A cc program consists of a store, which is a set
  of constraints, and a set of subprograms
  independently interacting with it. A subprogram
  can add a constraint to the store. It can also
  ask if a constraint is entailed by the store, and
  if so, reduce to a new set of subprograms. The
  output is the store when all subprograms are
  quiescent.
• Example program x 10, x 0, if x gt 0 then
  x -10

if y0 then z1
Basic combinators c add constraint c
to the store. if c then A if c is
entailed by the store, execute subprogram A,
otherwise wait. A, B execute subprograms
A and B in independently. unless c then A if c
will not be entailed in the current phase,
execute A (default cc).
x 0, x 10, x -10
Output
Gupta/Carlson
7
Default CCP

• A c gt A
• Unless c holds of the final store, run A
• (ask c) \/ A
• Leads to nondet behavior
• (c gt c)
• No behavior
• (c1 gt c2, c2 gt c1)
• gives c1 or c2
• (c gt d) gives d
• (c, cgtd) gives c

• A set S of pairs (c,d) st Sd c (c,d) in
  S denotes a clop.
• Operational implementation
• Backtracking search
• Open question
• compile-time analysis
• use negation as failure

8
Discrete Timed CCP

• Synchronicity principle
• System reacts instantaneously to the environment
• Semantic idea
• Run a default CCP program at each time point
• Add A next A
• No connection between the store at one point and
  the next.
• Future cannot affect past.

• Semantics
• Sets of sequences of (pairs of) constraints
• Non-empty
• Prefix-closed
• P after s d e s.e in P is denotation of a
  default CC program

9
Timed Default CCP basic results

• The usual combinators can be programmed
• always A
• do A watching c
• whenever c do A
• time A on c
• A general combinator can be defined
• time A on B the clock fed to A is determined by
  (agent) B

• Discrete timed synchronous programming language
  with the power of Esterel
• Present is translated using defaults
• Proof system
• Compilation to automata

10
jcc

• Implementation of Default Timed cc in Java
• Uses Java as a host language
• Full implementation of Timed Default CCP
• Promises, unification.
• More CSs can be added.
• Implements defaults via backtracking.
• Uses Java GC

• Very useful as a prototyping language.
• Currently only backend implemented.
• Available from sourceforge
• http//www.cse.psu.edu/saraswat/jcc.html
• LGPL

Saraswat,Jagadeesan, Gupta jcc Integrating
Timed Default CCP into Java
11
The Esterel stopwatch program
public void WatchAgent() watch WATCH
whenever (watchWATCH) unless
(changeMode) print("Watch Mode")
when (p UL) next
enterSetWatch
ENTER_SET_WATCH changeModeCHANGE_MODE
print("Set Watch Mode")
when (p LL) next
stopWatchSTOP_WATCH changeModeCHANGE_M
ODE print("Stop Watch Mode")
do always watch WATCH
watching ( changeMode ) unless
(changeMode) beep()
time.setAlarmBeeps(beeper)
12
Hybrid Systems

• Traditional Computer Science
• Discrete state, discrete change (assignment)
• E.g. Turing Machine
• Brittleness
• Small error ? major impact
• Devastating with large code!
• Traditional Mathematics
• Continuous variables (Reals)
• Smooth state change
• Mean-value theorem
• e.g. computing rocket trajectories
• Robustness in the face of change
• Stochastic systems (e.g. Brownian motion)

• Hybrid Systems combine both
• Discrete control
• Continuous state evolution
• Intuition Run program at every real value.
• Approximate by
• Discrete change at an instant
• Continuous change in an interval
• Primary application areas
• Engineering and Control systems
• Paper transport
• Autonomous vehicles
• And now.. Biological Computation.

Emerged in early 90s in the work of Nerode, Kohn,
Alur, Dill, Henzinger
13
Extending cc to Hybrid cc

• Basic assupmtion The evolution of a system is
  piecewise continuous. Thus, a system evolution
  can be modeled a sequence of alternating point
  and interval phases.
• Constraints will now include time-varying
  expressions e.g. ODEs.
• Execute a cc program in each phase to determine
  the output of that phase. This will also
  determine the cc program to be run in the next
  phase.
• In an interval phase, any constraints asked of
  the store are recorded as transition conditions.
  Integrate the ODEs in the store to evolve the
  time dependent variables, using the store in the
  previous point phase to determine the initial
  conditions. The phase ends when any transition
  condition changes status. The values of the
  variables at the end of the phase can be used by
  the next point phase.
• Example program x10,x0, hence if xgt0 then
  x-10, if prev(x)0 then x-0.5prev(x)

Gupta/Carlson
New combinator hence A execute a copy of A
in each phase (except the current point phase, if
any)
Gupta, Jagadeesan, Saraswat Computing with
continuous change, SCP 1998
14
Hybrid CC with interval constraints

• Arithmetic variables are interval valued.
  Arithmetic constraints are non-linear algebraic
  equations over these, using standard operators
  like , , , etc.
• Object-oriented system with methods and
  inheritance.
• Various combinators defined on the basic
  combinators e.g.
• do A watching c --- execute A, abort it
  when c becomes true
• when c do A --- start A at the first instant
  when c becomes true
• wait N do A --- start A after N time units
• forall C(X) do A(X) --- execute a copy of A
  for each object X of class C

• Methods and class definitions are constraints and
  can be changed during the course of a program.
  Recursive functions are allowed.
• Arithmetic expressions compiled to byte code and
  then machine code for efficiency. Common
  sub-expressions are recognized.
• Copying garbage collector speeds up execution,
  and allows taking snapshots of states.
• API from Java/C to use Hybrid cc as a library.
  System runs on Solaris, Linux, SGI and Windows
  NT.
• Users can easily add their own operators as C
  libraries (useful for connecting with external C
  tools, simulators etc.).

Carlson, Gupta Hybrid CC with Interval
Constraints
Gupta/Carlson
15
The arithmetic constraint system

• Constraints are used to narrow the intervals of
  their variables.
• For example, x2 y2 1 reduces the intervals
  for x and y to -1,1 each. Further adding x gt
  0.5 reduces the interval for x to 0.5, 1, and
  for y to -0.866, 0.866.
• Various interval pruning methods prune one
  variable at a time.

• Indexicals Given a constraint f(x,y) 0,
  rewrite it as x g(y). If x ? I and y ? J, then
  set x ? I ? g(J). Note y can be a vector of
  variables.
• Interval splitting If x ? a, b, do binary
  search to determine minimum c in a,b such that
  0 ? f(c,c, J), where y ? J. Similary determine
  maximum such d in a,b, and set x ? c,d.
• Newton Raphson Get minimum and maximum roots of
  f(x,J) 0, where y ? J. Set x as above.
• Simplex Given the constraints on x, find its
  minimum and maximum

These methods can be combined to increase
efficiency.
Gupta/Carlson
16
Integrating the differential equations

• Differential equations are just ordinary
  algebraic equations relating some variables and
  their derivatives, e.g. fma, xdxkx0
• We provide various integrators --- Euler, 4th
  order Runge Kutta, 4th order Runge Kutta with
  adaptive stepsize, Bulirsch-Stoer with polynomial
  extrapolation. Others can be added if necessary
• All integrators have been modified to integrate
  implicit differential equations, over interval
  valued variables.

• Exact determination of discrete changes (to
  determine the end of an interval phase) is done
  using cubic Hermite interpolation.
• For example, in the example program we need to
  check if x 0. We use the value of x and x at
  the beginning and end of an integration step to
  determine if x 0 anywhere in this step.
• If so, the step is rolled back, and a smaller
  step is taken based on the estimate of the time
  when x 0.
• This is repeated till the exact time when x 0
  is determined.

Gupta/Carlson
17
Example The solar system
Planet (m, initpx, initpy, initpz, initvx,
initvy,initvz) px, py, pz, mass px
initpx, py initpy, pz initpz, px'
initvx, py' initvy, pz'initvz, always
mass m, px'' sum(g P.mass
(P.px - px)/((P.px -px)2 (P.py -py)2 (P.pz
-pz)2)1.5, Planet(P), P ! Self), py''
sum(g P.mass (P.py - py)/((P.px -px)2
(P.py -py)2 (P.pz -pz)2)1.5, Planet(P), P !
Self), pz'' sum(g P.mass (P.pz -
pz)/((P.px -px)2 (P.py -py)2 (P.pz
-pz)2)1.5, Planet(P), P ! Self)
, always p Tg 8.88769408e-10, //Coordinates,
velocities on 1998-jul-01 0000 Planet(Sun,
332981.78652, -0.008348511782195148,
0.001967945668637627, 0.0002142251001467145,
-0.000001148114436325, -0.000008994958827348018,
0.00000006538635311283), Planet(Mercury,
0.0552765501, -0.4019000379850893,
-0.04633361689674035, 0.032392079927509,-0.0024238
75040887606, -0.02672168963230259,
-0.001959654820981497),
Planet(Venus, 0.8150026784,
0.6680247657009936, 0.2606201175567890,
-0.03529355196193388, -0.007293563117650372,
0.01879420958390879, 0.0006778739390714113), Plan
et(Earth, 1.0, 0.1508758612501242,
-1.002162526305211, 0.0002082851504420832,
0.01671098890724774, 0.002627047365383169,
-0.0000004771611907632339), / A fragment
of a model for the Solar system. The remaining
lines give the coordinates and velocities of the
other planets on July 1, 1998. The class planet
implements each planet as one of n bodies,
determining its acceleration to be the sum of the
accelerations due to all other bodies (this is
defined by the sum constraint). Units are
Earth-mass, Astronomical units and Earth
days./ Results of simulation Simulated time -
3321 units (9 years). CPU time 55 s. Accuracy
Mercury lt 4, Venus lt 1, other lt 0.0001 away
from actual positions after 9 years.
Gupta/Carlson
18
Programming in jcc
public class Furnace implements Plant const
Real heatR, coolR, initTemp public readOnly
FluentltRealgt temp public inputOnly FluentltBoolgt
switchOn public Furnace(Real heatR, Real coolR,
Real initT ) this.heatR heatR this.coolR
coolR this.initTemp initT public
void run() temp initT time always
tempheatR on switchOn time always
temp-coolR on switchOn
public class Controller Plant plant public
void setPlant(Plant p) this.plantp public
class ControlledFurnace Controller c Furnace
f public ControlledFurnace(Furnace f,
Controller c)
this.c c this.f f public void run()
c.run() c.setPlant(f) f.run()
19
Systems Biology
Goal To help the biologist model, simulate,
analyze, design and diagnose biological systems.

• Develop system-level understanding of biological
  systems
• Genomic DNA, Messenger RNA, proteins, information
  pathways, signaling networks
• Intra-cellular systems, Inter-cell regulation
• Cells, Organs, Organisms
• 12 orders of magnitude in space and time!

• Key question Function from Structure
• How do various components of a biological system
  interact in order to produce complex biological
  functions?
• How do you design systems with specific
  properties (e.g. organs from cells)?
• Share Formal Theories, Code, Models

Promises profound advances in Biology and
Computer Science
20
Chemical Reactions

• Cells host thousands of chemical reactions (e.g.
  citric acid cycle, glycolis)
• Chemical Reaction
• XY0 k0? XY0
• XY0 k-0 ? XY0
• Law of Mass Action
• Rate of reaction is proportional to product of
  conc of components
• X -k0XY k-0XY0
• YX
• XYk0XY-K-0XY0

• Conservation of Mass
• When multiple reactions, sum mass flows across
  all sources and sinks to get rate of change.
• Same analysis useful for enzyme-catalyzed
  reactions
• Michaelis-Menten kinetics
• May be simulated
• Using deterministic means.
• Using stochastic means (Gillespie algorithm).

At high concentration, species concentration can
be modeled as a continuous variable.
21
State dependent rate equations
/ Two constant signals fed on different
conditions generated from a sawtooth wave.
/ SAMPLE_INTERVAL_MAX 0.05 x0,y0, always
if (y lt 0.8) then x' -0.02x0.01, if (y gt
0.8) then x' -0.02x, y' 0.01x , sample(x,y
)

• Expression of gene x inhibits expression of gene
  y above a certain threshold, gene y inhibits
  expression of gene x

Bockmayr and Courtois Modeling biological
systems in hybrid concurrent constraint
programming
22
Cell division Delta-Notch signaling in X. Laevis

• Consider cell differentiation in a population of
  epidermic cells.
• Cells arranged in a hexagonal lattice.
• Each cell interacts concurrently with its
  neighbors.
• The concentration of Delta and Notch proteins in
  each cell varies continuously.
• Cell can be in one of four states Delta and
  Notch inhibited or expressed.

• Experimental Observations
• Delta (Notch) concentrations show typical spike
  at a threshold level.
• At equilibrium, cells are in only two states (D
  or N expressed other inhibited).

Ghosh, Tomlin Lateral inhibition through
Delta-Notch signaling A piece-wise affine hybrid
model, HSCC 2001
23
Delta-Notch Signaling

• Model
• VD, VN concentration of Delta and Notch protein
  in the cell.
• UD, UN Delta (Notch) production capacity of
  cell.
• UNsum_i VD_i (neighbors)
• UD -VN
• Parameters
• Threshold values HD,HN
• Degradation rates MD, MN
• Production rates RD, RN

• Model
• Cell in 1 of 4 states D,N x Expressed
  (above), Inhibited (below)
• Stochastic variables used to set random initial
  state.
• Model can be expressed directly in hcc.

if (UN(i,j) lt HN) VN -MNVN, if
(UN(I,j)gtHN)VNRN-MNVN, if
(UD(I,j)ltHD)VD-MDVD, if (UD(I,j)gtHD)VDRD-
MDVD,
Results Simulation confirms observations.
Tiwari/Lincoln prove that States 2 and 3 are
stable.
24
HCC2 Integration of Space

• Need to add continuous space.
• Need to add discrete space.
• Use same idea as when extending CCP to HCC
• Extend uniformly across space with an elsewhere
  ( spatial hence) operator.
• Think as if a default CC program is running
  simultaneously at each spatial point.

• Implementation
• Move to PDEs from ODEs.
• Much more complicated to solve.
• Need to generate meshes.
• Use Petsc (parallel ANL library).
• Uses MPI for parallel execution.

25
Generating code for parallel machines

• There is a large gap between a simple declarative
  representation and an efficient parallel
  implementation
• Cf Molecular dynamics
• Central challenge
• How can additional pieces of information (e.g.
  about target architecture, about mapping data to
  processors) be added compositionally so as to
  obtain efficient parallel algorithm?

• Need to support round-tripping.
• Identify patterns, integrate libraries of
  high-performance code (e.g. petsc).

26
The X10 language
Load input into array A striped into K locales in
parallel finish forall( i A) async Ai
int j f(Ai) async atomic Aj
Aj

• X10JavaThreadsLocales
• A locale is a region containing data and
  activities
• An activity may atomically execute statements
  that refer to data on current locale.
• Arrays are striped across locales.
• An activity may asynchronously spawn activities
  in other locales.
• Locales may be named through data.
• Barriers are used to detect termination.
• Supports SPMD computations.

The GUPS benchmark
A language for massively parallel non-uniform SMPs
27
The X10 Programming Model
Place
Place
Outbound activities
Inbound activities
Partitioned Global heap
Partitioned Global heap
Place-local heap
Place-local heap
Granularity of place can range from single
register file to an entire SMP system
. . .
Activities Activity-local storage
Activities Activity-local storage
. . .
. . .
Inbound activity replies
Outbound activity replies
Immutable Data

• A program is a collection of places, each
  containing resident data and a dynamic collection
  of activities.
• Program may distribute aggregate data (arrays)
  across places during allocation.
• Program may directly operate only on local data,
  using atomic sections.

• Program may spawn multiple (local or remote)
  activities in parallel.
• Program must use asynchronous operations to
  access/update remote data.
• Program may repeatedly detect quiescence of a
  programmer-specified, data-dependent, distributed
  set of activities.

Cluster Computing Common framework for P1 and P
gt 1
Shared Memory (P1)
MPI (P gt 1)
28
X10 Programming

• Subsumes MPI and OpenMP programming
• Closely related to Cilk

Load input into array A striped into K locales in
parallel finish forall( i A) async Ai
int j f(Ai) async atomic Aj
Aj
The GUPS benchmark
29
Integration of symbolic reasoning techniques

• Use state of the art constraint solvers
• ICS from SRI
• Shostak combination of theories (SAT, Herbrand,
  RCF, linear arithmetic over integers).
• Finite state analysis of hybrid systems
• Generate code for HAL

• Predicate abstraction techniques.
• Develop bounded model checking.
• Parameter search techniques.
• Use/Generate constraints on parameters to rule
  out portions of the space.
• Integrate QR work
• Qualitative simulation of hybrid systems

30
Conclusion

• We believe biological system modeling and
  analysis will be a very productive area for
  constraint programming and programming languages

• Handle continuous/discrete spacetime
• Handle stochastic descriptions
• Handle models varying over many orders of
  magnitude
• Handle symbolic analysis
• Handle parallel implementations

31
HCC references

• Gupta, Jagadeesan, Saraswat Computing with
  Continuous Change, Science of Computer
  Programming, Jan 1998, 30 (12), pp 3--49
• Saraswat, Jagadeesan, Gupta Timed Default
  Concurrent Constraint Programming, Journal of
  Symbolic Computation, Nov-Dec1996, 22 (56), pp
  475-520.
• Gupta, Jagadeesan, Saraswat Programming in
  Hybrid Constraint Languages, Nov 1995, Hybrid
  Systems II, LNCS 999.
• Alenius, Gupta Modeling an AERCam A case study
  in modeling with concurrent constraint
  languages, CP98 Workshop on Modeling and
  Constraints, Oct 1998.

32
Systems Biology

• Work subsumes past work on mathematical modeling
  in biology
• Hodgkin-Huxley model for neural firing
• Michaelis-Menten equation for Enzyme Kinetics
• Gillespie algorithm for Monte-Carlo simulation of
  stochastic systems.
• Bifurcation analysis for Xenopus cell cycle
• Flux balance analysis, metabolic control analysis

• Why Now?
• Exploiting genomic data
• Scale
• Across the internet, across space and time.
• Integration of computational tools
• Integration of new analysis techniques
• Collaboration using markup-based interlingua
  (SBML)
• Moores Law!

This is not the first time
33
Alternative splicing regulation

• Alternative splicing occurs in post
  transcriptional regulation of RNA
• Through selective elimination of introns, the
  same premessenger RNA can be used to generate
  many kinds of mature RNA
• The SR protein appears to control this process
  through activation and inhibition.

• Because of complexity, experimentation can focus
  on only one site at a time.
• Bockmayr et al use Hybrid CCP to model SR
  regulation at a single site.
• Michaelis-Menten model using 7 kinetic reactions
• This is used to create an n-site model by
  abstracting the action at one site via a splice
  efficiency function.

Results described in Alt, uses default
reasoning properties of HCC.
34
Controlling Cell division The p53-Mdm2 feedback
loop

• 1/ p53p530 p53Mdm2deg -dp53p53
• 2/ Mdm2p1p2max(In)/(KnIn)-dMdm2Mdm2
• I is some intermediary unknown mechanism
  induction of Mdm2 must be steep, n is usually gt
  10.
• May be better to use a discontinuous change?
• 3/ Iap53-kdelayI
• This introduces a time delay between the
  activation of p53 and the induction of Mdm2.
  There appears to be some hidden gearing up
  mechanism at work.
• 4/ ac1sig/(1c2Mdm2p53)
• 5/ sig-rsig(t)
• Models initial stimulus (signal) which decays
  rapidly, at a rate determined by repair.
• 6/ degdegbasal-kdegsig-thresh
• 7/ thresh-kdampthreshsig(t0)

Lev Bar-Or, Maya et al Generation of
oscillations by the p53-Mdm2 feedback loop..,2000
35
The p53-Mdm2 feedback loop

• Biologists are interested in
• Dependence of amplitude and width of first wave
  on different parameters
• Dependence of waveform on delay parameter.
• Constraint expressions on parameters that still
  lead to desired oscillatory waveform would be
  most useful!

• There is a more elaborate model of the kinetics
  of the G2 DNA damage checkpoint system.
• 23 species, rate equations
• Multiple interacting cycles/pathways/regulatory
  networks
• Signal transduction
• MPF
• Cdc25
• Wee1

Aguda A quantitative analysis of the kinetics of
the G2 DNA damage checkpoint system, 1999
36
The X10 Programming Model

• Support for productivity
• Axiom OO provides proven baseline productivity,
  maintenance, portability benefits.
• Axiom Design must rule out large classes of
  errors (Type safe, Memory safe, Pointer safe,
  Lock safe, Clock safe )
• Axiom Design must support incremental
  introduction of explicit place types/remote
  operations.
• Axiom PM must integrate with static tools
  (Eclipse) -- flag performance problems, refactor
  code, detect races.
• Axiom PM must support automatic static and
  dynamic optimization (CPO).

• Support for scalability
• Axiom Programmer must have explicit linguistic
  constructs to deal with non-uniformity of access.
• Axiom A program must be able to specify a large
  collection of asynchronous activities to allow
  for efficient overlap of computation and
  communication.
• Axiom A program must use scalable
  synchronization constructs.
• Axiom The runtime may implement aggregate
  operations more efficiently than user-specified
  iterations with index variables.
• Axiom The user may know more than the
  compiler/RTS.

Places / data distribution
async / future / force
intra-place atomic
scan / reduce / lift
Explicit Syntax