{  27/3/88  08:11   sim}

%include "consts.inc"
%include "formats.inc"
%include "graph1.5.inc"
%include "menus.inc"
%include "utils.inc"

%begin
%constinteger active=yellow
%real a                        {  Productivity}
%real b                        {  Productivity}
%constinteger background=black
%constinteger blob=white
%integer bottom
%conststring(6)%array capt(0:2)="Low   ","Medium","High  "
%integer ccx=500, ccy=150
%constintegerarray cx(0:3)=100, 200, 100, 200
%constintegerarray cy(0:3)=200, 200, 100, 100
%ownstring(255) datainstr= %c
"Move the cursor to the appropriate line. Press <RETURN> to alter the value.
Type in the new value and press <RETURN> again.
For help, type ? and press <RETURN>."
%externalstring(255)%fnspec date
%integer delay
%real depletion
%owninteger depth
%constinteger digit=black
%integerarray glhist(0:2)     {  Histogram of distribution of values of Gl}
%real gla, glb, gltot, glucose
%owninteger graphics=yes
%owninteger high=60           {  Switching value for glucose}
%constinteger hist colour=white
%integer i
%ownstring(31) initial=""
%integer interval, j
%integer la, lb                {  Latency}
%owninteger low=40             {  Switching value for glucose}
%integer maxtime               {Number of steps in simulation}
%constinteger on=1, off=0
%ownreal p0                    {Low attach probability}
%ownreal p1                    {High attach probability}
%ownreal p2                    {Low detach probability}
%ownreal p3                    {high detach probability}
%constinteger rblob=25
%constinteger simno=21
%record(dataf)%array simdata(0:simno)
%conststring(19)%array simitem(0:simno)= %c
"Start Simulation",   { 0}
"Productivity A",     { 1}
"Productivity B",     { 2}
"Production time A",  { 3}
"Production time B",  { 4}
"Depletion rate",     { 5}
"Inter-step delay",   { 6}
"Low",                { 7}
"High",               { 8}
"Low attach prob.",   { 9}
"High attach prob.",  {10}
"Low detach prob.",   {11}
"High detach prob.",  {12}
"Maxtime",            {13}
"Interval",           {14}
"Graphics",           {15}
"Stop",               {16}
"Time",               {17}
"Continue simulation",{18}
"Restart simulation", {19}
"Abandon Simulation", {20}
"Display clock",      {21}
""(*)
%integerarrayname simselect
%constintegerarray simselect1(0:simno+2)=
0, 17, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 21, 16, 20, -1, 0(*)
%constintegerarray simselect2(0:simno+2)=
18, 19, 17, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 21,
 16, 20, -1, 0(*)
%conststring(2)%array state(0:3)="0", "a", "b", "ab"
%integerarray statect(0:3)
%ownintegerarray strwidth(0:3)=0(*)
%integer t                     {  Number of time-steps}
%externalintegerfnspec test symbol
%externalstring(255)%fnspec time
%integer top
%integerarray v(-25:0)
%constinteger vduno=16
%ownrecord(dataf)%array vdudata(0:vduno)
%conststring(80) vduinstr= %c
"Type <RETURN> to be able to alter parameters."
%conststring(19)%array vduitem(0:vduno)= %c
"0",                 { 0}
"a",                 { 1}
"b",                 { 2}
"ab",                { 3}
"Steps",             { 4}
"Low",               { 5}
"Medium",            { 6}
"High",              { 7}
"Glucose from A",    { 8}
"Glucose from B",    { 9}
"Stop",              {10}
"Reset Data",        {11}
"Average glucose",   {12}
"Time",              {13}
"Parameters",        {14}
"Parameters",        {15}
""(*)
%constintegerarray vduselect(0:vduno+2)= %c
13, 0, 1, 2, 3, 4, 5, 6, 7, 12, 8, 9, 14, 15, -1, 0(*)
%constintegerarray vduselect2(0:vduno+2)= %c
11, 13, 0, 1, 2, 3, 4, 5, 6, 7, 12, 8, 9, 14, 15, 10, -1, 0(*)
%integer x, y, z               {  Used by the Random number generator}

%record(line80f)%map newline80
%record(line80f) pattern
%record(line80f)%name p
  p==new(pattern)
  p_next==nil; p_prev==nil
  p_line=""
  %result==p
%end

%routine read line(%string(*)%name a)
%integer j
  a=""
  read symbol(j) %until j#' '
  %cycle
    %exit %if j=nl
    a=a.tostring(j) %if j#' '
    read symbol(j)
  %repeat
%end

!  P0=low attach probability.
!  P1=high attach probability.
!  P2=low detach probability.
!  P3=high detach probability.

%realfn prob(%string(15) direction, %real glucose)
  %if direction="detach" %thenstart      ;!  Glucose detacheses.
    %if glucose<=low %then %result=p3 %else %result=p2
  %finishelsestart                     ;!  Glucose attaching.
    %if glucose<=high %then %result=p0 %else %result=p1
  %finish
%end

!  Productivity of proteins.

%realfn lambda(%integer state)
  %result=a %if state&1=1
  %result=0
%end

%realfn mu(%integer state)
  %result=b %if state&2=2
  %result=0
%end

%routine start random
%integer i,j
%string(15) t,u,v
  t=date
  t->t.("/").u
  u->u.("/").v
  i=171*dtoi(t)*dtoi(u)
  x=i-30269*(i//30269)
  t=time
  t->t.(".").u
  j=172*dtoi(u)
  y=j-30307*(j//30307)
  z=170*(i+j)
  z=z-30323*(z//30323)
%end

%realfn random
%owninteger first=yes
  start random %and first=no %if first=yes
  x=171*x
  x=x-30269*(x//30269)
  y=172*y
  y=y-30307*(y//30307)
  z=170*z
  z=z-30323*(z//30323)
  %result=fracpt(x/30269+y/30307+z/30323)
%end


%routine next state
%real ga, gb, p, q, r, s
  t=t+1
  %if lb>=la %then i=lb %else i=la
  v(i)=v(i+1) %for i=-i, 1, -1
!  P is the probability that  gene is switched on, q that it is switched off.
  p=prob("attach",glucose)
  q=prob("attach",glucose)
  r=random        {  In the range 0 to 1}
  s=random
  %if v(-1)&1=0 %thenstart
    v(0)=v(0)!1 %if 0<=r<=p
  %finishelsestart
    v(0)=v(0)&2 %if 0<=r<=q
  %finish
  %if v(-1)&2=0 %thenstart
    v(0)=v(0)!2 %if 0<=s<=p
  %finishelsestart
    v(0)=v(0)&1 %if 0<=s<=q
  %finish
  statect(v(0))=statect(v(0))+1
  ga=lambda(v(-la)); gb=mu(v(-lb))
  gla=gla+ga; glb=glb+gb
  glucose=glucose+ga+gb-depletion
  glucose=0 %if glucose<0
  %if glucose<low %then glhist(0)=glhist(0)+1 %else %c
  %if glucose<high %then glhist(1)=glhist(1)+1 %else %c
  glhist(2)=glhist(2)+1
  gltot=gltot+glucose
%end

%routine plant data
  select output(1)
  %for i=1, 1, 14 %cycle
    print string(simdata(i)_name."=".simdata(i)_val_head_line) %if i#12
    newline
  %repeat
  print string("Total   0     a    b    ab     Low  Medium  High")
  newline
%end

%routine plant values
%owninteger count=0
  select output(1)
  print string(itod(t)." ".itod(statect(0))." ".itod(statect(1)). %c
  " ".itod(statect(2))." ".itod(statect(3))."  ".itod(glhist(0)). %c
  " ".itod(glhist(1))." ".itod(glhist(2)))
  newline
  count=count+1
  %if count>25000 %thenstart
    close output
    select output(0)
    open output(1, "sim".seqno)
    select output(1)
    count=0
  %finish
%end

%routine set up
%integer i
  %for i=0, 1, simno %cycle
    simdata(i)_name=simitem(i)
    simdata(i)_val_head==nil; simdata(i)_val_tail==nil
  %repeat
  %for i=0, 1, vduno %cycle
    vdudata(i)_name=vduitem(i)
    vdudata(i)_val_head==nil; vdudata(i)_val_tail==nil
  %repeat
  %for i=0, 1, 3 %cycle
    strwidth(i)=int(string width(state(i))/2)
  %repeat
  altbase=512; base=0
  bottom=2; top=2
  font(0)
  depth=fontheight//2
  helpfile="simfile"
  interval=1000
  print string("File of initial values: ")
  prompt("Initial file: ")
  read line(initial)
  %if initial#"" %and exists(initial) %thenstart
    open input(1, initial)
    select input(1)
    read(a); read(b)
    read(la); read(lb)
    read(depletion)
    read(delay)
    read symbol(graphics) %until graphics=yes %or graphics=no
    read symbol(display clock) %until display clock=yes %or display clock=no
    display clock=no %if graphics=no
    read(maxtime); read(interval)
    read(p0); read(p1); read(p2); read(p3)
    close input
    select input(0)
  %finishelsestart
    a=0.5; b=0.5
    la=1; lb=15
    depletion=0.3
    delay=1
    display clock=yes
    graphics=yes
    maxtime=10000; interval=1000
    p0=0.1; p1=0.5; p2=0.1; p3=0.5
  %finish
  clear clock
  simselect==simselect1
  t=0
  start screen mode
%end

%routine set initial values
%integer i
  glhist(i)=0 %for i=0, 1, 2
  gla=0; glb=0         {Initial amount of glucose produced by each gene}
  glucose=0            {Initial concentration of glucose}
  gltot=0
  statect(i)=0 %for i=0, 1, 3
  t=0
  v(i)=0 %for i=-25, 1, 0
%end

%routine set val(%record(dataf)%arrayname data, %integer stage, %string(80) val)
%record(line80f)%name q
  q==newline80
  q_line=val
  delete list(data(stage)_val)
  append cell(q, data(stage)_val)
%end

%routine set display values
  set val(simdata, 1, rtod(a, 1, 3))
  set val(simdata, 2, rtod(b, 1, 3))
  set val(simdata, 3, itod(la))
  set val(simdata, 4, itod(lb))
  set val(simdata, 5, rtod(depletion, 1, 3))
  set val(simdata, 6, itod(delay))
  set val(simdata, 7, itod(low))
  set val(simdata, 8, itod(high))
  set val(simdata, 9, rtod(p0, 1, 3))
  set val(simdata, 10, rtod(p1, 1, 3))
  set val(simdata, 11, rtod(p2, 1, 3))
  set val(simdata, 12, rtod(p3, 1, 3))
  set val(simdata, 13, itod(maxtime))
  set val(simdata, 14, itod(interval))
  set val(simdata, 15, tostring(graphics))
  set val(simdata, 17, digital time)
  set val(simdata, 21, tostring(display clock))
%end

%routine display data
%integer i
  colour(hist colour)
  %for i=0, 512, 512 %cycle
    text at(10, i+10)
    show string(rtod(a, 1, 3)." ".rtod(b, 1, 3))
    show string("  ".itod(la)."  ".itod(lb))
    show string("  ".rtod(depletion, 1, 3))
    show string("  ".rtod(p0, 1, 3)."  ".rtod(p1, 1, 3))
    show string("  ".rtod(p2, 1, 3)."  ".rtod(p3, 1, 3))
    show string("  ".itod(low)."  ".itod(high))
    show string("  ".itod(maxtime)."  ".itod(interval))
  %repeat
  colour(background)
%end

%integerfn edit data
%integer stage
%string(80) s
%switch case(0:simno)

  %on %event 15 %start
    %result=no
  %finish

  menu instructions=datainstr
  restart screen
  %cycle
    select input(0); select output(0)
    advance time
    set val(simdata, 17, digital time)
    set up menu(simdata, simselect, top, bottom, "")
    clear screen
    write menu
    stage=cursor depth
    ->case(stage)

case(0):  ;!  Start simulation.
case(19):  ;!  Restart simulation.
    restart screen
    %if graphics=yes %thenstart
      clear
      colour(background)
      fill(0, 0, 1023, 1023)
      display data
    %finish
    %if stage=0 %then start clock(ccx, ccy) %else restart clock(ccx, ccy)
    set initial values
    set display values
    plant data
    %result=yes
case(1):  ;!  Productivity A.
    read screen line(s)
    a=dtor(s)
    ->cont
case(2):  ;!  Productivity B.
    read screen line(s)
    b=dtor(s)
    ->cont
case(3):  ;!  Production time A.
    read screen line(s)
    la=dtoi(s)
    ->cont
case(4):  ;!  Production time B.
    read screen line(s)
    lb=dtoi(s)
    ->cont
case(5):  ;!  Depletion rate.
    read screen line(s)
    depletion=dtor(s)
    ->cont
case(6):  ;!  Delay.
    read screen line(s)
    delay=dtoi(s)
    ->cont
case(7):  ;!  Low.
    read screen line(s)
    low=dtoi(s)
    ->cont
case(8):  ;!  High.
    read screen line(s)
    high=dtoi(s)
    ->cont
case(9):  ;!  Low attach probability.
    read screen line(s)
    p0=dtor(s)
    ->cont
case(10):  ;!  High attach probability.
    read screen line(s)
    p1=dtor(s)
    ->cont
case(11):  ;!  Low detach probability.
    read screen line(s)
    p2=dtor(s)
    ->cont
case(12):  ;!  High detach probability.
    read screen line(s)
    p3=dtor(s)
    ->cont
case(13):  ;!  Maxtime.
    read screen line(s)
    maxtime=dtoi(s)
    ->cont
case(14):  ;!  Interval.
    read screen line(s)
    interval=dtoi(s)
    ->cont
case(15):  ;!  Graphics.
    read screen line(s)
    lower(s)
    graphics=charno(s,1)
    %if yes#graphics#no %then s="The value must be  yes  or  no." %c
    %else s=tostring(graphics)
    display clock=no %if graphics=no
    ->cont
case(16):  ;!  Stop.
    stop screen mode
    %stop
case(17):  ;!  Time.
    read screen line(s)
    set time(s)
    ->cont
case(18):  ;!  Continue Simulation.
    %if graphics=yes %thenstart
      clear
      colour(background)
      fill(0, 0, 1023, 1023)
      display data
    %finish
    restart clock(ccx, ccy)
    set display values
    plant data
    %result=yes
case(20):  ;!  Abandon Simulation.
    stop screen mode
    %result=no
case(21):  ;!  Display Clock.
    read screen line(s)
    %if charno(s,1)=yes %or charno(s,1)=no %thenstart
      display clock=charno(s,1)
      s=tostring(display clock)
    %finishelse s="This requires either  yes  or  no"
    ->cont
cont:
    set val(simdata, stage, s)
  %repeat
%end

%string(80)%fn stars(%integer n)
%string(80) x
  x=""
  x=x."*" %for n=n, -1, 1
  %result=x
%end

%routine set vduvalues
%integer i, j,
%real sum
%string(80) s
  j=statect(0)+statect(1)+statect(2)+statect(3)
  j=1 %if j=0
  %for i=0, 1, 3 %cycle
    s=pad(itod(statect(i)),6,' ').pad(itod(int(100*statect(i)/j))."%",6,' '). %c
    ": ".stars(int(25*statect(i)/j))
    set val(vdudata, i, s)
  %repeat
  j=glhist(0)+glhist(1)+glhist(2)
  s=itod(j)
  set val(vdudata, 4, s)
  j=1 %if j=0
  %for i=0, 1, 2 %cycle
    s=pad(itod(glhist(i)),6,' ').pad(itod(int(100*glhist(i)/j))."%",6,' '). %c
    ": ".stars(int(25*glhist(i)/j))
    set val(vdudata, i+5, s)
  %repeat
  sum=gla+glb
  sum=1 %if sum<0.01
  s=pad(itod(int(gla)), 6, ' ').pad(itod(int(100*gla/sum)), 4, ' ')."%"
  set val(vdudata, 8, s)
  s=pad(itod(int(glb)), 6, ' ').pad(itod(int(100*glb/sum)), 4, ' ')."%"
  set val(vdudata, 9, s)
  j=t; j=1 %if j=0
  s=rtod(gltot/j,1,3)
  set val(vdudata, 12, s)
  set val(vdudata, 13, digital time)
  s=rtod(a, 1, 3)."  ".rtod(b, 1, 3)."  ".itod(la)."  ".itod(lb). %c
    "  ".rtod(depletion, 1, 3)."  ".rtod(p0, 1, 3)."  ".rtod(p1, 1, 3). %c
    "  ".rtod(p2, 1, 3)."  ".rtod(p3, 1, 3)
  set val(vdudata, 14, s)
  s=itod(low)."  ".itod(high)."  ".itod(maxtime)."  ".itod(interval)
  set val(vdudata, 15, s)
%end

%routine display
%integer i, j, n
%owninteger gap=0, oldsecs
%string(255) x
  advance time
  %if graphics=yes %thenstart
!  Show time at intervals.
    gap=gap+1
    %if gap>=9 %thenstart
      show time
      gap=0
    %finish
!  Show states.
    font(0)
    %for i=0, 1, 3 %cycle
      %continue %if i=v(0)
      colour(blob)
      disc(cx(i), cy(i)+base, rblob)
      colour(digit)
      text at(cx(i)-strwidth(i), cy(i)-depth+base)
      show string(state(i))
    %repeat
    colour(active)
    disc(cx(v(0)), cy(v(0))+base, rblob)
    colour(digit)
    text at(cx(v(0))-strwidth(v(0)), cy(v(0))-depth+base)
    show string(state(v(0)))
!  Show glucose concentration.
    text at(300, 100+base)
    %if glucose<=low %then colour(active) %else colour(blob)
    show string("low")
    text at(300, 150+base)
    %if low<=glucose<=high %then colour(active) %else colour(blob)
    show string("medium")
    text at(300,200+base)
    %if high<=glucose %then colour(active) %else colour(blob)
    show string("high")
!  Clear area for histograms.
    colour(background)
    fill(0, base+250, 700, base+511)
!  Plant Glhist.
    colour(histcolour)
    n=glhist(0)+glhist(1)+glhist(2)
    n=1 %if n=0
    %for i=0, 1, 2 %cycle
      text at(10, 500-25*i+base)
      show string(capt(i))
      text at(10+70, 500-25*i+base)
      j=int(100*glhist(i)/n)
      x=itod(j)."%"
      x=" ".x %while length(x)<5
      show string(x)
      j=int(50*glhist(i)/n)
      text at(10+130, 500-25*i+base)
      show string(stars(j))
    %repeat
    text at(10, 500-75+base)
    n=t; n=1 %if n=0
    show string(vduitem(12)."  ".rtod(gltot/n,1,3))
!  Plant Statect.
    n=statect(0)+statect(1)+statect(2)+statect(3)
    n=1 %if n=0
    %for i=0, 1, 3 %cycle
      text at(10, 400-25*i+base)
      show string(state(i))
      text at(10+70, 400-25*i+base)
      j=int(100*statect(i)/n)
      x=itod(j)."%"
      x=" ".x %while length(x)<5
      show string(x)
      text at(10+130, 400-25*i+base)
      j=int(50*statect(i)/n)
      show string(stars(j))
    %repeat
!  Show time.
    colour(histcolour)
    text at(10, base+256)
    show string("Steps: ".itod(t)."  gla=".itod(int(gla)). %c
    "  glb=".itod(int(glb))."  ".digital time)
!  Display new half-frame.
    offset(0, base)
!  Swop Y-bases of half-frames.
    n=altbase
    altbase=base
    base=n
    show time %if gap=0
  %finish
!  Display state at intervals.
  %if t=interval*(t//interval) %thenstart
    set vduvalues
    menu instructions=""
    set up menu(vdudata, vduselect, 0, 1, "")
    write menu
!  Plant values in file.
    plant values
  %finish
%end


{Main Program}

  open output(1,"Sim".seqno)
  set up
  set display values
  set initial values
  select input(0)
  %cycle
    %exit %if edit data=no
    display
    %cycle
      j=test symbol
      %exit %if j>0
      next state
      display
!  Delay before next step.
      j=0 %for i=1,1,delay
    %repeat %until t>=maxtime
    set vduvalues
    simselect==simselect2
    bottom=2; top=3
    menu instructions=vduinstr
    set up menu(vdudata, vduselect2, 1, 1, "")
    write menu
    j=cursor depth
    %exit %if j=10
  %repeat
  stop screen mode


%endofprogram