diff --git a/cma.py b/cma.py
index 01ddbfa..f1a5f4b 100755
--- a/cma.py
+++ b/cma.py
@@ -4686,6 +4686,7 @@ class _CMAParameters(object):
     def disp(self):
         pprint(self.__dict__)
 
+
 def fmin(objective_function, x0, sigma0,
          options=None,
          args=(),
@@ -4696,7 +4697,8 @@ def fmin(objective_function, x0, sigma0,
          eval_initial_x=False,
          noise_handler=None,
          noise_change_sigma_exponent=1,
-         noise_kappa_exponent=0  # TODO: add max kappa value as parameter
+         noise_kappa_exponent=0,  # TODO: add max kappa value as parameter
+         bipop=False
         ):
     """functional interface to the stochastic optimizer CMA-ES
     for non-convex function minimization.
@@ -4729,9 +4731,12 @@ def fmin(objective_function, x0, sigma0,
         `x0`
             list or `numpy.ndarray`, initial guess of minimum solution
             before the application of the geno-phenotype transformation
-            according to the ``transformation`` option. Otherwise
-            `x0` can also be a `cma.CMAEvolutionStrategy` object instance.
-            In the latter case `sigma0` can be ``None``.
+            according to the ``transformation`` option.  It can also be
+            a string that holds a Python expression that is evaluated
+            to yield the initial guess - this is important in case
+            restarts are performed so that they start from different
+            places.  Otherwise `x0` can also be a `cma.CMAEvolutionStrategy`
+            object instance - in that case `sigma0` can be ``None``.
         `sigma0`
             scalar, initial standard deviation in each coordinate.
             `sigma0` should be about 1/4th of the search domain width
@@ -4750,7 +4755,10 @@ def fmin(objective_function, x0, sigma0,
             `gradf` is called once in each iteration if
             ``gradf is not None``.
         ``restarts=0``
-            number of restarts
+            number of restarts; the restarts are performed with
+            increasing population size, with restarts=9 this function
+            is equivalent to the well-known IPOP-CMA-ES restart strategy;
+            to restart from different points, pass x0 as a string.
         ``restart_from_best=False``
             which point to restart from
         ``incpopsize=2``
@@ -4765,6 +4773,20 @@ def fmin(objective_function, x0, sigma0,
         ``noise_evaluations_as_kappa``
             instead of applying reevaluations, the "number of evaluations"
             is (ab)used as scaling factor kappa (experimental).
+        ``bipop``
+            if True, run as BIPOP-CMA-ES; BIPOP is a special restart
+            strategy switching between two population sizings - small
+            (like the default CMA, but with more focused search) and
+            large (progressively increased as in IPOP). This makes the
+            algorithm perform well both on functions with local optima
+            in irregular, wide basins of attraction (by frequently
+            restarting with small populations) and narrow local optima
+            following a general global structure (well suited for large
+            populations).  For the bipop parameter to actually
+            take effect, also select non-zero number of restarts;
+            the recommended setting is ``bipop=True, restarts=9``
+            and x0 passed as a string.  Note that small-population
+            restarts do not count into the total restart count.
 
     Optional Arguments
     ==================
@@ -4858,6 +4880,40 @@ def fmin(objective_function, x0, sigma0,
         es = cma.CMAEvolutionStrategy([0.1] * 10, 0.5,
                     options=options).optimize(cma.fcts.rosen)
 
+    The following example calls `fmin` optimizing the Rosenbrock function
+    in 3-D with random initial solution in [-2,2], initial step-size 0.5
+    and the BIPOP restart strategy (that progressively increases population).
+    The options are specified for the usage with the `doctest` module.
+
+    >>> import cma
+    >>> # cma.CMAOptions()  # returns all possible options
+    >>> options = {'seed':12345, 'verb_time':0, 'ftarget': 1e-8}
+    >>>
+    >>> res = cma.fmin(cma.fcts.rastrigin, '2. * np.random.rand(3) - 1', 0.5,
+    ...                options, restarts=9, bipop=True)
+    (3_w,7)-aCMA-ES (mu_w=2.3,w_1=58%) in dimension 3 (seed=12345)
+    Iterat #Fevals   function value    axis ratio  sigma  minstd maxstd min:sec
+        1       7 1.633489455763566e+01 1.0e+00 4.35e-01  4e-01  4e-01
+        2      14 9.762462950258016e+00 1.2e+00 4.12e-01  4e-01  4e-01
+        3      21 2.461107851413725e+01 1.4e+00 3.78e-01  3e-01  4e-01
+      100     700 9.949590571272680e-01 1.7e+00 5.07e-05  3e-07  5e-07
+      123     861 9.949590570932969e-01 1.3e+00 3.93e-06  9e-09  1e-08
+    termination on tolfun=1e-11
+    final/bestever f-value = 9.949591e-01 9.949591e-01
+    mean solution: [  9.94958638e-01  -7.19265205e-10   2.09294450e-10]
+    std deviation: [  8.71497860e-09   8.58994807e-09   9.85585654e-09]
+    ...
+    (4_w,9)-aCMA-ES (mu_w=2.8,w_1=49%) in dimension 3 (seed=12349)
+    Iterat #Fevals   function value    axis ratio  sigma  minstd maxstd min:sec
+        1  5342.0 2.114883315350800e+01 1.0e+00 3.42e-02  3e-02  4e-02
+        2  5351.0 1.810102940125502e+01 1.4e+00 3.79e-02  3e-02  4e-02
+        3  5360.0 1.340222457448063e+01 1.4e+00 4.58e-02  4e-02  6e-02
+       50  5783.0 8.631491965616078e-09 1.6e+00 2.01e-04  8e-06  1e-05
+    termination on ftarget=1e-08 after 4 restarts
+    final/bestever f-value = 8.316963e-09 8.316963e-09
+    mean solution: [ -3.10652459e-06   2.77935436e-06  -4.95444519e-06]
+    std deviation: [  1.02825265e-05   8.08348144e-06   8.47256408e-06]
+
     In either case, the method::
 
         cma.plot();
@@ -4905,9 +4961,57 @@ def fmin(objective_function, x0, sigma0,
         # checked that no options.ftarget =
         opts = CMAOptions(options.copy()).complement()
 
+        # BIPOP-related variables:
+        runs_with_small = 0
+        small_i = [0]
+        large_i = [0]
+        popsize0 = None  # to be evaluated after the first iteration
+        maxiter0 = None  # to be evaluated after the first iteration
+        base_evals = 0
+
         irun = 0
         best = BestSolution()
         while True:  # restart loop
+            sigma_factor = 1
+
+            # Adjust the population according to BIPOP after a restart.
+            if not bipop:
+                # BIPOP not in use, simply double the previous population
+                # on restart.
+                if irun > 0:
+                    popsize_multiplier = fmin_options['incpopsize'] ** (irun - runs_with_small)
+                    opts['popsize'] = popsize0 * popsize_multiplier
+
+            elif irun == 0:
+                # Initial run is with "normal" population size; it is
+                # the large population before first doubling, but its
+                # budget accounting is the same as in case of small
+                # population.
+                poptype = 'small'
+
+            elif sum(small_i) < sum(large_i):
+                # An interweaved run with small population size
+                poptype = 'small'
+
+                restarts += 1  # A small restart doesn't count in the total
+                runs_with_small += 1  # _Before_ it's used in popsize_lastlarge
+
+                sigma_factor = 0.01 ** np.random.uniform()  # Local search
+                popsize_multiplier = fmin_options['incpopsize'] ** (irun - runs_with_small)
+                opts['popsize'] = np.floor(popsize0 * popsize_multiplier ** (np.random.uniform() ** 2))
+                opts['maxiter'] = min(maxiter0, 0.5 * sum(large_i) / opts['popsize'])
+                # print('small basemul %s --> %s; maxiter %s' % (popsize_multiplier, opts['popsize'], opts['maxiter']))
+
+            else:
+                # A run with large population size; the population
+                # doubling is implicit with incpopsize.
+                poptype = 'large'
+
+                popsize_multiplier = fmin_options['incpopsize'] ** (irun - runs_with_small)
+                opts['popsize'] = popsize0 * popsize_multiplier
+                opts['maxiter'] = maxiter0
+                # print('large basemul %s --> %s; maxiter %s' % (popsize_multiplier, opts['popsize'], opts['maxiter']))
+
             # recover from a CMA object
             if irun == 0 and isinstance(x0, CMAEvolutionStrategy):
                 es = x0
@@ -4922,9 +5026,9 @@ def fmin(objective_function, x0, sigma0,
             else:  # default case
                 if irun and eval(str(fmin_options['restart_from_best'])):
                     print('CAVE: restart_from_best is often not useful')
-                    es = CMAEvolutionStrategy(best.x, sigma0, opts)
+                    es = CMAEvolutionStrategy(best.x, sigma_factor * sigma0, opts)
                 else:
-                    es = CMAEvolutionStrategy(x0, sigma0, opts)
+                    es = CMAEvolutionStrategy(x0, sigma_factor * sigma0, opts)
                 if eval_initial_x:
                     x = es.gp.pheno(es.mean,
                                     into_bounds=es.boundary_handler.repair,
@@ -4998,17 +5102,33 @@ def fmin(objective_function, x0, sigma0,
             best.update(es.best, es.sent_solutions)  # in restarted case
             # es.best.update(best)
 
+            this_evals = es.countevals - base_evals
+            base_evals = es.countevals
+
+            # BIPOP stats update
+
+            if irun == 0:
+                popsize0 = opts['popsize']
+                maxiter0 = opts['maxiter']
+                # XXX: This might be a bug? Reproduced from Matlab
+                small_i.append(this_evals)
+
+            if bipop:
+                if poptype == 'small':
+                    small_i.append(this_evals)
+                else:  # poptype == 'large'
+                    large_i.append(this_evals)
+
             # final message
             if opts['verb_disp']:
                 es.result_pretty(irun, time.asctime(time.localtime()),
                                  best.f)
 
             irun += 1
-            if irun > fmin_opts['restarts'] or 'ftarget' in es.stop() \
+            if irun > restarts or 'ftarget' in es.stop() \
                     or 'maxfevals' in es.stop(check=False):
                 break
             opts['verb_append'] = es.countevals
-            opts['popsize'] = fmin_opts['incpopsize'] * es.sp.popsize  # TODO: use rather options?
             opts['seed'] += 1
 
         # while irun
@@ -5035,6 +5155,7 @@ def fmin(objective_function, x0, sigma0,
             print(' in/outcomment ``raise`` in last line of cma.fmin to prevent/restore KeyboardInterrupt exception')
         raise  # cave: swallowing this exception can silently mess up experiments, if ctrl-C is hit
 
+
 # _____________________________________________________________________
 # _____________________________________________________________________
 #
