Realtime Covid-19 spread modelling

This is an upgrade to the previous project. Instead of taking hypothetical values for the SIR model parameters, this program retrieves the latest (daily) time series data for the Coronavirus pandemic from https://pomber.github.io/covid19/timeseries.json and simulates how the infection might play out in the future. The time series data is used to get an estimate of the SIR parameters β (the "contact rate" constant) and α (recovery constant). These vary enough through time to cause wild variations in the model outcome. So, I've included the option to use the latest, the minimum, the maximum and the mean values, out of which the mean seems to give the most sensible looking graph. The way the parameters are estimated is basically using two (shorter) of the three governing equations of the model:
delS/delT = -β*S*I/N
delI/delT = β*S*I/N - α*I
delR/delT = α*I
from where,
β = delS/(I*S/N)
α = delR/I
where delT = 1 day, since that is how frequent the available data is.
Now, the Susceptible population S, Infective population I and the Removed population R aren't available as such. The time series data source mentioned above gives the daily recorded data for each country since Feb. 22, 2020. This is converted into a global daily record time series data by summing over all the countries. The final form looks something like this, and this is for a particular date:

1. date: "2020-3-26"

1. confirmed: 529591

1. deaths: 23970

1. recovered: 122149

The 'confirmed' value corresponds to all the active infections (which meets the definition of the infective population in our model) + the 'deaths' + the 'recovered' values.

So, for the infective population, I = confirmed -  deaths - recovered

and for the susceptible population, that's the global population minus the people who are still in active infection and those who've recovered or died combined i.e. S = N - confirmed

I tried doing a simple power-law function fitting for the alpha and beta parameters using the Least Square Technique (LST) here but ultimately decided against using it in the model as I thought their time-evolution should be more logistic-like rather than monotonous. For that spreadsheet, I used data attributed to Mathemaniac after watching his video

Perhaps I could use the latest 5-day or 14-day averages for better estimating the model parameters. The model would do even better if I could figure out β(t) and α(t) for the entire duration of the epidemic since the parameters seem to be constantly fluctuating IRL.

Anyway, here is the code.

Edit: I've removed the (flawed) max/min parameter estimation and added last 5 and last 14 data based parameter estimation. Doesn't seem to do much good.
Edit: I've only now realized that web apps that don't require/have a back end i.e. web pages that don't require server-side coding can be hosted on GitHub itself! I have thus pushed this project to GitHub as a repository. The corresponding GitHub Pages site is here. I will be pushing all future JavaScript web apps directly to GitHub as repositories of their own; no more GitHub Gist for such programs. All update to the apps will be applied to their corresponding GitHub repositories, not their Gists. I will also be shortly porting all my past JavaScript web apps in this blog that I'd uploaded to GitHub Gist to GitHub as repositories.

Epidemic Spread Modelling with the SIR model

The recent national lockdown due to the Covid-19 global pandemic has afforded me all the free time in the world. So, I took this opportunity to learn more about how experts study epidemics to make inferences and devise policies. I inquired Google and YouTube on the topic and quickly discovered a wealth of relevant resources. Books, papers, articles and videos were aplenty. Herbert W. Hethcote (Three Basic Epidemiological Models), Aresh Dadlani et. al. (Deterministic Models in Epidemiology: from Modeling to Implementation), Fred Brauer and Carlos Castillo-Chávez (Epidemic Models) and Yiran Jing were the authors whose work I found really helpful in understanding how epidemic modelling works. YouTube videos from channels: Trefor Bazett, Tom Rocks Maths and VetenskapensHus were veritable gems on the subject owing to their platform-incentivized conciseness. Anyway, back to SIR - this model is arguably the simplest of all while still being able to yield useful insights into the workings of an epidemic.
There are two variants of the model - with and without vital dynamics (birth, death, migration). Here, I deal with the no vital dynamics variant which basically implies there's no significant change in the total population of the community under study, at least until the epidemic lasts. The idea is, the population is divided into three categories - Susceptible (normal, not yet infected individuals), Infective (infected individuals, who can infect others) and Recovered/Removed (individuals who have recovered from the disease and will not again go back to being infected, say because they've gained immunity). Populations of each compartment are constantly changing, moving from Susceptible(S) to Infective(I) to Recovered(R). The flux or rate of change of the population of each compartment (dS/dt, dI/dt, dR/dt) is assumed to be functions of one or multiple of such categorical populations. For example, the rate at which S changes into I is assumed to be proportional to the product of the Infective population and the Susceptible population. Such an ODE is written down for each of the compartments - giving a total of 3 ODEs that have to be solved simultaneously while satisfying the constant population assumption(S+I+R=N=constant). The system of ODEs happens to be impossible to solve explicitly owing to non-linearities in two of the three equations. So, only numeric solutions for I, S and R as functions of time are possible, via discretization as per the ole'faithful Finite Difference Method, which is what I've done here. The details of the model can be found on any of the resources listed above.
Here's how the graphs look:

Here's the code.
EDIT:
This is after I've realized what GitHub Pages is.
Live demo here
Corresponding repo here
Future updates, if any, will be made to the repo, not the Gist

Riverbed variation model

I haven't been able to manage enough time to code lately. It has got to do with me starting a job - out of a necessity of - I don't know what. But more on that on my other blog (that I hope I do write). Anyway, I thought I should put it here for my own sense of continuity rather than for the (dys)functionality slash (un)greatness of this code.

This is a problem that has to do with predicting the evolution of riverbed elevations in sediment carrying rivers. Yeah, also part of the (second semester of) M.Sc. Water Resources Engineering program's Water Induced Hazards::Riverbed Variation(RBV) module.
The problem of scraping together enough time to invest in writing a functioning program aside, this physical process inherently seems to me so very intractable to put into a working model, even now, when I believe I have a somewhat acceptable result it spit out a few minutes ago which was when I decided to do this write-up. What the program does is solve three PDEs describing the phenomenon using the much acclaimed MacCormack explicit finite differencing scheme.
When the equations were first presented in class, I don't believe there was much explanation regarding what each term meant  in a physical sense, and more importantly, how one interacted with  others. The presentation went something like this: write down the PDEs; discretize them to obtain equations for the next time level's riverbed elevation, water depth and specific discharge; write down the CFL condition; try and work out a numerical problem. So many questions were left unanswered: What should the boundary conditions be? Are they important? How sensitive is the output to boundary conditions? What does the sediment transport equation relating specific sediment discharge(qs) with specific discharge(q) mean? What does it mean to solve for specific discharge? Does specific discharge itself change as a result of riverbed variation? What output am I expecting? How should the computation procedure move in the finite difference grid? And as I put pen to paper, so to speak, there were so many possibilities to try - each requiring a great deal of effort - while most likely being wrong. One would easily be forgiven for giving up on even trying. Anyway, cut long story short, I tried to get it to work - a lot - and failed - a lot. Until I stumbled upon this article that referenced

Alam M (1998) Application of MacCormack scheme to the study of aggradation degradation in alluvial channels M.Sc. Engineering thesis work, Department of Water Resources Engineering, BUET, Dhaka.

 available here. The "Full thesis" was freely available for download at the second link. The document proved extremely helpful and answered most of my questions. It also detailed a set of controlled lab experiments on RBV and gave me a sense of what I should be expecting from my program.
So, here is what the program does:
1. Takes in all the relevant input parameters.
2. Populates the first element of an array of "profile"objects. The first object in this array has the initial condition values for h, q and z along with x and t. x is an array of downstream nodal points. h, q and z are arrays of water depths, specific discharges and bed elevations at their respective x positions. t is the time level of the profile.
3. Runs a loop from 0 to the stipulated maximum time of simulation. For each iteration of this loop, there's another loop for computing the next time level's h, q and z values using the MacCormack explicit formulae at all downstream nodes.
4. The sediment discharge at the first (zero) node is held at a constant value that is a user-specified multiple of the sediment carrying capacity(?) of the initial flow [i.e. qs=a(q/h)^b] for all time levels. So is the water discharge.
5. Displays the final profiles in a graph as an animation.
While in algorithm, the process seems pretty straightforward, there are subtleties involved to say the least in the actual implementation in code. The MacCormack scheme needs values of the unknown parameters one spatial node in excess to the up and downstream while time-marching from any node. This necessitates calculations of parameters in as many nodes downstream of the actual downstream boundary as the number of time levels desired. For nodes requiring info from upstream of the upstream boundary, the respective initial state values are assigned i.e. no fictitious nodes black magic is performed.
Also, keeping track of the profiles at different time levels and using the right ones in the scheme without screwing up in the slightest is easier said than done.
There is also the coordination between the main and the WebWorker thread and figuring out a good dynamic between the two.
Not all input parameters output a good result. In fact, most inputs lead to an explosion of the predicted parameters at some point in time. I've tried decreasing delT thinking it may be a manifestation of the scheme's instability because the CFL condition may have been violated. But even when there's a more than comfortable buffer between the used delT and the threshold between that predicted by the Courant condition, I've seen irrational predictions owing to negative q's, tiny q's and what not. While the aforementioned paper recommends using a new delT after each time march depending on the new q and h values using the CFL criterion, I'm skeptical it would offer any significant improvements.
Anyway, I'm done with this.
EDIT:
This is after I've realized what GitHub Pages is.
Live demo here
Corresponding repo here
Future updates, if any, will be made to the repo, not the Gist

I don't know why the first time march reduces the bed level to below the initial value here.

Convex hull

In the previously posted flownet generation program, there is a seemingly insurmountable problem. After equipotential points for a given head have been extracted from between the solved finite difference grid points (here, using linear interpolation) comes the hard part: properly joining them. I have outlined how I approached the problem in the first post about the program. There, I did point out that in equipotential points for a given head, discovered to the extreme left and right of the grid, a jumpy line springs into existence because of the limited space the grid actually represents which obviously clips the isobars. However, at the time, I hadn't realized that there was an even more sinister problem with the program and the 'connect the dots' logic - one which doesn't, at this point, seem solvable because of the weird ways in which these equipotential lines naturally seem to be able to bend. Take for example the following output with
contourInterval=5
numDys=13
numDxs=16
DxsTillDam=2
DamWidthInDxs=4
reservoirHead=100
downstreamBoundaryHead=50

The problem is visible for the 15m isobar. Starting from underneath the dam, the point near the 16.29m grid node should have been the one to be joined first rather than the point between the 11.06m and 19.26m grid nodes. It is easy enough for us to say via a simple visual inspection. But the 'order in descending Y coordinate and then by shortest consecutive distance' logic falls flat in this case since clearly, the program is doing exactly that. The point actually being joined is clearly closer to the equipotential point on dam than the 'right' point which happens to be a little too left.
I thought of trying all orderings of the points to minimize total length but quickly concluded that it would be computationally unviable - the number of equipotential points for a head regularly crossed 20 when I was testing the program - that's 20!, or about 10^18, number of combinations to try for drawing each isobar. What's more, I wasn't sure it would even produce the desired result. So, I scrapped the idea.

Then, I researched ways to do this online. Something called 'Convex hull' came up. Basically, it is a way of drawing a polygon around a bunch of points while passing through the right subset of the said points so that no point is left outside the polygon. The trick is to test each possible line segment (defined by two unique points) for convexity i.e. whether all other points lie towards the same side (of the line segment); if they do, the line segment being tested is part of the polygon - the 'convex hull'. The idea behind this test - whether a point lies toward one side of a line or the other - comes from checking the point's y-coordinate against that the line's, at the point's x-coordinate. If all the points lie above (the idea of a point being 'above' a line can be extended from the obvious case where the line is horizontal, right until it isn't completely vertical) the line, their y-coordinates are going to be greater than the y-coordinates of the line segment at their(the points') x-coordinates. Anyway, this thread's first answer is what I implemented - but not into the Fortran program itself. I was skeptical you see. I couldn't quite picture how this thing worked or what it would look like in action. So, I first made a small JavaScript program using p5.js so that I could have something concrete to look at and potentially tweak.

Here is the code for it.
Then, it became clear that no amount of tweaking the convex hull algorithm could manage to connect the equipotential points 'properly'. There would be too many exceptions - non-convex points that would need to be connected, but would not be. Leaving such points out would produce smooth and nice equipotential lines but that would also be misleading and far from accurate.
So, there's where I basically left the flownet program, without progress.
But at least I got to learn something new in the process.

Flownet update

This is a revised version of the last flownet generation program.

Changes:

• Head at dam heel is equal to the impounded water depth as should be(was assumed 0 before)
• Dam was fixed to be the size of a single grid box. Now, width of dam can be specified
• Dam's heel and toe are Dirichlet boundaries i.e. known heads but any node between are treated as unknowns
• Better code documentation during coefficient matrix construction.

The code file is 'Source1V2.f90'

My sketches for this revision.

An associated Excel spreadsheet for control.

Some more screenshots:

Flownet of the original problem with numDxs=4,numDys=3,DxsTillDam=2,DamWidthInDxs=2,reservoirHead=100,downstreamBoundaryHead=50 and contourInterval=10

Flownet of a problem with numDxs=5,numDys=5,DxsTillDam=2,DamWidthInDxs=3,reservoirHead=100,downstreamBoundaryHead=50 and contourInterval=10

Flownet of a problem with numDxs=15,numDys=15,DxsTillDam=5,DamWidthInDxs=3,reservoirHead=100,downstreamBoundaryHead=50 and contourInterval=5

The weird aggregation of isobars at the right(downstream) boundary between 0 and 50m is because of the type of boundary condition given. A constant head (50m) Dirichlet boundary in the vertical direction is unrealistic. An impermeable vertical plane (Neumann boundary) would have allowed a smooth and probably more realistic variation of head downstream.

Flownet Generation in Fortran

This was a 'challenge' assignment (at least that's how I perceived it) for us by our Fortran Simulation Lab professor. We were given a dam impounding 100m deep water. We had to solve the Laplace's equation in 2D using the 'implicit scheme' - basically solving a system of linear equations consisting of Laplace's finite difference equation for each unknown node. For the finite differencing domain, 4 grids along the horizontal and 3 grids along the vertical (each grid being a square of unspecified dimensions) were provided. We were also given Dirichlet boundary conditions on the nodes at the downstream end. When I first set out on manually solving the problem as a control for when actually developing code for it, I thought that the heel of the dam should have a head equal to the impounded water's depth, as all points upstream on the ground surface should - and still do think so. But that is how the problem was posed and so I ran with it.
The problem
I wanted to make the code 'general' in that it could work for other scenarios as well - not for just this specific combination of parameters. So, I decided to allow for 6 inputs: contourInterval, numDys, numDxs, DxsTillDam, reservoirHead, downstreamBoundaryHead. reservoirHead and contourInterval should be pretty obvious. numDxs and numDys are just that - total grids/boxes to work with along horizontal and vertical axes. DxsTillDam is the number of grid boxes to the left of dam's heel. downstreamBoundaryHead is the head at the downstream boundary. Based on these parameters, I had to think hard on what the coefficient matrix for the system of linear equations would look like. Since the matrix has a lot of structure in it, it was not that big of a problem. I just looked at the matrix I had prepared while working out the problem by pen and paper and coded in the logic to populate a 2D array with the exact same rule that the elements of the original problem's matrix followed. Another 2D array that functioned as a single column matrix i.e. vector representing the Dirichlet boundary conditions was created side by side while populating the coefficient matrix.
Matrix for reference and Manually solved grid for reference
Mazumder and Wen Shen were indispensable for understanding the matrix formulation.
The next step was to find the points on the grid lines that had the same head for each head. After some thinking, I settled on the idea of sweeping from the left boundary towards the right upto the right boundary(-1 to be precise), one column at a time, looking for a specific head to the immediate right and to the immediate above standing at each node moving up the column to top(-1 to be precise). Since the head being searched for rarely exist on the grid nodes themselves, I interpolated between the current node and the east node, and the current node and the north node to get the exact coordinates of the head - provided it existed between the nodes. I stored the interpolated coordinates off to two matrices. One matrix for the X coordinate, another for Y. In each matrix, each column corresponded to a head, viz. a column for 10m, another for 30m, etc.
Now that I had the equipotential coordinates/points for each head in question, all that was left was to join them. I first tried joining them in the order of their discovery i.e. according to their positions in the coordinate matrices. That didn't work very well as anyone might've guessed. Order matters when joining a bunch of equipotential points. You don't want jigjagging lines going back and forth for what should be a smooth isobaric curve. I then tried joining the points from the top down - basically sorting them in descending order of their Y coordinate. But that didn't work very well either. Since the isobars tend to be U-shaped, a lot of times, this method led to left-right-left jigjags. I then thought of joining the points from left to right. This worked for majority of the isobars. But it almost always messed up one isobar - right in the middle, when the curves shift from left facing to right facing. The isobars to the left and to the right are 'proper' functions i.e. unique Y for any given X. However, at the transition zone, between the two kinds - which happens to be right at the center of the base of dam(not here though, for reason described towards the top of this article) or more precisely at the ground level right between nodes of reservoirHead and 0 - the isobaric curve is somewhat vertical and hence exhibits 'ill' behavior as a function i.e. two or more Y's at some X. It is precisely due to the existence of such isobars that the left to right joining method fails towards the center of the grid. I finally ended up adopting a hybrid technique: joining from top to bottom but with shortest distance to the immediately next bottom point consideration. This involved sorting the equipotential points for a head in descending order of Y coordinate and then sorting by shortest distance between any two consecutive points. This gave better results but for some unknown reason, which may or may not be related to this technique, bottom most isobars starting making jump between left and right. I have not researched too much into this issue but it may just be how flownets are: loop-like and when a loop gets clipped due to our domain boundary, there's bound to be equipotential points to the right and to the left for a head with the bottom points clipped. In such a case, this 'jumpy' behavior is exactly what would be expected. Anyway, I just wrote a couple of anti-jump lines to prevent joining points farther apart than a specific X and Y distance.
My scribbles: 1 2 3
The viewport width and starting position(top left coordinates) are by default set to 600, 200, 8. Sometimes a larger viewport width may be desirable. The user may freely change these without fear of messing up the grid, irrespective of problem-centric parameters listed in the second paragraph.
Code

Flownet with reservoirHead=100, numDxs=4, numDys=3, contourInterval=10, downstreamBoundaryHead=50 and DxsTillDam=2

Flownet with reservoirHead=100, numDxs=20, numDys=20, contourInterval=10, downstreamBoundaryHead=50 and DxsTillDam=4

Dam-break simulation

This is a JavaScript based dam-break problem simulator. Something I am studying as part of Advanced Hydraulics curriculum. It's basically just one equation relating x, t and depth y.
Wanted to make something tangible. So, I did this. Just for fun. Makes use of Chart.js and Web Worker. Web Worker script requires to be kept in a separate .js file. Chrome doesn't allow html files to access local filesystem so, open using Firefox or use the flag --allow-file-access-from-files while starting Chrome and then open the html file. Code

EDIT:
This is after I've realized what GitHub Pages is.
Live demo here
Corresponding repo here
Future updates, if any, will be made to the repo, not the Gist