Hybrid Synchronous Languages

1
Hybrid Synchronous Languages

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

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
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

7
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

8
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

9
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, Dec 2003
10
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)
11
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
12
Background 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
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. Users can easily add their own
  operators as C libraries (useful for connecting
  with external C tools, simulators etc.).
• Object-oriented system with methods and
  inheritance. Methods and class definitions are
  constraints and can be changed during the course
  of a program. Recursive functions are allowed.
• 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
• Arithmetic expressions compiled to byte code and
  then machine code for efficiency. Common
  subexpressions 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.

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 values, and set it as above.
  Non-linear terms are treated as separate
  variables.
• These methods can be combined to increase
  efficiency. For example, we use Splitting only
  to reduce the size of the interval of x, then use
  Newton Raphson to get the root quickly.

Gupta/Carlson
16
Integrating the differential equations

• Differential equations are just ordinary
  algebraic equations relating some variables and
  their derivatives e.g. f m a, x dx
  kx 0.
• 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 pTg 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
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
21
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.
22
State dependent rate equations

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

if (y lt 0.8) x -0.02x 0.01, If (y gt 0.8)
x-0.02x, y0.01x
Bockmayr and Courtois Modeling biological
systems in hybrid concurrent constraint
programming
23
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
24
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.
25
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
26
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
27
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.

28
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).

29
The X10 language
Load input into array A striped into K locales in
parallel Barrier b new Barrier() forall( i
A) async Ai int j f(Ai) async
atomic Aj Aj before b await
b

• 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
SMMPs
30
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

31
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

32
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.

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.