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 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 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 cA 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 (ask c) \/ A Leads to nondet behavior 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 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 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. 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/~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 Traditional Computer Science 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. 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) 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 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 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 readOnly Fluent<Real> temp; 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 Degradation rates: MD, MN 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 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. 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). The X10 language 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; 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 Load input into array A 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 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 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.

Descargar
# The Biologist's Workbench