Me, blogging!

thoughts, science, code, jobs, and thoughts

Optimization problem solving with Pyomo

Next on this blog: An example about the optimization modeling language Pyomo.

 Tags: python, pyomo, optimization, models 2016/03/03
Exceedance probability curves

While working on a project proposal to estimate drought risks of some countries in the coming years, I learned something about risk analysis: It's rather simple! But, you need to be cautious if you want to use this information for guessing what might happen in the future. In the context of anthropogenic climate change that usually means that you can't rely on exceedance probabilities from observed (and hence past) data. Because the underlying climate is changing. And even if that's not the case, your record may be over- or underestimating the occurrences of specific events, just by pure chance.

However, let me show you how to quickly derive exceedance probabilities from some river discharge data.

From the monthly discharge data, we estimate the cumulative distribution function (CDF) by assigning each sorted monthluy value a number starting from 0 for the lowest up to 1 for the highest monthly discharge value. The exceedance probability is simply 1 - CDF. Plotting the exceedance probability against the sorted monthly data, et voilà, your exceedance probability curve. Here's the link to the CSV file elbe.csv and the Python code to process the data.

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from scipy import stats
import sys
import matplotlib
from matplotlib.ticker import ScalarFormatter, FormatStrFormatter, FixedLocator

matplotlib.style.use('bmh')

fnm  = 'elbe.csv'
raw = df['Elbe']
ser = raw.sort_values()

X = np.linspace(0.,max(ser.values)*1.1,100.)
s, loc, scale = stats.lognorm.fit(ser.values)

cum_dist = np.linspace(0.,1.,len(ser))
ser_cdf = pd.Series(cum_dist, index=ser)
ep = 1. - ser_cdf

fig = plt.figure(figsize=(11,8))

ax1 = plt.subplot(221)
ser_cdf.plot(ax=ax1,drawstyle='steps',label='data')
ax1.plot(X,stats.lognorm.cdf(X, s, loc, scale),label='lognormal')
ax1.set_xlabel('Discharge')
ax1.set_ylabel('CDF')
ax1.legend(loc=0,framealpha=0.5,fontsize=12)

ax2 = plt.subplot(223)
ax2.hist(ser.values, bins=21, normed=True, label='data')
ax2.plot(X,stats.lognorm.pdf(X, s, loc, scale),label='lognormal')
ax2.set_ylabel('Probability Density')
ax2.set_xlabel('Discharge')
ax2.legend(loc=0,framealpha=0.5,fontsize=12)

A = np.vstack([ep.values, np.ones(len(ep.values))]).T
[m, c], resid = np.linalg.lstsq(A, ser.values)[:2]
print m, c
r2 = 1 - resid / (len(ser.values) * np.var(ser.values))
print r2

ax3 = plt.subplot(222)
ax3.semilogx(100.*ep,ser.values,ls='',marker='o',label='data')
ax3.plot(100.*(1.-stats.lognorm.cdf(X, s, loc, scale)),X,label='lognormal')
minorLocator = FixedLocator([1,2,5,10,20,50,100])
ax3.xaxis.set_major_locator(minorLocator)
ax3.xaxis.set_major_formatter(ScalarFormatter())
ax3.xaxis.set_major_formatter(FormatStrFormatter("%d"))
ax3.set_xlim(1,100)
ax3.set_xlabel('Exceedance Probability (%)')
ax3.set_ylabel('Discharge')
ax3.invert_xaxis()
ax3.legend(loc=0,framealpha=0.5,fontsize=12)

ax4 = plt.subplot(224)
raw.plot(ax=ax4)
ax4.set_xlabel('time')
ax4.set_ylabel('Discharge')

plt.suptitle('Elbe discharge',y=1.01,fontsize=14)

plt.tight_layout()
plt.show()

fig.savefig('exceedance_curve_elbe.png',transparent=True,dpi=150,bbox_inches='tight')

 Tags: statistics, python, code, models 2016/02/20
Building a website generator

Hello everyone. This is my first posting on my blog where I want to provide some useful information about my job search. I started of as a physicist (TU-Berlin) about six years ago. Afterwards I went to Hamburg to do my graduate studies at the Hamburg University and the Max Planck Institute for Meteorology. After finishing my Ph.D. I went back to Berlin to work at the Potsdam Institute for Climate Impact Research as a postdoctoral researcher. My project about the melting of the Greenland Icesheet on different time scales ended two weeks ago. Now I'm taking my time to recollect what I've done so far and what I want to do next. This blog should help me and others to keep track about the expected and unexpected things to come, difficulties and the joy and pain of finding a new job outside of academia.

As of this week, I started writing a static blog engine in Python (pystable), which can be found on Github and cloned via

git clone https://github.com/mkrapp/pystable.git


Actually, this very web page is based on pystable.

 Tags: pystable, code, python, introduction 2014/11/28
f2py in action

Recently, I heavily made use of <code>f2py</code>—a tool that converts your FORTRAN code into a Python module. Here is how to use it:

\$ f2py -c -m module_name fortran_code.f90


In the current directory there is now a shared object of your module module_name.so. To use it in Python, just do the usual import. A docstring is generated for you as well, so you know how to call the methods, i.e., FORTRAN subroutines or functions

import module_name as mymodule
print mymodule.__doc__


Let's suppose you've written a subroutine that does some calculations that are faster in FORTRAN than in Python

! fortran_code.f90
subroutine some_calulations(a,b,c)
implicit none
real, intent(in) :: a, b
real, intent(out) :: c
real :: some_result

! start some calculations
...
c = some_result
end subroutine


To call that subroutine, our Python script has to look like this

import module_name as mymodule
print mymodule.__doc__

x1 = <some value>
x2 = <also some value>

result = mymodule.some_calculations(x1,x2)


That's it. Pretty easy, huh? Now you can use your fast code and let Python do the rest, like making a proper data analysis using numpy or scipy, or doing the graphics using matplotlib.

 Tags: f2py, fortran, python, tools 2013/08/02