ds_utils documentation

ds_utils

index /home/david/Dropbox/PythonCode/ds_utils.py

ds_utils.py
 
A library of miscellaneous utils that I find I use often. 
 
Version 1.0
 
dsinfield6@gmail.com/ Telegram or signal:omotforest 
 
Released under GPL
 
Use at your own risk. Subject to update or change at any time.

Modules
		getpass os platform

Classes
		ANSI Searches config   class ANSI     GetANSI(text,foreback,colour,style) text - text to print forecolour - [black|red|green|yellow|blue|magenta|cyan|white] backcolour - [black|red|green|yellow|blue|magenta|cyan|white] style - [normal|bold|italic|underline|inverse|strikethough]   Returns a string with ansi precursors to set foreground and background colours and style. Set to normal after printing.   Most terminals support 8 colours. Some will not understand bold.   Out of range parameters throw a dictionary error   Colours supported:     black     red     green     yellow     blue     magenta     cyan     white Styles supported:     normal     bold     italic     underline     inverse     strikethrough   Example: set foreground to red     print GetANSI('this is red text','red','none','none') Example: set foreground bold red, background cyan     print GetANSI('set properties','red','cyan','bold')     Methods defined here: GetANSI(self, text, forecolour='', backcolour='', style='') Constructs an ansi string to create multicoloured terminal display Parameters     text - the string to be output     forecolour - the colour to set the foreground. Optional. Default=normal     backcolour - the colour to set the background. Optional. Default=normal     style - one of the styles in dictStyle Returns      A string decorated with ansi codes __init__(self) Data and other attributes defined here: dictColour = {'black': 0, 'blue': 4, 'cyan': 6, 'green': 2, 'magenta': 5, 'red': 1, 'white': 7, 'yellow': 3} dictForeBack = {'back': 4, 'fore': 3} dictStyle = {'bold': 1, 'inverse': 7, 'italic': 3, 'normal': 0, 'strikethrough': 9, 'underline': 4} post = 'm' pre = '\x1b[' tonormal = '\x1b[0m'   class Searches     A class to contain search functions.   Initialisation __init__(head, tail,source,start=0,ignorcase=false,inclusivefind=false)     head - string to search from     tail - string to search to     start - integer position to start search. Optional. Default=0     ignorecase - set to true to make case insensitive. Optional. Default=false     inclusivefind - return the string from the start of head the end of tail. Optional.     overlapsearch - if True sets newstart to start+1          default is to return the string from the end of head to the beginning of tail  NB init calls get head to tail once so results of the first search are available after initialisation Reading results     newstart - set to the end of the head     result - set tp the string found or "" empty if not found     resultstatus - text description of result of search     resultsuccess - set to True if head and tail were found.     EOF - set to true if the search extends past end of file     listAll[] - set only if getAll has been called     Methods defined here: __init__(self, head, tail, source, start=0, ignorecase=False, inclusivefind=False, overlapsearch=False) getAll(self, head, tail) getHeadtoTail(self, head, tail, start=False) Data and other attributes defined here: EOF = False head = '' ignorecase = False inclusivefind = False listAll = [] newstart = 0 overlapsearch = False result = '' resultstatus = '' resultsuccess = False source = '' start = 0 tail = ''   class config     Config File Helper Class   Config files are in the format:     #start comments     ParameterName1 value1     ParameterName2 value2     .     .     .     #endcomments Assumes values are one word when reading. Anything after a second space is ignored.      Constructor takes 1 arguments     filename - the name of the file only the path will be derived depejding on the operating system                  Linux - will be a hidden file the user's home directory                  Windows will be a file in the same directory as the application                  Other will cause an exception     Methods defined here: __init__(self, filename) getConfigPath(self) Accessor function for file path property Call with object reference as parameter readParameter(self, parametername) Read the parameter from the file     Returns False if not found writeFile(self, dictParams, comments='', endcomments='') dictParams - a dictionary of parameters e.g {name1 : value1, name2 : value2} comments - a file header. Optional, default is empty string endcomments - a footnote. Optional, default is an empty string.            For neatness when doing a cat /n will be appended Data and other attributes defined here: config_file_name = '' config_file_path = '' platform_supported = False

Functions
		InputUntil(prompt, acceptable, case_sensitive=True, tries=0) Usage InputUntil(prompt,acceptable_input,[case sensitive=True],[tries=0]) Parameters:     prompt - Prompt to be displayed (the acceptable terms will be listed)     acceptable_input - String containing acceptable responses     [case_sensitive] if false ignores case.; Default is True     [tries] - Attempts before raisin an error. Keeps gettin input until one of the terms in acceptable is found. Works best for single characters See InputUntilList(prompt, acceptable, tries) for a word version   You should still check the returned input because in YNC, YN will match. InputUntilList(prompt, acceptable, case_sensitive=True, tries=0) Input until one of the items on the list is typed or the number of tries is reached. Returns the actual answer typed (irrespective of the case sensitive flag. Raises an error if tries exceeded catch it in your code to handle it.   For single letter (eg YN), the InputUntil in this module might be what you are looking for   Usage InputUntil(prompt,acceptable_input, [case_sensitive=True], [tries=0]) Parameters:     prompt - Prompt to be displayed (the acceptable terms will be listed)     acceptable_input - list containing acceptable responses       case_sensitve - David DAVID, daVid all match david if this is false     tries - Number of attempts until an error is raised. Set to <1 for infnite getParameterValue(parametername, args, unitary=False) Gets the command line parameter wih the name in parametername. Assumes the parameters are in -name value pairs. Longname parameters in the form --long_name value are supported but you must supply the -- Parameters with no value are not supported. (Coming soon: Optional argument to specify standalone which returns True/False) Usage: getParameterValue(parametername,args) Parameters:     parametername - the name of the parameter if no - then one is prepended     args - usually sys.argv. If you want to use another string you need a placeholder where the program name would be            e.g ["params",  "-p", "param1","-p2"] param2 works as expected     unitary - will not look for, a value just the presence or absence of a parameter makeDirectory(dirPath) Usage makeDirectory(directory_path)   Makes the directory in dirPath if it does not exist   Throws exception if it fails.

Data
		ARG_NOT_FOUND = False version = '1.0'

Robincam setup

Introduction

I wanted a remote webcam for watching the robin that lives around my apple tree. Trailcams are a bit pricey and I had an old android phone with a half decent  camera laying around doing nothing.

The Phone...

...s a very battered Blu Vivi 6. The screen is broken, the USB C port is loose and the battery lasts about 30 minutes because it has been frozen a few times..

The USB C port is fixed by gluing in a magnetic connector which also seals the port somewhat against the weather.

 I've installed IP Webcam on it from the play store. You can find it here https://play.google.com/store/apps/details?id=com.pas.webcam

I've gone for the Pro version but everything here works with the free one. The only limitations are Tasker integration and the ability to customise the UI. We are not going to use either of those here. I upgraded because I think the developer deserved a coffee (or Russian tea) on me.

Installing it is straightforward and I've also put it on a couple of old HTC mobiles that are being security cams.

To do everything in this article you are going to need a bit of memory. The Blu has 64GB and I'm using about half of it. That's about 1200 jpgs and several hours of mp4s.

I'll talk about how to get the images off to a pc and get the memory back later.

You'll also need a wifi signal where your camera is sited. You can use broadband and access your through the lvideon cloud see here https://www.ivideon.com/

I haven;t got a sim in the Blu so that's not an option.

I use this https://www.amazon.co.uk/NETGEAR-Wi-Fi-Range-Extender-EX2700/dp/B00NIUHAG6/ref=sr_1_5 network extender to get to the shed so I have coverage over all of my garden. Cheaper wifi extenders are available but this one works for me.

IP Webcam Setup

When you start IP Webcam it shows the setup screen. For plain vanilla you can just jump to the last item and it will work.

Before you set up anything you have to give IP Webcam some permissions to photos and sound.


For a start set up the Video Preferences I use the highest resolution possible and lock the camera in Landscape mode.

The other important settings are the motion and sound detection. I have motion detection on and sound detection off. I was getting a lot of false positives from the wind blowing across the camera. Getting the optimum sensitivity will take trial and error and you can play around with this in the web interface.

Starting the Webcam Server

IP Webcam comes with its own webserver so you can connect to it with a browser. When you start the server it tells you the web address that you need to type in to your browser.

The first thing you should do is configure your router to make this a address permanent. I have a BT Home Hub and the important bit of the settings screen looks like this:

 

If you find this doesn't work and the assigned ip address changes it could be that your phone is registered on the network more than once. Have a look for duplicate names or MAC addresses in your connected devices and remove the ones you are not using.

Power Management

Even with the all the power management options available the application is power hungry, especially if you opt for aggressive focusing which keeps the camera focused constantly so you can take a still image instantly. Note that if you want to record stills automatically you are going to need the IP Webcam add in. I'll blog about that later. If I'm out of reach of mains power, I use an external battery pack. A big RAV 280000mah battery pack like this https://www.amazon.co.uk/RAVPower-Portable-Charger-26800mAh-External-Black/dp/B012V9H3WA/ in a margarine tub with a hole for the power cable.

Weatherproofing

The Blu comes with a silicone case which is sealed except for the camera, headset and power sockets sealed. I've covered the headphone socket with tape. I wouldn't trust it in a raging storm but it has managed the odd shower without coming to any harm. I'm considering sealing it in a food box with only a hole for the camera lens which I can seal with elefant snot (that really is what Blue Tack is called in Danish)

All that is left is to site your camera in an appropriate place and wait for the motion sensors to pick stuff up. The video will record in a local file which you can watch on your phone. Until you get the motion sensitivity fine tuned, expect a lot of false positives. Especially if there is foliage in your camera field. You'll get a video every time the wind blows.

Uploading to cloud storage 

I use google photos to backup and manage the pictures. You can do this through the IP Webcam web interface but the google photos interface is a bit nicer and when they are on the cloud you can safely delete them from the phone.

First make sure you have recorded at least one video.

Setting up photos is pretty straightforward.go to Settings -> Backup and Sync --> Backup device folders and you should see a folder called 'modet.' Add that to the folders backed up and that's it. There's a bit of a lag between recording and the videos appearing in in Photos so be patient.

Expect to be clearing out the false positives frequently. The free allowance for original size photos is 50GB. You'll go through that quite quickly. The great thing is that deleting from Photos also deletes on your device so it keeps it cleaned up reasonably cleaned up.

There is a Dropbox plug in for IP Webcam but I seem to have broken it. I've emailed the developer to see what I've done wrong.

Coming Soon

Next article "Setting up IP Webcam to take still images" coming soon.

Underwhelmed

Data and Methodology

If he lockdown was intended to protect the NHS then it might have worked a bit too well. We were told we had to flatten the curve and we did. But by how much? We've seen the stories of empty hospitals and because elective and non urgent surgery cancelled.

Here's a chart generated from https://www.england.nhs.uk/statistics/statistical-work-areas/covid-19-daily-deaths/. I haven't done anything with it except cut and paste the totals of daily deaths by age into Libreoffice calc, Transform it from columns to rows and then copy it to Google sheets, which I've used for all the other manipulation.

And here's a chart I made from publically available data.

The blue bars are daily hospitalisations derived from 

• Daily deaths as reported on Worldometer (this snapshot was taken when only hospital deaths were reported)
• I've applied a 4 day smoothing to account to try to smooth out the weekend reporting systemic errors.
• I've multiplied the deaths by 3 as only one third of long term hospitalised patients survive. https://news.sky.com/story/coronavirus-a-third-of-all-covid-19-hospital-patients-die-with-dementia-and-obesity-carrying-the-highest-risk-11980894
• I've assumed that these are the patents that are occupying the critical care beds.

The red line is:

• On April 12th there were 2295 available NHS critical care beds. https://www.reuters.com/article/health-coronavirus-britain-nhs-idUSL9N26V04T
• I've added 2295 to the April 12th calculated hospitalisations
• this excludes the 4000 beds at the nightingale hospital which were not fully available at the time.

That's all data. I'm now going to step into the realm of models and speculation.

The big surprise came when I came to construct a tall thin bell curve with the same area which just failed to breach the capacity line.

The tall thin curve came later with a slower incremental part.

I know this can't be right. But as they say, "all models are wrong and some are useful."

This one is obviously wrong and not very useful. If you think you can improve it please do.

If I get the time and the inclination I will generate some numbers from my RCalc program to try to make a different shaped curve. I'll need to make a branch of the code so I can add add some fudge factors (something I  because the aim of the code was to find the bell curve, not fit it to real life) the results of R0=3 has 70 million people infected in 17 iterations and the opportunity factor doesn't get a chance to pull the curve down. The spreadsheet calculations are obviously crude and somehow wrong. My focus was to get a bell curve with the right shape and the same end total. The fact it had the same width is a mystery to me. It is not designed in to the formula. And I messed up the horizontal axis labels.

As ever all source data is on github here:

https://github.com/davidsinfield/RValue/blob/master/Underwhelmed.ods

It's freely available, do with it what you will.

Before you comment

• There is no doubt that this is a dangerous disease to the vulnerable. I am not trivialising it.
• The hospitals have done a grand job and the lack of masks and gowns was a fiasco. I'm not sure what could have been done about it. And I have no intention of being drawn in to a recriminations.
• I have no wish to kill you or your granny. I'm a grandad I'm well into the age danger zone with a highish BMI and I would be considered quite vulnerable. 
• The third chart is very definitely wrong. I don't need to be told I am open to suggestions as to how it can be improved.

Covid-19 Correlations

Where I started

I was thinking about why there was such a wide spread of values of deaths per million from or of Coronaviris across the 48 countries of Europe so I decided to look at some correlations. It's important to note these are not direct country to country comparisons, I think we all know by now they are very dangerous. Rather they are a statistical method to see what the possibility that a pair of sets of values are somehow linked. The strength of a correlation is measured in the P (for Pearson) value and they range form -1 negative correlation through 0 (no correlation) to +1 positive correlation. A negative correlation means that as one value gets bigger, the other value gets smaller. A positive correlation means that the values increase or decrease together. I used to do rank correlations by hand for my OPen University courses but that is hard work and we now have computers that can do that stuff for us. I used LibreOffice calc and it a matter of fractions of a second to get an answer that would certainly taken hours using the old pen and paper method.

"Correlation is not causation" is a very important concept. It might be that there is a spurious correlations. There are some great spurious correlations of time series here:

https://www.tylervigen.com/spurious-correlations

However given a strong correlation of P>.5 it is certainly legitimate to hypothesise for causation.

Population Density

My first instinct was that it was due to population densities. It seems obvious that if you cram people together the spread of disease will be faster and you may get a higher mortality per capita than more sparsely populated countries. So what does the correlation look like?

The P value is 0.234105007129706 which represents a weak positive correlation. Not very much above chance. I was surprised.

That is certainly not a good enough correlation warrant searching to hard for a causation

Life Expectancy

We generally think of Life expectancy as an indicator of the quality if health and social care and general health of the population. The longer the live expectancy the better we were doing. So you'd expect the longer the live expectancy, the better the health and social care available and there should be a negaitive correlation, that is as life expectancy goes up, mortality goes down.

The European values range from 71.72  to 85.42 a pretty wide range for what we would all expect to be developed industrial countries.

I tried life expectancy against deaths per million. P value is 0.529. That is certainly more significant. Yes, thats a positive correlation. as life expectancy goes up so does a country's mortality rate for the Chinese Cough.

It seems therefore that a high life expectancy is a disadvantage where it comes to mortality from this virus. So maybe it is not a good indicator of the quality of health, social care and fitness after all. Maybe it's just an indicator that we can keep sick people alive longer. I'm not that sure that's a good thing when a pandemic strikes.

There is an alternative explanation. It could be that the reporting in countries with a lower life expectancy is less accurate and they are underreporting deaths from or with to coronavirus.

I'll let you decide.

Data

I did this a while ago and the data is from Worldometer from Wednesday 29th April. The Uk has changed its reporting method since then to include non NHS deaths. I don't know whether this has a large effect but I'll do a revision if I get a chance.

The spreadsheet of values I have used is here 

https://github.com/davidsinfield/RValue/blob/master/Correlations1.ods

Risk takers and risk averse

Over the two hundred thousand years of the human race we have evolved a strategy for the survival of the species. It's worked pretty well. One of the cornerstones of the strategy finding a balance between the risk takers and the risk averse. The argument about whether this is nature or nurture is

Here's a very quickly drawn diagram (I have a love hate relationship with Inkscape.)

There's no scale and no attempt at quantifying it. The numbers might be distributed differently but it is what it is.

The important parts as is very often the case are the two edge cases.

The red zone are the risk takers, this is where innovation comes from. The first person to tame a horse, Montgolfier, the Wright Brothers, Edmund Hillary are all in the the red zone. So are all all the entrepreneurs and small businesses. This is the source of all innovation. Some fail. But the ones that succeed make the leaps forward that eventually filter up the other groups.

The green zone are the risk averse. They like the status quo. New ideas and things will only filter down to them after the early adopters and late adopter have proven them. It's very difficult to name them individually because their impact on the species is taken for granted. They take the absolutely safe path. They will survive. They will breed, they will continue the species. Luckily they will breed some a few more risk taker than risk averse because some of the risk takers fail in a spectacular fashion.

Over most of the 200,000 years of the human race this balance has been maintained successfully.  We have thrived. The risk takers have given us innovation and the risk averse have ensured we haven't entirely wiped ourselves out.

The reason for the balance is that we have been allowed to make our own risk assessments.

Consider what happens when we are forced by regulation to be come risk averse....

Playing with Reproduction Value

Background

When the chart of infection spread for CoronaVirus was published it intrigued me that both the tall thin curve representing unrestricted spread and the short fat curve representing controlled spread were normal distribution bell curves.

I wanted to know why the curves rose and fell and in particular what made the curves fall in such a symmetrical way. Particularly what caused them to fall. I decided the best way was to make my own model so I can play with it.

When I started reading I realised Basic Reproduction Value, R0 (are-nought) was the most important factor. It's pretty well the only property of the spread of the virus that at this time we are able to control.

The tall thin curve is the number of hospitalisations with time with no intervention. The native property of the virus.

The short fat curve is the number of hospitalisations if the R0 value is attenuated by reducing the opportunity to spread. But why the bell? Why not an ever increasing exponential curve until everyone who will be hospitalised is?

Assumptions

The charts use hospitalisations. That's so the dotted line (available hospital resources) is relatable to the numbers. I have assumed that infections follow that curve. There is a skewing factor that Is not taken into account in those curves and in my calculations and that is that there should be more hospitalisations (though not infections) in as the more vulnerable will be hospitalised early.

I have used a simple model for the reproduction value iteration time. In reality the virus is infectious for a long period. I have seen number like 7-10 days in the incubation period and up to several weeks in the symptomatic period. I haven't modelled for this. For convenience I've called an iteration period a day.

The model uses a chance or opportunity factor. I have assumed that the opportunity factor is linear with the number of people who have been infected. There has been some talk of a small number of reinfections in South Korea. If infection does not confer immunity in the vast number of cases then all bets are off and we are all dead.

Population density will have a huge effect on R value (more specifically on oportunity factor.) I haven't even attempted to model how high density areas within the same population affect the overall shape of the curve. My speculation is that there will be a much longer tail as more isolated communities become infected over time.

The model has no 'fudge factors' to aim to give a preconceived result. The R based calculations are taken directly from the definition of R

Method

The model is very simple. Day zero there are n infectious cases. Which means on day two there will be n*R. The total infected is n+(n*R)* The next day the infections are the (total infected on the previous day) * RO. They are added to the running total. 

I then apply an opportunity factor which is the chance of encountering RO number of people who are do not have acquired immunity. Which is the population/the running total of infected people.

The model runs until the number of new infections is zero or given number of iteraterations or a target number of infections has been reached

.

I have pasted the code after the footnotes. It runs under Python 3. I use version 3.6 on various versions of Ubuntu. This has been tested on 18.04.

Results

I'm presenting the results in a standard form:

• Parameters
• Raw output
• csv
• chart
• analysis

The charts are created using LibreOffice Calc

Run one

I wanted a picture of how the tall thin curve might look so I used the profile for infecting the whole population at R0=3.

The parameters for my first run were:

population = 70m

day zero infections = 1

iterations 100

R value = 3

I chose the R0 value because Niall Ferguson had revised his estimate of 200-500k deaths to 20k based on a R) value of "probably greater than 3."

Here are the results. 

If you want the comma separated values to load into a spreadsheet you can cut and paste these:

=================================================================

day,infected,todaysinfections, opportunityfactor

1,4.0,3.0,1.0

2,12.999999614285715,8.999999614285715,0.9999999571428572

3,39.99999382857178,26.999994214286065,0.9999998285714341

4,120.99993134287536,80.99993751430358,0.999999442857231

5,363.99932731491725,242.9993959720419,0.9999982857152666

6,1092.9937348617311,728.9944075468138,0.999994814295324

7,3279.942840759639,2186.9491058979074,0.999984400089502

8,9840.48283497705,6560.539994217413,0.9999531579594177

9,29519.33628960095,19678.8534546239,0.9998594359595003

10,88531.0014957269,59011.66520612595,0.9995783094815771

11,265342.09842229204,176811.09692656514,0.9987352856929181

12,793764.7356017698,528422.6371794777,0.9962094128796816

13,2361056.530291741,1567291.794689971,0.9886605180628318

14,6904340.644671642,4543284.114379901,0.9662706352815467

15,19189833.989471916,12285493.344800275,0.9013665765046909

16,45942461.21758762,26752627.228115708,0.7258595287218298

17,73525420.95292085,27582959.73533324,0.3436791397487483

==============================================================

The chart of iterations v infections is:

Definitely starting to get a bell curve. The exponential part works. And the opportunity factor starts to kick in but doesn't have chance to pull the curve down before the target infections is reached.

The code seems to work as expected although I was surprise the opportunity factor didn't get a chance to kick in.

Run 2

For this run I wanted to try to get the full bell curve. And I did with these parameters.

population = 70m

day zero infections = 1

iterations 100

R value = 2

Raw output:

I was quite pleased to see the bell in the results.

If you want to play with them in a spreadsheet, here are the comma separated values:

=====================================================================
day,infected,todaysinfections, opportunityfactor
1,3.0,2.0,1.0
2,6.999999885714286,3.9999998857142858,0.9999999714285714
3,14.999998971428605,7.999999085714319,0.999999914285716
4,30.99999394285784,15.999994971429238,0.9999998000000146
5,62.99997017143768,31.999976228579836,0.9999995714286578
6,126.99986594295245,63.999895771514765,0.9999991142861404
7,254.99942708660234,127.99956114364987,0.999998200001915
8,510.99762046489644,255.9981933782941,0.999996371436756
9,1022.9902769794999,511.9926565146035,0.9999927143197077
10,2046.9606399653685,1023.9703629858686,0.9999854001389004
11,4094.8415087068443,2047.8808687414755,0.9999707719908576
12,8190.363711912501,4095.5222032056563,0.9999415165498757
13,16380.449840583464,8190.086128670962,0.9998830090898299
14,32756.789266355827,16376.339425772363,0.9997660078594202
15,65494.14183438098,32737.35256802515,0.9995320601533376
16,130907.58776828786,65413.44593390688,0.9990643836880803
17,261489.82103606238,130582.23326777453,0.9981299058890245
18,521678.6934508844,260188.87241482205,0.9962644454137707
19,1038178.3031141586,516499.6096632741,0.9925474615221302
20,2055857.003245857,1017678.7001316986,0.9851688956697977
21,4031437.2359366487,1975580.2326907916,0.9706306285250591
22,7755042.680259095,3723605.444322446,0.9424080537723336
23,14377204.556859408,6622161.876600313,0.8892136902820129
24,24901294.90187841,10524090.345019003,0.7946113777591514
25,38461947.971023865,13560653.069145458,0.6442672299731655
26,50681279.27394551,12219331.302921642,0.4505436146996591
27,57425903.87736294,6744624.60341743,0.2759817389436356
28,59848977.1579086,2423073.2805456645,0.1796299588948152
29,60551739.29053585,702762.1326272473,0.14501462631559137
30,60741450.163351946,189710.87281609673,0.13497516727805925
31,60791634.38507288,50184.22172093443,0.13226501195211507
32,60804837.662555486,13203.277482604371,0.13154809449895885
33,60808306.413785264,3468.751229775726,0.13135947624920735
34,60809217.37669668,910.9629114187214,0.13130992266021052
35,60809456.589925475,239.21322879153985,0.1312969089043331
36,60809519.40420555,62.81428007674185,0.1312934915724932
37,60809535.898305126,16.494099571382108,0.13129259422563497
38,60809540.2294036,4.331098471276219,0.13129235859564106
39,60809541.366683334,1.137279731252975,0.1312922967228057
40,60809541.665315434,0.29863209891056264,0.13129228047595237

=================================================================

As you can see the run ended because there were no new infections. There were still ~10 million uninfected.

I had acheived my bell curve:

Run 3

For run 3 I wanted to find an R0 value that gave me bell curve and infected every one.

I could have added some goal seeking to my code but decided to do it by hand. I knew the value was greater than 2 and less than 3 so I did increments of 0.1. The sensitivity is greater than I expected so I had to go to a second decimal place.

I won't post all the results, just the the one that got me closest:

population = 70m

day zero infections = 1

iterations 100

R value = 2.56

csv:

========================================================

day,infected,todaysinfections, opportunityfactor

1,3.56,2.56,1.0

2,10.113599760325487,6.553599760325486,0.9999999634285714

3,26.89081296246121,16.777213202135727,0.9999998698057176

4,69.84046287418919,42.949649911727974,0.9999996301312434

5,179.79145851842884,109.95099564423965,0.999999016564816

6,461.26528843560914,281.4738299171803,0.9999974458363069

7,1181.8335550987078,720.5682666630987,0.9999934247815937

8,3026.457200181365,1844.6236450826568,0.9999831309492129

9,7748.489632671324,4722.032432489959,0.9999567791828545

10,19835.55473445807,12087.065101786746,0.9998893215766762

11,50769.67370390198,30934.11896944391,0.9997166492180792

12,129903.58341521795,79133.90971131597,0.9992747332328012

13,332110.4489872361,202206.86557201814,0.998144248808354

14,847304.0774886198,515193.6285013837,0.9952555793001824

15,2150235.418221515,1302931.3407328953,0.9878956703215912

16,5383280.8504639575,3233045.4322424424,0.9692823654539783

17,13023373.812327115,7640092.961863157,0.9230960021362292

18,28943162.735471915,15919788.9231448,0.8139518169667554

19,52846840.87116953,23903678.135697614,0.5865262609218297

20,67841990.3561394,14995149.484969873,0.24504514469757815

21,69025430.52783662,1183440.1716972135,0.030828723483722762

22,69067610.03309497,42179.50525835205,0.013922435316619703

23,69069048.30654654,1438.2734515715865,0.01331987095578611

24,69069097.27439271,48.96784616825324,0.0132993241921923

25,69069098.94147752,1.667084815835123,0.013298624651532727

26,69069098.99823245,0.05675493254555362,0.013298600836035396

=======================================================

This is quite a long post so I'll end the results there. 

Conclusions

Whether you call it 'herd immunity' or 'widespread acquired infection' the opportunity factor serves to form the expected bell curve in all cases where R0 is between 1.2 (results not shown here for brevity. No new cases iteration 158, 22602160 infected) and 2.56.

R0 of 2 or less leaves a large population of potential infections if the R0 is increased by removing restrictions unless there are no or very few infections when restrictions are relieved.**

Even with this simple model the theoretical expected and possibly the real world results for cases/deaths (we wont know until at least one country has fully played out their infection.)

Any discussion of comparison with actual infections is futile because there is no way of telling what they are except other models.

A more valid result might be to take the deaths and compare them to infections calculated by this model to get an estimate of mortality. This might become another blog post at some time maybe.

Comments

Most comments are welcome either here or by mentioning me on  twitter. I will answer as many and as soon as I can.

I will ignore stuff like "my old aunty died of CoronaVirus so you are wrong." Your anecdotes do not interest me.

Or anything which refers to something I didn't actually say.

If you are commenting please be specific. "Cobblers" or any other such stuff will be either ignored or published to show what an idiot I am dealing with.

If you have a set of parameters you wish me to run please send them. Only the parameters show in the results section please.

If you want me to tailor the model or submit change requests then you'd better be ready with some paypal. I might or might not improve, change, remove it at any time. My choice.

You may take what additional inferences from the results. They are not my responsibility.

You use the information here at your own risk. Any decisions you base on it are your problem not mibe.

Footnotes

*The brackets are extraneous but I know a lot of people have problems with operated precedence and I didn't want to use RPN although I'm probably one of the last people alive who can.

** I may revise my model so the R0 value can be changed in the course of the curve. That might  show any subsequent peak(s).

Code:

Also available on GitHub here

If anyone requires the results files or images I can zip them up and make them available.

=========================================================
#---- License ----
"""
This code is the property of David Sinfield (dsinfield6@gmail.com) It is  released under The GNU public license the short version of which is you may copy, 
distribute, amend but this section of header should remain intact.

The long version is at https://www.gnu.org/licenses/gpl-3.0.html
"""

#------------------------------------------------------------------------------------------------------------------------------
# This is a Transmissibility generator V1.0
# The parameters are prompted for on the command line. It takes no command line parameters.
#blog post description is at https://davidsdiamondsandrust.blogspot.com/2020/04/playing-with-reproduceability-value.html
#------------------------------------------------------------------------------------------------------------------------------

#usage
print ("This program models the number of people infected in over a number of days for a given R value.")
print ()
print ("Day Zero are the number of infections on Day 0")
print("The R value is the number of people each infected person can transmit to.")
print(" R may be a decimal value. The result on day N will be rounded up to an integer.")
print("The program assumes that R declines linearly with the number of infections as the chance of contacting R uninfected declines")
print ("the chance parameter calculation may be a little simplistic but is calculated as (population-infected)/population ")
print ("The definition of R does not include a time element so this assumes one iteration a day.")
#functions
#all the get functions require validity checking for input type numeric.
def getpopulation():
    ipopulation=input("Population(default 70m) - ")
    if len(ipopulation) > 0:
        population=int(ipopulation)
    else:
        print("70m used")
        population=70000000
        return population
    
def getdayzero():
    idayzero=input("Day Zero infections (default =1) - ")
    if len(idayzero)>0:
        dayzero=int(idayzero)
    else:
        print ("1 used")
        dayzero=int(1)
    return dayzero
    #ToDo validty check/default
def getdays():
    return int(input("Days - "))
      #ToDo validty check/default
def getrvalue():
    return float(input("R Value - "))
def gettarget():
    itarget=input("Target value (leave blank for population=target) - ")
    if len(itarget)>0:
        target=int(itarget)
    else:
        print(str(population)+" used")
        target=population
    return target
#end functions

#start procedural
population=getpopulation()
dayzero=getdayzero()
days=getdays()
rvalue=getrvalue()
target=gettarget()

infected=dayzero
newcases=dayzero
csv=["day,infected,todaysinfections, opportunityfactor\n"]
for day in range (1, days):
    # calculations
    chance=((population+1)-infected)/population
    todayscases=float((newcases*rvalue)*chance)
    infected=infected+todayscases
    newcases=float(todayscases)
    #show results
    print ("Chance: "+str(chance)+" Day number: " + str(day) + " - Infected: " + str(int(infected))+" Newcases today:"+str(int(todayscases+0.9)))
    csv.append(str(day)+","+str(infected)+","+str(todayscases)+","+str(chance)+"\n")
    #exit conditions
    if target > 0:
        if infected > target: 
            break
    if newcases<1:
        print("No new cases")
        break
print("Program end")

#Option to write a file of data
dofile=input("Enter y to write a result file: ")
if dofile=="y" or dofile== "Y":
    filename=input("Filename (exension .csv will be added): ")
    if len(filename)==0:
        print ("No Filename results.csv will be used")
        f=open("results.csv","a")
    else:
        f=open(filename+".csv", "a")
    print ("Writing "+ f.name)
    for x in csv:
        f.write(x)
    f.close
# end procedural
 

===============================================

Optimists, rationalists and pessimists

Thing change. It's the only constant.

And as things change we like to guess outcomes.

Optimists like to guess things will be fine.

Pessimists guess at bad outcomes.

Rationalists recognise that of the outcomes some are better than others.

On this continuum I sit between rational and optimistic. What steers me from pure rationalism?

1. Faith in human ingenuity to steer the outcomes towards the good.

2. The fact that from a historical perspective optimists are more often proved right. Unless you believe a nasty Brutus and short life is better than we have now, overall things improve.

Tell me I'm wrong.