- 64 Views
- Uploaded on
- Presentation posted in: General

Hybrid Synchronous Languages

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -

Hybrid Synchronous Languages

Vijay Saraswat

IBM TJ Watson Research Center

March 2004

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

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

A,B,D ::= atomic formulae | A&B |X^A

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

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

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

A ::= c ~~> A

Unless c holds of the final store, run A

ask c \/ A

Leads to nondet behavior

(c ~~> c)

No behavior

(c1 ~~> c2, c2 ~~> c1)

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

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

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

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

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

}

}

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…

if x>0 then

x’’=-10

x’’=-10

if x>0 then

x’’=-10

if y=0 then

z=1

if y=0 then

z=1

X = 10

X’ = 0

if y=0 then

z=1

X = 10

x’ = 0

x’ = 0, x = 10

x’ = 0, x = 10

if x>0 then

x’’=-10

store

store

store

if y=0 then

z=1

store

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 > 0 then x’’ = -10

if y=0 then

z=1

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

x = 10,

x’ = 0

if x>0 thenx’’ = -10

if prev(x)=0 then x’=-0.5*prev(x’)

if x>0 thenx’’ = -10

if prev(x)=0 then x’=-0.5*prev(x’)

cc prog

x = 10,

x’ = 0

x = 0,

x’ = 7.07

output

x > 0,

prev(x) = 0

x’’=-10

x’’=-10

t = 0+

t = 1.414

t = 0

x=0,x’=-14.14

x=10,x’=0

x > 0,

prev(x) = 0

t = 1.414-

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’)}

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

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

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

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

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

public class Furnace implements Plant {

const Real heatR, coolR, initTemp;

public readOnly Fluent<Real> temp;

public inputOnly Fluent<Bool> 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 { temp’=heatR} on switchOn;

time always {temp’=-coolR} on ~switchOn;

}}

public class Controller {

Plant plant;

public void setPlant(Plant p) { this.plant=p;}

…

}

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();

}

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 …

Goal: To help the biologist model, simulate, analyze,design and diagnosebiological systems.

Promises profound advances in Biology and Computer Science

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)

Moore’s Law!

This is not the first time…

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.

- Expression of gene x inhibits expression of gene y; above a certain threshold, gene y inhibits expression of gene x:

if (y < 0.8) {x’= -0.02*x + 0.01},

If (y >= 0.8) {x’=-0.02*x, y’=0.01*x}

Bockmayr and Courtois: Modeling biological systems in hybrid concurrent constraint programming

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

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

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) < 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},

Results: Simulation confirms observations. Tiwari/Lincoln prove that States 2 and 3 are stable.

- 1/ [p53]’=[p53]0 –[p53]*[Mdm2]*deg -dp53*[p53]
- 2/ [Mdm2]’=p1+p2max*(I^n)/(K^n+I^n)-dMdm2*[Mdm2]
- I is some intermediary unknown mechanism; induction of [Mdm2] must be steep, n is usually > 10.
- May be better to use a discontinuous change?

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

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

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.

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

X10=Java–Threads+Locales

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.

Load input into array A striped into K locales in parallel;

Barrier b = new Barrier();

forall( i : A) {

async A[i] {

int j = f(A[i]);

async atomic A[j]{

A[j]++;

} before b;

}

await b;

The GUPS benchmark

A language for massively parallel non-uniform SMMPs

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

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

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

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.