[From Bill Powers (941117.2115 MST)]
E Coli Modelers: Ecol4bwp.pas
Here is my last stab before we go on (I hope) to talking about the
operant conditioning models we are going to use for Bruce's experiments.
In the model that follows, there are two substances diffusing from two
centers, and the organism has two higher-level control systems each
seeking a certain reference concentration of one substance. Only the
reference level for the left-hand substance is variable now by pressing
the + or - keys. When you start the program, the spot wiggles its way to
a more or less final position. At that point you can press the 'z' key to
clear the screen and redraw it, with circles indicating the places where
the reference concentration for each substance can be found. You will
find that both higher-level systems can satisfy their reference
conditions: the spot seeks an intersection of the two circles.
Now you can raise the reference level for one system using the + key. The
spot will move to a new location and hang around it. Hit the 'z' key to
draw the circles showing where the new reference concentrations are to be
found. Again you will find the spot near the intersection of the circles.
When there are no intersections of the circles, the spot just seeks the
best compromise position that minimizes both errors.
Enjoy.
Best,
Bill P.
P.S. I found that in the original NutConcen procedure, the concentration
was actually falling off as the fourth power of distance. I have fixed
that so it falls off as the inverse square, and slightly simplified the
routine. I've separated out some initializations for purposes of
redrawing the screen without reinitializing some variables. And I've
added the routine "findradius" which is entered with a reference
concentration and finds, by using a little control system, the radius at
which that concentration occurs. The resulting radius is returned as an
integer for use by the circle procedure.
···
======================================================================
{***********************************************************************
* E. COLI "REINFORCEMENT" SIMULATION 4awp *
* *
* *
* Modification of Bruce Abbott's Ecoli4a to use continuous control *
* and to seek a specific level of concentration by varying the gain *
* of the loop that controls rate of change of concentration. *
* Change ref level of concentration with + and - keys in steps of *
* 500. *
* WTP, 941116 *
* further modified for two sources of concentration 941117
*
***********************************************************************}
program Ecol4bwp;
uses
CRT, Graph, GrUtils;
const
TWOPI = PI * 2;
ENDSESSION = 50000;
var
MaxX, MaxY: integer;
NutrX, NutrY: array[1..2] of integer;
X, Y: integer;
NutMag,OldCon,NewCon,RefCon: array[1..2] of real;
EcoliX, EcoliY, timer,timelim, dNut,RefdNut, gain, error,
Speed, Angle: real;
Ch: char;
Clock: longint;
i,t: integer;
procedure InitScreen;
begin
MaxX := GetMaxX; MaxY := GetMaxY;
Rectangle(0, 0, MaxX, MaxY);
OutTextXY(MaxX div 2 - 220, 5,
'E. COLI SIMULATION # 4AWP: 2-LEVEL CONTROL MODEL');
OutTextXY(MaxX - 200, 50, 'C ref level');
OutTextXY(MaxX - 200, 60, 'C perceived ');
OutTextXY(20, MaxY-60, 'Press ESC or q to Quit...');
OutTextXY(20, MaxY-45, 'Press + or - to change left ref
concentration...');
OutTextXY(20, MaxY-30, 'Press z to draw circles at ref
concentrations');
OutTextXY(20, MaxY-15, 'Press space to pause or continue...');
end;
procedure ShowReal(x,y: integer; v: real);
var s: string;
begin
str(v:6:1, s);
setfillstyle(0,0);
bar (x,y,x+textwidth(s),y+textheight(s));
outtextxy(x,y,s);
end;
procedure placetargets;
var i: integer;
begin
for i := 1 to 2 do
begin
Rectangle(NutrX[i]-2, NutrY[i]-2, NutrX[i]+2, NutrY[i]+2);
NutMag[i] := 100.0; { max concentration }
end;
end;
procedure InitSim;
begin
Randomize;
Speed := 1.0;
for i := 1 to 2 do
begin
NutrX[i] := MAxx div 2 - MaxX div 10 + (i-1)*maxx div 5;
NutrY[i] := MaxY div 2;
RefCon[i] := 2000;
end;
EcoliX := maxx div 8 + 10 + random(100);
EcoliY := maxy div 8 - 10 + random(100);
Angle := TwoPi * Random;
EcoliX := EcoliX + Speed * cos(Angle);
EcoliY := EcoliY + Speed * sin(Angle);
X := Round(EcoliX);
Y := Round(EcoliY);
PlaceTargets;
Clock := 0;
timelim := 20.0;
timer := 0;
RefdNut := 1;
gain := -0.01;
end;
function NutConcen(X, Y: real; i: integer): real;
{ Nutient concentration at point X, Y: environment function }
begin
NutConcen := 7E2*NutMag[i] / (10.0 + 0.001*(Sqr(X) + Sqr(Y)) );
end;
function FindRadius(ref: real): integer;
var conc: real;
X, OldX: real;
begin
X := 100.0;
repeat
OldX := X;
Conc := 7E2*NutMag[i] / (10.0 + 0.001*(X*X) );
X := X - 0.01*(ref - conc)
until abs(Oldx - X) < 0.001*X;
FindRadius := round(X);
end;
procedure StepEColi;
var
NewNut: real;
begin
EcoliX := EcoliX + Speed * cos(Angle);
EcoliY := EcoliY + Speed * sin(Angle);
X := Round(EcoliX);
Y := Round(EcoliY);
PutPixel(X, Y, white);
gain := 0;
dNut := 0;
for i := 1 to 2 do
begin
NewCon[i] := NutConcen(EcoliX - NutrX[i],EcoliY - NutrY[i],i);
dNut := dNut + NewCon[i] - OldCon[i];
OldCon[i] := NewCon[i];
gain := gain + 0.01*(RefCon[i] - NewCon[i]);
timer := timer + 2 + gain * (RefdNut - dNut);
end;
if timer < 0.0 then timer := 0.0;
if timer > timelim then
begin
Angle := TwoPi * Random;
timer := 0.0;
end;
for i := 1 to 2 do
begin
showreal(MaxX - 100, 50 + 30*(i-1),RefCon[i]);
showreal(MaxX - 100, 60 + 30*(i-1),NewCon[i]);
end;
end;
begin
ClrScr;
InitGraphics;
InitScreen;
InitSim;
randomize;
repeat
inc(Clock);
ch := chr(0);
StepEcoli;
if Keypressed then Ch := ReadKey;
if ch = '+' then if RefCon[1] < 50000.0 then RefCon[1] :=
RefCon[1]*1.1;
if ch = '-' then if Refcon[1] > 500.0 then RefCon[1] :=
RefCon[1]/1.1;
if ch = ' ' then ch := readkey;
if ch = 'z' then
begin
clearviewport;
InitScreen;
PlaceTargets;
for i := 1 to 2 do
circle(Nutrx[i],NutrY[i],FindRadius(RefCon[i]));
end;
until (Ch = #27) or (Clock >= ENDSESSION) or (ch = 'q');
RestoreCRTMode;
CloseGraph;
end.