```Hybrid Synchronous
Languages
Vijay Saraswat
IBM TJ Watson Research Center
March 2005
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.
Constraint systems
A,B,D ::= atomic formulae | A&B |X^A


Any (intuitionistic,
classical) system of
partial information
formulae, the basic
relationship is:



A1,…, An |- A
Read as “If each of the
A1,…, An hold, then A
holds”
Require conjunction,
existential quantification
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, A&B |- D
(&-L) G |- A G|- B  G |- A&B
(^-R) G |- A[t/X]  G |- X^A
(^-L,*) G,A |- D  G,X^A |- D
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
Concurrent Constraint Programming


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.
(Agents) A::= c
cA
A,B
X^A
(Config) G ::= A,…, A
G, A&B  G,A,B
G, X^A  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
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 > 0 then x’’ = -10
X = 10
X’ = 0
X = 10
if x>0 then
x’’=-10
x’ = 0
store
if y=0 then
z=1
store
if x>0 then
x’’=-10
if x>0 then
x’’=-10
if y=0 then
z=1
if y=0 then
z=1
x’’=-10
if y=0 then
z=1
x’ = 0, x = 10
store
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
store
if y=0 then
z=1
x’ = 0, x =10,
x’’ = -10
Output
Gupta/Carlson
Default CCP

A ::= c ~~> A




(c ~~> c)



No behavior
(c1 ~~> c2, c2 ~~> c1)


Unless c holds of the
final store, run A
gives c1 or c2
(c ~~> d): gives d
(c, c~~>d): 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
Discrete Timed CCP

Synchronicity principle


System reacts
instantaneously to the
environment
Semantic idea




Run a default CCP program
at each time point
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
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
jcc

Implementation of
Default Timed cc in
Java


Uses Java as a host
language
Full implementation of
Timed Default CCP




Promises, unification.
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/~sar
aswat/jcc.html
LGPL
Saraswat,Jagadeesan, Gupta “jcc: Integrating Timed Default CCP into Java”
The Esterel stopwatch program
public void WatchAgent() {
watch = WATCH;
whenever (watch==WATCH) {
unless (changeMode) {
print("Watch Mode");
}
when (p == UL) {
next
enterSetWatch =
ENTER_SET_WATCH;
changeMode=CHANGE_MODE;
print("Set Watch Mode");
}
when (p == LL) {
next stopWatch=STOP_WATCH;
changeMode=CHANGE_MODE;
print("Stop Watch Mode");
}
do {
always watch = WATCH;
} watching ( changeMode )
unless (changeMode) {
beep();
}
time.setAlarmBeeps(beeper);
}
}
Hybrid Systems




Discrete state, discrete change
(assignment)
E.g. Turing Machine
Brittleness:




Hybrid Systems combine both



Small error  major impact
Devastating with large code!
Discrete control
Continuous state evolution
Intuition: Run program at every
real value.





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)
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…
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. ODE’s.
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 ODE’s 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: x=10,x’=0, hence {if x>0 then x”=-10, if prev(x)=0 then x’=-0.5*prev(x’)}
cc prog
x = 10,
x’ = 0
output
x = 10,
x’ = 0
t=0
if x>0 then x’’ = -10
if prev(x)=0 then x’=-0.5*prev(x’)
x’’=-10
x=10,x’=0
t = 0+
x’’=-10
x=0,x’=-14.14
t = 1.414-
New combinator: hence A
x > 0,
prev(x) = 0
if x>0 then x’’ = -10
if prev(x)=0 then x’=-0.5*prev(x’)
x = 0,
x’ = 7.07
t = 1.414
x > 0,
prev(x) = 0
Gupta/Carlson
execute a copy of A in each phase (except the current point phase, if any)
Gupta, Jagadeesan, Saraswat “Computing with continuous change”, SCP 1998
Hybrid CC with interval constraints



Arithmetic variables are interval
valued. Arithmetic constraints are nonlinear 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
The arithmetic constraint system

Constraints are used to narrow the
intervals of their variables.

For example, x^2 + y^2 = 1
reduces the intervals for x and y to
[-1,1] each. Further adding x >=
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
Integrating the differential equations


Differential equations are just ordinary
algebraic equations relating some
variables and their derivatives, e.g.
f=m*a’’, x’’+d*x’+k*x=0
We provide various integrators --Euler, 4th order Runge Kutta, 4th order
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
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 00:00
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.002423875040887606, 0.02672168963230259, -0.001959654820981497),
Planet(Venus, 0.8150026784,
0.6680247657009936, 0.2606201175567890,
-0.03529355196193388, -0.007293563117650372,
0.01879420958390879, 0.0006778739390714113),
Planet(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 < 4°, Venus < 1°, other < 0.0001° away
from actual positions after 9 years.
Gupta/Carlson
Programming in jcc
public class Furnace implements Plant {
public class Controller {
const Real heatR, coolR, initTemp;
Plant plant;
public void setPlant(Plant p) { this.plant=p;}
public inputOnly Fluent<Bool> switchOn;
…
public Furnace(Real heatR, Real coolR, Real initT )
}
{
public class ControlledFurnace {
this.heatR = heatR; this.coolR = coolR;
Controller c;
this.initTemp = initT;
Furnace f;
}
public ControlledFurnace(Furnace f,
public void run() {
Controller c) {
temp = initT;
time always { temp’=heatR} on switchOn;
this.c = c; this.f = f;}
public void run() {
time always {temp’=-coolR} on ~switchOn;
}}
c.run(); c.setPlant(f); f.run();
}
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, Intercell 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
Chemical Reactions


Cells host thousands of
chemical reactions (e.g. citric
acid cycle, glycolis…)
Chemical Reaction



X+Y0 –k0 XY0
XY0 –k-0  X+Y0
Law of Mass Action




Rate of reaction is proportional
to product of conc of
components
[X]’= -k0[X][Y] + k-0[XY0]
[Y]’=[X]’
[XY]’=k0[X][Y]-K-0[XY0]




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.
State dependent rate equations

Expression of gene x
inhibits expression of
gene y; above a certain
threshold, gene y
inhibits expression of
gene x:
/*
Two constant signals fed on different
conditions generated from a
sawtooth wave.
*/
#SAMPLE_INTERVAL_MAX 0.05
x=0,y=0,
always {
if (y < 0.8) then x' = -0.02*x+0.01,
if (y >= 0.8) then x' = -0.02*x,
y' = 0.01*x
},
sample(x,y)
Bockmayr and Courtois: Modeling biological systems in hybrid concurrent
constraint programming
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 piecewise affine hybrid model”, HSCC 2001
Delta-Notch Signaling

Model:






VD, VN: concentration of
Delta and Notch protein in
the cell.
UD, UN: Delta (Notch)
production capacity of cell.
UN=sum_i VD_i (neighbors)
UD = -VN
Parameters:



Threshold values: HD,HN
Production rates: RD, RN
Model:

Cell in 1 of 4 states: {D,N} x
{Expressed (above),
Inhibited (below)}
if (UN(i,j) < HN) {VN’= -MN*VN},
if (UN(I,j)>=HN){VN’=RN-MN*VN},
if (UD(I,j)<HD){VD’=-MD*VD},
if (UD(I,j)>=HD){VD’=RD-MD*VD},


Stochastic variables used to
set random initial state.
Model can be expressed
directly in hcc.
Results: Simulation confirms observations. Tiwari/Lincoln prove that
States 2 and 3 are stable.
HCC2: Integration of Space



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.
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
data to processors) be
to obtain efficient parallel
algorithm?

Need to support round-tripping.
 Identify patterns, integrate
libraries of high-performance
code (e.g. petsc).
The X10 language








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.
striped into K locales in
parallel;
finish forall( i : A) {
async A[i] {
int j = f(A[i]);
async atomic A[j]{
A[j]++;
}
}
The GUPS benchmark
A language for massively parallel non-uniform SMPs
The X10 Programming Model
Place
Place
Partitioned Global heap
Granularity of
place can range
from single
register file to an
entire SMP system
Outbound
activities
Inbound
activities
Place-local heap
Place-local heap
Activities &
Activity-local storage
heap
stack
control
...
Activities &
Activity-local storage
heap
...
stack
control
Partitioned Global heap
heap
Inbound
activity
replies
stack
Outbound
activity
replies
heap
...
control
stack
control
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 P=1 and P > 1
Shared Memory (P=1)
MPI (P > 1)
X10 Programming


Subsumes MPI and
OpenMP programming
Closely related to Cilk
striped into K locales in
parallel;
finish forall( i : A) {
async A[i] {
int j = f(A[i]);
async atomic A[j]{
A[j]++;
}
}
The GUPS benchmark
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
Conclusion

We believe biological
system modeling and
analysis will be a very
productive area for
constraint programming and
programming languages





Handle continuous/discrete
space+time
Handle stochastic
descriptions
Handle models varying over
many orders of magnitude
Handle symbolic analysis
Handle parallel
implementations
HCC references




Gupta, Jagadeesan, Saraswat “Computing with Continuous
Change”, Science of Computer Programming, Jan 1998, 30
(1—2), pp 3--49
Saraswat, Jagadeesan, Gupta “Timed Default Concurrent
Constraint Programming”, Journal of Symbolic
Computation, Nov-Dec1996, 22 (5—6), 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”, CP’98
Workshop on Modeling and Constraints, Oct 1998.
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…
This is not the first time…

Why Now?


Exploiting genomic data
Scale





Across the internet, across
space and time.
Integration of computational
tools
Integration of new analysis
techniques
Collaboration using markupbased interlingua (SBML)
Moore’s Law!
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 nsite model by abstracting
the action at one site via a
splice efficiency function.
Results described in [Alt], uses default reasoning properties of HCC.
Controlling Cell division:
The p53-Mdm2 feedback loop


1/ [p53]’=[p53]0 –[p53]*[Mdm2]*deg -dp53*[p53]
2/ [Mdm2]’=p1+p2max*(I^n)/(K^n+I^n)-dMdm2*[Mdm2]



3/ [I]’=a*[p53]-kdelay*I




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/ a=c1*sig/(1+c2*[Mdm2]*[p53])
5/ sig’=-r*sig(t)


I is some intermediary unknown mechanism; induction of [Mdm2] must
be steep, n is usually > 10.
May be better to use a discontinuous change?
Models initial stimulus (signal) which decays rapidly, at a rate determined
by repair.
6/ deg=degbasal-[kdeg*sig-thresh]
7/ thresh’=-kdamp*thresh*sig(t=0)
Lev Bar-Or, Maya et al “Generation of oscillations by the p53-Mdm2 feedback loop..”,2000
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
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
The X10 Programming Model

Support for productivity

Support for scalability


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).
Places /
data
distribut
ion
async /
future
/
force
intraplace
atomic

scan /
reduce /
lift

Explicit
Syntax

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 userspecified iterations with index
variables.
Axiom: The user may know more
than the compiler/RTS.
```