Showing posts with label software. Show all posts
Showing posts with label software. Show all posts

Sunday, May 3, 2020

Playing with Virus data

Thanks
to the  Sage for the mathematical framework,
to the ECDC for the worldwide pandemic data,
to the China for the virus,

we can now play with the statistics and create some interesting graphs,
  • like the current number of days to double number of victims [should be high]
  • or the exponent of the current cases grows (if considered exponential) [should be low].
Copy the program, download the data (in CSV) and play .....


import csv
from datetime import date

# ------------------------------------------------------------------------------

date0 = date(2019, 12, 31)

# ------------------------------------------------------------------------------

# Functions defining graphs

# Cases / Deaths

def cases(v):
  return anything(v, 'cases')
  
def deaths(v):
  return anything(v, 'deaths')
 
def anything(v, c):
  v1 = {}
  #for x in v:
  #  v1[x] = v[x][c]
  for k in sorted(v.keys()):
      v1[k] = v[k][c]
  return v1
  
# Daily increase [%]
  
def delta_cases(v):
  return delta_any(v, 'cases')

def delta_deaths(v):
  return delta_any(v, 'deaths')
  
def delta_any(v, c):
  v1 = {}
  sum = 0
  for k in sorted(v.keys()):
    sum += v[k][c]
    if sum != 0:
      v1[k] = 100 * v[k][c] / sum
  return v1
  
# Number of days to double numbers 
  
def double_cases(v):
  return double_any(v, 'cases')
  
def double_deaths(v):
  return double_any(v, 'deaths')
   
def double_any(v, c):
  v1 = {}
  sum = 0
  for k in sorted(v.keys()):
    sum += v[k][c]
    if v[k][c] != 0:
      v1[k] = sum / v[k][c]
  return v1
 
# Exponet for exponential fit of latest days 
 
def exp_cases(v):
  return exp_any(v, 'cases')
 
def exp_deaths(v):
  return exp_any(v, 'deaths')
  
def exp_any(v, c):
  v1 = {}
  y0 = -1
  for k in sorted(v.keys()):
    if y0 > 0 and v[k][c] > 0:
      v1[k] = log(float(v[k][c])) - y0
    y0 = log(float(v[k][c]))
  return v1
 
# Mortality 
 
def mortality(v):
  v1 = {}
  sum = 0
  for k in sorted(v.keys()):
    sum += v[k]['deaths']
    if v[k]['cases'] != 0:
      v1[k] = 100 * v[k]['deaths'] / v[k]['cases']
  return v1
  
# Sigmoid fit  

def sigmoidFit(name, values):
  data = []
  sum = 0
  for v in values:
    sum += v[1]
    data += [(v[0], sum)]
  var('L, k, x0, x')
  model(x) = L / (1 + exp(- k * (x - x0)))
  sol = find_fit(data, model)
  print(sol)
  f(x) = model(L = sol[0].rhs(), k = sol[1].rhs(), x0 = sol[2].rhs())
  h = plot([])
  h += text(name,                                                                      ( 30, sol[0].rhs() / 2      ), fontsize = 10, color = 'black')
  h += text('today is ' + str(len(data) - round(sol[2].rhs())) + ' days after peak', (130, sol[0].rhs() / 2 * 0.9), fontsize = 10, color = 'black')
  h += plot(f(x), x, [0, len(data) * 2])
  h += list_plot(data, color = 'red')
  show(h)

# Reading and interpreting country data

def readCountry(country, f, normalised, integrated, offset, smear):
  val = {}
  with open('virus.csv', newline='') as csvfile:
    reader = csv.DictReader(csvfile)
    for row in reader:
      date1 = date(int(row['year']), int(row['month']), int(row['day']))
      days = (date1 - date0).days
      if row['countriesAndTerritories'] == country:
        v = {}
        if normalised:
          v['cases']  = float(row['cases'])  / float(row['popData2018'])
          v['deaths'] = float(row['deaths']) / float(row['popData2018'])
        else:
          v['cases']  = int(row['cases'])
          v['deaths'] = int(row['deaths'])
        val[int(days)] = v
  if integrated:
    valI = {}
    sumC = 0
    sumD = 0
    for k in sorted(val.keys()):
      sumC += val[k]['cases']
      sumD += val[k]['deaths']
      valI[k] = {'cases':sumC, 'deaths':sumD}
    val = valI
  if smear > 0:
    valS = {}
    for k in sorted(val.keys()):
      sumC = 0
      sumD = 0
      n = 0
      for i in [0..smear]:
        if k + i in val.keys():
          n += 1
          sumC += val[k + i]['cases']
          sumD += val[k + i]['deaths']
      valS[k] = {'cases':(sumC / n), 'deaths':(sumD / n)}
    val = valS
  if offset > 0:
    valO = {}
    d    = 0
    for k in sorted(val.keys()):
      d += 1
      if d > offset:
        valO[k] = val[k]
    val = valO
  return f(val)

# Creating histogram

def plotHisto(values, logarithmic, col, country, doSpline):
  data = []
  maxx = 0
  for val in values:
    data += [(val, values[val])]
    if val > maxx:
      maxx = val
  sc = 'linear'
  if logarithmic:
    sc = 'semilogy'
  p = list_plot(data, scale = sc, color = col)
  if (doSpline):
    s = spline(data)
    p += plot(s, 0, maxx,  scale = sc, color = col, legend_label = country)
  return p

# ------------------------------------------------------------------------------
       
# Define group of countries to analyse            
countries = [
             ('France',                   'red'),
             ('Italy',                    'green'),
             ('Germany',                  'grey'),
             ('Belgium',                  'pink'),
             ('Russia',                   'purple'),
             ('Spain',                    'orange'),
             ('Czechia',                  'blue'), 
             ('Sweden',                   'brown'), 
             ('Japan',                    'yellow'), 
             ('United_Kingdom',           'violet'), 
             ('United_States_of_America', 'black')
             ]
countries0 = [
              ('France',                   'red'),
              ('Italy',                    'green'),
              ('Czechia',                  'blue')
              ]
countriesF = [
              ('France',                   'red'),
              ]
       
# Set graphs to create       
# - for analytical graphs              
optionsList = [
               {'normalised':true,   'integrated':false, 'logarithmic':false, 'smear':5,  'offset':60, 'doSpline':true, 'countries':countries,  'f':[cases, deaths]},
               {'normalised':true,   'integrated':false, 'logarithmic':false, 'smear':0,  'offset':60, 'doSpline':true, 'countries':countries,  'f':[cases, deaths]},
               {'normalised':true,   'integrated':true,  'logarithmic':false, 'smear':0,  'offset':60, 'doSpline':true, 'countries':countries,  'f':[cases, deaths]},
               {'normalised':false,  'integrated':false, 'logarithmic':false, 'smear':0,  'offset':60, 'doSpline':true, 'countries':countries0, 'f':[cases, deaths]},
               {'normalised':false,  'integrated':true,  'logarithmic':false, 'smear':0,  'offset':60, 'doSpline':true, 'countries':countries0, 'f':[cases, deaths]},
               {'normalised':true,  'integrated':true,  'logarithmic':true,  'smear':5,  'offset':60, 'doSpline':true, 'countries':countries,  'f':[cases, deaths]},
               {'normalised':false, 'integrated':false, 'logarithmic':false, 'smear':5,  'offset':60, 'doSpline':true, 'countries':countries,  'f':[delta_cases, delta_deaths]},
               {'normalised':false, 'integrated':false, 'logarithmic':false, 'smear':5,  'offset':60, 'doSpline':true, 'countries':countries,  'f':[double_cases, double_deaths]},
               {'normalised':false, 'integrated':false, 'logarithmic':false, 'smear':5,  'offset':60, 'doSpline':true, 'countries':countries,  'f':[exp_cases, exp_deaths]},
               {'normalised':false, 'integrated':false, 'logarithmic':false, 'smear':10, 'offset':60, 'doSpline':true, 'countries':countries,  'f':[mortality]},
               ]
# - for sigmoid fit
countries4fit = countriesF
  
# ------------------------------------------------------------------------------
   
# Main loop for analytical graphs
for options in optionsList:
  normalised  = options['normalised']
  integrated  = options['integrated']
  logarithmic = options['logarithmic']
  smear       = options['smear']
  offset      = options['offset']
  doSpline    = options['doSpline']
  countries    = options['countries']
  for f in options['f']:
    if f == delta_cases or f == delta_deaths or f == double_cases or f == double_deaths:
      integrated = false
    if f == exp_cases or f == exp_deaths:
      integrated = true
      normalised = false
    #if integrated:
    #  offset = 0
    h = plot([])
    for (country, color) in countries:
      h += plotHisto(readCountry(country, f, normalised, integrated, offset, smear), logarithmic, color, country, doSpline)
    x = (3 * h.get_minmax_data()['xmin'] + h.get_minmax_data()['xmax']) / 4
    y = (3 * h.get_minmax_data()['ymin'] + h.get_minmax_data()['ymax']) / 4
    t = f.__name__
    if normalised:
      t += ",\nnormalised"
    if integrated:
      t += ",\nintegrated"
    if logarithmic:
      t += ",\nlogarithmic"
    if smear > 0:
      t += ",\nsmear = " + str(smear)
    h += text(t, (x, y), fontsize = 10, color = 'black')
    show(h)

# Main loop for sigmoid fit
for country in countries4fit:
  for h in [cases, deaths]:
    values = readCountry(country[0], h, false, false, 0, 0)
    a = []
    for v in values:
      a += [(v, values[v])] 
    sigmoidFit(country[0] + "\n" + h.__name__, a)

# ------------------------------------------------------------------------------

Thursday, April 23, 2020

SIR with Sage

Very simple implementation of the SIR epidemic model using Sage
mathematical framework.
# Infected
IInit = 1
# Susceptible
SInit = N - IInit
# Resistant
RInit = 0
# Transmission rate
betaInit = 0.5
# Transmission rate attenuation
delta = 0.01
# Recovery rate
gamma = 0.1
# Population size
N = 10000
# End time
tMax = 100
# Number of points
points = 1000

# Standard SIR model
def ODE_RHS(t, Y):
  (S, I, R, beta) = Y    
  dS = - beta * S * I / N
  dI = beta * S * I / N - gamma * I
  dR = gamma * I
  dbeta = - beta * delta
  return (dS, dI, dR, dbeta)

# Set up numerical solution of ODE
solver = ode_solver(function = ODE_RHS,
                    y_0 = (SInit, IInit, RInit, betaInit),
                    t_span = (0, tMax),
                    algorithm = 'rk8pd')

# Numerically solve
solver.ode_solve(num_points = points)

# Plot solution
show(plot(solver.interpolate_solution(i = 0), 0, tMax, legend_label = 'S(t)', color = 'green')
   + plot(solver.interpolate_solution(i = 1), 0, tMax, legend_label = 'I(t)', color = 'red')
   + plot(solver.interpolate_solution(i = 2), 0, tMax, legend_label = 'R(t)', color = 'blue')) 

Result:

Wednesday, April 22, 2020

Sigmoidal Fit

Interesting.
When I do extremely naive sigmoidal fit to the evolution of the number of virus victims in France, I get about the same prevision as prof.Raoult - that the epidemie will end very soon.
This can mean two things:
- Either I'm as clever as prof.Raoult.
- Or the other way.

Thursday, February 28, 2019

30th/26th Anniversary of the World Wide Web


On March, 12, it will be 30 years from the start of the World Wide Web in CERN.


And 26 years from the start of the World Wide Web in the Central Europe - in the Prague Institute of Physics, in spring 1993.

and
(24 years old, in czech).

Friday, October 14, 2016

HEP SW = Blind Prophet

In HEP software, we often hit a problem before others. We find our solution to it.

Then, the problem becomes widespread. Better solutions are found elsewhere.

But we are stuck with our home-grown hacked-together solution.

Examples ?
CW ntuples, reflective language, Map/Reduce data analyses, Big Data, 3D graphics, projection histograms,...

Friday, June 6, 2014

C++ is Fun

C++ is a great language. It allows to solve all problems ..... which don't exist in other languages.

And fixing C++ deficiencies is fun - more fun than doing something useful. 

Friday, October 18, 2013

CHEP2013 Observations


C++11: So in far future C++ may have some of the features of JVM (Java++) had long time ago. Impressive.  Not really.

Recurrent pattern: Root inability to interface to X is interpreted as a weakness of X. X is bad because Root doesn't like it. Compatibility with Root is taken for the measure of quality.

Many presentations & posters about Hadoop: o o o o o o o o o o

Friday, March 22, 2013

Android Applications for Particle Physics


 AWork by S.Chekanov

AWork is an interactive environment for scientific computation, data analysis and data visualization designed for scientists, engineers and students. AWork is an Android port of the jHepWork (jWork) Java project (for data analysis) and jMathLab Java project for numerical and symbolic calculations.

The AWork project supports symbolic and numeric calculations using Octave-like high-level interpreted language. A similar desktop application is called jMathLab. For data-analysis using Java numerical libraries, it uses BeanShell scripting language, with the full access to Java API and jHepWork collection of Java numeric libraries.

 LHSee by Ch.Boddy

LHSee: Discover what's happening at the LHC.

You want to find out how to hunt the Higgs boson with her cell phone? Ever wondered how the experiments work at the Large Hadron Collider, and the collisions look like?

Scientists at the world's largest scientific experiment - the Large Hadron Collider at CERN in Geneva - try to answer basic questions about the universe, the origin of mass and the structure of space and time. LHSee makes the hadron collider for all who have an Android smartphone or tablet PC, accessible. Written by scientists at Oxford University in collaboration with the ATLAS experiment at CERN, this app is for use by experts and laymen alike conceived.

For the first time, you can now get transferred from the collision events in Geneva underground detectors directly in real time on your device. The app includes a variety of educational resources, and allows you to interact with the 3D shown in collision events. You can find out how to work the different parts of the detector and learn how to recognize different types of collisions. With the game "The hunt for the Higgs boson," Your new skills will be put to the test.

 arXiv mobile by Nonlinear Apps

Keep up with the latest science research: physics, math, computers science arXiv mobile is science on your phone! Browse daily science articles (Physics, Astronomy, Math...) at arXiv.org and search the entire arXiv collection. This is a must have app for every scientist and physicist. You can read PDFs, save them for later reading and share articles via email/twitter etc... A homescreen widget lets you know of any new articles in your favorite science categories (dozens of physics, math, computer science categories) or custom searches.

 PDG Particle Listings by Ch.H.Christensen

An Android version of the famous PDG booklet (almost). This app shows the basic properties of various particles, as well as their decay modes.

 ATLAS Status by Atlas

An expert monitoring app for ATLAS on-call experts. This app requires CERN login rights in order to be used. It presents the current operational status of the ATLAS partition, including Run Status (i.e. Running, Configured, Initial, etc.), a summary of the busy status, the included hardware, and the LHC page 1 message.

 Particle Physics Booklet by I.Kreslo

Mobile version of PDG Particle Physics Booklet 2010.

 Particle Properties by I.Kreslo

Offline database for properties of elementary particles. Based on the table of masses, widths, and PDG Monte Carlo particle ID numbers published by Particle Database Group, issue of 2008.

 The Standard Model by Appocalypse

The standard model of particle physics reference app. As complete as they come. Will update multiple times.

 ioda by G.Barrand

ioda should be read "IO-DA", IO for Input/Output and DA for Data Analysis. It is a scientific application that permits to read files at various formats as FITS used in astronomy, AIDA and ROOT used in high energy physics (HEP), JPG to store images and FOG developed at CEA/Saclay (France) to describe the LHC/ATLAS detector geometry. ioda permits to browse these files and visualize some of their data.

For AIDA files, the histograms 1D, 2D, profiles 1D, 2D can be plotted. For ROOT files, ioda can plot TH[1,2][D,F] objects along than TProfile ones. For JPG files, the image is visualized. For FITS files, the "HDUs" can be listed and their keys can be seen. If the HDU is an IMAGE_HDU type, ioda attempts to visualize it as an image (which is one of the primary goal of the FITS format : store astronomical pictures). If the HDU is a BINARY_TBL, ioda shows a description of the columns (name, type) and proposes to histogram and plot a choosen column. Files at the FOG CEA/Saclay format permit to visualize LHC/ATLAS sub detectors

 Agora by G.Barrand

This is a first version of a smartphone/tablet like event display for the LHC/ATLAS experiment. This first version is more a proof of concept (a "concept app") than an application ready to do physics. Anyway it can give to someone novice in high energy physics (HEP) a glance at what a HEP detector looks like.

Part of the geometry comes from an atlas.root file found on the CERN/ROOT web site (courtesy of the CERN/EP/SFT group) but also from FOG files created some time ago by L.Chevalier, C.Baudry, Marc Virchaux of CEA/Saclay (France). The su4small.dst file coming with the app and containing some events is a courtesy of the ATLAS LAL group. Special thanks here to L.Duflot for having prepared this very light event file. Else, the "way of doing" is the same than for the ioda application and the agora app is developed by using the inlib/exlib as explained in the ioda web pages.

 PMX by G.Barrand

This is a first version of a smartphone/tablet like event display for the LHC/LHCb experiment. This first version is more a proof of concept (a "concept app") than an application ready to do physics. Anyway it can give to someone novice in high energy physics (HEP) a glance at what a HEP detector looks like.

The geometry comes from an lhcbfull.root file found on the CERN/ROOT web site (courtesy of the CERN/EP/SFT group). The bdpi_ghost.dst file coming with the app and containing two events is the default event file, created long time ago, for the LHCb Panoramix event display (courtesy of the LHCb software team). Else, the "way of doing" is the same than for the ioda application and the pmx app is developed by using the inlib/exlib as explained in the ioda web pages.

 g4view by G.Barrand

g4view is a scientific application that permits to read geometry files at the Geant4 GDML format describing particle physics detectors. g4view permits to visualize the geometry of a loaded detector. It can read also "scenarios" files that permit to choose the volumes seen from GDML and strongly customize their graphical attributes as the coloring or the wire frame/solid rendering. A scenario permits also to customize an adapted particle source. From the main menu, you can then start a "particle through matter" simulation of the detector which is done by using the Geant4 toolkit. g4view can also read and execute "mac" files containing Geant4 scripting commands. Under the examples menu, there are predefined typical setups, as a piece of an electromagnetic calorimeter, intended to be used for particle physics detector teaching. We recall that Geant4 is a particle through matter simulation toolkit heavily used in high energy physics (for example at the CERN LHC experiments), but also for spatial and medical science.

 LHC Status by G.Barrand

LHCStatus permits to view the beam status of the CERN LHC accelerator.

 PACLC by J.Hrivnac

An Android client to Atlas Athenaeum (remote Athena) framework.

 Alidroid by Androworks

Android interface to MonALISA repository for ALICE.

 Radiation to Electronics by Ch.Daviet

Sunday, March 25, 2012

How to improve Samsung Galaxy SII (i9100)

  1. Install Ice Cream Sandwich (Android 4) (instructions in french)
  2. Get Root Access
    1. Install CF-Root-SGS2_XX_XEO_LPQ-PROPER-v5.4-CWM5 kernel (via CWM)
  3. Install AdWay to block Ads
How to mount SD card to USB-connected PC
  1. Settings
  2. Wireless and network / More...
  3. USB utilities
  4. Connect storage to PC
  5. Connect USB cable from phone to PC
  6. Turn on USB storage
How to save screenshot

Volume Down + Power Buttons
The screenshot will be added to /sdcard/Pictures/Screenshots folder then.

Friday, March 2, 2012

How to import Facebook Friends in Google+

Mr.Facebook tries hard to prevent people from migrating to Google+. At the beginning, it was possible to export Facebook directly or via Yahoo. Facebook has stopped those possibilities. But there is still one way, how to do it, in several minutes.
  1. You need an Android smartphone. If you don't have a physical one, you may try the emulator for PC or Mac (not tested).
  2. Enable Facebook on your smartphone with complete synchronization. You Facebook friends will appear in your Android Contacs.
  3. Export you Contacts into sdcard. They will be written into file 00001.vcf (or so).
  4. Copy that file into your PC/Mac and import it into your Google Mail.
  5. All your Facebook friends will then appear also in your Google Find People / Gmail list.

Wednesday, February 29, 2012

Google vs Facebook

ComScore a publié l'étude qui dit que les gens passent beaucoup plus de temps sur Facebook que sur Google +. Mais cette étude est gravement biaisée.

Si vous comparez FB avec Google, vous devez inclure la vidéo (FB vidéo vs Google / YouTube), photos (FB Photos vs Google Picasa), ... aussi Google Docs, Google Apps, Google Android, Gmail, Google Maps, Google Talk, ... - Puis le solde total est différent.

Dans le FB, tout (photos, vidéo, mail, chat, ...) est sous la même couverture. Dans Google, vous disposez d'un ensemble d'applications différentes.

Friday, February 17, 2012

New Swarm

New version (1.1.0) of Swarm Parallel Framework has been released.
You can try it as a WebStart application (Java 6 is required).



New features:
  1. BusMembers (Producers, Consumers, Algorithms) are associated to BusItems (data) by their types.
  2. Threads are dynamicaly balanced to provide the optimal data throughput and performance.
  3. Barious bug fixes, GUI cleaning and small improvements.

Monday, February 13, 2012

Clean Facebook

I've come to the conclusion, that there is no reason, why Facebook should own all my Internet history (as he states in his Contract). So I've cleaned up my Facebook Timeline and I will only keep entries not older than 1 month.

How to delete (all) your entries from the Facebook Wall, but keep your Friends.
  1. Use Firefox.
  2. Install iMacros for Firefox.
  3. Go to your Facebook Activity Log.
  4. In iMacros Edit/Options, change the Timeout Limit to 4 seconds.
  5. Go through all activities and for each:
    1. Start recording by Rec/Record.
    2. Perform task (deleting/hiding/unliking/...).
    3. Stop recording by Rec/Stop.
    4. Add following line to the script via Edit/Edit Macro: SET !ERRORIGNORE YES
    5. Execute corresponding script by Play (Loop).
    6. Sometimes, the script gets stuck and you will have to stop it, remove the current item by hand, or even reload the whole page.
  6. Even after all that, you will have to go through your Timeline and remove some last entries by hand.
Clearly, Mr.Facebook doesn't like you doing that. So he will try to create additional obstacles by modifying the interface etc.

Then, to keep Facebook clean, you can use Exfoliate Android application. It will remove all your Facebook entries older than same period. If you don't have an Android gadget, you can run Exfoliate in an emulator. (Exfoliate fails for very old Facebook entries so you can't use it for the initial cleanup.) 

It's also a good idea to remove (almost) all Facebook applications.

All my important posts&links will stay in my Blog and Google+ account, where I'm their owner.

Thursday, June 9, 2011

Athenaeum CLI JobConfig modifs via Groovy Slurp

I've implemented a possibility to modify JobConfig.xml file from the
Athenaeum CLI. There is a new argument -modifs which contains modification instruction
script string.
The syntax is very intuitive, it comes from the Groovy Slurp XML syntax.
athenaeum variable represents JobConfig.xml file to be modified.

Some examples:

// change email element
athenaeum.email="someone@cern.ch"

// change extract id attribute
athenaeum.extract.@id="newid"

// change logical collection physics condition
athenaeum.extract.logical_collections.collection_ref.physics_cond.@sql=
  'NLooseElectron>1 AND abs(LooseElectronPt1)>1000.'

// add a step element into profile element to describe this processing step
athenaeum.profile.appendNode {step (processor:"CLI-GroovyXMLSlurper",
                                    version:"1.2.3",
                                    ip:"123.123.123.123",
                                    time:"0")}

// more complex: for each collection with specified name: add RunNumber query to SQL
athenaeum.collections.collection.find
  {it.@name == "data10_7TeV_physics_MinBias_r1239_p134_READ"}.each
  {c -> c.query.@sql = (c.query.@sql + "RunNumber = 12345")}

Tuesday, June 7, 2011

Athenaeum Remote Worker Management

Athenaeum Workers can be managed remotely from any Athenaeum client provided there is at least one Worker alive.
Five actions are available:
  • status gives a status of all related Workers on the Worker machine
  • stop stops a Worker (for obvious reasons this command cannot be undone)
  • restart restarts the Worker
  • update updates the Worker to the newest version taken from the related Manager
The command result is not given immediately but will be sent to the provided E-mail address.

Linux CLI


Portable Client


Web Portal


Android Client