## Corona virus simulator

General FreeBASIC programming questions.
Tourist Trap
Posts: 2901
Joined: Jun 02, 2015 16:24

### Re: Corona virus simulator

I can't say what exactly we are supposed to see. Impressive anyway :)
UEZ
Posts: 614
Joined: May 05, 2017 19:59
Location: Germany

### Re: Corona virus simulator

@badidea: lol, looks like an infection of a borg collective.

@BasicScience: is there a special candidate you are asking?
BasicScience
Posts: 489
Joined: Apr 18, 2008 4:09
Location: Los Angeles, CA
Contact:

### Re: Corona virus simulator

@UEZ the changes I suggested could be added to any of the posted simulations. Would be interesting with yours, for example.
BasicScience
Posts: 489
Joined: Apr 18, 2008 4:09
Location: Los Angeles, CA
Contact:

### Re: Corona virus simulator

I see you are planning to implement a maximum radius of motion for each person. I guess the initial position should be a random location in the personal circle. Otherwise the motion vector will be along a radius line, and the reflection at the circle boarder will just be an oscillation back and forth along a diameter. I started looking into coding this but... calculating where a ray intersects the circle would require solving quadratic (computationally expensive). Of course this would be done only for those steps when the distance from origin > radius. Perhaps a reflection could be done at the end of this step just outside the circle. The boundary would be fuzzy (but real life is fuzzy too), and it would save having to solve a quadratic.
BasicScience
Posts: 489
Joined: Apr 18, 2008 4:09
Location: Los Angeles, CA
Contact:

### Re: Corona virus simulator

Pseudo-random movement within a circle for "safer at home" simulation

Code: Select all

`Type __person    dim as single x_init, y_init, x , y, sqradius    dim as single  Vx, Vy    dim as ulong   cEND TYPEDIM as integer N_people = 2, iPsize = 5ScreenRes 640, 480, 32DIM Shared as __person person(N_people-1)Randomize'InitializePerson(0).x_init = 100person(0).y_init = 100person(0).Vx = 2*(rnd-.5)person(0).Vy = 2*(rnd-.5)person(0).c      = &hFFFFFFFFperson(0).sqradius = 2500person(1).x_init = 300person(1).y_init = 200Person(1).Vx = (rnd-.5)Person(1).Vy = (rnd-.5)Person(1).c      = &hFFFFFFFFPerson(1).sqradius = 6400For i as integer = 0 to N_People -1    Person(i).x = Person(i).x_init + sqr(person(i).sqradius)*(rnd-.5)    person(i).y = person(i).y_init + sqr(person(i).sqradius)*(rnd-.5)    Line (Person(i).x, Person(i).y)-(Person(i).x + iPSize, Person(i).y + iPSize), Person(i).c, BF     circle (person(i).x_init, person(i).y_init),sqr(person(i).sqradius), rgb(255,0,0)Next ido until inkey = "" : loopDIM as integer cnt = 0Do    CLS        For i as integer = 0 to N_people-1        IF CNT > 200 then            Person(i).Vx = Person(i).Vx * sgn(rnd-.5)            Person(i).Vy = Person(i).Vy * sgn(rnd-.5)            Cnt = 0        End IF        Cnt += 1                Person(i).x += Person(i).Vx         Person(i).y += Person(i).Vy                IF (Person(i).x - Person(i).x_init) * (person(i).x - person(i).x_init) + (person(i).y - Person(i).y_init) * (person(i).y -person(i).y_init) > person(i).sqradius then                        'Get back inside circle            Person(i).x -= Person(i).Vx            Person(i).y -= Person(i).Vy                        'Assume current x,y is on circle (not really true)            'Computer new velocities base on reflection to tangent                        DIM as single Nx, Ny, scale  ' normal            Nx = -(Person(i).x - Person(i).X_init)            Ny = -(Person(i).y - Person(i).Y_init)            scale = 2*(Person(i).Vx * Nx + Person(i).Vy * Ny)/(Nx*Nx + Ny*Ny)            Person(i).Vx = Person(i).Vx - scale * Nx            Person(i).Vy = Person(i).Vy - Scale * Ny                        Person(i).x += Person(i).Vx            Person(i).y += Person(i).Vy        End if        Line (Person(i).x, Person(i).y)-(Person(i).x + iPSize, Person(i).y + iPSize), Person(i).c, BF         circle (person(i).x_init, person(i).y_init),sqr(person(i).sqradius), rgb(255,0,0)    next i    Sleep 5Loop while inkey = ""sleep`
Posts: 2134
Joined: May 24, 2007 22:10
Location: The Netherlands

### Re: Corona virus simulator

Virus outbreak form the patient zero perspective. <space> to trigger. Uses a semi-random walk.
To do:
* Sectioning to handle more persons
* Social distancing control

Code: Select all

`#Include "fbgfx.bi"#Include "string.bi"'===============================================================================Type int2d   As Integer x, y   Declare Constructor   Declare Constructor(x As Integer, y As Integer)   Declare Operator Cast () As StringEnd TypeConstructor int2dEnd ConstructorConstructor int2d(x As Integer, y As Integer)   This.x = x : This.y = yEnd Constructor' "x, y"Operator int2d.cast () As String  Return Str(x) & "," & Str(y)End OperatorOperator = (a As int2d, b As int2d) As boolean   If a.x <> b.x Then Return false   If a.y <> b.y Then Return false   Return trueEnd OperatorOperator <> (a As int2d, b As int2d) As boolean   If a.x = b.x And a.y = b.y Then Return false   Return trueEnd Operator' a + b Operator + (a As int2d, b As int2d) As int2d   Return Type(a.x + b.x, a.y + b.y)End Operator' a - bOperator - (a As int2d, b As int2d) As int2d   Return Type(a.x - b.x, a.y - b.y)End Operator' -aOperator - (a As int2d) As int2d   Return Type(-a.x, -a.y)End Operator' a * bOperator * (a As int2d, b As int2d) As int2d   Return Type(a.x * b.x, a.y * b.y)End Operator' a * mulOperator * (a As int2d, mul As Integer) As int2d   Return Type(a.x * mul, a.y * mul)End Operator' a \ bOperator \ (a As int2d, b As int2d) As int2d   Return Type(a.x \ b.x, a.y \ b.y)End Operator' a \ divOperator \ (a As int2d, divider As Integer) As int2d   Return Type(a.x \ divider, a.y \ divider)End Operator'===============================================================================Type sgl2d   As Single x, y   Declare Constructor   Declare Constructor(x As Single, y As Single)   Declare Operator Cast() As String   Declare Function cross(b As sgl2d) As Single   Declare Function lengthSqrd() As Single   Declare Function dist(b As sgl2d) As Single   Declare Function distSqrd(b As sgl2d) As Single   Declare Function normalise() As sgl2dEnd TypeConstructor sgl2dEnd ConstructorConstructor sgl2d(x As Single, y As Single)   This.x = x : This.y = yEnd ConstructorFunction sgl2d.cross(b As sgl2d) As Single   Return This.x * b.y - This.y * b.xEnd FunctionFunction sgl2d.lengthSqrd() As Single   Return (This.x * This.x) + (This.y * This.y) End FunctionFunction sgl2d.dist(b As sgl2d) As Single   Dim As Single dx = This.x - b.x   Dim As Single dy = This.y - b.y   Return Sqr((dx * dx) + (dy * dy)) End FunctionFunction sgl2d.distSqrd(b As sgl2d) As Single   Dim As Single dx = This.x - b.x   Dim As Single dy = This.y - b.y   Return (dx * dx) + (dy * dy) End FunctionFunction sgl2d.normalise() As sgl2d   Dim As Single length = Sqr((This.x * This.x) + (This.y * This.y))   Return sgl2d(This.x / length, This.y / length)End Function' "x, y"Operator sgl2d.cast() As String   Return Str(x) & "," & Str(y)End Operator'---- operators ---' distance / lenthOperator Len (a As sgl2d) As Single   Return Sqr(a.x * a.x + a.y * a.y)End Operator' a = b ?Operator = (a As sgl2d, b As sgl2d) As boolean   If a.x <> b.x Then Return false   If a.y <> b.y Then Return false   Return trueEnd Operator' a != b ?Operator <> (a As sgl2d, b As sgl2d) As boolean   If a.x = b.x And a.y = b.y Then Return false   Return trueEnd Operator' a + b Operator + (a As sgl2d, b As sgl2d) As sgl2d   Return Type(a.x + b.x, a.y + b.y)End Operator' a - bOperator - (a As sgl2d, b As sgl2d) As sgl2d   Return Type(a.x - b.x, a.y - b.y)End Operator' -aOperator - (a As sgl2d) As sgl2d   Return Type(-a.x, -a.y)End Operator' a * bOperator * (a As sgl2d, b As sgl2d) As sgl2d   Return Type(a.x * b.x, a.y * b.y)End Operator' a * mulOperator * (a As sgl2d, mul As Single) As sgl2d   Return Type(a.x * mul, a.y * mul)End Operator' a / divOperator / (a As sgl2d, div As Single) As sgl2d   Return Type(a.x / div, a.y / div)End Operator'===============================================================================#Define MOUSE_IDLE 0#Define MOUSE_POS_CHANGED 1#Define MOUSE_LB_PRESSED 2#Define MOUSE_LB_RELEASED 3#Define MOUSE_RB_PRESSED 4#Define MOUSE_RB_RELEASED 5#Define MOUSE_MB_PRESSED 6#Define MOUSE_MB_RELEASED 7#Define MOUSE_WHEEL_UP 8#Define MOUSE_WHEEL_DOWN 9Type mouse_type   Pos As int2d   posChange As int2d   wheel As Integer   buttons As Integer   lb As Integer 'left button   rb As Integer 'right button   mb As Integer 'middle buttonEnd TypeFunction handleMouse(ByRef mouse As mouse_type) As Integer   Static previous As mouse_type   Dim As Integer change = MOUSE_IDLE   GetMouse mouse.pos.x, mouse.pos.y, mouse.wheel, mouse.buttons   If (mouse.buttons = -1) Then      mouse.lb = 0      mouse.rb = 0      mouse.mb = 0      mouse.posChange.x = 0      mouse.posChange.y = 0   Else      mouse.lb = (mouse.buttons And 1)      mouse.rb = (mouse.buttons ShR 1) And 1      mouse.mb = (mouse.buttons ShR 2) And 1      If (previous.pos.x <> mouse.pos.x Or previous.pos.y <> mouse.pos.y) Then         change = MOUSE_POS_CHANGED      End If      mouse.posChange.x = mouse.pos.x - previous.pos.x      mouse.posChange.y = mouse.pos.y - previous.pos.y      If (previous.buttons <> mouse.buttons) Then         If (previous.lb = 0 And mouse.lb = 1) Then change = MOUSE_LB_PRESSED         If (previous.lb = 1 And mouse.lb = 0) Then change = MOUSE_LB_RELEASED         If (previous.rb = 0 And mouse.rb = 1) Then change = MOUSE_RB_PRESSED         If (previous.rb = 1 And mouse.rb = 0) Then change = MOUSE_RB_RELEASED         If (previous.mb = 0 And mouse.mb = 1) Then change = MOUSE_MB_PRESSED         If (previous.mb = 1 And mouse.mb = 0) Then change = MOUSE_MB_RELEASED      End If      If (mouse.wheel > previous.wheel) Then change = MOUSE_WHEEl_UP      If (mouse.wheel < previous.wheel) Then change = MOUSE_WHEEl_DOWN      previous = mouse   End If   Return changeEnd Function'==============================================================================='Note: y+ = up, x+ = right, (0,0) = centerType scaled_graphics_type   Dim As Single scale = 1 ' = 1 / pixel_size 'pixels / meter   'dim as int2d offset' = (scrn_w \ 2, h \ 2) 'offset in pixels   Dim As sgl2d offset   Dim As Integer w = -1, h = -1   Dim As Integer wc = -1, hc = -1 'center x,y   Declare Sub setScreen(w As Integer, h As Integer)   Declare Sub setScaling(scale As Single, offset As sgl2d)   Declare Sub clearScreen(c As ULong)   Declare Function pos2screen(p As sgl2d) As int2d   Declare Sub drawPixel(p As sgl2d, c As ULong)   Declare Sub drawCircle(p As sgl2d, r As Single, c As ULong)   'declare sub drawCircleFilled(p as sgl2d, r as single, c as ulong, cFill as ulong)   Declare Sub drawCircleFilled(p As sgl2d, r As Single, c As ULong)   Declare Sub drawElipse(p As sgl2d, r As Single, aspect As Single, c As ULong)   Declare Sub drawLine(p1 As sgl2d, p2 As sgl2d, c As ULong)End TypeSub scaled_graphics_type.setScreen(w As Integer, h As Integer)   This.w = w 'width   This.h = h 'height   wc = w \ 2   hc = h \ 2   ScreenRes w, h, 32   Width w \ 8, h \ 16 'bigger fontEnd SubSub scaled_graphics_type.setScaling(scale As Single, offset As sgl2d)   This.scale = scale   This.offset = offsetEnd SubSub scaled_graphics_type.clearScreen(c As ULong)   Line(0, 0)-(w - 1, h - 1), c, bfEnd SubFunction scaled_graphics_type.pos2screen(p As sgl2d) As int2d   Return int2d(Int(wc + (p.x - offset.x) * scale), h - Int(hc + (p.y - offset.y) * scale))End FunctionSub scaled_graphics_type.drawPixel(p As sgl2d, c As ULong)   Dim As int2d posScrn = pos2screen(p)   PSet(posScrn.x, posScrn.y), cEnd SubSub scaled_graphics_type.drawCircle(p As sgl2d, r As Single, c As ULong)   Dim As int2d posScrn = pos2screen(p)   Circle(posScrn.x, posScrn.y), r * scale, cEnd SubSub scaled_graphics_type.drawCircleFilled(p As sgl2d, r As Single, c As ULong)   Dim As int2d posScrn = pos2screen(p)   Circle(posScrn.x, posScrn.y), r * scale, c,,,,fEnd SubSub scaled_graphics_type.drawElipse(p As sgl2d, r As Single, aspect As Single, c As ULong)   Dim As int2d posScrn = pos2screen(p)   Circle(posScrn.x, posScrn.y), r * scale, c, , , aspectEnd SubSub scaled_graphics_type.drawLine(p1 As sgl2d, p2 As sgl2d, c As ULong)   Dim As int2d posScrn1 = pos2screen(p1)   Dim As int2d posScrn2 = pos2screen(p2)   Line(posScrn1.x, posScrn1.y)-(posScrn2.x, posScrn2.y), cEnd Sub'===============================================================================Const As Single INFEC_DIST = 0.20 'mConst As Single INFEC_DIST_SQRD = INFEC_DIST ^ 2Const As Single SICK_DAYS_MIN = 7, SICK_DAYS_MAX = 21Const As Single MORTALITY = 3 / 100'const as single SICK_SPEED_FACTOR = 0.1 'disable due to some bugConst SEC_PER_DAY = 24 * 60 * 60Const As Single DAY_PER_SEC = 1 / SEC_PER_DAY Const PERSONS = 500Const As Single V_MAX = 1 / 3600 '1 m/hrConst As Single MAP_X_MAX = +50, MAP_X_MIN = -MAP_X_MAXConst As Single MAP_Y_MAX = +35, MAP_Y_MIN = -MAP_Y_MAXConst DRAW_INTERVAL = 1 '1 to 10 Const FOLLOW_PATIENT_0 = TRUE 'FALSE'screen stuffConst As Single PPM = 20.0 'pixels per meter (set zoom level)Const SCRN_W = 1050, SCRN_H = 750Dim Shared As scaled_graphics_type sgsg.setScreen(SCRN_W, SCRN_H)sg.setScaling(PPM, sgl2d(0, 0))Enum HEALTH_STATE   ST_INIT '0   ST_INFEC '1   ST_RECOV '2   ST_DEAD '3   ST_LAST 'number of statesEnd EnumDim As Integer statCounters(ST_LAST - 1)Dim As ULong stateColor(ST_LAST - 1) = _   {RGB(0, 150, 0), RGB(150, 0, 0), RGB(150, 150, 0), RGB(0, 0, 0)}Type person_type   p As sgl2d 'position [m]   v As sgl2d 'velocity [m/s]   r As Single 'radius [m]   sickEndTime As Single 'not double!   state As Integer 'healthEnd TypeFunction drawUpdate(interval As Integer) As boolean   Static As Integer counter = 0   counter += 1   If counter >= interval Then      counter = 0      Return TRUE   Else      Return FALSE   End IfEnd FunctionFunction rndRangeSgl(min As Single, max As Single) As Single   Return Rnd() * (max - min) + minEnd FunctionDim As Integer mouseEvent, mouseDragDim As mouse_type mouseDim As person_type person(0 To PERSONS-1)Randomize Timer'initialise personsFor i As Integer = 0 To UBound(person)   person(i).p.x = rndRangeSgl(MAP_X_MIN, MAP_X_MAX) 'm   person(i).p.y = rndRangeSgl(MAP_Y_MIN, MAP_Y_MAX) 'm   person(i).v.x = rndRangeSgl(-V_MAX, +V_MAX) 'm/s   person(i).v.y = rndRangeSgl(-V_MAX, +V_MAX) 'm/s   person(i).r = 0.25   'person(i).state = int(rndRangeSgl(0, csng(ST_LAST)))   person(i).state = ST_INITNext'time step such that a max speed, position change = 10 % of infection distanceDim As Single timeStep = (INFEC_DIST * 0.1) / V_MAX 'secDim As Single simulDays, simulTime = 0Dim As Integer numInit, numDead, numRecovWhile Not MultiKey(FB.SC_ESCAPE)   'trigger outbreak   If MultiKey(FB.SC_SPACE) Then      person(0).state = ST_INFEC      person(0).sickEndTime = simulTime _         + rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY   End If   'zoom / drag view by mouse   mouseEvent = handleMouse(mouse)      If (mouse.buttons <> -1) Then      If (mouseEvent = MOUSE_LB_PRESSED) Then mouseDrag = 1      If (mouseEvent = MOUSE_LB_RELEASED) Then mouseDrag = 0      If (mouseEvent = MOUSE_WHEEl_UP) Then sg.scale *= 1.1      If (mouseEvent = MOUSE_WHEEl_DOWN) Then sg.scale /= 1.1   End If   If (mouseDrag) Then      sg.offset.x -= (mouse.posChange.x / PPM)      sg.offset.y += (mouse.posChange.y / PPM)   End If   'patient 0 perpective   If FOLLOW_PATIENT_0 = TRUE Then      sg.offset.x = (person(0).p.x)      sg.offset.y = (person(0).p.y)   End If   'update positions   For i As Integer = 0 To UBound(person)      person(i).p += person(i).v * timeStep   Next   'random direction change   For i As Integer = 1 To 5      Dim As Integer iPerson = Int(rndRangeSgl(0, PERSONS))      If Rnd() > 0.5 Then         person(iPerson).v.x = -person(iPerson).v.x      Else         person(iPerson).v.y = -person(iPerson).v.y      End If   Next   'check wall collisions   For i As Integer = 0 To UBound(person)      If person(i).p.x < MAP_X_MIN Then person(i).v.x = +Abs(person(i).v.x)      If person(i).p.x > MAP_X_MAX Then person(i).v.x = -abs(person(i).v.x)      If person(i).p.y < MAP_Y_MIN Then person(i).v.y = +Abs(person(i).v.y)      If person(i).p.y > MAP_Y_MAX Then person(i).v.y = -abs(person(i).v.y)   Next   'check end of sickness   For i As Integer = 0 To UBound(person)      If person(i).state = ST_INFEC Then         If simulTime > person(i).sickEndTime Then            If Rnd() < MORTALITY Then               person(i).state = ST_DEAD               person(i).v = sgl2d(0, 0)            Else               person(i).state = ST_RECOV               'person(i).v /= SICK_SPEED_FACTOR            End If         End If      End If   Next   'check for disease transmission (iSrc = source, iTar = target)   For iSrc As Integer = 0 To UBound(person)      If person(iSrc).state = ST_INFEC Then         For iTar As Integer = 0 To UBound(person)            If person(iTar).state = ST_INIT Then               If person(iSrc).p.distSqrd(person(iTar).p) < INFEC_DIST_SQRD Then                  person(iTar).state = ST_INFEC                  person(iTar).sickEndTime = simulTime _                     + rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY                  'person(iTar).v *= SICK_SPEED_FACTOR               End If            End If         Next      End If   Next   'clear stats   For i As Integer = 0 To UBound(statCounters)      statCounters(i) = 0   Next   'update stats   For i As Integer = 0 To UBound(person)      statCounters(person(i).state) += 1   Next   'draw world   If (drawUpdate(DRAW_INTERVAL) = TRUE) Then      ScreenLock      sg.clearScreen(RGB(250, 250, 250))      sg.drawCircle(person(0).p, person(0).r * 3, stateColor(person(0).state))      For i As Integer = 0 To UBound(person)         Dim As ULong c = stateColor(person(i).state)         sg.drawCircleFilled(person(i).p, person(i).r, c)      Next      Draw String (10, 10), "Day: " & Format(simulDays, "0.#0"), RGB(0, 0, 0)      Draw String (10, 10 + 16), "Initial: " & statCounters(ST_INIT), stateColor(ST_INIT)      Draw String (10, 10 + 32), "Infected: " & statCounters(ST_INFEC), stateColor(ST_INFEC)      Draw String (10, 10 + 48), "Recovered: " & statCounters(ST_RECOV), stateColor(ST_RECOV)      Draw String (10, 10 + 64), "Dead: " & statCounters(ST_DEAD), stateColor(ST_DEAD)      ScreenUnLock      Sleep 1,1 'don't sleep every loop   End If   simulTime += timeStep   simulDays = simulTime * DAY_PER_SECWendEnd`
albert
Posts: 5889
Joined: Sep 28, 2006 2:41
Location: California, USA

### Re: Corona virus simulator

I came up with a cure for Covid-19

You take the white blood cells from persons who have recovered the fastest..
And transplant them into persons who aren't recovering well.. ( people in the hospitals )

It would also work for the regular influenza that kills far more people each year..
Posts: 2134
Joined: May 24, 2007 22:10
Location: The Netherlands

### Re: Corona virus simulator

A version with grid based spacial partitioning. I hope that I did not make an error, with 6 nested for-loops and a bunch of if-statements, a bit nasty.
It is a bit heavy on the drawing, so I only display 1 in 10 frames. Can be speed up further with macros, ditching the classes and testing the optimal spatial sectioning grid size.

Code: Select all

`#Include "fbgfx.bi"#Include "string.bi"Type int2d   As Integer x, y   Declare Constructor   Declare Constructor(x As Integer, y As Integer)   Declare Operator Cast () As StringEnd TypeConstructor int2dEnd ConstructorConstructor int2d(x As Integer, y As Integer)   This.x = x : This.y = yEnd Constructor' "x, y"Operator int2d.cast () As String  Return Str(x) & "," & Str(y)End OperatorOperator = (a As int2d, b As int2d) As boolean   If a.x <> b.x Then Return false   If a.y <> b.y Then Return false   Return trueEnd OperatorOperator <> (a As int2d, b As int2d) As boolean   If a.x = b.x And a.y = b.y Then Return false   Return trueEnd Operator' a + b Operator + (a As int2d, b As int2d) As int2d   Return Type(a.x + b.x, a.y + b.y)End Operator' a - bOperator - (a As int2d, b As int2d) As int2d   Return Type(a.x - b.x, a.y - b.y)End Operator' -aOperator - (a As int2d) As int2d   Return Type(-a.x, -a.y)End Operator' a * bOperator * (a As int2d, b As int2d) As int2d   Return Type(a.x * b.x, a.y * b.y)End Operator' a * mulOperator * (a As int2d, mul As Integer) As int2d   Return Type(a.x * mul, a.y * mul)End Operator' a \ bOperator \ (a As int2d, b As int2d) As int2d   Return Type(a.x \ b.x, a.y \ b.y)End Operator' a \ divOperator \ (a As int2d, divider As Integer) As int2d   Return Type(a.x \ divider, a.y \ divider)End Operator'===============================================================================Type sgl2d   As Single x, y   Declare Constructor   Declare Constructor(x As Single, y As Single)   Declare Operator Cast() As String   Declare Function cross(b As sgl2d) As Single   Declare Function lengthSqrd() As Single   Declare Function dist(b As sgl2d) As Single   Declare Function distSqrd(b As sgl2d) As Single   Declare Function normalise() As sgl2dEnd TypeConstructor sgl2dEnd ConstructorConstructor sgl2d(x As Single, y As Single)   This.x = x : This.y = yEnd ConstructorFunction sgl2d.cross(b As sgl2d) As Single   Return This.x * b.y - This.y * b.xEnd FunctionFunction sgl2d.lengthSqrd() As Single   Return (This.x * This.x) + (This.y * This.y) End FunctionFunction sgl2d.dist(b As sgl2d) As Single   Dim As Single dx = This.x - b.x   Dim As Single dy = This.y - b.y   Return Sqr((dx * dx) + (dy * dy)) End FunctionFunction sgl2d.distSqrd(b As sgl2d) As Single   Dim As Single dx = This.x - b.x   Dim As Single dy = This.y - b.y   Return (dx * dx) + (dy * dy) End FunctionFunction sgl2d.normalise() As sgl2d   Dim As Single length = Sqr((This.x * This.x) + (This.y * This.y))   Return sgl2d(This.x / length, This.y / length)End Function' "x, y"Operator sgl2d.cast() As String   Return Str(x) & "," & Str(y)End Operator'---- operators ---' distance / lenthOperator Len (a As sgl2d) As Single   Return Sqr(a.x * a.x + a.y * a.y)End Operator' a = b ?Operator = (a As sgl2d, b As sgl2d) As boolean   If a.x <> b.x Then Return false   If a.y <> b.y Then Return false   Return trueEnd Operator' a != b ?Operator <> (a As sgl2d, b As sgl2d) As boolean   If a.x = b.x And a.y = b.y Then Return false   Return trueEnd Operator' a + b Operator + (a As sgl2d, b As sgl2d) As sgl2d   Return Type(a.x + b.x, a.y + b.y)End Operator' a - bOperator - (a As sgl2d, b As sgl2d) As sgl2d   Return Type(a.x - b.x, a.y - b.y)End Operator' -aOperator - (a As sgl2d) As sgl2d   Return Type(-a.x, -a.y)End Operator' a * bOperator * (a As sgl2d, b As sgl2d) As sgl2d   Return Type(a.x * b.x, a.y * b.y)End Operator' a * mulOperator * (a As sgl2d, mul As Single) As sgl2d   Return Type(a.x * mul, a.y * mul)End Operator' a / divOperator / (a As sgl2d, div As Single) As sgl2d   Return Type(a.x / div, a.y / div)End Operator'==============================================================================='list can grow, but never shrink, for performance, non-sortedType grow_list   Dim As Integer index(Any)   Private:   Dim As Integer numItems   Public:   Declare Constructor(startSize As Integer)   Declare Constructor()   Declare Destructor()   Declare Sub Add(newIndex As Integer)   Declare Sub empty()   Declare Function numAlloc() As Integer   Declare Function numInUse() As Integer   Declare Sub show()   'non-list methodsEnd TypeConstructor grow_list(startSize As Integer)   If startSize > 0 Then      ReDim index(startSize - 1)   End IfEnd ConstructorConstructor grow_list()   This.constructor(0)End ConstructorDestructor grow_list()   Erase(index)End DestructorSub grow_list.add(newIndex As Integer)   Dim As Integer ub = UBound(index)   'if list is full, increase list size by 1   If numItems = ub + 1 Then      ReDim Preserve index(numItems)   End If   index(numItems) = newIndex   numItems += 1End SubSub grow_list.empty()   numItems = 0End SubFunction grow_list.numAlloc() As Integer   Return UBound(index) + 1End FunctionFunction grow_list.numInUse() As Integer   Return numItemsEnd Function'for debuggingSub grow_list.show()   Print "--- " & numInUse() & " / " & numAlloc() & " ---"   For i As Integer = 0 To numItems - 1      Print i, index(i)   NextEnd Sub'===============================================================================#Define MOUSE_IDLE 0#Define MOUSE_POS_CHANGED 1#Define MOUSE_LB_PRESSED 2#Define MOUSE_LB_RELEASED 3#Define MOUSE_RB_PRESSED 4#Define MOUSE_RB_RELEASED 5#Define MOUSE_MB_PRESSED 6#Define MOUSE_MB_RELEASED 7#Define MOUSE_WHEEL_UP 8#Define MOUSE_WHEEL_DOWN 9Type mouse_type   Pos As int2d   posChange As int2d   wheel As Integer   buttons As Integer   lb As Integer 'left button   rb As Integer 'right button   mb As Integer 'middle buttonEnd TypeFunction handleMouse(ByRef mouse As mouse_type) As Integer   Static previous As mouse_type   Dim As Integer change = MOUSE_IDLE   GetMouse mouse.pos.x, mouse.pos.y, mouse.wheel, mouse.buttons   If (mouse.buttons = -1) Then      mouse.lb = 0      mouse.rb = 0      mouse.mb = 0      mouse.posChange.x = 0      mouse.posChange.y = 0   Else      mouse.lb = (mouse.buttons And 1)      mouse.rb = (mouse.buttons ShR 1) And 1      mouse.mb = (mouse.buttons ShR 2) And 1      If (previous.pos.x <> mouse.pos.x Or previous.pos.y <> mouse.pos.y) Then         change = MOUSE_POS_CHANGED      End If      mouse.posChange.x = mouse.pos.x - previous.pos.x      mouse.posChange.y = mouse.pos.y - previous.pos.y      If (previous.buttons <> mouse.buttons) Then         If (previous.lb = 0 And mouse.lb = 1) Then change = MOUSE_LB_PRESSED         If (previous.lb = 1 And mouse.lb = 0) Then change = MOUSE_LB_RELEASED         If (previous.rb = 0 And mouse.rb = 1) Then change = MOUSE_RB_PRESSED         If (previous.rb = 1 And mouse.rb = 0) Then change = MOUSE_RB_RELEASED         If (previous.mb = 0 And mouse.mb = 1) Then change = MOUSE_MB_PRESSED         If (previous.mb = 1 And mouse.mb = 0) Then change = MOUSE_MB_RELEASED      End If      If (mouse.wheel > previous.wheel) Then change = MOUSE_WHEEl_UP      If (mouse.wheel < previous.wheel) Then change = MOUSE_WHEEl_DOWN      previous = mouse   End If   Return changeEnd Function'==============================================================================='Note: y+ = up, x+ = right, (0,0) = centerType scaled_graphics_type   Dim As Single scale = 1 ' = 1 / pixel_size 'pixels / meter   'dim as int2d offset' = (scrn_w \ 2, h \ 2) 'offset in pixels   Dim As sgl2d offset   Dim As Integer w = -1, h = -1   Dim As Integer wc = -1, hc = -1 'center x,y   Declare Sub setScreen(w As Integer, h As Integer)   Declare Sub setScaling(scale As Single, offset As sgl2d)   Declare Sub clearScreen(c As ULong)   Declare Function pos2screen(p As sgl2d) As int2d   Declare Sub drawPixel(p As sgl2d, c As ULong)   Declare Sub drawCircle(p As sgl2d, r As Single, c As ULong)   'declare sub drawCircleFilled(p as sgl2d, r as single, c as ulong, cFill as ulong)   Declare Sub drawCircleFilled(p As sgl2d, r As Single, c As ULong)   Declare Sub drawElipse(p As sgl2d, r As Single, aspect As Single, c As ULong)   Declare Sub drawLine(p1 As sgl2d, p2 As sgl2d, c As ULong)   Declare Sub drawRect(p1 As sgl2d, p2 As sgl2d, c As ULong)   Declare Sub drawRectFilled(p1 As sgl2d, p2 As sgl2d, c As ULong)End TypeSub scaled_graphics_type.setScreen(w As Integer, h As Integer)   This.w = w 'width   This.h = h 'height   wc = w \ 2   hc = h \ 2   ScreenRes w, h, 32   Width w \ 8, h \ 16 'bigger fontEnd SubSub scaled_graphics_type.setScaling(scale As Single, offset As sgl2d)   This.scale = scale   This.offset = offsetEnd SubSub scaled_graphics_type.clearScreen(c As ULong)   Line(0, 0)-(w - 1, h - 1), c, bfEnd SubFunction scaled_graphics_type.pos2screen(p As sgl2d) As int2d   Return int2d(Int(wc + (p.x - offset.x) * scale), h - Int(hc + (p.y - offset.y) * scale))End FunctionSub scaled_graphics_type.drawPixel(p As sgl2d, c As ULong)   Dim As int2d posScrn = pos2screen(p)   PSet(posScrn.x, posScrn.y), cEnd SubSub scaled_graphics_type.drawCircle(p As sgl2d, r As Single, c As ULong)   Dim As int2d posScrn = pos2screen(p)   Circle(posScrn.x, posScrn.y), r * scale, cEnd SubSub scaled_graphics_type.drawCircleFilled(p As sgl2d, r As Single, c As ULong)   Dim As int2d posScrn = pos2screen(p)   Circle(posScrn.x, posScrn.y), r * scale, c,,,,fEnd SubSub scaled_graphics_type.drawElipse(p As sgl2d, r As Single, aspect As Single, c As ULong)   Dim As int2d posScrn = pos2screen(p)   Circle(posScrn.x, posScrn.y), r * scale, c, , , aspectEnd SubSub scaled_graphics_type.drawLine(p1 As sgl2d, p2 As sgl2d, c As ULong)   Dim As int2d posScrn1 = pos2screen(p1)   Dim As int2d posScrn2 = pos2screen(p2)   Line(posScrn1.x, posScrn1.y)-(posScrn2.x, posScrn2.y), cEnd SubSub scaled_graphics_type.drawRect(p1 As sgl2d, p2 As sgl2d, c As ULong)   Dim As int2d posScrn1 = pos2screen(p1)   Dim As int2d posScrn2 = pos2screen(p2)   Line(posScrn1.x, posScrn1.y)-(posScrn2.x, posScrn2.y), c, bEnd SubSub scaled_graphics_type.drawRectFilled(p1 As sgl2d, p2 As sgl2d, c As ULong)   Dim As int2d posScrn1 = pos2screen(p1)   Dim As int2d posScrn2 = pos2screen(p2)   Line(posScrn1.x, posScrn1.y)-(posScrn2.x, posScrn2.y), c, bfEnd Sub'===============================================================================Const As Single INFEC_DIST = 0.25 'mConst As Single INFEC_DIST_SQRD = INFEC_DIST ^ 2Const As Single SICK_DAYS_MIN = 7, SICK_DAYS_MAX = 21Const As Single MORTALITY = 3 / 100'const as single SICK_SPEED_FACTOR = 0.1 'disable due to some bugConst SEC_PER_DAY = 24 * 60 * 60Const As Single DAY_PER_SEC = 1 / SEC_PER_DAY Const PERSONS = 10000Const As Single V_MAX = 1 / 3600 '1 m/hrConst As Single MAP_X_MAX = +50, MAP_X_MIN = -MAP_X_MAXConst As Single MAP_Y_MAX = +35, MAP_Y_MIN = -MAP_Y_MAXConst As Single MAP_W = MAP_X_MAX - MAP_X_MINConst As Single MAP_H = MAP_Y_MAX - MAP_Y_MINConst DRAW_INTERVAL = 10 '1 to 10 Const FOLLOW_PATIENT_0 = FALSEConst NUM_SECT_X = 16 * 3, NUM_SECT_Y = 12 * 3Dim As grow_list sector(NUM_SECT_X - 1, NUM_SECT_Y -1)'var pSector = new grow_list[NUM_SECT_X, NUM_SECT_Y]'screen stuffConst As Single PPM = 10.0 'pixels per meter (set zoom level)Const SCRN_W = 1050, SCRN_H = 750Dim Shared As scaled_graphics_type sgsg.setScreen(SCRN_W, SCRN_H)sg.setScaling(PPM, sgl2d(0, 0))Enum HEALTH_STATE   ST_INIT '0   ST_INFEC '1   ST_RECOV '2   ST_DEAD '3   ST_LAST 'number of statesEnd EnumDim As Integer statCounters(ST_LAST - 1)Dim As ULong stateColor(ST_LAST - 1) = _   {RGB(0, 150, 0), RGB(150, 0, 0), RGB(150, 150, 0), RGB(0, 0, 0)}Type person_type   p As sgl2d 'position [m]   v As sgl2d 'velocity [m/s]   r As Single 'radius [m]   sickEndTime As Single 'not double!   state As Integer 'healthEnd TypeFunction drawUpdate(interval As Integer) As boolean   Static As Integer counter = 0   counter += 1   If counter >= interval Then      counter = 0      Return TRUE   Else      Return FALSE   End IfEnd FunctionFunction rndRangeSgl(min As Single, max As Single) As Single   Return Rnd() * (max - min) + minEnd FunctionSub limitInt(ByRef value As Integer, min As Integer, max As Integer)   If value < min Then value = min   If value > max Then value = maxEnd SubDim As Integer mouseEvent, mouseDragDim As mouse_type mouseDim As person_type person(0 To PERSONS-1)Randomize Timer'initialise personsFor i As Integer = 0 To UBound(person)   person(i).p.x = rndRangeSgl(MAP_X_MIN, MAP_X_MAX) 'm   person(i).p.y = rndRangeSgl(MAP_Y_MIN, MAP_Y_MAX) 'm   person(i).v.x = rndRangeSgl(-V_MAX, +V_MAX) 'm/s   person(i).v.y = rndRangeSgl(-V_MAX, +V_MAX) 'm/s   person(i).r = 0.25   'person(i).state = int(rndRangeSgl(0, csng(ST_LAST)))   person(i).state = ST_INITNext'time step such that a max speed, position change = 10 % of infection distanceDim As Single timeStep = (INFEC_DIST * 0.1) / V_MAX 'secDim As Single simulDays, simulTime = 0Dim As Integer numInit, numDead, numRecovDim As boolean sectioning = TRUEWhile Not MultiKey(FB.SC_ESCAPE)   'trigger outbreak   If MultiKey(FB.SC_SPACE) Then      person(0).state = ST_INFEC      person(0).sickEndTime = simulTime _         + rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY   End If   'zoom / drag view by mouse   mouseEvent = handleMouse(mouse)      If (mouse.buttons <> -1) Then      If (mouseEvent = MOUSE_LB_PRESSED) Then mouseDrag = 1      If (mouseEvent = MOUSE_LB_RELEASED) Then mouseDrag = 0      If (mouseEvent = MOUSE_WHEEl_UP) Then sg.scale *= 1.1      If (mouseEvent = MOUSE_WHEEl_DOWN) Then sg.scale /= 1.1   End If   If (mouseDrag) Then      sg.offset.x -= (mouse.posChange.x / PPM)      sg.offset.y += (mouse.posChange.y / PPM)   End If   'patient 0 perpective   If FOLLOW_PATIENT_0 = TRUE Then      sg.offset.x = (person(0).p.x)      sg.offset.y = (person(0).p.y)   End If   'update positions   For i As Integer = 0 To UBound(person)      person(i).p += person(i).v * timeStep   Next   'random direction change   '~ for i as integer = 1 to 5      '~ dim as integer iPerson = int(rndRangeSgl(0, PERSONS))      '~ if rnd() > 0.5 then         '~ person(iPerson).v.x = -person(iPerson).v.x      '~ else         '~ person(iPerson).v.y = -person(iPerson).v.y      '~ end if   '~ next   'check wall collisions   For i As Integer = 0 To UBound(person)      If person(i).p.x < MAP_X_MIN Then person(i).v.x = +Abs(person(i).v.x)      If person(i).p.x > MAP_X_MAX Then person(i).v.x = -abs(person(i).v.x)      If person(i).p.y < MAP_Y_MIN Then person(i).v.y = +Abs(person(i).v.y)      If person(i).p.y > MAP_Y_MAX Then person(i).v.y = -abs(person(i).v.y)   Next   'check end of sickness   For i As Integer = 0 To UBound(person)      If person(i).state = ST_INFEC Then         If simulTime > person(i).sickEndTime Then            If Rnd() < MORTALITY Then               person(i).state = ST_DEAD               person(i).v = sgl2d(0, 0)            Else               person(i).state = ST_RECOV               'person(i).v /= SICK_SPEED_FACTOR            End If         End If      End If   Next   'clear sectors   For xi As Integer = 0 To NUM_SECT_X - 1      For yi As Integer = 0 To NUM_SECT_Y - 1         sector(xi, yi).empty()      Next   Next   'assign persons to sectors   For i As Integer = 0 To UBound(person)      Dim As Integer xi = Int((person(i).p.x - MAP_X_MIN) / (MAP_W / NUM_SECT_X))      Dim As Integer yi = Int((person(i).p.y - MAP_Y_MIN) / (MAP_H / NUM_SECT_Y))      limitInt(xi, 0, NUM_SECT_X - 1)      limitInt(yi, 0, NUM_SECT_Y - 1)      sector(xi, yi).add(i)   Next   If sectioning = TRUE Then      For xi As Integer = 0 To NUM_SECT_X - 1         For yi As Integer = 0 To NUM_SECT_Y - 1            'loop source in 1 sector            For iiSrc As Integer = 0 To sector(xi, yi).numInUse - 1               Dim As Integer iSrc = sector(xi, yi).index(iiSrc)               If person(iSrc).state = ST_INFEC Then                  'check against targets in near sectors, including this sector                  For xdi As Integer = -1 To +1                     If xi + xdi < 0 Then Continue For                     If xi + xdi >= NUM_SECT_X Then Continue For                     For ydi As Integer = -1 To +1                        If yi + ydi < 0 Then Continue For                        If yi + ydi >= NUM_SECT_Y Then Continue For                        'loop targets with 1 (near) sector                        For iiTar As Integer = 0 To sector(xi + xdi, yi + ydi).numInUse - 1                           Dim As Integer iTar = sector(xi + xdi, yi + ydi).index(iiTar)                           If person(iTar).state = ST_INIT Then                              If person(iSrc).p.distSqrd(person(iTar).p) < INFEC_DIST_SQRD Then                                 person(iTar).state = ST_INFEC                                 person(iTar).sickEndTime = simulTime _                                    + rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY                                 'person(iTar).v *= SICK_SPEED_FACTOR                              End If                           End If                        Next                     Next                  Next               End If            Next         Next      Next   Else      'case: sectioning = FALSE      'check for disease transmission (iSrc = source, iTar = target)      For iSrc As Integer = 0 To UBound(person)         If person(iSrc).state = ST_INFEC Then            For iTar As Integer = 0 To UBound(person)               If person(iTar).state = ST_INIT Then                  If person(iSrc).p.distSqrd(person(iTar).p) < INFEC_DIST_SQRD Then                     person(iTar).state = ST_INFEC                     person(iTar).sickEndTime = simulTime _                        + rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY                     'person(iTar).v *= SICK_SPEED_FACTOR                  End If               End If            Next         End If      Next   End If   'clear stats   For i As Integer = 0 To UBound(statCounters)      statCounters(i) = 0   Next   'update stats   For i As Integer = 0 To UBound(person)      statCounters(person(i).state) += 1   Next   'draw world   If (drawUpdate(DRAW_INTERVAL) = TRUE) Then      ScreenLock      sg.clearScreen(RGB(250, 250, 250))      For xi As Integer = 0 To NUM_SECT_X - 1         For yi As Integer = 0 To NUM_SECT_Y - 1            Dim As Single x1 = MAP_W * (xi / NUM_SECT_X) + MAP_X_MIN            Dim As Single y1 = MAP_H * (yi / NUM_SECT_Y) + MAP_Y_MIN            Dim As Single x2 = MAP_W * ((xi + 1) / NUM_SECT_X) + MAP_X_MIN            Dim As Single y2 = MAP_H * ((yi + 1) / NUM_SECT_Y) + MAP_Y_MIN            sg.drawRect(sgl2d(x1, y1), sgl2d(x2, y2), RGB(0, 0, 0))            Dim As int2d scrnPos = sg.pos2screen(sgl2d(x1, y1))            Draw String (scrnPos.x + 2, scrnPos.y - 16), Str(sector(xi, yi).numInUse), 0            'draw string (scrnPos.x + 2, scrnPos.y - 32), str(sector(xi, yi).numAlloc), 0         Next      Next      sg.drawCircle(person(0).p, person(0).r * 3, stateColor(person(0).state))      For i As Integer = 0 To UBound(person)         Dim As ULong c = stateColor(person(i).state)         sg.drawCircleFilled(person(i).p, person(i).r, c)      Next      Draw String (10, 10), "Day: " & Format(simulDays, "0.#0"), RGB(0, 0, 0)      Draw String (10, 10 + 16), "Initial: " & statCounters(ST_INIT), stateColor(ST_INIT)      Draw String (10, 10 + 32), "Infected: " & statCounters(ST_INFEC), stateColor(ST_INFEC)      Draw String (10, 10 + 48), "Recovered: " & statCounters(ST_RECOV), stateColor(ST_RECOV)      Draw String (10, 10 + 64), "Dead: " & statCounters(ST_DEAD), stateColor(ST_DEAD)      ScreenUnLock      Sleep 1,1 'don't sleep every loop   End If   simulTime += timeStep   simulDays = simulTime * DAY_PER_SECWendEnd`
UEZ
Posts: 614
Joined: May 05, 2017 19:59
Location: Germany

### Re: Corona virus simulator

badidea wrote:A version with grid based spacial partitioning. I hope that I did not make an error, with 6 nested for-loops and a bunch of if-statements, a bit nasty.
It is a bit heavy on the drawing, so I only display 1 in 10 frames. Can be speed up further with macros, ditching the classes and testing the optimal spatial sectioning grid size.

Code: Select all

`#Include "fbgfx.bi"#Include "string.bi"Type int2d   As Integer x, y   Declare Constructor   Declare Constructor(x As Integer, y As Integer)   Declare Operator Cast () As StringEnd TypeConstructor int2dEnd ConstructorConstructor int2d(x As Integer, y As Integer)   This.x = x : This.y = yEnd Constructor' "x, y"Operator int2d.cast () As String  Return Str(x) & "," & Str(y)End OperatorOperator = (a As int2d, b As int2d) As boolean   If a.x <> b.x Then Return false   If a.y <> b.y Then Return false   Return trueEnd OperatorOperator <> (a As int2d, b As int2d) As boolean   If a.x = b.x And a.y = b.y Then Return false   Return trueEnd Operator' a + b Operator + (a As int2d, b As int2d) As int2d   Return Type(a.x + b.x, a.y + b.y)End Operator' a - bOperator - (a As int2d, b As int2d) As int2d   Return Type(a.x - b.x, a.y - b.y)End Operator' -aOperator - (a As int2d) As int2d   Return Type(-a.x, -a.y)End Operator' a * bOperator * (a As int2d, b As int2d) As int2d   Return Type(a.x * b.x, a.y * b.y)End Operator' a * mulOperator * (a As int2d, mul As Integer) As int2d   Return Type(a.x * mul, a.y * mul)End Operator' a \ bOperator \ (a As int2d, b As int2d) As int2d   Return Type(a.x \ b.x, a.y \ b.y)End Operator' a \ divOperator \ (a As int2d, divider As Integer) As int2d   Return Type(a.x \ divider, a.y \ divider)End Operator'===============================================================================Type sgl2d   As Single x, y   Declare Constructor   Declare Constructor(x As Single, y As Single)   Declare Operator Cast() As String   Declare Function cross(b As sgl2d) As Single   Declare Function lengthSqrd() As Single   Declare Function dist(b As sgl2d) As Single   Declare Function distSqrd(b As sgl2d) As Single   Declare Function normalise() As sgl2dEnd TypeConstructor sgl2dEnd ConstructorConstructor sgl2d(x As Single, y As Single)   This.x = x : This.y = yEnd ConstructorFunction sgl2d.cross(b As sgl2d) As Single   Return This.x * b.y - This.y * b.xEnd FunctionFunction sgl2d.lengthSqrd() As Single   Return (This.x * This.x) + (This.y * This.y) End FunctionFunction sgl2d.dist(b As sgl2d) As Single   Dim As Single dx = This.x - b.x   Dim As Single dy = This.y - b.y   Return Sqr((dx * dx) + (dy * dy)) End FunctionFunction sgl2d.distSqrd(b As sgl2d) As Single   Dim As Single dx = This.x - b.x   Dim As Single dy = This.y - b.y   Return (dx * dx) + (dy * dy) End FunctionFunction sgl2d.normalise() As sgl2d   Dim As Single length = Sqr((This.x * This.x) + (This.y * This.y))   Return sgl2d(This.x / length, This.y / length)End Function' "x, y"Operator sgl2d.cast() As String   Return Str(x) & "," & Str(y)End Operator'---- operators ---' distance / lenthOperator Len (a As sgl2d) As Single   Return Sqr(a.x * a.x + a.y * a.y)End Operator' a = b ?Operator = (a As sgl2d, b As sgl2d) As boolean   If a.x <> b.x Then Return false   If a.y <> b.y Then Return false   Return trueEnd Operator' a != b ?Operator <> (a As sgl2d, b As sgl2d) As boolean   If a.x = b.x And a.y = b.y Then Return false   Return trueEnd Operator' a + b Operator + (a As sgl2d, b As sgl2d) As sgl2d   Return Type(a.x + b.x, a.y + b.y)End Operator' a - bOperator - (a As sgl2d, b As sgl2d) As sgl2d   Return Type(a.x - b.x, a.y - b.y)End Operator' -aOperator - (a As sgl2d) As sgl2d   Return Type(-a.x, -a.y)End Operator' a * bOperator * (a As sgl2d, b As sgl2d) As sgl2d   Return Type(a.x * b.x, a.y * b.y)End Operator' a * mulOperator * (a As sgl2d, mul As Single) As sgl2d   Return Type(a.x * mul, a.y * mul)End Operator' a / divOperator / (a As sgl2d, div As Single) As sgl2d   Return Type(a.x / div, a.y / div)End Operator'==============================================================================='list can grow, but never shrink, for performance, non-sortedType grow_list   Dim As Integer index(Any)   Private:   Dim As Integer numItems   Public:   Declare Constructor(startSize As Integer)   Declare Constructor()   Declare Destructor()   Declare Sub Add(newIndex As Integer)   Declare Sub empty()   Declare Function numAlloc() As Integer   Declare Function numInUse() As Integer   Declare Sub show()   'non-list methodsEnd TypeConstructor grow_list(startSize As Integer)   If startSize > 0 Then      ReDim index(startSize - 1)   End IfEnd ConstructorConstructor grow_list()   This.constructor(0)End ConstructorDestructor grow_list()   Erase(index)End DestructorSub grow_list.add(newIndex As Integer)   Dim As Integer ub = UBound(index)   'if list is full, increase list size by 1   If numItems = ub + 1 Then      ReDim Preserve index(numItems)   End If   index(numItems) = newIndex   numItems += 1End SubSub grow_list.empty()   numItems = 0End SubFunction grow_list.numAlloc() As Integer   Return UBound(index) + 1End FunctionFunction grow_list.numInUse() As Integer   Return numItemsEnd Function'for debuggingSub grow_list.show()   Print "--- " & numInUse() & " / " & numAlloc() & " ---"   For i As Integer = 0 To numItems - 1      Print i, index(i)   NextEnd Sub'===============================================================================#Define MOUSE_IDLE 0#Define MOUSE_POS_CHANGED 1#Define MOUSE_LB_PRESSED 2#Define MOUSE_LB_RELEASED 3#Define MOUSE_RB_PRESSED 4#Define MOUSE_RB_RELEASED 5#Define MOUSE_MB_PRESSED 6#Define MOUSE_MB_RELEASED 7#Define MOUSE_WHEEL_UP 8#Define MOUSE_WHEEL_DOWN 9Type mouse_type   Pos As int2d   posChange As int2d   wheel As Integer   buttons As Integer   lb As Integer 'left button   rb As Integer 'right button   mb As Integer 'middle buttonEnd TypeFunction handleMouse(ByRef mouse As mouse_type) As Integer   Static previous As mouse_type   Dim As Integer change = MOUSE_IDLE   GetMouse mouse.pos.x, mouse.pos.y, mouse.wheel, mouse.buttons   If (mouse.buttons = -1) Then      mouse.lb = 0      mouse.rb = 0      mouse.mb = 0      mouse.posChange.x = 0      mouse.posChange.y = 0   Else      mouse.lb = (mouse.buttons And 1)      mouse.rb = (mouse.buttons ShR 1) And 1      mouse.mb = (mouse.buttons ShR 2) And 1      If (previous.pos.x <> mouse.pos.x Or previous.pos.y <> mouse.pos.y) Then         change = MOUSE_POS_CHANGED      End If      mouse.posChange.x = mouse.pos.x - previous.pos.x      mouse.posChange.y = mouse.pos.y - previous.pos.y      If (previous.buttons <> mouse.buttons) Then         If (previous.lb = 0 And mouse.lb = 1) Then change = MOUSE_LB_PRESSED         If (previous.lb = 1 And mouse.lb = 0) Then change = MOUSE_LB_RELEASED         If (previous.rb = 0 And mouse.rb = 1) Then change = MOUSE_RB_PRESSED         If (previous.rb = 1 And mouse.rb = 0) Then change = MOUSE_RB_RELEASED         If (previous.mb = 0 And mouse.mb = 1) Then change = MOUSE_MB_PRESSED         If (previous.mb = 1 And mouse.mb = 0) Then change = MOUSE_MB_RELEASED      End If      If (mouse.wheel > previous.wheel) Then change = MOUSE_WHEEl_UP      If (mouse.wheel < previous.wheel) Then change = MOUSE_WHEEl_DOWN      previous = mouse   End If   Return changeEnd Function'==============================================================================='Note: y+ = up, x+ = right, (0,0) = centerType scaled_graphics_type   Dim As Single scale = 1 ' = 1 / pixel_size 'pixels / meter   'dim as int2d offset' = (scrn_w \ 2, h \ 2) 'offset in pixels   Dim As sgl2d offset   Dim As Integer w = -1, h = -1   Dim As Integer wc = -1, hc = -1 'center x,y   Declare Sub setScreen(w As Integer, h As Integer)   Declare Sub setScaling(scale As Single, offset As sgl2d)   Declare Sub clearScreen(c As ULong)   Declare Function pos2screen(p As sgl2d) As int2d   Declare Sub drawPixel(p As sgl2d, c As ULong)   Declare Sub drawCircle(p As sgl2d, r As Single, c As ULong)   'declare sub drawCircleFilled(p as sgl2d, r as single, c as ulong, cFill as ulong)   Declare Sub drawCircleFilled(p As sgl2d, r As Single, c As ULong)   Declare Sub drawElipse(p As sgl2d, r As Single, aspect As Single, c As ULong)   Declare Sub drawLine(p1 As sgl2d, p2 As sgl2d, c As ULong)   Declare Sub drawRect(p1 As sgl2d, p2 As sgl2d, c As ULong)   Declare Sub drawRectFilled(p1 As sgl2d, p2 As sgl2d, c As ULong)End TypeSub scaled_graphics_type.setScreen(w As Integer, h As Integer)   This.w = w 'width   This.h = h 'height   wc = w \ 2   hc = h \ 2   ScreenRes w, h, 32   Width w \ 8, h \ 16 'bigger fontEnd SubSub scaled_graphics_type.setScaling(scale As Single, offset As sgl2d)   This.scale = scale   This.offset = offsetEnd SubSub scaled_graphics_type.clearScreen(c As ULong)   Line(0, 0)-(w - 1, h - 1), c, bfEnd SubFunction scaled_graphics_type.pos2screen(p As sgl2d) As int2d   Return int2d(Int(wc + (p.x - offset.x) * scale), h - Int(hc + (p.y - offset.y) * scale))End FunctionSub scaled_graphics_type.drawPixel(p As sgl2d, c As ULong)   Dim As int2d posScrn = pos2screen(p)   PSet(posScrn.x, posScrn.y), cEnd SubSub scaled_graphics_type.drawCircle(p As sgl2d, r As Single, c As ULong)   Dim As int2d posScrn = pos2screen(p)   Circle(posScrn.x, posScrn.y), r * scale, cEnd SubSub scaled_graphics_type.drawCircleFilled(p As sgl2d, r As Single, c As ULong)   Dim As int2d posScrn = pos2screen(p)   Circle(posScrn.x, posScrn.y), r * scale, c,,,,fEnd SubSub scaled_graphics_type.drawElipse(p As sgl2d, r As Single, aspect As Single, c As ULong)   Dim As int2d posScrn = pos2screen(p)   Circle(posScrn.x, posScrn.y), r * scale, c, , , aspectEnd SubSub scaled_graphics_type.drawLine(p1 As sgl2d, p2 As sgl2d, c As ULong)   Dim As int2d posScrn1 = pos2screen(p1)   Dim As int2d posScrn2 = pos2screen(p2)   Line(posScrn1.x, posScrn1.y)-(posScrn2.x, posScrn2.y), cEnd SubSub scaled_graphics_type.drawRect(p1 As sgl2d, p2 As sgl2d, c As ULong)   Dim As int2d posScrn1 = pos2screen(p1)   Dim As int2d posScrn2 = pos2screen(p2)   Line(posScrn1.x, posScrn1.y)-(posScrn2.x, posScrn2.y), c, bEnd SubSub scaled_graphics_type.drawRectFilled(p1 As sgl2d, p2 As sgl2d, c As ULong)   Dim As int2d posScrn1 = pos2screen(p1)   Dim As int2d posScrn2 = pos2screen(p2)   Line(posScrn1.x, posScrn1.y)-(posScrn2.x, posScrn2.y), c, bfEnd Sub'===============================================================================Const As Single INFEC_DIST = 0.25 'mConst As Single INFEC_DIST_SQRD = INFEC_DIST ^ 2Const As Single SICK_DAYS_MIN = 7, SICK_DAYS_MAX = 21Const As Single MORTALITY = 3 / 100'const as single SICK_SPEED_FACTOR = 0.1 'disable due to some bugConst SEC_PER_DAY = 24 * 60 * 60Const As Single DAY_PER_SEC = 1 / SEC_PER_DAY Const PERSONS = 10000Const As Single V_MAX = 1 / 3600 '1 m/hrConst As Single MAP_X_MAX = +50, MAP_X_MIN = -MAP_X_MAXConst As Single MAP_Y_MAX = +35, MAP_Y_MIN = -MAP_Y_MAXConst As Single MAP_W = MAP_X_MAX - MAP_X_MINConst As Single MAP_H = MAP_Y_MAX - MAP_Y_MINConst DRAW_INTERVAL = 10 '1 to 10 Const FOLLOW_PATIENT_0 = FALSEConst NUM_SECT_X = 16 * 3, NUM_SECT_Y = 12 * 3Dim As grow_list sector(NUM_SECT_X - 1, NUM_SECT_Y -1)'var pSector = new grow_list[NUM_SECT_X, NUM_SECT_Y]'screen stuffConst As Single PPM = 10.0 'pixels per meter (set zoom level)Const SCRN_W = 1050, SCRN_H = 750Dim Shared As scaled_graphics_type sgsg.setScreen(SCRN_W, SCRN_H)sg.setScaling(PPM, sgl2d(0, 0))Enum HEALTH_STATE   ST_INIT '0   ST_INFEC '1   ST_RECOV '2   ST_DEAD '3   ST_LAST 'number of statesEnd EnumDim As Integer statCounters(ST_LAST - 1)Dim As ULong stateColor(ST_LAST - 1) = _   {RGB(0, 150, 0), RGB(150, 0, 0), RGB(150, 150, 0), RGB(0, 0, 0)}Type person_type   p As sgl2d 'position [m]   v As sgl2d 'velocity [m/s]   r As Single 'radius [m]   sickEndTime As Single 'not double!   state As Integer 'healthEnd TypeFunction drawUpdate(interval As Integer) As boolean   Static As Integer counter = 0   counter += 1   If counter >= interval Then      counter = 0      Return TRUE   Else      Return FALSE   End IfEnd FunctionFunction rndRangeSgl(min As Single, max As Single) As Single   Return Rnd() * (max - min) + minEnd FunctionSub limitInt(ByRef value As Integer, min As Integer, max As Integer)   If value < min Then value = min   If value > max Then value = maxEnd SubDim As Integer mouseEvent, mouseDragDim As mouse_type mouseDim As person_type person(0 To PERSONS-1)Randomize Timer'initialise personsFor i As Integer = 0 To UBound(person)   person(i).p.x = rndRangeSgl(MAP_X_MIN, MAP_X_MAX) 'm   person(i).p.y = rndRangeSgl(MAP_Y_MIN, MAP_Y_MAX) 'm   person(i).v.x = rndRangeSgl(-V_MAX, +V_MAX) 'm/s   person(i).v.y = rndRangeSgl(-V_MAX, +V_MAX) 'm/s   person(i).r = 0.25   'person(i).state = int(rndRangeSgl(0, csng(ST_LAST)))   person(i).state = ST_INITNext'time step such that a max speed, position change = 10 % of infection distanceDim As Single timeStep = (INFEC_DIST * 0.1) / V_MAX 'secDim As Single simulDays, simulTime = 0Dim As Integer numInit, numDead, numRecovDim As boolean sectioning = TRUEWhile Not MultiKey(FB.SC_ESCAPE)   'trigger outbreak   If MultiKey(FB.SC_SPACE) Then      person(0).state = ST_INFEC      person(0).sickEndTime = simulTime _         + rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY   End If   'zoom / drag view by mouse   mouseEvent = handleMouse(mouse)      If (mouse.buttons <> -1) Then      If (mouseEvent = MOUSE_LB_PRESSED) Then mouseDrag = 1      If (mouseEvent = MOUSE_LB_RELEASED) Then mouseDrag = 0      If (mouseEvent = MOUSE_WHEEl_UP) Then sg.scale *= 1.1      If (mouseEvent = MOUSE_WHEEl_DOWN) Then sg.scale /= 1.1   End If   If (mouseDrag) Then      sg.offset.x -= (mouse.posChange.x / PPM)      sg.offset.y += (mouse.posChange.y / PPM)   End If   'patient 0 perpective   If FOLLOW_PATIENT_0 = TRUE Then      sg.offset.x = (person(0).p.x)      sg.offset.y = (person(0).p.y)   End If   'update positions   For i As Integer = 0 To UBound(person)      person(i).p += person(i).v * timeStep   Next   'random direction change   '~ for i as integer = 1 to 5      '~ dim as integer iPerson = int(rndRangeSgl(0, PERSONS))      '~ if rnd() > 0.5 then         '~ person(iPerson).v.x = -person(iPerson).v.x      '~ else         '~ person(iPerson).v.y = -person(iPerson).v.y      '~ end if   '~ next   'check wall collisions   For i As Integer = 0 To UBound(person)      If person(i).p.x < MAP_X_MIN Then person(i).v.x = +Abs(person(i).v.x)      If person(i).p.x > MAP_X_MAX Then person(i).v.x = -abs(person(i).v.x)      If person(i).p.y < MAP_Y_MIN Then person(i).v.y = +Abs(person(i).v.y)      If person(i).p.y > MAP_Y_MAX Then person(i).v.y = -abs(person(i).v.y)   Next   'check end of sickness   For i As Integer = 0 To UBound(person)      If person(i).state = ST_INFEC Then         If simulTime > person(i).sickEndTime Then            If Rnd() < MORTALITY Then               person(i).state = ST_DEAD               person(i).v = sgl2d(0, 0)            Else               person(i).state = ST_RECOV               'person(i).v /= SICK_SPEED_FACTOR            End If         End If      End If   Next   'clear sectors   For xi As Integer = 0 To NUM_SECT_X - 1      For yi As Integer = 0 To NUM_SECT_Y - 1         sector(xi, yi).empty()      Next   Next   'assign persons to sectors   For i As Integer = 0 To UBound(person)      Dim As Integer xi = Int((person(i).p.x - MAP_X_MIN) / (MAP_W / NUM_SECT_X))      Dim As Integer yi = Int((person(i).p.y - MAP_Y_MIN) / (MAP_H / NUM_SECT_Y))      limitInt(xi, 0, NUM_SECT_X - 1)      limitInt(yi, 0, NUM_SECT_Y - 1)      sector(xi, yi).add(i)   Next   If sectioning = TRUE Then      For xi As Integer = 0 To NUM_SECT_X - 1         For yi As Integer = 0 To NUM_SECT_Y - 1            'loop source in 1 sector            For iiSrc As Integer = 0 To sector(xi, yi).numInUse - 1               Dim As Integer iSrc = sector(xi, yi).index(iiSrc)               If person(iSrc).state = ST_INFEC Then                  'check against targets in near sectors, including this sector                  For xdi As Integer = -1 To +1                     If xi + xdi < 0 Then Continue For                     If xi + xdi >= NUM_SECT_X Then Continue For                     For ydi As Integer = -1 To +1                        If yi + ydi < 0 Then Continue For                        If yi + ydi >= NUM_SECT_Y Then Continue For                        'loop targets with 1 (near) sector                        For iiTar As Integer = 0 To sector(xi + xdi, yi + ydi).numInUse - 1                           Dim As Integer iTar = sector(xi + xdi, yi + ydi).index(iiTar)                           If person(iTar).state = ST_INIT Then                              If person(iSrc).p.distSqrd(person(iTar).p) < INFEC_DIST_SQRD Then                                 person(iTar).state = ST_INFEC                                 person(iTar).sickEndTime = simulTime _                                    + rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY                                 'person(iTar).v *= SICK_SPEED_FACTOR                              End If                           End If                        Next                     Next                  Next               End If            Next         Next      Next   Else      'case: sectioning = FALSE      'check for disease transmission (iSrc = source, iTar = target)      For iSrc As Integer = 0 To UBound(person)         If person(iSrc).state = ST_INFEC Then            For iTar As Integer = 0 To UBound(person)               If person(iTar).state = ST_INIT Then                  If person(iSrc).p.distSqrd(person(iTar).p) < INFEC_DIST_SQRD Then                     person(iTar).state = ST_INFEC                     person(iTar).sickEndTime = simulTime _                        + rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY                     'person(iTar).v *= SICK_SPEED_FACTOR                  End If               End If            Next         End If      Next   End If   'clear stats   For i As Integer = 0 To UBound(statCounters)      statCounters(i) = 0   Next   'update stats   For i As Integer = 0 To UBound(person)      statCounters(person(i).state) += 1   Next   'draw world   If (drawUpdate(DRAW_INTERVAL) = TRUE) Then      ScreenLock      sg.clearScreen(RGB(250, 250, 250))      For xi As Integer = 0 To NUM_SECT_X - 1         For yi As Integer = 0 To NUM_SECT_Y - 1            Dim As Single x1 = MAP_W * (xi / NUM_SECT_X) + MAP_X_MIN            Dim As Single y1 = MAP_H * (yi / NUM_SECT_Y) + MAP_Y_MIN            Dim As Single x2 = MAP_W * ((xi + 1) / NUM_SECT_X) + MAP_X_MIN            Dim As Single y2 = MAP_H * ((yi + 1) / NUM_SECT_Y) + MAP_Y_MIN            sg.drawRect(sgl2d(x1, y1), sgl2d(x2, y2), RGB(0, 0, 0))            Dim As int2d scrnPos = sg.pos2screen(sgl2d(x1, y1))            Draw String (scrnPos.x + 2, scrnPos.y - 16), Str(sector(xi, yi).numInUse), 0            'draw string (scrnPos.x + 2, scrnPos.y - 32), str(sector(xi, yi).numAlloc), 0         Next      Next      sg.drawCircle(person(0).p, person(0).r * 3, stateColor(person(0).state))      For i As Integer = 0 To UBound(person)         Dim As ULong c = stateColor(person(i).state)         sg.drawCircleFilled(person(i).p, person(i).r, c)      Next      Draw String (10, 10), "Day: " & Format(simulDays, "0.#0"), RGB(0, 0, 0)      Draw String (10, 10 + 16), "Initial: " & statCounters(ST_INIT), stateColor(ST_INIT)      Draw String (10, 10 + 32), "Infected: " & statCounters(ST_INFEC), stateColor(ST_INFEC)      Draw String (10, 10 + 48), "Recovered: " & statCounters(ST_RECOV), stateColor(ST_RECOV)      Draw String (10, 10 + 64), "Dead: " & statCounters(ST_DEAD), stateColor(ST_DEAD)      ScreenUnLock      Sleep 1,1 'don't sleep every loop   End If   simulTime += timeStep   simulDays = simulTime * DAY_PER_SECWendEnd`

Pretty cool.

You and dodicat using often the operator within you code. That's apparently a good way to make things easier. I need to understand it first...
Posts: 2134
Joined: May 24, 2007 22:10
Location: The Netherlands

### Re: Corona virus simulator

I am running a comparison With and Without spacial partitioning (sectors in code).

I have added accumulative timers for simulation steps (excluding rendering):
step 1 = irrelevant (mouse and keyboard handling)
step 2 = update positions
step 3 = wall collisions
step 4 = end of sickness check
step 5 = irrelevant (clear sectors)
step 6 = assign persons to sectors
step 7 = virus transmission check
step 8 = irrelevant (update statistics)
step 9 = not a step, sum of timers above

Result with with spacial partitioning:

Result with without spacial partitioning:
(I let you know in an hour or more, going very slow, now at day 4 after 1400 seconds)
Done, that took some time (step 5 and 6 were not needed, but not significant):

The time from first outbreak to last patient recovered (or dead) is the same. This make me confident the the 'fast code' is correct.

Code: Select all

`#Include "fbgfx.bi"#Include "string.bi"Type int2d   As Integer x, y   Declare Constructor   Declare Constructor(x As Integer, y As Integer)   Declare Operator Cast () As StringEnd TypeConstructor int2dEnd ConstructorConstructor int2d(x As Integer, y As Integer)   This.x = x : This.y = yEnd Constructor' "x, y"Operator int2d.cast () As String  Return Str(x) & "," & Str(y)End OperatorOperator = (a As int2d, b As int2d) As boolean   If a.x <> b.x Then Return false   If a.y <> b.y Then Return false   Return trueEnd OperatorOperator <> (a As int2d, b As int2d) As boolean   If a.x = b.x And a.y = b.y Then Return false   Return trueEnd Operator' a + b Operator + (a As int2d, b As int2d) As int2d   Return Type(a.x + b.x, a.y + b.y)End Operator' a - bOperator - (a As int2d, b As int2d) As int2d   Return Type(a.x - b.x, a.y - b.y)End Operator' -aOperator - (a As int2d) As int2d   Return Type(-a.x, -a.y)End Operator' a * bOperator * (a As int2d, b As int2d) As int2d   Return Type(a.x * b.x, a.y * b.y)End Operator' a * mulOperator * (a As int2d, mul As Integer) As int2d   Return Type(a.x * mul, a.y * mul)End Operator' a \ bOperator \ (a As int2d, b As int2d) As int2d   Return Type(a.x \ b.x, a.y \ b.y)End Operator' a \ divOperator \ (a As int2d, divider As Integer) As int2d   Return Type(a.x \ divider, a.y \ divider)End Operator'===============================================================================Type sgl2d   As Single x, y   Declare Constructor   Declare Constructor(x As Single, y As Single)   Declare Operator Cast() As String   Declare Function cross(b As sgl2d) As Single   Declare Function lengthSqrd() As Single   Declare Function dist(b As sgl2d) As Single   Declare Function distSqrd(b As sgl2d) As Single   Declare Function normalise() As sgl2dEnd TypeConstructor sgl2dEnd ConstructorConstructor sgl2d(x As Single, y As Single)   This.x = x : This.y = yEnd ConstructorFunction sgl2d.cross(b As sgl2d) As Single   Return This.x * b.y - This.y * b.xEnd FunctionFunction sgl2d.lengthSqrd() As Single   Return (This.x * This.x) + (This.y * This.y) End FunctionFunction sgl2d.dist(b As sgl2d) As Single   Dim As Single dx = This.x - b.x   Dim As Single dy = This.y - b.y   Return Sqr((dx * dx) + (dy * dy)) End FunctionFunction sgl2d.distSqrd(b As sgl2d) As Single   Dim As Single dx = This.x - b.x   Dim As Single dy = This.y - b.y   Return (dx * dx) + (dy * dy) End FunctionFunction sgl2d.normalise() As sgl2d   Dim As Single length = Sqr((This.x * This.x) + (This.y * This.y))   Return sgl2d(This.x / length, This.y / length)End Function' "x, y"Operator sgl2d.cast() As String   Return Str(x) & "," & Str(y)End Operator'---- operators ---' distance / lenthOperator Len (a As sgl2d) As Single   Return Sqr(a.x * a.x + a.y * a.y)End Operator' a = b ?Operator = (a As sgl2d, b As sgl2d) As boolean   If a.x <> b.x Then Return false   If a.y <> b.y Then Return false   Return trueEnd Operator' a != b ?Operator <> (a As sgl2d, b As sgl2d) As boolean   If a.x = b.x And a.y = b.y Then Return false   Return trueEnd Operator' a + b Operator + (a As sgl2d, b As sgl2d) As sgl2d   Return Type(a.x + b.x, a.y + b.y)End Operator' a - bOperator - (a As sgl2d, b As sgl2d) As sgl2d   Return Type(a.x - b.x, a.y - b.y)End Operator' -aOperator - (a As sgl2d) As sgl2d   Return Type(-a.x, -a.y)End Operator' a * bOperator * (a As sgl2d, b As sgl2d) As sgl2d   Return Type(a.x * b.x, a.y * b.y)End Operator' a * mulOperator * (a As sgl2d, mul As Single) As sgl2d   Return Type(a.x * mul, a.y * mul)End Operator' a / divOperator / (a As sgl2d, div As Single) As sgl2d   Return Type(a.x / div, a.y / div)End Operator'==============================================================================='list can grow, but never shrink, for performance, non-sortedType grow_list   Dim As Integer index(Any)   Private:   Dim As Integer numItems   Public:   Declare Constructor(startSize As Integer)   Declare Constructor()   Declare Destructor()   Declare Sub Add(newIndex As Integer)   Declare Sub empty()   Declare Function numAlloc() As Integer   Declare Function numInUse() As Integer   Declare Sub show()   'non-list methodsEnd TypeConstructor grow_list(startSize As Integer)   If startSize > 0 Then      ReDim index(startSize - 1)   End IfEnd ConstructorConstructor grow_list()   This.constructor(0)End ConstructorDestructor grow_list()   Erase(index)End DestructorSub grow_list.add(newIndex As Integer)   Dim As Integer ub = UBound(index)   'if list is full, increase list size by 1   If numItems = ub + 1 Then      ReDim Preserve index(numItems)   End If   index(numItems) = newIndex   numItems += 1End SubSub grow_list.empty()   numItems = 0End SubFunction grow_list.numAlloc() As Integer   Return UBound(index) + 1End FunctionFunction grow_list.numInUse() As Integer   Return numItemsEnd Function'for debuggingSub grow_list.show()   Print "--- " & numInUse() & " / " & numAlloc() & " ---"   For i As Integer = 0 To numItems - 1      Print i, index(i)   NextEnd Sub'===============================================================================#Define MOUSE_IDLE 0#Define MOUSE_POS_CHANGED 1#Define MOUSE_LB_PRESSED 2#Define MOUSE_LB_RELEASED 3#Define MOUSE_RB_PRESSED 4#Define MOUSE_RB_RELEASED 5#Define MOUSE_MB_PRESSED 6#Define MOUSE_MB_RELEASED 7#Define MOUSE_WHEEL_UP 8#Define MOUSE_WHEEL_DOWN 9Type mouse_type   Pos As int2d   posChange As int2d   wheel As Integer   buttons As Integer   lb As Integer 'left button   rb As Integer 'right button   mb As Integer 'middle buttonEnd TypeFunction handleMouse(ByRef mouse As mouse_type) As Integer   Static previous As mouse_type   Dim As Integer change = MOUSE_IDLE   GetMouse mouse.pos.x, mouse.pos.y, mouse.wheel, mouse.buttons   If (mouse.buttons = -1) Then      mouse.lb = 0      mouse.rb = 0      mouse.mb = 0      mouse.posChange.x = 0      mouse.posChange.y = 0   Else      mouse.lb = (mouse.buttons And 1)      mouse.rb = (mouse.buttons ShR 1) And 1      mouse.mb = (mouse.buttons ShR 2) And 1      If (previous.pos.x <> mouse.pos.x Or previous.pos.y <> mouse.pos.y) Then         change = MOUSE_POS_CHANGED      End If      mouse.posChange.x = mouse.pos.x - previous.pos.x      mouse.posChange.y = mouse.pos.y - previous.pos.y      If (previous.buttons <> mouse.buttons) Then         If (previous.lb = 0 And mouse.lb = 1) Then change = MOUSE_LB_PRESSED         If (previous.lb = 1 And mouse.lb = 0) Then change = MOUSE_LB_RELEASED         If (previous.rb = 0 And mouse.rb = 1) Then change = MOUSE_RB_PRESSED         If (previous.rb = 1 And mouse.rb = 0) Then change = MOUSE_RB_RELEASED         If (previous.mb = 0 And mouse.mb = 1) Then change = MOUSE_MB_PRESSED         If (previous.mb = 1 And mouse.mb = 0) Then change = MOUSE_MB_RELEASED      End If      If (mouse.wheel > previous.wheel) Then change = MOUSE_WHEEl_UP      If (mouse.wheel < previous.wheel) Then change = MOUSE_WHEEl_DOWN      previous = mouse   End If   Return changeEnd Function'==============================================================================='Note: y+ = up, x+ = right, (0,0) = centerType scaled_graphics_type   Dim As Single scale = 1 ' = 1 / pixel_size 'pixels / meter   'dim as int2d offset' = (scrn_w \ 2, h \ 2) 'offset in pixels   Dim As sgl2d offset   Dim As Integer w = -1, h = -1   Dim As Integer wc = -1, hc = -1 'center x,y   Declare Sub setScreen(w As Integer, h As Integer)   Declare Sub setScaling(scale As Single, offset As sgl2d)   Declare Sub clearScreen(c As ULong)   Declare Function pos2screen(p As sgl2d) As int2d   Declare Sub drawPixel(p As sgl2d, c As ULong)   Declare Sub drawCircle(p As sgl2d, r As Single, c As ULong)   'declare sub drawCircleFilled(p as sgl2d, r as single, c as ulong, cFill as ulong)   Declare Sub drawCircleFilled(p As sgl2d, r As Single, c As ULong)   Declare Sub drawElipse(p As sgl2d, r As Single, aspect As Single, c As ULong)   Declare Sub drawLine(p1 As sgl2d, p2 As sgl2d, c As ULong)   Declare Sub drawRect(p1 As sgl2d, p2 As sgl2d, c As ULong)   Declare Sub drawRectFilled(p1 As sgl2d, p2 As sgl2d, c As ULong)End TypeSub scaled_graphics_type.setScreen(w As Integer, h As Integer)   This.w = w 'width   This.h = h 'height   wc = w \ 2   hc = h \ 2   ScreenRes w, h, 32   Width w \ 8, h \ 16 'bigger fontEnd SubSub scaled_graphics_type.setScaling(scale As Single, offset As sgl2d)   This.scale = scale   This.offset = offsetEnd SubSub scaled_graphics_type.clearScreen(c As ULong)   Line(0, 0)-(w - 1, h - 1), c, bfEnd SubFunction scaled_graphics_type.pos2screen(p As sgl2d) As int2d   Return int2d(Int(wc + (p.x - offset.x) * scale), h - Int(hc + (p.y - offset.y) * scale))End FunctionSub scaled_graphics_type.drawPixel(p As sgl2d, c As ULong)   Dim As int2d posScrn = pos2screen(p)   PSet(posScrn.x, posScrn.y), cEnd SubSub scaled_graphics_type.drawCircle(p As sgl2d, r As Single, c As ULong)   Dim As int2d posScrn = pos2screen(p)   Circle(posScrn.x, posScrn.y), r * scale, cEnd SubSub scaled_graphics_type.drawCircleFilled(p As sgl2d, r As Single, c As ULong)   Dim As int2d posScrn = pos2screen(p)   Circle(posScrn.x, posScrn.y), r * scale, c,,,,fEnd SubSub scaled_graphics_type.drawElipse(p As sgl2d, r As Single, aspect As Single, c As ULong)   Dim As int2d posScrn = pos2screen(p)   Circle(posScrn.x, posScrn.y), r * scale, c, , , aspectEnd SubSub scaled_graphics_type.drawLine(p1 As sgl2d, p2 As sgl2d, c As ULong)   Dim As int2d posScrn1 = pos2screen(p1)   Dim As int2d posScrn2 = pos2screen(p2)   Line(posScrn1.x, posScrn1.y)-(posScrn2.x, posScrn2.y), cEnd SubSub scaled_graphics_type.drawRect(p1 As sgl2d, p2 As sgl2d, c As ULong)   Dim As int2d posScrn1 = pos2screen(p1)   Dim As int2d posScrn2 = pos2screen(p2)   Line(posScrn1.x, posScrn1.y)-(posScrn2.x, posScrn2.y), c, bEnd SubSub scaled_graphics_type.drawRectFilled(p1 As sgl2d, p2 As sgl2d, c As ULong)   Dim As int2d posScrn1 = pos2screen(p1)   Dim As int2d posScrn2 = pos2screen(p2)   Line(posScrn1.x, posScrn1.y)-(posScrn2.x, posScrn2.y), c, bfEnd Sub'===============================================================================Const As Single INFEC_DIST = 0.25 'mConst As Single INFEC_DIST_SQRD = INFEC_DIST ^ 2Const As Single SICK_DAYS_MIN = 7, SICK_DAYS_MAX = 21Const As Single MORTALITY = 3 / 100'const as single SICK_SPEED_FACTOR = 0.1 'disable due to some bugConst SEC_PER_DAY = 24 * 60 * 60Const As Single DAY_PER_SEC = 1 / SEC_PER_DAY Const PERSONS = 10000Const As Single V_MAX = 1 / 3600 '1 m/hrConst As Single MAP_X_MAX = +50, MAP_X_MIN = -MAP_X_MAXConst As Single MAP_Y_MAX = +35, MAP_Y_MIN = -MAP_Y_MAXConst As Single MAP_W = MAP_X_MAX - MAP_X_MINConst As Single MAP_H = MAP_Y_MAX - MAP_Y_MINConst DRAW_INTERVAL = 10 '1 to 10 Const FOLLOW_PATIENT_0 = FALSEConst NUM_SECT_X = 16 * 3, NUM_SECT_Y = 12 * 3Dim As grow_list sector(NUM_SECT_X - 1, NUM_SECT_Y -1)'var pSector = new grow_list[NUM_SECT_X, NUM_SECT_Y]'screen stuffConst As Single PPM = 10.0 'pixels per meter (set zoom level)Const SCRN_W = 1050, SCRN_H = 750Dim Shared As scaled_graphics_type sgsg.setScreen(SCRN_W, SCRN_H)sg.setScaling(PPM, sgl2d(0, 0))Enum HEALTH_STATE   ST_INIT '0   ST_INFEC '1   ST_RECOV '2   ST_DEAD '3   ST_LAST 'number of statesEnd EnumDim As Integer statCounters(ST_LAST - 1)Dim As ULong stateColor(ST_LAST - 1) = _   {RGB(0, 150, 0), RGB(150, 0, 0), RGB(150, 150, 0), RGB(0, 0, 0)}Type person_type   p As sgl2d 'position [m]   v As sgl2d 'velocity [m/s]   r As Single 'radius [m]   sickEndTime As Single 'not double!   state As Integer 'healthEnd TypeFunction drawUpdate(interval As Integer) As boolean   Static As Integer counter = 0   counter += 1   If counter >= interval Then      counter = 0      Return TRUE   Else      Return FALSE   End IfEnd FunctionFunction rndRangeSgl(min As Single, max As Single) As Single   Return Rnd() * (max - min) + minEnd FunctionSub limitInt(ByRef value As Integer, min As Integer, max As Integer)   If value < min Then value = min   If value > max Then value = maxEnd SubDim As Integer mouseEvent, mouseDragDim As mouse_type mouseDim As person_type person(0 To PERSONS-1)Randomize 1234'initialise personsFor i As Integer = 0 To UBound(person)   person(i).p.x = rndRangeSgl(MAP_X_MIN, MAP_X_MAX) 'm   person(i).p.y = rndRangeSgl(MAP_Y_MIN, MAP_Y_MAX) 'm   person(i).v.x = rndRangeSgl(-V_MAX, +V_MAX) 'm/s   person(i).v.y = rndRangeSgl(-V_MAX, +V_MAX) 'm/s   person(i).r = 0.25   'person(i).state = int(rndRangeSgl(0, csng(ST_LAST)))   person(i).state = ST_INITNext'time step such that a max speed, position change = 10 % of infection distanceDim As Single timeStep = (INFEC_DIST * 0.1) / V_MAX 'secDim As Single simulDays, simulTime = 0Dim As Integer numInit, numDead, numRecovDim As boolean sectioning = TRUEsg.offset.x = -15 'move view for numberDim As Double tNow(0 To 9), tAcc(0 To 9)Dim As Integer trigger = 1While Not MultiKey(FB.SC_ESCAPE)   tNow(0) = Timer   'trigger outbreak   If MultiKey(FB.SC_SPACE) Then trigger = 1   If trigger = 1 Then      trigger = 0      person(0).state = ST_INFEC      person(0).sickEndTime = simulTime _         + rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY   End If   'zoom / drag view by mouse   mouseEvent = handleMouse(mouse)      If (mouse.buttons <> -1) Then      If (mouseEvent = MOUSE_LB_PRESSED) Then mouseDrag = 1      If (mouseEvent = MOUSE_LB_RELEASED) Then mouseDrag = 0      If (mouseEvent = MOUSE_WHEEl_UP) Then sg.scale *= 1.1      If (mouseEvent = MOUSE_WHEEl_DOWN) Then sg.scale /= 1.1   End If   If (mouseDrag) Then      sg.offset.x -= (mouse.posChange.x / PPM)      sg.offset.y += (mouse.posChange.y / PPM)   End If   'patient 0 perpective   If FOLLOW_PATIENT_0 = TRUE Then      sg.offset.x = (person(0).p.x)      sg.offset.y = (person(0).p.y)   End If   tNow(1) = Timer   'update positions   For i As Integer = 0 To UBound(person)      person(i).p += person(i).v * timeStep   Next   'random direction change   '~ for i as integer = 1 to 5      '~ dim as integer iPerson = int(rndRangeSgl(0, PERSONS))      '~ if rnd() > 0.5 then         '~ person(iPerson).v.x = -person(iPerson).v.x      '~ else         '~ person(iPerson).v.y = -person(iPerson).v.y      '~ end if   '~ next   tNow(2) = Timer   'check wall collisions   For i As Integer = 0 To UBound(person)      If person(i).p.x < MAP_X_MIN Then person(i).v.x = +Abs(person(i).v.x)      If person(i).p.x > MAP_X_MAX Then person(i).v.x = -abs(person(i).v.x)      If person(i).p.y < MAP_Y_MIN Then person(i).v.y = +Abs(person(i).v.y)      If person(i).p.y > MAP_Y_MAX Then person(i).v.y = -abs(person(i).v.y)   Next   tNow(3) = Timer   'check end of sickness   For i As Integer = 0 To UBound(person)      If person(i).state = ST_INFEC Then         If simulTime > person(i).sickEndTime Then            If Rnd() < MORTALITY Then               person(i).state = ST_DEAD               person(i).v = sgl2d(0, 0)            Else               person(i).state = ST_RECOV               'person(i).v /= SICK_SPEED_FACTOR            End If         End If      End If   Next      tNow(4) = Timer   'clear sectors   For xi As Integer = 0 To NUM_SECT_X - 1      For yi As Integer = 0 To NUM_SECT_Y - 1         sector(xi, yi).empty()      Next   Next      tNow(5) = Timer   'assign persons to sectors   For i As Integer = 0 To UBound(person)      Dim As Integer xi = Int((person(i).p.x - MAP_X_MIN) / (MAP_W / NUM_SECT_X))      Dim As Integer yi = Int((person(i).p.y - MAP_Y_MIN) / (MAP_H / NUM_SECT_Y))      limitInt(xi, 0, NUM_SECT_X - 1)      limitInt(yi, 0, NUM_SECT_Y - 1)      sector(xi, yi).add(i)   Next   tNow(6) = Timer   'spread the virus   If sectioning = TRUE Then      For xi As Integer = 0 To NUM_SECT_X - 1         For yi As Integer = 0 To NUM_SECT_Y - 1            'loop source in 1 sector            For iiSrc As Integer = 0 To sector(xi, yi).numInUse - 1               Dim As Integer iSrc = sector(xi, yi).index(iiSrc)               If person(iSrc).state = ST_INFEC Then                  'check against targets in near sectors, including this sector                  For xdi As Integer = -1 To +1                     If xi + xdi < 0 Then Continue For                     If xi + xdi >= NUM_SECT_X Then Continue For                     For ydi As Integer = -1 To +1                        If yi + ydi < 0 Then Continue For                        If yi + ydi >= NUM_SECT_Y Then Continue For                        'loop targets with 1 (near) sector                        For iiTar As Integer = 0 To sector(xi + xdi, yi + ydi).numInUse - 1                           Dim As Integer iTar = sector(xi + xdi, yi + ydi).index(iiTar)                           If person(iTar).state = ST_INIT Then                              If person(iSrc).p.distSqrd(person(iTar).p) < INFEC_DIST_SQRD Then                                 person(iTar).state = ST_INFEC                                 person(iTar).sickEndTime = simulTime _                                    + rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY                                 'person(iTar).v *= SICK_SPEED_FACTOR                              End If                           End If                        Next                     Next                  Next               End If            Next         Next      Next   Else      'case: sectioning = FALSE      'check for disease transmission (iSrc = source, iTar = target)      For iSrc As Integer = 0 To UBound(person)         If person(iSrc).state = ST_INFEC Then            For iTar As Integer = 0 To UBound(person)               If person(iTar).state = ST_INIT Then                  If person(iSrc).p.distSqrd(person(iTar).p) < INFEC_DIST_SQRD Then                     person(iTar).state = ST_INFEC                     person(iTar).sickEndTime = simulTime _                        + rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY                     'person(iTar).v *= SICK_SPEED_FACTOR                  End If               End If            Next         End If      Next   End If   tNow(7) = Timer   'clear stats   For i As Integer = 0 To UBound(statCounters)      statCounters(i) = 0   Next   'update stats   For i As Integer = 0 To UBound(person)      statCounters(person(i).state) += 1   Next      tNow(8) = Timer   'draw world   If (drawUpdate(DRAW_INTERVAL) = TRUE) Then      ScreenLock      sg.clearScreen(RGB(250, 250, 250))      For xi As Integer = 0 To NUM_SECT_X - 1         For yi As Integer = 0 To NUM_SECT_Y - 1            Dim As Single x1 = MAP_W * (xi / NUM_SECT_X) + MAP_X_MIN            Dim As Single y1 = MAP_H * (yi / NUM_SECT_Y) + MAP_Y_MIN            Dim As Single x2 = MAP_W * ((xi + 1) / NUM_SECT_X) + MAP_X_MIN            Dim As Single y2 = MAP_H * ((yi + 1) / NUM_SECT_Y) + MAP_Y_MIN            sg.drawRect(sgl2d(x1, y1), sgl2d(x2, y2), RGB(0, 0, 0))            Dim As int2d scrnPos = sg.pos2screen(sgl2d(x1, y1))            Draw String (scrnPos.x + 2, scrnPos.y - 16), Str(sector(xi, yi).numInUse), 0            'draw string (scrnPos.x + 2, scrnPos.y - 32), str(sector(xi, yi).numAlloc), 0         Next      Next      sg.drawCircle(person(0).p, person(0).r * 3, stateColor(person(0).state))      For i As Integer = 0 To UBound(person)         Dim As ULong c = stateColor(person(i).state)         sg.drawCircleFilled(person(i).p, person(i).r, c)      Next      Draw String (10, 10), "Day: " & Format(simulDays, "0.#0"), RGB(0, 0, 0)      Draw String (10, 10 + 16), "Initial: " & statCounters(ST_INIT), stateColor(ST_INIT)      Draw String (10, 10 + 32), "Infected: " & statCounters(ST_INFEC), stateColor(ST_INFEC)      Draw String (10, 10 + 48), "Recovered: " & statCounters(ST_RECOV), stateColor(ST_RECOV)      Draw String (10, 10 + 64), "Dead: " & statCounters(ST_DEAD), stateColor(ST_DEAD)      For i As Integer = 1 To 9         Draw String (10, 90 + i * 16), "Step: " & i & " = " &  Format(tAcc(i), "0.#0"), RGB(0, 0, 0)      Next      ScreenUnLock      Sleep 1,1 'don't sleep every loop   End If   simulTime += timeStep   simulDays = simulTime * DAY_PER_SEC   tAcc(9) = 0 'sum other timers   For i As Integer = 1 To 8      tAcc(i) += (tNow(i) - tNow(i - 1))      tAcc(9) += tAcc(i) 'sum   Next   If statCounters(ST_INFEC) = 0 Then Exit WhileWendLocate 1,1: Print "End of simulation!" While InKey = "" : Sleep 1 : WendEnd`

Next steps:
* Find optimal sector / grid size
* Make time graphs
* Apply social distancing measures

Edit: Current grid of 48 x 36 seems close to optimal. Larger grid is really slower, smaller grid does not give much improvement.
Posts: 2134
Joined: May 24, 2007 22:10
Location: The Netherlands

### Re: Corona virus simulator

Simulation of virus containment failure:

Same conditions without borders:

Code: Select all

`#Include "fbgfx.bi"#Include "string.bi"Type int2d   As Integer x, y   Declare Constructor   Declare Constructor(x As Integer, y As Integer)   Declare Operator Cast () As StringEnd TypeConstructor int2dEnd ConstructorConstructor int2d(x As Integer, y As Integer)   This.x = x : This.y = yEnd Constructor' "x, y"Operator int2d.cast () As String  Return Str(x) & "," & Str(y)End OperatorOperator = (a As int2d, b As int2d) As boolean   If a.x <> b.x Then Return false   If a.y <> b.y Then Return false   Return trueEnd OperatorOperator <> (a As int2d, b As int2d) As boolean   If a.x = b.x And a.y = b.y Then Return false   Return trueEnd Operator' a + b Operator + (a As int2d, b As int2d) As int2d   Return Type(a.x + b.x, a.y + b.y)End Operator' a - bOperator - (a As int2d, b As int2d) As int2d   Return Type(a.x - b.x, a.y - b.y)End Operator' -aOperator - (a As int2d) As int2d   Return Type(-a.x, -a.y)End Operator' a * bOperator * (a As int2d, b As int2d) As int2d   Return Type(a.x * b.x, a.y * b.y)End Operator' a * mulOperator * (a As int2d, mul As Integer) As int2d   Return Type(a.x * mul, a.y * mul)End Operator' a \ bOperator \ (a As int2d, b As int2d) As int2d   Return Type(a.x \ b.x, a.y \ b.y)End Operator' a \ divOperator \ (a As int2d, divider As Integer) As int2d   Return Type(a.x \ divider, a.y \ divider)End Operator'===============================================================================Type sgl2d   As Single x, y   Declare Constructor   Declare Constructor(x As Single, y As Single)   Declare Operator Cast() As String   Declare Function cross(b As sgl2d) As Single   Declare Function lengthSqrd() As Single   Declare Function dist(b As sgl2d) As Single   Declare Function distSqrd(b As sgl2d) As Single   Declare Function normalise() As sgl2dEnd TypeConstructor sgl2dEnd ConstructorConstructor sgl2d(x As Single, y As Single)   This.x = x : This.y = yEnd ConstructorFunction sgl2d.cross(b As sgl2d) As Single   Return This.x * b.y - This.y * b.xEnd FunctionFunction sgl2d.lengthSqrd() As Single   Return (This.x * This.x) + (This.y * This.y) End FunctionFunction sgl2d.dist(b As sgl2d) As Single   Dim As Single dx = This.x - b.x   Dim As Single dy = This.y - b.y   Return Sqr((dx * dx) + (dy * dy)) End FunctionFunction sgl2d.distSqrd(b As sgl2d) As Single   Dim As Single dx = This.x - b.x   Dim As Single dy = This.y - b.y   Return (dx * dx) + (dy * dy) End FunctionFunction sgl2d.normalise() As sgl2d   Dim As Single length = Sqr((This.x * This.x) + (This.y * This.y))   Return sgl2d(This.x / length, This.y / length)End Function' "x, y"Operator sgl2d.cast() As String   Return Str(x) & "," & Str(y)End Operator'---- operators ---' distance / lenthOperator Len (a As sgl2d) As Single   Return Sqr(a.x * a.x + a.y * a.y)End Operator' a = b ?Operator = (a As sgl2d, b As sgl2d) As boolean   If a.x <> b.x Then Return false   If a.y <> b.y Then Return false   Return trueEnd Operator' a != b ?Operator <> (a As sgl2d, b As sgl2d) As boolean   If a.x = b.x And a.y = b.y Then Return false   Return trueEnd Operator' a + b Operator + (a As sgl2d, b As sgl2d) As sgl2d   Return Type(a.x + b.x, a.y + b.y)End Operator' a - bOperator - (a As sgl2d, b As sgl2d) As sgl2d   Return Type(a.x - b.x, a.y - b.y)End Operator' -aOperator - (a As sgl2d) As sgl2d   Return Type(-a.x, -a.y)End Operator' a * bOperator * (a As sgl2d, b As sgl2d) As sgl2d   Return Type(a.x * b.x, a.y * b.y)End Operator' a * mulOperator * (a As sgl2d, mul As Single) As sgl2d   Return Type(a.x * mul, a.y * mul)End Operator' a / divOperator / (a As sgl2d, div As Single) As sgl2d   Return Type(a.x / div, a.y / div)End Operator'==============================================================================='list can grow, but never shrink, for performance, non-sortedType grow_list   Dim As Integer index(Any)   Private:   Dim As Integer numItems   Public:   Declare Constructor(startSize As Integer)   Declare Constructor()   Declare Destructor()   Declare Sub Add(newIndex As Integer)   Declare Sub empty()   Declare Function numAlloc() As Integer   Declare Function numInUse() As Integer   Declare Sub show()   'non-list methodsEnd TypeConstructor grow_list(startSize As Integer)   If startSize > 0 Then      ReDim index(startSize - 1)   End IfEnd ConstructorConstructor grow_list()   This.constructor(0)End ConstructorDestructor grow_list()   Erase(index)End DestructorSub grow_list.add(newIndex As Integer)   Dim As Integer ub = UBound(index)   'if list is full, increase list size by 1   If numItems = ub + 1 Then      ReDim Preserve index(numItems)   End If   index(numItems) = newIndex   numItems += 1End SubSub grow_list.empty()   numItems = 0End SubFunction grow_list.numAlloc() As Integer   Return UBound(index) + 1End FunctionFunction grow_list.numInUse() As Integer   Return numItemsEnd Function'for debuggingSub grow_list.show()   Print "--- " & numInUse() & " / " & numAlloc() & " ---"   For i As Integer = 0 To numItems - 1      Print i, index(i)   NextEnd Sub'===============================================================================#Define MOUSE_IDLE 0#Define MOUSE_POS_CHANGED 1#Define MOUSE_LB_PRESSED 2#Define MOUSE_LB_RELEASED 3#Define MOUSE_RB_PRESSED 4#Define MOUSE_RB_RELEASED 5#Define MOUSE_MB_PRESSED 6#Define MOUSE_MB_RELEASED 7#Define MOUSE_WHEEL_UP 8#Define MOUSE_WHEEL_DOWN 9Type mouse_type   Pos As int2d   posChange As int2d   wheel As Integer   buttons As Integer   lb As Integer 'left button   rb As Integer 'right button   mb As Integer 'middle buttonEnd TypeFunction handleMouse(ByRef mouse As mouse_type) As Integer   Static previous As mouse_type   Dim As Integer change = MOUSE_IDLE   GetMouse mouse.pos.x, mouse.pos.y, mouse.wheel, mouse.buttons   If (mouse.buttons = -1) Then      mouse.lb = 0      mouse.rb = 0      mouse.mb = 0      mouse.posChange.x = 0      mouse.posChange.y = 0   Else      mouse.lb = (mouse.buttons And 1)      mouse.rb = (mouse.buttons ShR 1) And 1      mouse.mb = (mouse.buttons ShR 2) And 1      If (previous.pos.x <> mouse.pos.x Or previous.pos.y <> mouse.pos.y) Then         change = MOUSE_POS_CHANGED      End If      mouse.posChange.x = mouse.pos.x - previous.pos.x      mouse.posChange.y = mouse.pos.y - previous.pos.y      If (previous.buttons <> mouse.buttons) Then         If (previous.lb = 0 And mouse.lb = 1) Then change = MOUSE_LB_PRESSED         If (previous.lb = 1 And mouse.lb = 0) Then change = MOUSE_LB_RELEASED         If (previous.rb = 0 And mouse.rb = 1) Then change = MOUSE_RB_PRESSED         If (previous.rb = 1 And mouse.rb = 0) Then change = MOUSE_RB_RELEASED         If (previous.mb = 0 And mouse.mb = 1) Then change = MOUSE_MB_PRESSED         If (previous.mb = 1 And mouse.mb = 0) Then change = MOUSE_MB_RELEASED      End If      If (mouse.wheel > previous.wheel) Then change = MOUSE_WHEEl_UP      If (mouse.wheel < previous.wheel) Then change = MOUSE_WHEEl_DOWN      previous = mouse   End If   Return changeEnd Function'==============================================================================='Note: y+ = up, x+ = right, (0,0) = centerType scaled_graphics_type   Dim As Single scale = 1 ' = 1 / pixel_size 'pixels / meter   'dim as int2d offset' = (scrn_w \ 2, h \ 2) 'offset in pixels   Dim As sgl2d offset   Dim As Integer w = -1, h = -1   Dim As Integer wc = -1, hc = -1 'center x,y   Declare Sub setScreen(w As Integer, h As Integer)   Declare Sub setScaling(scale As Single, offset As sgl2d)   Declare Sub clearScreen(c As ULong)   Declare Function pos2screen(p As sgl2d) As int2d   Declare Sub drawPixel(p As sgl2d, c As ULong)   Declare Sub drawCircle(p As sgl2d, r As Single, c As ULong)   'declare sub drawCircleFilled(p as sgl2d, r as single, c as ulong, cFill as ulong)   Declare Sub drawCircleFilled(p As sgl2d, r As Single, c As ULong)   Declare Sub drawElipse(p As sgl2d, r As Single, aspect As Single, c As ULong)   Declare Sub drawLine(p1 As sgl2d, p2 As sgl2d, c As ULong)   Declare Sub drawRect(p1 As sgl2d, p2 As sgl2d, c As ULong)   Declare Sub drawRectFilled(p1 As sgl2d, p2 As sgl2d, c As ULong)End TypeSub scaled_graphics_type.setScreen(w As Integer, h As Integer)   This.w = w 'width   This.h = h 'height   wc = w \ 2   hc = h \ 2   ScreenRes w, h, 32   Width w \ 8, h \ 16 'bigger fontEnd SubSub scaled_graphics_type.setScaling(scale As Single, offset As sgl2d)   This.scale = scale   This.offset = offsetEnd SubSub scaled_graphics_type.clearScreen(c As ULong)   Line(0, 0)-(w - 1, h - 1), c, bfEnd SubFunction scaled_graphics_type.pos2screen(p As sgl2d) As int2d   Return int2d(Int(wc + (p.x - offset.x) * scale), h - Int(hc + (p.y - offset.y) * scale))End FunctionSub scaled_graphics_type.drawPixel(p As sgl2d, c As ULong)   Dim As int2d posScrn = pos2screen(p)   PSet(posScrn.x, posScrn.y), cEnd SubSub scaled_graphics_type.drawCircle(p As sgl2d, r As Single, c As ULong)   Dim As int2d posScrn = pos2screen(p)   Circle(posScrn.x, posScrn.y), r * scale, cEnd SubSub scaled_graphics_type.drawCircleFilled(p As sgl2d, r As Single, c As ULong)   Dim As int2d posScrn = pos2screen(p)   Circle(posScrn.x, posScrn.y), r * scale, c,,,,fEnd SubSub scaled_graphics_type.drawElipse(p As sgl2d, r As Single, aspect As Single, c As ULong)   Dim As int2d posScrn = pos2screen(p)   Circle(posScrn.x, posScrn.y), r * scale, c, , , aspectEnd SubSub scaled_graphics_type.drawLine(p1 As sgl2d, p2 As sgl2d, c As ULong)   Dim As int2d posScrn1 = pos2screen(p1)   Dim As int2d posScrn2 = pos2screen(p2)   Line(posScrn1.x, posScrn1.y)-(posScrn2.x, posScrn2.y), cEnd SubSub scaled_graphics_type.drawRect(p1 As sgl2d, p2 As sgl2d, c As ULong)   Dim As int2d posScrn1 = pos2screen(p1)   Dim As int2d posScrn2 = pos2screen(p2)   Line(posScrn1.x, posScrn1.y)-(posScrn2.x, posScrn2.y), c, bEnd SubSub scaled_graphics_type.drawRectFilled(p1 As sgl2d, p2 As sgl2d, c As ULong)   Dim As int2d posScrn1 = pos2screen(p1)   Dim As int2d posScrn2 = pos2screen(p2)   Line(posScrn1.x, posScrn1.y)-(posScrn2.x, posScrn2.y), c, bfEnd Sub'===============================================================================Const As Single INFEC_DIST = 0.05 'mConst As Single INFEC_DIST_SQRD = INFEC_DIST ^ 2Const As Single SICK_DAYS_MIN = 7, SICK_DAYS_MAX = 21Const As Single MORTALITY = 3 / 100'const as single SICK_SPEED_FACTOR = 0.25 'disable due to some bugConst SEC_PER_DAY = 24 * 60 * 60Const As Single DAY_PER_SEC = 1 / SEC_PER_DAY Const PERSONS = 5000Const As Single V_MAX = 0.1 / 3600 '1 m/hrConst As Single MAP_X_MAX = +50, MAP_X_MIN = -MAP_X_MAXConst As Single MAP_Y_MAX = +35, MAP_Y_MIN = -MAP_Y_MAXConst As Single MAP_W = MAP_X_MAX - MAP_X_MINConst As Single MAP_H = MAP_Y_MAX - MAP_Y_MINConst DRAW_INTERVAL = 10 '1 to 10 Const FOLLOW_PATIENT_0 = FALSEConst NUM_SECT_X = 16 * 3, NUM_SECT_Y = 12 * 3Dim As grow_list sector(NUM_SECT_X - 1, NUM_SECT_Y -1)'screen stuffConst As Single PPM = 6.0 'pixels per meter (set zoom level)Const SCRN_W = 1050, SCRN_H = 750Dim Shared As scaled_graphics_type sgsg.setScreen(SCRN_W, SCRN_H)sg.setScaling(PPM, sgl2d(0, 0))Enum HEALTH_STATE   ST_INIT '0   ST_INFEC '1   ST_RECOV '2   ST_DEAD '3   ST_LAST 'number of statesEnd EnumDim As Integer statCounters(ST_LAST - 1)Dim As ULong stateColor(ST_LAST - 1) = _   {RGB(0, 150, 0), RGB(150, 0, 0), RGB(150, 150, 0), RGB(0, 0, 0)}Dim As Single statPercentage(ST_LAST - 1)Const As ULong WHITE = RGB(250, 250, 250)Const As ULong BLACK = RGB(0, 0, 0)Type person_type   p As sgl2d 'position [m]   pPrev As sgl2d 'position [m]   v As sgl2d 'velocity [m/s]   r As Single 'radius [m]   sickEndTime As Single 'not double!   state As Integer 'healthEnd TypeFunction drawUpdate(interval As Integer) As boolean   Static As Integer counter = 0   counter += 1   If counter >= interval Then      counter = 0      Return TRUE   Else      Return FALSE   End IfEnd FunctionFunction rndRangeSgl(min As Single, max As Single) As Single   Return Rnd() * (max - min) + minEnd FunctionSub limitInt(ByRef value As Integer, min As Integer, max As Integer)   If value < min Then value = min   If value > max Then value = maxEnd SubDim As Integer mouseEvent, mouseDragDim As mouse_type mouseDim As person_type person(0 To PERSONS-1)Randomize 1234'initialise personsFor i As Integer = 0 To UBound(person)   person(i).p.x = rndRangeSgl(MAP_X_MIN, MAP_X_MAX) 'm   person(i).p.y = rndRangeSgl(MAP_Y_MIN, MAP_Y_MAX) 'm   person(i).pPrev = person(i).p   person(i).v.x = rndRangeSgl(-V_MAX, +V_MAX) 'm/s   person(i).v.y = rndRangeSgl(-V_MAX, +V_MAX) 'm/s   person(i).r = 0.35   'person(i).state = int(rndRangeSgl(0, csng(ST_LAST)))   person(i).state = ST_INITNext'time step such that a max speed, position change = 20 % of infection distanceDim As Single timeStep = (INFEC_DIST * 0.2) / V_MAX 'secDim As Single simulDays, simulTime = 0Dim As Integer numInit, numDead, numRecovDim As boolean sectioning = TRUE'move viewsg.offset.x = -35sg.offset.y = -25Dim As Double tNow(0 To 9), tAcc(0 To 9)Dim As Integer trigger = 1sg.clearScreen(WHITE)While Not MultiKey(FB.SC_ESCAPE)   tNow(0) = Timer   'trigger outbreak   If MultiKey(FB.SC_SPACE) Then trigger = 1   If trigger = 1 Then      trigger = 0      person(0).state = ST_INFEC      person(0).sickEndTime = simulTime _         + rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY   End If   'zoom / drag view by mouse   mouseEvent = handleMouse(mouse)      If (mouse.buttons <> -1) Then      If (mouseEvent = MOUSE_LB_PRESSED) Then mouseDrag = 1      If (mouseEvent = MOUSE_LB_RELEASED) Then mouseDrag = 0      If (mouseEvent = MOUSE_WHEEl_UP) Then sg.scale *= 1.1      If (mouseEvent = MOUSE_WHEEl_DOWN) Then sg.scale /= 1.1   End If   If (mouseDrag) Then      sg.offset.x -= (mouse.posChange.x / PPM)      sg.offset.y += (mouse.posChange.y / PPM)   End If   'patient 0 perpective   If FOLLOW_PATIENT_0 = TRUE Then      sg.offset.x = (person(0).p.x)      sg.offset.y = (person(0).p.y)   End If   tNow(1) = Timer   'update positions   For i As Integer = 0 To UBound(person)      person(i).p += person(i).v * timeStep   Next   'random direction change   For i As Integer = 1 To 1      Dim As Integer iPerson = Int(rndRangeSgl(0, PERSONS))      If Rnd() > 0.5 Then         person(iPerson).v.x = -person(iPerson).v.x      Else         person(iPerson).v.y = -person(iPerson).v.y      End If   Next   tNow(2) = Timer   'check wall collisions   For i As Integer = 0 To UBound(person)      With person(i)         If .p.x < MAP_X_MIN Then .v.x = +Abs(.v.x)         If .p.x > MAP_X_MAX Then .v.x = -abs(.v.x)         If .p.y < MAP_Y_MIN Then .v.y = +Abs(.v.y)         If .p.y > MAP_Y_MAX Then .v.y = -abs(.v.y)         If .pPrev.x > 0 And .p.x < 0 Then .v.x = +Abs(.v.x)         If .pPrev.x < 0 And .p.x > 0 Then .v.x = -abs(.v.x)         If .pPrev.y > 0 And .p.y < 0 Then .v.y = +Abs(.v.y)         If .pPrev.y < 0 And .p.y > 0 Then .v.y = -abs(.v.y)      End With   Next   tNow(3) = Timer   'check end of sickness   For i As Integer = 0 To UBound(person)      If person(i).state = ST_INFEC Then         If simulTime > person(i).sickEndTime Then            If Rnd() < MORTALITY Then               person(i).state = ST_DEAD               person(i).v = sgl2d(0, 0)            Else               person(i).state = ST_RECOV               'person(i).v /= SICK_SPEED_FACTOR 'somthing wrong with this               'person(i).v.x = rndRangeSgl(-V_MAX, +V_MAX) 'm/s               'person(i).v.y = rndRangeSgl(-V_MAX, +V_MAX) 'm/s            End If         End If      End If   Next      tNow(4) = Timer   'clear sectors   For xi As Integer = 0 To NUM_SECT_X - 1      For yi As Integer = 0 To NUM_SECT_Y - 1         sector(xi, yi).empty()      Next   Next      tNow(5) = Timer   'assign persons to sectors   For i As Integer = 0 To UBound(person)      Dim As Integer xi = Int((person(i).p.x - MAP_X_MIN) / (MAP_W / NUM_SECT_X))      Dim As Integer yi = Int((person(i).p.y - MAP_Y_MIN) / (MAP_H / NUM_SECT_Y))      limitInt(xi, 0, NUM_SECT_X - 1)      limitInt(yi, 0, NUM_SECT_Y - 1)      sector(xi, yi).add(i)   Next   tNow(6) = Timer   'spread the virus   If sectioning = TRUE Then      For xi As Integer = 0 To NUM_SECT_X - 1         For yi As Integer = 0 To NUM_SECT_Y - 1            'loop source in 1 sector            For iiSrc As Integer = 0 To sector(xi, yi).numInUse - 1               Dim As Integer iSrc = sector(xi, yi).index(iiSrc)               If person(iSrc).state = ST_INFEC Then                  'check against targets in near sectors, including this sector                  For xdi As Integer = -1 To +1                     If xi + xdi < 0 Then Continue For                     If xi + xdi >= NUM_SECT_X Then Continue For                     For ydi As Integer = -1 To +1                        If yi + ydi < 0 Then Continue For                        If yi + ydi >= NUM_SECT_Y Then Continue For                        'loop targets with 1 (near) sector                        For iiTar As Integer = 0 To sector(xi + xdi, yi + ydi).numInUse - 1                           Dim As Integer iTar = sector(xi + xdi, yi + ydi).index(iiTar)                           If person(iTar).state = ST_INIT Then                              If person(iSrc).p.distSqrd(person(iTar).p) < INFEC_DIST_SQRD Then                                 person(iTar).state = ST_INFEC                                 person(iTar).sickEndTime = simulTime _                                    + rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY                                 'person(iTar).v *= SICK_SPEED_FACTOR                              End If                           End If                        Next                     Next                  Next               End If            Next         Next      Next   Else      'case: sectioning = FALSE      'check for disease transmission (iSrc = source, iTar = target)      For iSrc As Integer = 0 To UBound(person)         If person(iSrc).state = ST_INFEC Then            For iTar As Integer = 0 To UBound(person)               If person(iTar).state = ST_INIT Then                  If person(iSrc).p.distSqrd(person(iTar).p) < INFEC_DIST_SQRD Then                     person(iTar).state = ST_INFEC                     person(iTar).sickEndTime = simulTime _                        + rndRangeSgl(SICK_DAYS_MIN, SICK_DAYS_MAX) * SEC_PER_DAY                     'person(iTar).v *= SICK_SPEED_FACTOR                  End If               End If            Next         End If      Next   End If   tNow(7) = Timer   'clear stats   For i As Integer = 0 To UBound(statCounters)      statCounters(i) = 0   Next   'update stats   For i As Integer = 0 To UBound(person)      statCounters(person(i).state) += 1   Next      tNow(8) = Timer   'draw world   If (drawUpdate(DRAW_INTERVAL) = TRUE) Then      ScreenLock      'sg.clearScreen(WHITE)      'clear part of screen      sg.drawRectFilled(sgl2d(MAP_X_MIN - 2, MAP_Y_MIN - 2), _         sgl2d(MAP_X_MAX + 2, MAP_Y_MAX + 2), WHITE)      sg.drawLine(sgl2d(MAP_X_MIN, 0), sgl2d(MAP_X_MAX, 0), BLACK)      sg.drawLine(sgl2d(0, MAP_Y_MIN), sgl2d(0, MAP_Y_MAX), BLACK)      For xi As Integer = 0 To NUM_SECT_X - 1         For yi As Integer = 0 To NUM_SECT_Y - 1            Dim As Single x1 = MAP_W * (xi / NUM_SECT_X) + MAP_X_MIN            Dim As Single y1 = MAP_H * (yi / NUM_SECT_Y) + MAP_Y_MIN            Dim As Single x2 = MAP_W * ((xi + 1) / NUM_SECT_X) + MAP_X_MIN            Dim As Single y2 = MAP_H * ((yi + 1) / NUM_SECT_Y) + MAP_Y_MIN            'sg.drawRect(sgl2d(x1, y1), sgl2d(x2, y2), BLACK) 'grid            'dim as int2d scrnPos = sg.pos2screen(sgl2d(x1, y1))            'draw string (scrnPos.x + 2, scrnPos.y - 16), str(sector(xi, yi).numInUse), 0            'draw string (scrnPos.x + 2, scrnPos.y - 32), str(sector(xi, yi).numAlloc), 0         Next      Next      sg.drawCircle(person(0).p, person(0).r * 3, stateColor(person(0).state))      For i As Integer = 0 To UBound(person)         Dim As ULong c = stateColor(person(i).state)         sg.drawCircleFilled(person(i).p, person(i).r, c)      Next      Line(0, 0)-(150, 100), WHITE, bf      Draw String (10, 10 + 00), "Day: " & Format(simulDays, "0.#0"), BLACK      Draw String (10, 10 + 16), "Initial: " & statCounters(ST_INIT), stateColor(ST_INIT)      Draw String (10, 10 + 32), "Infected: " & statCounters(ST_INFEC), stateColor(ST_INFEC)      Draw String (10, 10 + 48), "Recovered: " & statCounters(ST_RECOV), stateColor(ST_RECOV)      Draw String (10, 10 + 64), "Dead: " & statCounters(ST_DEAD), stateColor(ST_DEAD)      '~ for i as integer = 1 to 9         '~ draw string (10, 90 + i * 16), "Step: " & i & " = " &  format(tAcc(i), "0.#0"), BLACK      '~ next      For i As Integer = 0 To ST_LAST - 1         statPercentage(i) = 100 * (statCounters(i) / PERSONS)         Circle(10 + 4 * simulDays, (SCRN_H - 10) - 3 * statPercentage(i)), 1, stateColor(i),,,,f      Next      ScreenUnLock      Sleep 1,1 'don't sleep every loop   End If   simulTime += timeStep   simulDays = simulTime * DAY_PER_SEC   tAcc(9) = 0 'sum other timers   For i As Integer = 1 To 8      tAcc(i) += (tNow(i) - tNow(i - 1))      tAcc(9) += tAcc(i) 'sum   Next   If statCounters(ST_INFEC) = 0 Then Exit WhileWendDraw String (10, SCRN_W - 26), "End of simulation", BLACKWhile InKey = "" : Sleep 1 : WendEnd`
Posts: 2134
Joined: May 24, 2007 22:10
Location: The Netherlands

### Re: Corona virus simulator

A different version. No more fancy math, all integers and random walk. And no collision detection. I set a location on the map as infected when an infected person is at that spot. When a healthy (not immune) person is also at this spot, the disease is transmitted. Map is cleaned each cycle. Simulation is with 2500 persons, but 10 times more is no problem (switch display to pixel instead of circle). The daily movement of persons per day is very limited. I think that this shows that although 1 person is 'only' sick for 14 days, the total duration of the outbreak lasts much longer.

Code: Select all

`const SCRN_W = 1000, SCRN_H = 700screenres SCRN_W, SCRN_H, 32width SCRN_W \ 18, SCRN_H \ 16const ST_INIT = 0, ST_INFEC = 1, ST_RECOV = 2, ST_DEAD = 3, ST_INVALID = 4const SICK_DAYS = 14const as single MORTALITY = 3 / 100const as ulong WHITE = rgb(200, 200, 200), BLACK = rgb(0, 0, 0)dim as ulong stateColor(ST_INIT to ST_DEAD) = _   {rgb(0, 200, 0), rgb(250, 0, 0), rgb(200, 200, 0), rgb(200, 200, 200)}dim as integer statCounters(ST_INIT to ST_DEAD)dim as string statNames(ST_INIT to ST_DEAD) = _   {"Initial", "Infected", "Recovered", "Dead"}const SEC_MIN = 60, SEC_HR = SEC_MIN * 60, SEC_DAY = SEC_HR * 24type person_type   as integer x, y   as integer state   as integer sickEndTimeend typeconst PERSONS = 2500dim as person_type person(any)redim person(PERSONS - 1)'map array dynamic, else too large for stackconst MAP_X = 700, MAP_Y = MAP_Xdim as integer map(any, any)redim map(MAP_X - 1, MAP_Y - 1)sub clearMap(map() as integer)   clear map(0), 0, MAP_Y * MAP_X   'set map borders   for xi as integer = 0 to MAP_X - 1      map(xi, 0) = ST_INVALID      map(xi, MAP_Y - 1) = ST_INVALID   next   for yi as integer = 0 to MAP_Y - 1      map(0, yi) = ST_INVALID      map(MAP_X - 1, yi) = ST_INVALID   nextend subclearMap(map())const RND_DEF = 0, RND_CRT = 1, RND_FST = 2, RND_MTW = 3, RND_QBC = 4, RND_SYS = 5randomize 222, RND_FSTfor i as integer = 0 to PERSONS - 1   person(i).x = int(rnd * (MAP_X - 2)) + 1   person(i).y = int(rnd * (MAP_Y - 2)) + 1nextperson(0).state = ST_INFECperson(0).sickEndTime = SICK_DAYS * SEC_DAYdim as integer quit = 0dim as integer simTime = 0, simTimeStep = SEC_MIN * 10while quit = 0   if inkey = chr(27) then quit = 1   'move persons   for i as integer = 0 to PERSONS - 1      'dead don't walk      if person(i).state < ST_DEAD then         select case int(rnd * 4)         case 0 'try right            if map(person(i).x + 1, person(i).y) <> ST_INVALID then person(i).x += 1         case 1 'try left            if map(person(i).x - 1, person(i).y) <> ST_INVALID then person(i).x -= 1         case 2 'try down            if map(person(i).x, person(i).y + 1) <> ST_INVALID then person(i).y += 1         case 3 'try up            if map(person(i).x, person(i).y - 1) <> ST_INVALID then person(i).y -= 1         end select         'set infected locations         if person(i).state = ST_INFEC then map(person(i).x, person(i).y) = ST_INFEC      end if   next   'set new infections   for i as integer = 0 to PERSONS - 1      if person(i).state = ST_INIT then         if map(person(i).x, person(i).y) = ST_INFEC then            person(i).state = ST_INFEC            person(i).sickEndTime = simTime + SICK_DAYS * SEC_DAY         end if      end if   next   'reset map & check end of sickness   for i as integer = 0 to PERSONS - 1      if person(i).state = ST_INFEC then         map(person(i).x, person(i).y) = 0         if simTime > person(i).sickEndTime then            person(i).state = iif(rnd() < MORTALITY, ST_DEAD, ST_RECOV)         end if      end if   next   'clearMap(map())   'clear stats   for i as integer = 0 to ubound(statCounters)      statCounters(i) = 0   next   'update stats   for i as integer = 0 to ubound(person)      statCounters(person(i).state) += 1   next   'draw   screenlock   line(0, 0) - (SCRN_W - 1, SCRN_H - 1), BLACK, bf   for i as integer = 0 to PERSONS - 1      'line(150 + person(i).x * 2, person(i).y * 2)-step(1, 1), stateColor(person(i).state), bf      'pset(150 + person(i).x, person(i).y), stateColor(person(i).state)      circle(150 + person(i).x, person(i).y), 3, stateColor(person(i).state),,,,f   next   locate 1, 1: print "days: " & (simTime \ SEC_DAY)   for i as integer = 0 to ubound(statCounters)      locate 2 + i, 1: print statNames(i) & ": " & statCounters(i)   next   screenunlock      simTime += simTimeStep   sleep 1   if statCounters(ST_INFEC) = 0 then quit = 1wendsleep`
ShawnLG
Posts: 137
Joined: Dec 25, 2008 20:21

### Re: Corona virus simulator

I have updated my CoronaSim with long distance vacation travel, simulate communal movement and use hashing to improve the simulator's time complexity from O(n^2) to O(n). Code is updated in the first post. It is like watching the Game of Life but with virus VS humans.
Posts: 2134
Joined: May 24, 2007 22:10
Location: The Netherlands

### Re: Corona virus simulator

ShawnLG wrote:I have updated my CoronaSim with long distance vacation travel, simulate communal movement and use hashing to improve the simulator's time complexity from O(n^2) to O(n). Code is updated in the first post. It is like watching the Game of Life but with virus VS humans.

I have some weird lines flashing on the screen, what are those? Oh, that is the long distance travel.
ShawnLG
Posts: 137
Joined: Dec 25, 2008 20:21

### Re: Corona virus simulator

badidea wrote:I have some weird lines flashing on the screen, what are those? Oh, that is the long distance travel.

Yes they are.