Cacher is the code snippet organizer for pro developers

from __future__ import division
import numpy as np
from scipy.optimize import brentq as root
from borg import *

nvars = 100 # decision variables: pollution for each timestep
nobjs = 4
nconstr = 1
nsamples = 100 # 100 scenarios of natural inflows
b = 0.42 # decay rate for P in lake (0.42 = irreversible)
q = 2.0 # recycling exponent
alpha = 0.4 # utility from pollution
delta = 0.98 # (1-r) discount rate
Pcrit = root(lambda x: x**q/(1+x**q) - b*x, 0.01, 1.5)

def lake(*decisions):
  objs = [0.0]*nobjs
  constr = [0.0]
  X = np.zeros((nvars,))
  average_daily_P = np.zeros((nvars,))
  decisions = np.array(decisions)

  for s in xrange(nsamples):
    X[0] = 0
    natural_inflows = np.random.lognormal(-3.52, 0.105, size=nvars)
    
    for t in xrange(1,nvars):
      X[t] = (1-b)*X[t-1] + X[t-1]**q/(1+X[t-1]**q) + decisions[t-1] + natural_inflows[t-1]
      average_daily_P[t] += X[t]/nsamples

    objs[3] -= np.sum(X < Pcrit)/(nsamples*nvars) # Maximize timesteps satisfying X < Pcrit
  
  objs[0] = np.max(average_daily_P) # Minimize the maximum daily P in lake
  objs[1] = -1*np.sum(alpha*decisions*np.power(delta,np.arange(nvars))) # Maximize the average sum of discounted benefits
  objs[2] = -1*np.sum(np.diff(decisions) > -0.02)/(nvars-1) # Maximize timesteps meeting inertia threshold

  constr[0] = Constraint.greaterThan(-1*objs[3], 0.95)
  return (objs, constr)

borg = Borg(nvars, nobjs, nconstr, lake)
borg.setBounds(*[[0.0, 0.1]]*nvars)
borg.setEpsilons(*[0.01]*nobjs)

result = borg.solve({"maxEvaluations":10000})

for solution in result:
  print solution.getObjectives() + solution.getConstraints()