RE: law of effect versus control.
The code follows.
···
----------------------------------------------------------------------
Best,
Bill P.
{***********************************************************************
* E. COLI "REINFORCEMENT" SIMULATION 4a *
* *
* Programmer: Dr. Bruce B. Abbott *
* Psychological Sciences *
* Indiana U. - Purdue U. *
* Fort Wayne, IN 46805-1499 *
* (219) 481-6399 *
* Language: Borland/Turbo Pascal 7.0 *
* Created: 11/14/94 *
* *
* This program implements a discriminated operant simulation of the *
* "tumble-and-swim" behavior of an "e. coli" capable of learning from *
* experience. If a tumble in the presence of S+ (rising nutrient *
* level) results in a more positive rate of nutrient change, *
* p(Tumble|S+) increases, otherwise it decreases. If a tumble in the *
* presence of S- (decreasing or steady nutrient level) results in a *
* more positive rate of nutrient increase, then p(Tumble|S- increases, *
* otherwise it decreases. Thus the effect of the change in the rate *
* of nutrient change following a tumble on p(Tumble|S+) and *
* p(Tumble|S-) is symmetrical. Experience with the consequences *
* of tumbling gradually "shapes" e. coli's behavior so as to maximize *
* nutrient levels. *
* *
* MODIFIED 941128 (WTP) SO ENVIRONMENT RULE IS: *
* IF ANGLE BETWEEN DIRECTION OF TRAVEL AND DIRECTION TO TARGET *
* IS GREATER THAN 0.3*PI (54 DEGREES) ON EITHER SIDE, SET DNUT TO *
* -1, ELSE SET DNUT TO +1. THUS THE PROBABILITY OF AN UNFAVORABLE *
* RESULT IS INCREASED. *
* *
* THE 'c' KEY TOGGLES BETWEEN THE CONTROL MODEL AND THE LAW OF *
* EFFECT MODEL. THE CONTROL MODEL SIMPLY SAYS "TUMBLE IF DNUT IS *
* LESS THAN ZERO' AND MAKES NO USE OF INFORMATION ABOUT PAST VALUES *
* OF DNUT OR ABOUT IMPROVEMENTS IN DNUT AFTER A TUMBLE. THE LAW OF *
* EFFECT MODEL USES THE SAME CONDITIONAL STATEMENTS AS BEFORE. THE *
* OLD COMPUTATIONS OF PROBABILITIES ARE STILL IN PLACE BUT ARE NOT *
* USED WHEN THE CONTROL MODE IS IN EFFECT. *
* *
* *
***********************************************************************}
program Ecoli4a;
uses
CRT, Graph, GrUtils;
const
TWOPI = PI * 2;
ENDSESSION = 50000;
var
MaxX, MaxY: integer;
NutrX, NutrY, X, Y: integer;
NutMag, NutCon, dNut, NutSave: real;
EcoliX, EcoliY,
Speed, Angle, LearnRate,
pTumbleGivenSplus, pTumbleGivenSminus,
pMax, pMin: real;
JustTumbled: boolean;
Ch: char;
Clock: longint;
t,alter: integer;
control: boolean;
procedure InitScreen;
begin
ClrScr;
InitGraphics;
MaxX := GetMaxX; MaxY := GetMaxY;
Rectangle(0, 0, MaxX, MaxY);
OutTextXY(MaxX div 2 - 220, Y+5,
'E. COLI SIMULATION # 4BBA');
OutTextXY(MaxX - 200, 50, ' dNutrient');
OutTextXY(MaxX - 200, 60, 'p(Tumble|S+)');
OutTextXy(MaxX - 200, 70, 'p(Tumble|S-)');
OutTextXY(20, MaxY-50, 'Press ESC to Quit...');
end;
procedure ShowReal(x,y: integer; v: real);
var s: string;
begin
str(v:8:4, s);
setfillstyle(0,0);
bar (x,y,x+textwidth(s),y+textheight(s));
outtextxy(x,y,s);
end;
procedure ShowInt(x,y: integer; v: integer);
var s: string;
begin
str(v:8, s);
setfillstyle(0,0);
bar (x,y,x+textwidth(s),y+textheight(s));
outtextxy(x,y,s);
end;
function findangle(x1,y1,x2,y2: integer): real;
var u,v,x,y,theta: real;
begin
x := x2 - x1; y := y2 - y1;
u := abs(x); v := abs(y);
if (u = 0) and (v = 0) then
begin findangle := 0.0; exit; end;
if v <= u then theta := arctan(v/u)
else theta := Pi/2.0 - arctan(u/v);
if x <= 0.0 then theta := pi - theta;
if y <= 0.0 then theta := - theta;
findangle := theta;
end;
procedure Tumble(var Angle: real);
begin
Angle := TwoPi * Random;
end;
procedure InitSim;
begin
Randomize;
NutrX := MaxX div 2;
NutrY := MaxY div 2;
EcoliX := 50.0;
EcoliY := 50.0;
X := Round(EcoliX);
Y := Round(EcoliY);
Rectangle(NutrX-2, NutrY-2, NutrX+2, NutrY+2);
Speed := 1.0;
LearnRate := 0.01;
JustTumbled := false;
NutSave := 0;
pMax := 1.00; { Maximum tumble rate }
pMin := 0.005; { Minimum tumble rate }
pTumbleGivenSplus := 0.50; { initial tumble rates }
pTumbleGivenSminus := 0.50;
NutMag := 100.0; { max concentration }
repeat Tumble(Angle) until (Angle < PI/2);
Clock := 0;
end;
function NutConcen(X, Y: real): real;
{ Nutient concentration at point X, Y: environment function }
var
Dist: real;
begin
Dist := Sqrt(Sqr(X - NutrX) + Sqr(Y - NutrY));
NutConcen := NutMag / (1 + 0.001*(Sqr(Dist)));
end;
procedure StepEColi;
var
NewNut: real;
Theta1,Theta2,Theta3: real;
procedure ReinforceOrPunish;
var
DeltaNutRate: real;
begin
DeltaNutRate := dNut - NutSave; { Change in the rate of change in }
{ nutrient following a tumble. }
{ + = improvement = reinforcement }
{ - = deterioration = punishment }
If DeltaNutRate > 0 then { Nutrient rate increased by tumble: reinforce }
begin { tumbling }
If NutSave > 0 then { If S+ present during tumble then }
begin { increase probability of tumble given S+ }
pTumbleGivenSplus := pTumbleGivenSplus + LearnRate;
if pTumbleGivenSplus > pMax then pTumbleGivenSplus := pMax;
end
else { S- present when last tumbled then }
begin { increase probability of tumble given S- }
pTumbleGivenSminus := pTumbleGivenSminus + LearnRate;
if pTumbleGivenSminus > pMax then pTumbleGivenSminus := pMax;
end
end
else
if DeltaNutrate <= 0 then { Nutrient rate decreased by tumble: punish
}
If NutSave > 0 then { If S+ present when last tumbled then }
begin { decrease probability of tumble given S+ }
pTumbleGivenSplus := pTumbleGivenSplus - LearnRate;
if pTumbleGivenSplus < pMin then pTumbleGivenSplus := pMin;
end
else { If S- present when last tumbled then }
begin { decrease probability of tumble given S- }
pTumbleGivenSminus := pTumbleGivenSminus - LearnRate;
if pTumbleGivenSminus < pMin then pTumbleGivenSminus := pMin;
end;
end;
procedure DoTumble;
begin
Tumble(Angle);
JustTumbled := true;
NutSave := dNut; { NutSave is nutrient rate of change immediately }
end; { after a tumble }
begin
EcoliX := EcoliX + Speed * cos(Angle);
EcoliY := EcoliY + Speed * sin(Angle);
X := Round(EcoliX);
Y := Round(EcoliY);
if not control then PutPixel(X, Y, white)
else PutPixel(X, Y, lightred);
NewNut := NutConcen(EcoliX, EcoliY);
{New way of making dNut depend on angle of travel}
Theta1 := Angle;
Theta2 := Findangle(X,Y,NutrX,NutrY);
Theta3 := Theta2 - Theta1;
while Theta3 > Pi do Theta3 := Theta3 - TwoPi;
while Theta3 < -Pi do Theta3 := Theta3 + TwoPi;
if abs(Theta3) > 0.3*Pi then
dNut := -1.0 else dNut := 1.0;
{ dNut := (NewNut - NutCon);} {Old way of determining dNut}
if JustTumbled then ReinforceOrPunish;
if control then
if dNut < 0 then DoTumble else JustTumbled := false
else
begin
if dNut > 0 then { S+ present; tumble probability determined by S+ }
begin
if (Random < pTumbleGivenSplus) then DoTumble
else JustTumbled := false;
end
else { S- present; tumble probability determined by S- }
begin
if (Random < pTumbleGivenSminus) then DoTumble
else JustTumbled := false;
end;
end;
NutCon := NewNut;
showreal(maxx - 100, 50, dNut);
showreal(maxx - 100, 60, pTumbleGivenSplus);
showreal(maxx - 100, 70, pTumbleGivenSminus);
end;
begin
InitScreen;
InitSim;
control := false;
setfillstyle(0,0);
setcolor(white);
outtextxy(maxx div 2,maxy - 20,'LAW OF EFFECT');
repeat
ch := chr(0);
inc(Clock);
StepEcoli;
if Keypressed then Ch := ReadKey;
if ch = ' ' then ch := readkey;
if (ch = 'c') or (ch = 'C') then
begin
control := not control;
bar(maxx div 2,maxy - 20,maxx div 2 + 110, maxy - 10);
if control then
begin
setcolor(lightred);
outtextxy(maxx div 2,maxy - 20,'CONTROL MODE')
end
else
outtextxy(maxx div 2,maxy - 20,'LAW OF EFFECT');
end;
setcolor(white);
until (Ch = #27) or (Clock >= ENDSESSION);
if Ch <> #27 then Ch := ReadKey;
RestoreCRTMode;
CloseGraph;
end.