Corona virus simulator

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

Re: Corona virus simulator

Postby Tourist Trap » Mar 29, 2020 0:39

badidea wrote:After exactly 350 seconds it return to initial state.

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

Re: Corona virus simulator

Postby UEZ » Mar 29, 2020 0:41

@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

Postby BasicScience » Mar 29, 2020 1:06

@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

Postby BasicScience » Mar 29, 2020 1:13

@Badidea - mesmerizing.

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

Postby BasicScience » Mar 29, 2020 5:54

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   c
END TYPE

DIM as integer N_people = 2, iPsize = 5

ScreenRes 640, 480, 32

DIM Shared as __person person(N_people-1)
Randomize

'Initialize
Person(0).x_init = 100
person(0).y_init = 100
person(0).Vx = 2*(rnd-.5)
person(0).Vy = 2*(rnd-.5)
person(0).c      = &hFFFFFFFF
person(0).sqradius = 2500

person(1).x_init = 300
person(1).y_init = 200
Person(1).Vx = (rnd-.5)
Person(1).Vy = (rnd-.5)
Person(1).c      = &hFFFFFFFF
Person(1).sqradius = 6400

For 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 i

do until inkey = "" : loop

DIM as integer cnt = 0

Do
    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 5
Loop while inkey = ""

sleep
badidea
Posts: 2000
Joined: May 24, 2007 22:10
Location: The Netherlands

Re: Corona virus simulator

Postby badidea » Mar 29, 2020 23:03

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 String
End Type

Constructor int2d
End Constructor

Constructor int2d(x As Integer, y As Integer)
   This.x = x : This.y = y
End Constructor

' "x, y"
Operator int2d.cast () As String
  Return Str(x) & "," & Str(y)
End Operator

Operator = (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 true
End Operator

Operator <> (a As int2d, b As int2d) As boolean
   If a.x = b.x And a.y = b.y Then Return false
   Return true
End 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 - b
Operator - (a As int2d, b As int2d) As int2d
   Return Type(a.x - b.x, a.y - b.y)
End Operator

' -a
Operator - (a As int2d) As int2d
   Return Type(-a.x, -a.y)
End 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 * mul
Operator * (a As int2d, mul As Integer) As int2d
   Return Type(a.x * mul, a.y * mul)
End 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 \ div
Operator \ (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 sgl2d
End Type

Constructor sgl2d
End Constructor

Constructor sgl2d(x As Single, y As Single)
   This.x = x : This.y = y
End Constructor

Function sgl2d.cross(b As sgl2d) As Single
   Return This.x * b.y - This.y * b.x
End Function

Function sgl2d.lengthSqrd() As Single
   Return (This.x * This.x) + (This.y * This.y)
End Function

Function 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 Function

Function 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 Function

Function 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 / lenth
Operator 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 true
End 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 true
End 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 - b
Operator - (a As sgl2d, b As sgl2d) As sgl2d
   Return Type(a.x - b.x, a.y - b.y)
End Operator

' -a
Operator - (a As sgl2d) As sgl2d
   Return Type(-a.x, -a.y)
End 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 * mul
Operator * (a As sgl2d, mul As Single) As sgl2d
   Return Type(a.x * mul, a.y * mul)
End Operator

' a / div
Operator / (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 9

Type 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 button
End Type

Function 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 change
End Function

'===============================================================================

'Note: y+ = up, x+ = right, (0,0) = center
Type 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 Type

Sub 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 font
End Sub

Sub scaled_graphics_type.setScaling(scale As Single, offset As sgl2d)
   This.scale = scale
   This.offset = offset
End Sub

Sub scaled_graphics_type.clearScreen(c As ULong)
   Line(0, 0)-(w - 1, h - 1), c, bf
End Sub

Function 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 Function

Sub scaled_graphics_type.drawPixel(p As sgl2d, c As ULong)
   Dim As int2d posScrn = pos2screen(p)
   PSet(posScrn.x, posScrn.y), c
End Sub

Sub 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, c
End Sub

Sub 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,,,,f
End Sub

Sub 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, , , aspect
End Sub

Sub 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), c
End Sub

'===============================================================================

Const As Single INFEC_DIST = 0.20 'm
Const As Single INFEC_DIST_SQRD = INFEC_DIST ^ 2
Const As Single SICK_DAYS_MIN = 7, SICK_DAYS_MAX = 21
Const As Single MORTALITY = 3 / 100
'const as single SICK_SPEED_FACTOR = 0.1 'disable due to some bug

Const SEC_PER_DAY = 24 * 60 * 60
Const As Single DAY_PER_SEC = 1 / SEC_PER_DAY

Const PERSONS = 500
Const As Single V_MAX = 1 / 3600 '1 m/hr
Const As Single MAP_X_MAX = +50, MAP_X_MIN = -MAP_X_MAX
Const As Single MAP_Y_MAX = +35, MAP_Y_MIN = -MAP_Y_MAX

Const DRAW_INTERVAL = 1 '1 to 10
Const FOLLOW_PATIENT_0 = TRUE 'FALSE

'screen stuff
Const As Single PPM = 20.0 'pixels per meter (set zoom level)
Const SCRN_W = 1050, SCRN_H = 750
Dim Shared As scaled_graphics_type sg
sg.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 states
End Enum

Dim 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 'health
End Type

Function 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 If
End Function

Function rndRangeSgl(min As Single, max As Single) As Single
   Return Rnd() * (max - min) + min
End Function

Dim As Integer mouseEvent, mouseDrag
Dim As mouse_type mouse
Dim As person_type person(0 To PERSONS-1)

Randomize Timer
'initialise persons
For 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_INIT
Next

'time step such that a max speed, position change = 10 % of infection distance
Dim As Single timeStep = (INFEC_DIST * 0.1) / V_MAX 'sec
Dim As Single simulDays, simulTime = 0
Dim As Integer numInit, numDead, numRecov

While 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_SEC
Wend
End
albert
Posts: 5634
Joined: Sep 28, 2006 2:41
Location: California, USA

Re: Corona virus simulator

Postby albert » Mar 29, 2020 23:15

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..
badidea
Posts: 2000
Joined: May 24, 2007 22:10
Location: The Netherlands

Re: Corona virus simulator

Postby badidea » Mar 30, 2020 23:36

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 String
End Type

Constructor int2d
End Constructor

Constructor int2d(x As Integer, y As Integer)
   This.x = x : This.y = y
End Constructor

' "x, y"
Operator int2d.cast () As String
  Return Str(x) & "," & Str(y)
End Operator

Operator = (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 true
End Operator

Operator <> (a As int2d, b As int2d) As boolean
   If a.x = b.x And a.y = b.y Then Return false
   Return true
End 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 - b
Operator - (a As int2d, b As int2d) As int2d
   Return Type(a.x - b.x, a.y - b.y)
End Operator

' -a
Operator - (a As int2d) As int2d
   Return Type(-a.x, -a.y)
End 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 * mul
Operator * (a As int2d, mul As Integer) As int2d
   Return Type(a.x * mul, a.y * mul)
End 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 \ div
Operator \ (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 sgl2d
End Type

Constructor sgl2d
End Constructor

Constructor sgl2d(x As Single, y As Single)
   This.x = x : This.y = y
End Constructor

Function sgl2d.cross(b As sgl2d) As Single
   Return This.x * b.y - This.y * b.x
End Function

Function sgl2d.lengthSqrd() As Single
   Return (This.x * This.x) + (This.y * This.y)
End Function

Function 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 Function

Function 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 Function

Function 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 / lenth
Operator 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 true
End 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 true
End 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 - b
Operator - (a As sgl2d, b As sgl2d) As sgl2d
   Return Type(a.x - b.x, a.y - b.y)
End Operator

' -a
Operator - (a As sgl2d) As sgl2d
   Return Type(-a.x, -a.y)
End 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 * mul
Operator * (a As sgl2d, mul As Single) As sgl2d
   Return Type(a.x * mul, a.y * mul)
End Operator

' a / div
Operator / (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-sorted
Type 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 methods
End Type

Constructor grow_list(startSize As Integer)
   If startSize > 0 Then
      ReDim index(startSize - 1)
   End If
End Constructor

Constructor grow_list()
   This.constructor(0)
End Constructor

Destructor grow_list()
   Erase(index)
End Destructor

Sub 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 += 1
End Sub

Sub grow_list.empty()
   numItems = 0
End Sub

Function grow_list.numAlloc() As Integer
   Return UBound(index) + 1
End Function

Function grow_list.numInUse() As Integer
   Return numItems
End Function

'for debugging
Sub grow_list.show()
   Print "--- " & numInUse() & " / " & numAlloc() & " ---"
   For i As Integer = 0 To numItems - 1
      Print i, index(i)
   Next
End 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 9

Type 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 button
End Type

Function 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 change
End Function

'===============================================================================

'Note: y+ = up, x+ = right, (0,0) = center
Type 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 Type

Sub 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 font
End Sub

Sub scaled_graphics_type.setScaling(scale As Single, offset As sgl2d)
   This.scale = scale
   This.offset = offset
End Sub

Sub scaled_graphics_type.clearScreen(c As ULong)
   Line(0, 0)-(w - 1, h - 1), c, bf
End Sub

Function 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 Function

Sub scaled_graphics_type.drawPixel(p As sgl2d, c As ULong)
   Dim As int2d posScrn = pos2screen(p)
   PSet(posScrn.x, posScrn.y), c
End Sub

Sub 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, c
End Sub

Sub 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,,,,f
End Sub

Sub 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, , , aspect
End Sub

Sub 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), c
End Sub

Sub 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, b
End Sub

Sub 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, bf
End Sub

'===============================================================================

Const As Single INFEC_DIST = 0.25 'm
Const As Single INFEC_DIST_SQRD = INFEC_DIST ^ 2
Const As Single SICK_DAYS_MIN = 7, SICK_DAYS_MAX = 21
Const As Single MORTALITY = 3 / 100
'const as single SICK_SPEED_FACTOR = 0.1 'disable due to some bug

Const SEC_PER_DAY = 24 * 60 * 60
Const As Single DAY_PER_SEC = 1 / SEC_PER_DAY

Const PERSONS = 10000
Const As Single V_MAX = 1 / 3600 '1 m/hr
Const As Single MAP_X_MAX = +50, MAP_X_MIN = -MAP_X_MAX
Const As Single MAP_Y_MAX = +35, MAP_Y_MIN = -MAP_Y_MAX
Const As Single MAP_W = MAP_X_MAX - MAP_X_MIN
Const As Single MAP_H = MAP_Y_MAX - MAP_Y_MIN

Const DRAW_INTERVAL = 10 '1 to 10
Const FOLLOW_PATIENT_0 = FALSE

Const NUM_SECT_X = 16 * 3, NUM_SECT_Y = 12 * 3
Dim As grow_list sector(NUM_SECT_X - 1, NUM_SECT_Y -1)
'var pSector = new grow_list[NUM_SECT_X, NUM_SECT_Y]

'screen stuff
Const As Single PPM = 10.0 'pixels per meter (set zoom level)
Const SCRN_W = 1050, SCRN_H = 750
Dim Shared As scaled_graphics_type sg
sg.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 states
End Enum

Dim 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 'health
End Type

Function 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 If
End Function

Function rndRangeSgl(min As Single, max As Single) As Single
   Return Rnd() * (max - min) + min
End Function

Sub limitInt(ByRef value As Integer, min As Integer, max As Integer)
   If value < min Then value = min
   If value > max Then value = max
End Sub

Dim As Integer mouseEvent, mouseDrag
Dim As mouse_type mouse
Dim As person_type person(0 To PERSONS-1)

Randomize Timer
'initialise persons
For 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_INIT
Next

'time step such that a max speed, position change = 10 % of infection distance
Dim As Single timeStep = (INFEC_DIST * 0.1) / V_MAX 'sec
Dim As Single simulDays, simulTime = 0
Dim As Integer numInit, numDead, numRecov
Dim As boolean sectioning = TRUE

While 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_SEC
Wend
End
UEZ
Posts: 519
Joined: May 05, 2017 19:59
Location: Germany

Re: Corona virus simulator

Postby UEZ » Mar 31, 2020 6:14

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 String
End Type

Constructor int2d
End Constructor

Constructor int2d(x As Integer, y As Integer)
   This.x = x : This.y = y
End Constructor

' "x, y"
Operator int2d.cast () As String
  Return Str(x) & "," & Str(y)
End Operator

Operator = (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 true
End Operator

Operator <> (a As int2d, b As int2d) As boolean
   If a.x = b.x And a.y = b.y Then Return false
   Return true
End 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 - b
Operator - (a As int2d, b As int2d) As int2d
   Return Type(a.x - b.x, a.y - b.y)
End Operator

' -a
Operator - (a As int2d) As int2d
   Return Type(-a.x, -a.y)
End 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 * mul
Operator * (a As int2d, mul As Integer) As int2d
   Return Type(a.x * mul, a.y * mul)
End 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 \ div
Operator \ (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 sgl2d
End Type

Constructor sgl2d
End Constructor

Constructor sgl2d(x As Single, y As Single)
   This.x = x : This.y = y
End Constructor

Function sgl2d.cross(b As sgl2d) As Single
   Return This.x * b.y - This.y * b.x
End Function

Function sgl2d.lengthSqrd() As Single
   Return (This.x * This.x) + (This.y * This.y)
End Function

Function 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 Function

Function 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 Function

Function 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 / lenth
Operator 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 true
End 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 true
End 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 - b
Operator - (a As sgl2d, b As sgl2d) As sgl2d
   Return Type(a.x - b.x, a.y - b.y)
End Operator

' -a
Operator - (a As sgl2d) As sgl2d
   Return Type(-a.x, -a.y)
End 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 * mul
Operator * (a As sgl2d, mul As Single) As sgl2d
   Return Type(a.x * mul, a.y * mul)
End Operator

' a / div
Operator / (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-sorted
Type 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 methods
End Type

Constructor grow_list(startSize As Integer)
   If startSize > 0 Then
      ReDim index(startSize - 1)
   End If
End Constructor

Constructor grow_list()
   This.constructor(0)
End Constructor

Destructor grow_list()
   Erase(index)
End Destructor

Sub 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 += 1
End Sub

Sub grow_list.empty()
   numItems = 0
End Sub

Function grow_list.numAlloc() As Integer
   Return UBound(index) + 1
End Function

Function grow_list.numInUse() As Integer
   Return numItems
End Function

'for debugging
Sub grow_list.show()
   Print "--- " & numInUse() & " / " & numAlloc() & " ---"
   For i As Integer = 0 To numItems - 1
      Print i, index(i)
   Next
End 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 9

Type 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 button
End Type

Function 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 change
End Function

'===============================================================================

'Note: y+ = up, x+ = right, (0,0) = center
Type 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 Type

Sub 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 font
End Sub

Sub scaled_graphics_type.setScaling(scale As Single, offset As sgl2d)
   This.scale = scale
   This.offset = offset
End Sub

Sub scaled_graphics_type.clearScreen(c As ULong)
   Line(0, 0)-(w - 1, h - 1), c, bf
End Sub

Function 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 Function

Sub scaled_graphics_type.drawPixel(p As sgl2d, c As ULong)
   Dim As int2d posScrn = pos2screen(p)
   PSet(posScrn.x, posScrn.y), c
End Sub

Sub 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, c
End Sub

Sub 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,,,,f
End Sub

Sub 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, , , aspect
End Sub

Sub 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), c
End Sub

Sub 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, b
End Sub

Sub 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, bf
End Sub

'===============================================================================

Const As Single INFEC_DIST = 0.25 'm
Const As Single INFEC_DIST_SQRD = INFEC_DIST ^ 2
Const As Single SICK_DAYS_MIN = 7, SICK_DAYS_MAX = 21
Const As Single MORTALITY = 3 / 100
'const as single SICK_SPEED_FACTOR = 0.1 'disable due to some bug

Const SEC_PER_DAY = 24 * 60 * 60
Const As Single DAY_PER_SEC = 1 / SEC_PER_DAY

Const PERSONS = 10000
Const As Single V_MAX = 1 / 3600 '1 m/hr
Const As Single MAP_X_MAX = +50, MAP_X_MIN = -MAP_X_MAX
Const As Single MAP_Y_MAX = +35, MAP_Y_MIN = -MAP_Y_MAX
Const As Single MAP_W = MAP_X_MAX - MAP_X_MIN
Const As Single MAP_H = MAP_Y_MAX - MAP_Y_MIN

Const DRAW_INTERVAL = 10 '1 to 10
Const FOLLOW_PATIENT_0 = FALSE

Const NUM_SECT_X = 16 * 3, NUM_SECT_Y = 12 * 3
Dim As grow_list sector(NUM_SECT_X - 1, NUM_SECT_Y -1)
'var pSector = new grow_list[NUM_SECT_X, NUM_SECT_Y]

'screen stuff
Const As Single PPM = 10.0 'pixels per meter (set zoom level)
Const SCRN_W = 1050, SCRN_H = 750
Dim Shared As scaled_graphics_type sg
sg.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 states
End Enum

Dim 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 'health
End Type

Function 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 If
End Function

Function rndRangeSgl(min As Single, max As Single) As Single
   Return Rnd() * (max - min) + min
End Function

Sub limitInt(ByRef value As Integer, min As Integer, max As Integer)
   If value < min Then value = min
   If value > max Then value = max
End Sub

Dim As Integer mouseEvent, mouseDrag
Dim As mouse_type mouse
Dim As person_type person(0 To PERSONS-1)

Randomize Timer
'initialise persons
For 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_INIT
Next

'time step such that a max speed, position change = 10 % of infection distance
Dim As Single timeStep = (INFEC_DIST * 0.1) / V_MAX 'sec
Dim As Single simulDays, simulTime = 0
Dim As Integer numInit, numDead, numRecov
Dim As boolean sectioning = TRUE

While 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_SEC
Wend
End


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...
badidea
Posts: 2000
Joined: May 24, 2007 22:10
Location: The Netherlands

Re: Corona virus simulator

Postby badidea » Mar 31, 2020 16:11

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:
Image

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):
Image

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 String
End Type

Constructor int2d
End Constructor

Constructor int2d(x As Integer, y As Integer)
   This.x = x : This.y = y
End Constructor

' "x, y"
Operator int2d.cast () As String
  Return Str(x) & "," & Str(y)
End Operator

Operator = (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 true
End Operator

Operator <> (a As int2d, b As int2d) As boolean
   If a.x = b.x And a.y = b.y Then Return false
   Return true
End 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 - b
Operator - (a As int2d, b As int2d) As int2d
   Return Type(a.x - b.x, a.y - b.y)
End Operator

' -a
Operator - (a As int2d) As int2d
   Return Type(-a.x, -a.y)
End 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 * mul
Operator * (a As int2d, mul As Integer) As int2d
   Return Type(a.x * mul, a.y * mul)
End 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 \ div
Operator \ (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 sgl2d
End Type

Constructor sgl2d
End Constructor

Constructor sgl2d(x As Single, y As Single)
   This.x = x : This.y = y
End Constructor

Function sgl2d.cross(b As sgl2d) As Single
   Return This.x * b.y - This.y * b.x
End Function

Function sgl2d.lengthSqrd() As Single
   Return (This.x * This.x) + (This.y * This.y)
End Function

Function 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 Function

Function 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 Function

Function 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 / lenth
Operator 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 true
End 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 true
End 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 - b
Operator - (a As sgl2d, b As sgl2d) As sgl2d
   Return Type(a.x - b.x, a.y - b.y)
End Operator

' -a
Operator - (a As sgl2d) As sgl2d
   Return Type(-a.x, -a.y)
End 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 * mul
Operator * (a As sgl2d, mul As Single) As sgl2d
   Return Type(a.x * mul, a.y * mul)
End Operator

' a / div
Operator / (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-sorted
Type 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 methods
End Type

Constructor grow_list(startSize As Integer)
   If startSize > 0 Then
      ReDim index(startSize - 1)
   End If
End Constructor

Constructor grow_list()
   This.constructor(0)
End Constructor

Destructor grow_list()
   Erase(index)
End Destructor

Sub 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 += 1
End Sub

Sub grow_list.empty()
   numItems = 0
End Sub

Function grow_list.numAlloc() As Integer
   Return UBound(index) + 1
End Function

Function grow_list.numInUse() As Integer
   Return numItems
End Function

'for debugging
Sub grow_list.show()
   Print "--- " & numInUse() & " / " & numAlloc() & " ---"
   For i As Integer = 0 To numItems - 1
      Print i, index(i)
   Next
End 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 9

Type 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 button
End Type

Function 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 change
End Function

'===============================================================================

'Note: y+ = up, x+ = right, (0,0) = center
Type 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 Type

Sub 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 font
End Sub

Sub scaled_graphics_type.setScaling(scale As Single, offset As sgl2d)
   This.scale = scale
   This.offset = offset
End Sub

Sub scaled_graphics_type.clearScreen(c As ULong)
   Line(0, 0)-(w - 1, h - 1), c, bf
End Sub

Function 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 Function

Sub scaled_graphics_type.drawPixel(p As sgl2d, c As ULong)
   Dim As int2d posScrn = pos2screen(p)
   PSet(posScrn.x, posScrn.y), c
End Sub

Sub 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, c
End Sub

Sub 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,,,,f
End Sub

Sub 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, , , aspect
End Sub

Sub 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), c
End Sub

Sub 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, b
End Sub

Sub 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, bf
End Sub

'===============================================================================

Const As Single INFEC_DIST = 0.25 'm
Const As Single INFEC_DIST_SQRD = INFEC_DIST ^ 2
Const As Single SICK_DAYS_MIN = 7, SICK_DAYS_MAX = 21
Const As Single MORTALITY = 3 / 100
'const as single SICK_SPEED_FACTOR = 0.1 'disable due to some bug

Const SEC_PER_DAY = 24 * 60 * 60
Const As Single DAY_PER_SEC = 1 / SEC_PER_DAY

Const PERSONS = 10000
Const As Single V_MAX = 1 / 3600 '1 m/hr
Const As Single MAP_X_MAX = +50, MAP_X_MIN = -MAP_X_MAX
Const As Single MAP_Y_MAX = +35, MAP_Y_MIN = -MAP_Y_MAX
Const As Single MAP_W = MAP_X_MAX - MAP_X_MIN
Const As Single MAP_H = MAP_Y_MAX - MAP_Y_MIN

Const DRAW_INTERVAL = 10 '1 to 10
Const FOLLOW_PATIENT_0 = FALSE

Const NUM_SECT_X = 16 * 3, NUM_SECT_Y = 12 * 3
Dim As grow_list sector(NUM_SECT_X - 1, NUM_SECT_Y -1)
'var pSector = new grow_list[NUM_SECT_X, NUM_SECT_Y]

'screen stuff
Const As Single PPM = 10.0 'pixels per meter (set zoom level)
Const SCRN_W = 1050, SCRN_H = 750
Dim Shared As scaled_graphics_type sg
sg.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 states
End Enum

Dim 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 'health
End Type

Function 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 If
End Function

Function rndRangeSgl(min As Single, max As Single) As Single
   Return Rnd() * (max - min) + min
End Function

Sub limitInt(ByRef value As Integer, min As Integer, max As Integer)
   If value < min Then value = min
   If value > max Then value = max
End Sub

Dim As Integer mouseEvent, mouseDrag
Dim As mouse_type mouse
Dim As person_type person(0 To PERSONS-1)

Randomize 1234
'initialise persons
For 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_INIT
Next

'time step such that a max speed, position change = 10 % of infection distance
Dim As Single timeStep = (INFEC_DIST * 0.1) / V_MAX 'sec
Dim As Single simulDays, simulTime = 0
Dim As Integer numInit, numDead, numRecov
Dim As boolean sectioning = TRUE

sg.offset.x = -15 'move view for number
Dim As Double tNow(0 To 9), tAcc(0 To 9)
Dim As Integer trigger = 1

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 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 While
Wend
Locate 1,1: Print "End of simulation!"
While InKey = "" : Sleep 1 : Wend
End

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.
badidea
Posts: 2000
Joined: May 24, 2007 22:10
Location: The Netherlands

Re: Corona virus simulator

Postby badidea » Apr 01, 2020 22:23

Simulation of virus containment failure:
Image
Same conditions without borders:
Image

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 String
End Type

Constructor int2d
End Constructor

Constructor int2d(x As Integer, y As Integer)
   This.x = x : This.y = y
End Constructor

' "x, y"
Operator int2d.cast () As String
  Return Str(x) & "," & Str(y)
End Operator

Operator = (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 true
End Operator

Operator <> (a As int2d, b As int2d) As boolean
   If a.x = b.x And a.y = b.y Then Return false
   Return true
End 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 - b
Operator - (a As int2d, b As int2d) As int2d
   Return Type(a.x - b.x, a.y - b.y)
End Operator

' -a
Operator - (a As int2d) As int2d
   Return Type(-a.x, -a.y)
End 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 * mul
Operator * (a As int2d, mul As Integer) As int2d
   Return Type(a.x * mul, a.y * mul)
End 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 \ div
Operator \ (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 sgl2d
End Type

Constructor sgl2d
End Constructor

Constructor sgl2d(x As Single, y As Single)
   This.x = x : This.y = y
End Constructor

Function sgl2d.cross(b As sgl2d) As Single
   Return This.x * b.y - This.y * b.x
End Function

Function sgl2d.lengthSqrd() As Single
   Return (This.x * This.x) + (This.y * This.y)
End Function

Function 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 Function

Function 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 Function

Function 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 / lenth
Operator 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 true
End 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 true
End 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 - b
Operator - (a As sgl2d, b As sgl2d) As sgl2d
   Return Type(a.x - b.x, a.y - b.y)
End Operator

' -a
Operator - (a As sgl2d) As sgl2d
   Return Type(-a.x, -a.y)
End 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 * mul
Operator * (a As sgl2d, mul As Single) As sgl2d
   Return Type(a.x * mul, a.y * mul)
End Operator

' a / div
Operator / (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-sorted
Type 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 methods
End Type

Constructor grow_list(startSize As Integer)
   If startSize > 0 Then
      ReDim index(startSize - 1)
   End If
End Constructor

Constructor grow_list()
   This.constructor(0)
End Constructor

Destructor grow_list()
   Erase(index)
End Destructor

Sub 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 += 1
End Sub

Sub grow_list.empty()
   numItems = 0
End Sub

Function grow_list.numAlloc() As Integer
   Return UBound(index) + 1
End Function

Function grow_list.numInUse() As Integer
   Return numItems
End Function

'for debugging
Sub grow_list.show()
   Print "--- " & numInUse() & " / " & numAlloc() & " ---"
   For i As Integer = 0 To numItems - 1
      Print i, index(i)
   Next
End 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 9

Type 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 button
End Type

Function 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 change
End Function

'===============================================================================

'Note: y+ = up, x+ = right, (0,0) = center
Type 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 Type

Sub 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 font
End Sub

Sub scaled_graphics_type.setScaling(scale As Single, offset As sgl2d)
   This.scale = scale
   This.offset = offset
End Sub

Sub scaled_graphics_type.clearScreen(c As ULong)
   Line(0, 0)-(w - 1, h - 1), c, bf
End Sub

Function 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 Function

Sub scaled_graphics_type.drawPixel(p As sgl2d, c As ULong)
   Dim As int2d posScrn = pos2screen(p)
   PSet(posScrn.x, posScrn.y), c
End Sub

Sub 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, c
End Sub

Sub 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,,,,f
End Sub

Sub 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, , , aspect
End Sub

Sub 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), c
End Sub

Sub 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, b
End Sub

Sub 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, bf
End Sub

'===============================================================================

Const As Single INFEC_DIST = 0.05 'm
Const As Single INFEC_DIST_SQRD = INFEC_DIST ^ 2
Const As Single SICK_DAYS_MIN = 7, SICK_DAYS_MAX = 21
Const As Single MORTALITY = 3 / 100
'const as single SICK_SPEED_FACTOR = 0.25 'disable due to some bug

Const SEC_PER_DAY = 24 * 60 * 60
Const As Single DAY_PER_SEC = 1 / SEC_PER_DAY

Const PERSONS = 5000
Const As Single V_MAX = 0.1 / 3600 '1 m/hr
Const As Single MAP_X_MAX = +50, MAP_X_MIN = -MAP_X_MAX
Const As Single MAP_Y_MAX = +35, MAP_Y_MIN = -MAP_Y_MAX
Const As Single MAP_W = MAP_X_MAX - MAP_X_MIN
Const As Single MAP_H = MAP_Y_MAX - MAP_Y_MIN

Const DRAW_INTERVAL = 10 '1 to 10
Const FOLLOW_PATIENT_0 = FALSE

Const NUM_SECT_X = 16 * 3, NUM_SECT_Y = 12 * 3
Dim As grow_list sector(NUM_SECT_X - 1, NUM_SECT_Y -1)

'screen stuff
Const As Single PPM = 6.0 'pixels per meter (set zoom level)
Const SCRN_W = 1050, SCRN_H = 750
Dim Shared As scaled_graphics_type sg
sg.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 states
End Enum

Dim 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 'health
End Type

Function 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 If
End Function

Function rndRangeSgl(min As Single, max As Single) As Single
   Return Rnd() * (max - min) + min
End Function

Sub limitInt(ByRef value As Integer, min As Integer, max As Integer)
   If value < min Then value = min
   If value > max Then value = max
End Sub

Dim As Integer mouseEvent, mouseDrag
Dim As mouse_type mouse
Dim As person_type person(0 To PERSONS-1)

Randomize 1234
'initialise persons
For 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_INIT
Next

'time step such that a max speed, position change = 20 % of infection distance
Dim As Single timeStep = (INFEC_DIST * 0.2) / V_MAX 'sec
Dim As Single simulDays, simulTime = 0
Dim As Integer numInit, numDead, numRecov
Dim As boolean sectioning = TRUE

'move view
sg.offset.x = -35
sg.offset.y = -25
Dim As Double tNow(0 To 9), tAcc(0 To 9)
Dim As Integer trigger = 1

sg.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 While
Wend
Draw String (10, SCRN_W - 26), "End of simulation", BLACK
While InKey = "" : Sleep 1 : Wend
End
badidea
Posts: 2000
Joined: May 24, 2007 22:10
Location: The Netherlands

Re: Corona virus simulator

Postby badidea » Apr 03, 2020 10:31

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 = 700

screenres SCRN_W, SCRN_H, 32
width SCRN_W \ 18, SCRN_H \ 16

const ST_INIT = 0, ST_INFEC = 1, ST_RECOV = 2, ST_DEAD = 3, ST_INVALID = 4
const SICK_DAYS = 14
const as single MORTALITY = 3 / 100

const 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 * 24

type person_type
   as integer x, y
   as integer state
   as integer sickEndTime
end type

const PERSONS = 2500
dim as person_type person(any)
redim person(PERSONS - 1)

'map array dynamic, else too large for stack
const MAP_X = 700, MAP_Y = MAP_X
dim 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
   next
end sub

clearMap(map())

const RND_DEF = 0, RND_CRT = 1, RND_FST = 2, RND_MTW = 3, RND_QBC = 4, RND_SYS = 5
randomize 222, RND_FST

for i as integer = 0 to PERSONS - 1
   person(i).x = int(rnd * (MAP_X - 2)) + 1
   person(i).y = int(rnd * (MAP_Y - 2)) + 1
next
person(0).state = ST_INFEC
person(0).sickEndTime = SICK_DAYS * SEC_DAY

dim as integer quit = 0
dim as integer simTime = 0, simTimeStep = SEC_MIN * 10
while 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 = 1
wend
sleep
ShawnLG
Posts: 127
Joined: Dec 25, 2008 20:21

Re: Corona virus simulator

Postby ShawnLG » Apr 14, 2020 16:42

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.
badidea
Posts: 2000
Joined: May 24, 2007 22:10
Location: The Netherlands

Re: Corona virus simulator

Postby badidea » Apr 14, 2020 19:18

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: 127
Joined: Dec 25, 2008 20:21

Re: Corona virus simulator

Postby ShawnLG » Apr 15, 2020 0:16

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


Yes they are.

Return to “General”

Who is online

Users browsing this forum: Google [Bot] and 3 guests