new_contents: ######################################################################
  # Copyright (C) 2015 Jaakko Luttinen
  #
  # This file is licensed under Version 3.0 of the GNU General Public
  # License. See LICENSE for a text of the license.
  ######################################################################
  
  ######################################################################
  # This file is part of BayesPy.
  #
  # BayesPy is free software: you can redistribute it and/or modify it
  # under the terms of the GNU General Public License version 3 as
  # published by the Free Software Foundation.
  #
  # BayesPy is distributed in the hope that it will be useful, but
  # WITHOUT ANY WARRANTY; without even the implied warranty of
  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  # General Public License for more details.
  #
  # You should have received a copy of the GNU General Public License
  # along with BayesPy.  If not, see <http://www.gnu.org/licenses/>.
  ######################################################################
  
  """
  Demonstration of deterministic annealing.
  
  Deterministic annealing aims at avoiding convergence to local optima and
  finding the global optimum :cite:`Katahira:2008`.
  """
  
  
  import numpy as np
  import scipy
  
  import matplotlib.pyplot as plt
  import bayespy.plot as myplt
  
  from bayespy.utils import misc
  from bayespy.utils import random
  from bayespy.nodes import GaussianARD, Categorical, Mixture
  
  from bayespy.inference.vmp.vmp import VB
  from bayespy.inference.vmp import transformations
  
  import bayespy.plot as bpplt
  
  from bayespy.demos import pca
  
  
  def run(N=500, seed=42, maxiter=100, plot=True):
      """
      Run deterministic annealing demo for 1-D Gaussian mixture.
      """
  
      if seed is not None:
          np.random.seed(seed)
  
      mu = GaussianARD(0, 1,
                       plates=(2,),
                       name='means')
      Z = Categorical([0.3, 0.7],
                      plates=(N,),
                      name='classes')
      Y = Mixture(Z, GaussianARD, mu, 1,
                  name='observations')
  
      # Generate data
      z = Z.random()
      data = np.empty(N)
      for n in range(N):
          data[n] = [4, -4][z[n]]
  
      Y.observe(data)
  
      # Initialize means closer to the inferior local optimum
      mu.initialize_from_value([0, 6])
      #mu.initialize_from_value([6, 0])
  
      Q = VB(Y, Z, mu)
      Q.save()
  
      #
      # Standard VB-EM algorithm
      #
      Q.update(repeat=maxiter)
      ## if plot:
      ##     bpplt.pyplot.plot(Q.L, 'k-')
  
      mu_vbem = mu.u[0].copy()
      L_vbem = Q.compute_lowerbound()
  
      #
      # VB-EM with deterministic annealing
      #
      Q.load()
      beta = 0.03
      while True:
          Q.set_annealing(beta)
          print("Set annealing to %.2f" % beta)
          Q.update(repeat=maxiter, tol=1e-8)
          if beta == 1:
              break
          else:
              beta = min(beta*2, 1.0)
  
      ## if plot:
      ##     bpplt.pyplot.plot(Q.L, 'r:')
  
      ##     bpplt.pyplot.xlabel('Iterations')
      ##     bpplt.pyplot.ylabel('VB lower bound')
      ##     bpplt.pyplot.legend(['VB-EM', 'Deterministic annealing'],
      ##                         loc='lower right')
  
      mu_anneal = mu.u[0].copy()
      L_anneal = Q.compute_lowerbound()
  
      print("True component probabilities:", np.array([0.3, 0.7]))
      print("True component means:", np.array([4, -4]))
      print("VB-EM component means:", mu_vbem)
      print("VB-EM lower bound:", L_vbem)
      print("Annealed VB-EM component means:", mu_anneal)
      print("Annealed VB-EM lower bound:", L_anneal)
      
      return
  
  
  if __name__ == '__main__':
      import sys, getopt, os
      try:
          opts, args = getopt.getopt(sys.argv[1:],
                                     "",
                                     ["n=",
                                      "seed=",
                                      "maxiter="])
      except getopt.GetoptError:
          print('python annealing.py <options>')
          print('--n=<INT>        Number of data points')
          print('--maxiter=<INT>  Maximum number of VB iterations')
          print('--seed=<INT>     Seed (integer) for the random number generator')
          sys.exit(2)
  
      kwargs = {}
      for opt, arg in opts:
          if opt == "--maxiter":
              kwargs["maxiter"] = int(arg)
          elif opt == "--seed":
              kwargs["seed"] = int(arg)
          elif opt in ("--n",):
              kwargs["N"] = int(arg)
  
      run(**kwargs)
      plt.show()
  
  
subject: Add deterministic annealing demo