/beam-3dveglab-vlab/src/main/resources/auxdata/VLabImpl.py

Large files files are truncated, but you can click here to view the full file

#
# Copyright (C) 2010-2014 Netcetera Switzerland (info@netcetera.com)
# 
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the Free
# Software Foundation; either version 3 of the License, or (at your option)
# any later version.
# This program 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 this program; if not, see http://www.gnu.org/licenses/
# 
# @(#) $Id: $
#
# Authors: Cyrill Schenkel, Daniel Kueckenbrink, Joshy Cyriac, Marcel Kessler, Jason Brazile
#

##############################################################################
# Three ways to run it:
# 
# 1. embedded within the BEAM processor (normal case)
#
# 2. standalone tests - headless (either jython or python)
#
#    jython -Dpython.path=jcommon-1.0.16.jar:jfreechart-1.0.13.jar VLabImpl.py
#    python VLabImpl.py
#
# 3. standalone with a 'fake' swing-based gui
# 
#    java -jar ${HOME}/beam-4.11/lib/jython-2.5.2.jar -Dvlab.fakebeam=1 -Dpython.path=${HOME}/beam-4.11/lib/jcommon-1.0.16.jar:${HOME}/beam-4.11/lib/jfreechart-1.0.13.jar VLabImpl.py
#
# Those BEAM-supplied jars can also be obtained like this:
#    wget -U "Mozilla/5.0" http://repo1.maven.org/maven2/jfree/jfreechart/1.0.13/jfreechart-1.0.13.jar
#    wget -U "Mozilla/5.0" http://repo1.maven.org/maven2/jfree/jcommon/1.0.16/jcommon-1.0.16.jar
#
#

import sys, math, operator, array, time, struct
from array import array

class VLAB:
  """VLAB contains conf. constants, static utility methods, and test methods"""

  COPYRIGHT_INFO     = 'Copyright (C) 2010-2014 Netcetera Switzerland (info@netcetera.com)'
  PROCESSOR_NAME     = 'BEAM VLab Processor'
  PROCESSOR_SNAME    = 'beam-vlab'
  REQUEST_TYPE       = 'VLAB'
  UI_TITLE           = 'VLab - Processor'
  VERSION_STRING     = '1.0 (20 July 2014)'
  DEFAULT_LOG_PREFIX = 'vlab'
  LOGGER_NAME        = 'beam.processor.vlab'

  D_PRODNAME         = 'vlab_out.dim'
  P_CONDITION        = '.condition'
  P_EXPRESSION       = '.expression'
  P_OUTPUT           = '.output'

  # NOTE: Once released, random number generation should NOT be reproducible
  CONF_RND_REPRODUCE = True

  JCB                = 'JComboBox'
  JTF                = 'JTextField'
  JTD                = (JTF, False)

  K_LIBRAT           = 'librat'
  K_DART             = 'DART'
  K_DUMMY            = 'dummy'

  K_LAEGERN          = 'Laegern'
  K_THARANDT         = 'Tharandt'
  K_RAMI             = 'RAMI'
  K_USER_DEFINED     = 'User defined (UserDefined.obj)'

  K_YES              = 'Yes'
  K_NO               = 'No'

  K_SENTINEL2        = 'MSI (Sentinel 2)'
  K_SENTINEL3        = 'OLCI (Sentinel 3)'
  K_MODIS            = 'MODIS'
  K_MERIS            = 'MERIS'
  K_LANDSAT_OLI      = 'Landsat (OLI)'
  K_LANDSAT_ETM      = 'Landsat (ETM)'

  K_RURAL            = 'Rural'
  K_MARITIME         = 'Maritime'
  K_URBAN            = 'Urban'
  K_TROPOSPHERIC     = 'Tropospheric'

  DBG                = 'debug' 
  INF                = 'info'
  ERR                = 'error'

  plst               = []

  if sys.platform.startswith('java'):
    from java.util.logging import Logger
    from java.util.logging import FileHandler
    from java.util.logging import Level
    logger = Logger.getLogger(LOGGER_NAME)
    # comment these out
    # logfh = FileHandler('%s.log' % LOGGER_NAME)
    # logfh.setLevel(Level.FINEST)
    # logger.addHandler(logfh)
  else:
    import logging
    logger = logging.getLogger(LOGGER_NAME)
    logger.setLevel(logging.DEBUG)
    logch = logging.StreamHandler()
    logch.setLevel(logging.DEBUG)
    logger.addHandler(logch)
    # comment these out
    # logfh = logging.FileHandler('%s.log' % LOGGER_NAME)
    # logfh.setLevel(logging.DEBUG)
    # logger.addHander(logfh)

  model = (
{'Forward Modeling': (
{'Model Selection': (
 ('3D Scene',          '3dScene',             JCB, (K_RAMI, K_LAEGERN, K_THARANDT, K_USER_DEFINED)),
 ('RT Processor',      'RTProcessor',         JCB, (K_DUMMY, K_LIBRAT, K_DART)))},
{'Spectral Characteristics': (
 ('Sensor',            'Sensor',              JCB, (K_SENTINEL2, K_SENTINEL3, K_MODIS, K_MERIS, K_LANDSAT_OLI, K_LANDSAT_ETM)),
 ('Bands',             'Bands',               JTF, '1, 2, 3, 4, 5, 6, 7, 8, 9, 10'))},
{'Viewing Characteristics': (
 ('Zenith',            'ViewingZenith',       JTF, '20.0'),
 ('Azimuth',           'ViewingAzimuth',      JTF, '0.0'))},
{'Illumination Characteristics':(
 ('Zenith',            'IlluminationZenith',  JTF, '20.0'),
 ('Azimuth',           'IlluminationAzimuth', JTF, '0.0'))},
{'Scene Parameters': (
 ('Pixel Size',        'ScenePixel',          JTD, 'not available yet'),
 ('(Alt A) Filename',  'SceneLocFile',        JTD, 'not available yet'),
 ('(Alt B) XC',        'SceneXC',             JTD, 'not available yet'),
 ('(Alt B) XW',        'SceneXW',             JTD, 'not available yet'),
 ('(Alt B) YC',        'SceneYC',             JTD, 'not available yet'),
 ('(Alt B) YW',        'SceneYW',             JTD, 'not available yet'))},
{'Atmospheric Parameters': (
 ('CO2 Mixing Ratio',  'AtmosphereCO2',       JTF, '1.6'),
 ('Aerosol Profile',   'AtmosphereAerosol',   JCB, (K_RURAL, K_MARITIME, K_URBAN, K_TROPOSPHERIC)),
 ('Water Vapor',       'AtmosphereWater',     JTF, '0.0'),
 ('Ozone Column',      'AtmosphereOzone',     JTF, '300'))},
{'Output Parameters': (
 ('Result file prefix','OutputPrefix',        JTF, 'HET01_DIS_UNI_NIR_20.obj'),
 ('Result Directory',  'OutputDirectory',     JTF, 'dart.rpv.rami.2'),
 ('Image file',        'ImageFile',           JCB, (K_YES, K_NO)),
 ('Ascii file',        'AsciiFile',           JCB, (K_YES, K_NO)))}
)},
{'DHP Simulation': (
{'Model Selection': (
 ('3D Scene',          'DHP_3dScene',         JCB, (K_RAMI, K_LAEGERN, K_THARANDT)),
 (),
 ('Resolution',        'DHP_Resolution',      JTF, '4000000'),
 ())},
{'DHP Location': (
 ('X',                 'DHP_X',               JTD, 'not available yet'),
 ('Y',                 'DHP_Y',               JTD, 'not available yet'))},
{'DHP Properties': (
 ('Zenith',            'DHP_Zenith',          JTF, '0.0'),
 ('Azimuth',           'DHP_Azimuth',         JTF, '0.0'))},
{'DHP Imaging Plane': (
 ('Orientation',       'DHP_Orientation',     JTD, 'not available yet'),
 ('Height(z)',         'DHP_Height',          JTD, 'not available yet'))},
{'Output Parameters': (
 ('Result file prefix','DHP_OutputPrefix',    JTF, 'RAMI00_'),
 ('Result Directory',  'DHP_OutputDirectory', JTF, ''),
 ('Image file',        'DHP_ImageFile',       JCB, (K_YES, K_NO)),
 ('Ascii file',        'DHP_AsciiFile',       JCB, (K_YES, K_NO)))}
)},
)
  # set parameter names
  for tabGroups in model:
    for tabName in tabGroups:
      for group in tabGroups[tabName]:
        for groupName in group:
          for groupTuple in group[groupName]:
            if len(groupTuple) == 4:
              (lbl, nm, typ, vals) = groupTuple
              exec('P_' + nm + ' = nm')
              exec('L_' + nm + ' = lbl')
              plst.append(nm)

  def __init__(self):
    self.cmap = {}
    self.vmap = {}

  def me():
    """Returns name of currently executing method"""
    nm = ''
    try:
      raise ZeroDivisionError
    except ZeroDivisionError:
      nm = sys.exc_info()[2].tb_frame.f_back.f_code.co_name
    return nm+'()'
  me = staticmethod(me)
  def lineSeparator():
    """Return the OS line separator"""
    if sys.platform.startswith('java'):
      from java.lang import System
      if sys.platform.startswith('java1.6'):
        return System.getProperty('line.separator')
      elif sys.platform.startswith('java1.7'):
        return System.lineSeparator()
    else:
      import os
      os.linesep
  lineSeparator = staticmethod(lineSeparator)
  def listdir(path):
    """list files in the directory given by path"""
    if sys.platform.startswith('java'):
      from java.io import File
      array = File(path).list()
      listFile = []
      if array != None:
        for i in xrange(len(array)):
          listFile.append(array[i])
      return listFile
    else:
      import os
      return os.listdir(path)
  listdir = staticmethod(listdir)
  def getenv(key, default=None):
    if sys.platform.startswith('java'):
      from java.lang.System import getenv
      return getenv(key)
    else:
      import os
      return os.getenv(key)
  getenv = staticmethod(getenv)
  def checkFile(fname):
    """open a file if it exists, otherwise die"""
    try:
      fp = open(fname, 'r+')
      return fp
    except IOError, e:
      raise RuntimeError(e)
  checkFile = staticmethod(checkFile)
  def fileExists(fname):
    """check if fname exists as a file"""
    if sys.platform.startswith('java'):
      from java.io import File
      return File(fname).exists()
    else:
      import os
      return os.path.exists(fname)
  fileExists = staticmethod(fileExists)
  def getFullPath(fname):
    """return canonical/absolute path of given fname"""
    if sys.platform.startswith('java'):
      from java.io import File
      return File(fname).getCanonicalPath()
    else:
      import os
      return os.path.abspath(fname)
  getFullPath = staticmethod(getFullPath)
  def copyDir(dname, target):
    """Recursively copy 'dname' directory to 'target'.
       /!\If 'target' already exists it will be removed or overwrited.
    """
    if sys.platform.startswith('java'):
      # import java module
      from java.io import File
      from java.io import FileInputStream
      from java.io import FileOutputStream
      from java.util import Scanner
      from java.lang import String
      dnameFile = File(dname)
      targetFile = File(target)
      # recursive copy
      if dnameFile.isDirectory():
        # Create folder if not exists
        if not targetFile.exists():
          targetFile.mkdir()
        # Copy all content recursively
        for fname in dnameFile.list().tolist():
          VLAB.copyDir(dname + File.separator + fname, target + File.separator + fname)
      else:
        # Read dname file
        istream = FileInputStream(dname)
        scan = Scanner(istream).useDelimiter("\\Z")
        # Test if file is empty
        if scan.hasNextLine():
          content = String(scan.next())
        else:
          content = String("")
        scan.close()
        # Create and write target
        if not targetFile.exists():
          targetFile.createNewFile()
        ostream = FileOutputStream(target)
        ostream.write(content.getBytes())
        ostream.flush()
        ostream.close()
    else:
      import shutil, os
      # remove exisiting target
      if os.path.isdir(target) or os.path.isfile(target):
        shutil.rmtree(target)
      # recursive copy of dnma to target
      shutil.copytree(dname, target)
  copyDir = staticmethod(copyDir)
  def XMLEditNode(fname, nodeName, attributName, value, multiple=False):
    """ Edit a given node (nodeName) in a given XML file (fname)
    attributName and value could be either a list of string or a string
    """
    if sys.platform.startswith('java'):
      from javax.xml.parsers import DocumentBuilderFactory
      from javax.xml.transform import TransformerFactory
      from javax.xml.transform import OutputKeys
      from javax.xml.transform.dom import DOMSource
      from javax.xml.transform.stream import StreamResult
      from java.io import File
      # Get whole tree
      tree = DocumentBuilderFactory.newInstance().newDocumentBuilder().parse(fname)
      # Get nodeName
      node = tree.getElementsByTagName(nodeName)
      # Check if we get only one node (as expected)
      if node.getLength() == 0:
        raise IOError("Cannot found '%s' in file '%s'" % (nodeName, fname))
      elif node.getLength() == 1:
        nodes = [node.item(0)]
      else:
        if not multiple:
          raise IOError("Get multiple nodes for '%s' in file '%s'" % (nodeName, fname))
        else:
          nodes = [ node.item(i) for i in xrange(node.getLength()) ]
      for node in nodes:
        # Modify the node attribute
        if isinstance(attributName, list) and isinstance(value, list):
          for att, val in zip(attributName, value):
            node.setAttribute(att, val)
        elif isinstance(attributName, str) and isinstance(value, str):
          node.setAttribute(attributName, value)
        else:
          raise ValueError("Wrong parameter used: attributName and value should be both either a list of string or a string")
      # Write new XML tree in fname
      transformer = TransformerFactory.newInstance().newTransformer()
      transformer.setOutputProperty(OutputKeys.INDENT, "yes")
      source = DOMSource(tree)
      result = StreamResult(File(fname))
      transformer.transform(source, result)
    else:
      import xml.etree.ElementTree as ET
      # Get whole tree from xml
      tree = ET.parse(fname)
      # Get nodeName
      #nodes = tree.findall(".//*%s" % nodeName) # This line seems to not work for root child node!! Bug?
      nodes = tree.findall(".//*../%s" % nodeName)
      # Check if we get only one node (as expected)
      if len(nodes) == 0:
        raise IOError("Cannot found '%s' in file '%s'" % (nodeName, fname))
      elif len(nodes) > 1 and not multiple:
        raise IOError("Get multiple nodes for '%s' in file '%s'" % (nodeName, fname))
      for node in nodes:
        # Modify the node attribute
        if isinstance(attributName, list) and isinstance(value, list):
          for att, val in zip(attributName, value):
            node.set(att, val)
        elif isinstance(attributName, str) and isinstance(value, str):
          node.set(attributName, value)
        else:
          raise ValueError("Wrong parameter used: attributName and value should be both either a list of string or a string")
      # Write new XML tree in fname
      tree.write(fname)
  XMLEditNode = staticmethod(XMLEditNode)
  def XMLReplaceNodeContent(fname, parent, subnode, attributName, value, spectralBands=False):
    """ Edit an XML file (fname) and replace the content of a node with subnode(s) (subnode) within attribute(s) and value(s).
    attributName and value could be either a list of string or a string
    """
    if sys.platform.startswith('java'):
      from javax.xml.parsers import DocumentBuilderFactory
      from javax.xml.transform import TransformerFactory
      from javax.xml.transform import OutputKeys
      from javax.xml.transform.dom import DOMSource
      from javax.xml.transform.stream import StreamResult
      from java.io import File
      # Get whole tree
      tree = DocumentBuilderFactory.newInstance().newDocumentBuilder().parse(fname)
      # Get nodeName
      node = tree.getElementsByTagName(parent)
      # Check if we get only one node (as expected)
      if node.getLength() > 1:
        raise IOError("Get multiple nodes for '%s' in file '%s'" % (parent, fname))
      elif node.getLength() == 0:
        raise IOError("Cannot found '%s' in file '%s'" % (parent, fname))
      else:
        node = node.item(0)
      # Remove content node
      while node.hasChildNodes():
        node.removeChild(node.getFirstChild())
      # Modify the node attribute
      elem = tree.createElement(subnode)
      if spectralBands:
        elem.setAttribute("bandNumber", "0")
        elem.setAttribute("spectralDartMode", "0")
      if isinstance(attributName, list) and isinstance(value, list):
        if isinstance(value[0], list):
          for bandNumber, val in enumerate(value):
            if spectralBands:
              elem.setAttribute("bandNumber", str(bandNumber))
            for atr, v in zip(attributName, val):
              elem.setAttribute(atr, v)
            node.appendChild(elem.cloneNode(True))
        else:
          elem = tree.createElement(subnode)
          for atr, v in zip(attributName, value):
            elem.setAttribute(atr, v)
          node.appendChild(elem)
      elif isinstance(attributName, str) and isinstance(value, str):
        elem = tree.createElement(subnode)
        elem.setAttribute(attributName, value)
        node.appendChild(elem)
      else:
        raise ValueError("Wrong parameter used: attributName and value should be both either a list of string or a string")
      # Write new XML tree in fname
      transformer = TransformerFactory.newInstance().newTransformer()
      transformer.setOutputProperty(OutputKeys.INDENT, "yes")
      source = DOMSource(tree)
      result = StreamResult(File(fname))
      transformer.transform(source, result)
    else:
      import xml.etree.ElementTree as ET
      # Get whole tree from xml
      tree = ET.parse(fname)
      # Get nodeName
      #nodes = tree.findall(".//*%s" % nodeName) # This line seems to not work for root child node!! Bug?
      node = tree.findall(".//*../%s" % parent)
      # Check if we get only one node (as expected)
      if len(node) > 1:
        raise IOError("Get multiple nodes for '%s' in file '%s'" % (parent, fname))
      elif len(node) == 0:
        raise IOError("Cannot found '%s' in file '%s'" % (parent, fname))
      else:
        node = node[0]
      # Remove content of node
      node.clear()
      # Modify the node attribute
      if spectralBands:
        attrib = {"bandNumber":"0", "spectralDartMode":"0"}
      else:
        attrib = {}
      if isinstance(attributName, list) and isinstance(value, list):
        if isinstance(value[0], list):
          for bandNumber, val in enumerate(value):
            if spectralBands:
              attrib["bandNumber"] = str(bandNumber)
            for atr, v in zip(attributName, val):
              attrib[atr] = v
            node.append(ET.Element(subnode, attrib=attrib))
        else:
          for atr, v in zip(attributName, value):
            attrib[atr] = v
          node.append(ET.Element(subnode, attrib=attrib))
      elif isinstance(attributName, str) and isinstance(value, str):
        attrib[attributName] = value
        node.append(ET.Element(subnode, attrib=attrib))
      else:
        raise ValueError("Wrong parameter used: attributName and value should be both either a list of string or a string")
      # Write new XML tree in fname
      tree.write(fname)
  XMLReplaceNodeContent = staticmethod(XMLReplaceNodeContent)
  def XMLAddNode(fname, parent, treeNodes, nodesSetup):
    """ Add a node (subnode) to a given parent in the provided fname file.
    """
    if sys.platform.startswith('java'):
      from javax.xml.parsers import DocumentBuilderFactory
      from javax.xml.transform import TransformerFactory
      from javax.xml.transform import OutputKeys
      from javax.xml.transform.dom import DOMSource
      from javax.xml.transform.stream import StreamResult
      from java.io import File
      # Get whole tree
      tree = DocumentBuilderFactory.newInstance().newDocumentBuilder().parse(fname)
      # Get nodeName
      node = tree.getElementsByTagName(parent)
      # Check if we get only one node (as expected)
      if node.getLength() > 1:
        raise IOError("Get multiple nodes for '%s' in file '%s'" % (parent, fname))
      elif node.getLength() == 0:
        raise IOError("Cannot found '%s' in file '%s'" % (parent, fname))
      else:
        node = node.item(0)
      # Create all the new nodes
      nodes = {parent:node}
      for name in treeNodes:
        nodes[name] = (tree.createElement(name))
        for atr, val in zip(*[ nodesSetup[name][key] for key in nodesSetup[name].iterkeys() if key in ["attribute", "value"] ]):
          nodes[name].setAttribute(atr, val)
      # Add all created node to its parent (previous node in the list)
      for name, elem in nodes.iteritems():
        if name != parent:
          nodes[nodesSetup[name]["parent"]].appendChild(elem)
      # Write new XML tree in fname
      transformer = TransformerFactory.newInstance().newTransformer()
      transformer.setOutputProperty(OutputKeys.INDENT, "yes")
      source = DOMSource(tree)
      result = StreamResult(File(fname))
      transformer.transform(source, result)
    else:
      import xml.etree.ElementTree as ET
      # Get whole tree from xml
      tree = ET.parse(fname)
      # Get nodeName
      #nodes = tree.findall(".//*%s" % nodeName) # This line seems to not work for root child node!! Bug?
      node = tree.findall(".//*../%s" % parent)
      # Check if we get only one node (as expected)
      if len(node) > 1:
        raise IOError("Get multiple nodes for '%s' in file '%s'" % (parent, fname))
      elif len(node) == 0:
        raise IOError("Cannot found '%s' in file '%s'" % (parent, fname))
      else:
        node = node[0]
      # Create all the new nodes
      nodes = {parent:node}
      for name in treeNodes:
        attrib = {}
        for atr, val in zip(*[ nodesSetup[name][key] for key in nodesSetup[name].iterkeys() if key in ["attribute", "value"] ]):
          attrib[atr] = val
        nodes[name] = ET.Element(name, attrib)
      # Add all created node to its parent (previous node in the list)
      for name, elem in nodes.iteritems():
        if name != parent:
          nodes[nodesSetup[name]["parent"]].append(elem)
      # Write new XML tree in fname
      tree.write(fname)
  XMLAddNode = staticmethod(XMLAddNode)
  def XMLGetNodeAttributeValue(fname, nodeName, attributName):
    """ Return the value of the given attributName for the given nodeName
    """
    if sys.platform.startswith('java'):
      from javax.xml.parsers import DocumentBuilderFactory
      from org.w3c.dom import Element
      # Get whole tree
      tree = DocumentBuilderFactory.newInstance().newDocumentBuilder().parse(fname)
      # Get nodeName
      node = tree.getElementsByTagName(nodeName)
      # Check if we get only one node (as expected)
      if node.getLength() > 1:
        raise IOError("Get multiple nodes for '%s' in file '%s'" % (nodeName, fname))
      elif node.getLength() == 0:
        raise IOError("Cannot found '%s' in file '%s'" % (nodeName, fname))
      else:
        node = node.item(0)
      # Check if that node has attributName
      if node.hasAttribute(attributName):
        return node.getAttributes().getNamedItem(attributName).getNodeValue()
      else:
        raise IOError("Attribute name '%s' not found for node '%s' in file '%s'" % (attributName, nodeName, fname))
    else:
      import xml.etree.ElementTree as ET
      # Get whole tree from xml
      tree = ET.parse(fname)
      # Get nodeName
      #nodes = tree.findall(".//*%s" % nodeName) # This line seems to not work for root child node!! Bug?
      node = tree.findall(".//*../%s[@%s]" % (nodeName, attributName))
      # Check if we get only one node (as expected)
      if len(node) > 1:
        raise IOError("Get multiple nodes for '%s' in file '%s'" % (nodeName, fname))
      elif len(node) == 0:
        raise IOError("Cannot found '%s' in file '%s'" % (nodeName, fname))
      else:
        node = node[0]
      # Check if that node has an attributName
      if attributName in node.keys():
        return node.get(attributName)
      else:
        raise IOError("Attribute name '%s' not found for node '%s' in file '%s'" % (attributName, nodeName, fname))
  XMLGetNodeAttributeValue = staticmethod(XMLGetNodeAttributeValue)
  def getBandsFromGUI(bands):
    """ Return a DART spectral bands list: ["deltaLambda", "meanLambda"] in micro meter
    In case of several bands the result should be a list of list:
    [["deltaLambda0", "meanLambda0"], ["deltaLambda1", "meanLambda1"], ["deltaLambda2", "meanLambda2"], ...]
    e.g.: [["0.02", "0.56"], ["0.02", "0.58"], ["0.02", "0.60"], ["0.02", "0.62"]]
    """
    # TODO: You should write the spectral band converter here!
    return [["0.02", "0.56"], ["0.02", "0.58"], ["0.02", "0.60"], ["0.02", "0.62"]]
  getBandsFromGUI = staticmethod(getBandsFromGUI)
  class path:
    def exists(path):
      if sys.platform.startswith('java'):
        from java.io import File
        return File(path).exists()
      else:
        import os
        return os.path.exists(path)
    exists = staticmethod(exists)
    def isabs(path):
      if sys.platform.startswith('java'):
        from java.io import File
        return File(path).isAbsolute()
      else:
        import os
        return os.path.isabs(path)
    isabs = staticmethod(isabs)
    def isdir(path):
      if sys.platform.startswith('java'):
        from java.io import File
        return File(path).isDirectory()
      else:
        import os
        return os.path.isdir(path)
    isdir = staticmethod(isdir)
    def join(path, *args):
      if sys.platform.startswith('java'):
        from java.io import File
        f = File(path)
        for a in args:
          g = File(a)
          if g.isAbsolute() or len(f.getPath()) == 0:
            f = g
          else:
            f = File(f, a)
        return f.getPath()
      else:
        import os
        return os.path.join(path, *args)
    join = staticmethod(join)
    def isfile(path):
      if sys.platform.startswith('java'):
        from java.io import File
        return File(path).isfile()
      else:
        import os
        return os.path.isfile(path)
    isfile = staticmethod(isfile)
    def normcase(path):
      if sys.platform.startswith('java'):
        from java.io import File
        return File(path).getPath()
      else:
        import os
        return os.path.normcase(path)
    normcase = staticmethod(normcase)
    def normpath(path):
      if sys.platform.startswith('java'):
        from java.io import File
        return File(path).getCanonicalPath()
      else:
        import os
        return os.path.normpath(path)
    normpath = staticmethod(normpath)
    def split(path):
      if sys.platform.startswith('java'):
        from java.io import File
        f=File(path)
        d=f.getParent()
        if not d:
          if f.isAbsolute():
            d = path
          else:
            d = ""
        return (d, f.getName())
      else:
        import os
        return os.path.split(path)
    split = staticmethod(split)
  def frange(end, start=0, inc=0):
    """a jython-compatible range function with float increments"""
    import math
    if not start:
      start = end + 0.0
      end = 0.0
    else: end += 0.0
    if not inc:
      inc = 1.0
    count = int(math.ceil((start - end) / inc))
    L = [None] * count
    L[0] = end
    for i in (xrange(1,count)):
      L[i] = L[i-1] + inc
    return L
  frange = staticmethod(frange)
  def openFileIfNotExists(filename):
    """open file exclusively"""
    if sys.platform.startswith('java'):
      from java.io import File
      File(filename).getParentFile().mkdirs()
      if File(filename).createNewFile():
        return open(filename, 'w')
      else:
        return None
    else:
      import os
      try:
        # os.O_EXCL => open exclusive => acquire a lock on the file
        fd = os.open(filename, os.O_CREAT|os.O_EXCL|os.O_WRONLY|os.O_TRUNC)
      except:
        return None
      fobj = os.fdopen(fd,'w')
      return fobj
  openFileIfNotExists = staticmethod(openFileIfNotExists)
  def rndInit(seed = None):
    """initialize the randon number generator"""
    if sys.platform.startswith('java'):
      from java.util import Random
      if seed == None:
        return Random()
      else:
        return Random(seed)
    else:
      import random
      random.seed(seed)
      return None
  rndInit = staticmethod(rndInit)
  def rndNextFloat(randState):
    """return the next pseudo-random floating point number in the sequence"""
    if sys.platform.startswith('java'):
      from java.util import Random
      return randState.nextFloat()
    else:
      import random
      return random.random()
  rndNextFloat = staticmethod(rndNextFloat)
  def r2d(v):
    """jython-compatible conversion of radians to degrees"""
    if sys.platform.startswith('java'):
      from java.lang import Math
      return Math.toDegrees(v)
    else:
      return math.degrees(v)
  r2d = staticmethod(r2d)
  def d2r(v):
    """jython-compatible conversion of degrees to radians"""
    if sys.platform.startswith('java'):
      from java.lang import Math
      return Math.toRadians(float(v))
    else:
      return math.radians(v)
  d2r = staticmethod(d2r)
  def mkDirPath(path):
    """create directory (including non-existing parents) for the given path"""
    if sys.platform.startswith('java'):
      from java.io import File
      if not File(path).isDirectory():
        if not File(path).mkdirs():
          raise RuntimeError('failed to mkdir %s' % path)
    else:
      import os
      try:
        os.stat(path)
      except:
        os.makedirs(path)
  mkDirPath = staticmethod(mkDirPath)
  def fPath(d,n):
    """get file path of the file defined by directory d and file name n"""
    if sys.platform.startswith('java'):
      from java.io import File
      return File(d, n).getPath()
    else:
      import os
      return os.path.join(d, n)
  fPath = staticmethod(fPath)
  def savetxt(a, b, fmt=False):
    """save text b in a text file named a"""
    fh = open(a, 'w')
    if not fmt:
      fmt = '%s'
    for row in b:
      for element in row:
        fh.write(fmt % element + ' ')
      fh.write('\n')
    fh.close()
  savetxt = staticmethod(savetxt)
  def osName():
    """return which OS we are currently running on"""
    if sys.platform.startswith('java'):
      from java.lang import System
      oname = System.getProperty('os.name')
    else:
      import os
      oname = os.name
      if not oname.endswith('x'): oname = 'Windows'
    return oname
  osName = staticmethod(osName)
  def expandEnv(instr):
    """potentially expand environment variables in the given string"""
    outstr = instr
    m = {'$HOME':'HOME','%HOMEDRIVE%':'HOMEDRIVE','%HOMEPATH%':'HOMEPATH'}
    for e in m:
      if outstr.find(e) != -1:
        if sys.platform.startswith('java'):
          from java.lang import System
          repl = System.getenv(m[e])
        else:
          import os
          repl = os.getenv(m[e])
        if repl != None:
          outstr = outstr.replace(e, repl)
    return outstr
  expandEnv = staticmethod(expandEnv)

  if sys.platform.startswith('java'):
    from java.lang import Runnable
    class Helper(Runnable):
      def __init__(self, nm, strm, fName):
        self.nm=nm; self.strm=strm; self.fp=None
        if fName != None:
          if VLAB.fileExists(fName):
            self.fp = open(fName, 'w')
          else:
            self.fp = VLAB.openFileIfNotExists(fName)
      def run(self):
        """helper class for slurping up child streams"""
        from java.io import BufferedReader
        from java.io import InputStreamReader
        line = None; br = BufferedReader(InputStreamReader(self.strm))
        line = br.readLine()
        while (line != None):
          if self.fp != None:
            self.fp.write(line + VLAB.lineSeparator())
            self.fp.flush()
          VLAB.logger.info('%s %s' %(self.nm, line.rstrip()))
          line = br.readLine()
        br.close()
        if self.fp != None:
          self.fp.close()
  else:
    def helper(nm, strm):
      """helper method for slurping up child streams"""
      for line in strm: VLAB.logger.info('%s %s' %(nm, line.rstrip()))
      if not strm.closed: strm.close()
    helper = staticmethod(helper)
  def doExec(cmdrec):
    """run the specified external program under windows or unix"""
    cmdLine = []
    osName = VLAB.osName()
    if osName.startswith('Windows'):
      cmd=cmdrec['windows']
      cmdLine = ['cmd', '/c']
    else:
      cmd=cmdrec['linux']
    exe = VLAB.expandEnv(cmd['exe'])
    if not (VLAB.fileExists(exe) or osName.startswith('Windows')):
      raise RuntimeError('Cannot find exe "%s"' % exe)
    cmdLine.append(exe)
    for i in cmd['cmdline']:
      cmdLine.append(VLAB.expandEnv(i))

    VLAB.logger.info('cmdLine is [%s]' % cmdLine)
    if sys.platform.startswith('java'):
      from java.lang import ProcessBuilder
      from java.lang import Thread
      from java.io import BufferedWriter
      from java.io import OutputStreamWriter
      from java.io import File
      pb = ProcessBuilder(cmdLine)
      if 'cwd' in cmd and cmd['cwd'] != None:
        pb.directory(File(VLAB.expandEnv(cmd['cwd'])))
      if 'env' in cmd and cmd['env'] != None:
        env = pb.environment()
        cmdenv = cmd['env']
        for e in cmdenv:
          env[e] = VLAB.expandEnv(cmdenv[e])
      proc = pb.start()
      stdoutfName = None
      if 'stdout' in cmd and cmd['stdout'] != None:
        stdoutfName = VLAB.expandEnv(cmd['stdout'])
      stderrfName = None
      if 'stderr' in cmd and cmd['stderr'] != None:
        stderrfName = VLAB.expandEnv(cmd['stderr'])
      outs = Thread(VLAB.Helper('out', proc.getInputStream(), stdoutfName))
      errs = Thread(VLAB.Helper('err', proc.getErrorStream(), stderrfName))
      outs.start(); errs.start()
      bw = BufferedWriter(OutputStreamWriter(proc.getOutputStream()))
      if 'stdin' in cmd and cmd['stdin'] != None:
        inFile = VLAB.expandEnv(cmd['stdin'])
        VLAB.logger.info('stdin is [%s]' % inFile)
        if 'cwd' in cmd and cmd['cwd'] != None:
          if not VLAB.fileExists(inFile):
            # try pre-pending the cwd
            inFile = VLAB.fPath(VLAB.expandEnv(cmd['cwd']), inFile)
            if not VLAB.fileExists(inFile):
              raise RuntimeError('Cannot find stdin "%s"' % cmd['cwd'])
        fp = open(inFile, 'r')
        for line in fp:
          bw.write(line)
        bw.close()
        fp.close()
      exitCode = proc.waitFor()
      outs.join(); errs.join()
    else:
      import threading, subprocess, os
      if 'cwd' in cmd and cmd['cwd'] != None:
        os.chdir(VLAB.expandEnv(cmd['cwd']))
      if 'env' in cmd and cmd['env'] != None:
        cmdenv = cmd['env']
        for e in cmdenv:
          os.putenv(e, VLAB.expandEnv(cmdenv[e]))
      if 'stdin' in cmd and cmd['stdin'] != None:
        proc = subprocess.Popen(cmdLine, stdin=open(VLAB.expandEnv(cmd['stdin']),'r'),stdout=subprocess.PIPE, stderr=subprocess.PIPE)
      else:
        proc = subprocess.Popen(cmdLine, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
      t1 = threading.Thread(target=VLAB.helper, args=('out', proc.stdout))
      t2 = threading.Thread(target=VLAB.helper, args=('err', proc.stderr))
      t1.start(); t2.start()
      exitCode = proc.wait()
      t1.join(); t2.join()
    VLAB.logger.info('exitCode=%d' % exitCode)
  doExec = staticmethod(doExec)

  def valuesfromfile(path, transpose=False):
    """Returns a 2D array with the values of the csv file at 'path'.

    Keyword arguments:
    transpose -- transpose the matrix with the values

    """
    fp = open(path)
    values = [line.strip().split() for line in fp
              if not line.startswith('#')]
    fp.close()
    values = [[float(value) for value in row] for row in values]
    if transpose:
      values = zip(*values)
    return values
  valuesfromfile = staticmethod(valuesfromfile)
  def fabsa(x):
    """Return the element-wise abs() of the given array."""
    return map(lambda x : math.fabs(x), x)
  fabsa = staticmethod(fabsa)
  def cosa(arcs):
    """Return the element-wise cos() of 'arcs' array given in degrees."""
    return map(lambda x : math.cos(VLAB.d2r(x)), arcs)
  cosa = staticmethod(cosa)
  def replacerectinarray(array, replacement, xul, yul, xlr, ylr):
    """Replace the array with the specified rectangle substituted with
    replacement.
    (array[xul:xlr,yul:ylr])
     +---------------+
     |(xul, yul)     |
     |               |
     |     (xlr, ylr)|
     +---------------+
    """
    ri = 0
    for x in xrange(xul, xlr):
      array[x][yul:ylr] = replacement[ri]
      ri += 1
    return array
  replacerectinarray = staticmethod(replacerectinarray)
  def replaceinarray(haystack, predicate, replacement):
    """Return 'haystack' with 'predicate' matches replaced by 'replacement'"""
    return map(lambda item : {True: replacement, False: item}[predicate(item)],
               haystack)
  replaceinarray = staticmethod(replaceinarray)
  def sqrta(values):
    """Return the element-wise sqrt() of the given array."""
    return map(math.sqrt, values)
  sqrta = staticmethod(sqrta)
  def suba(lista, listb):
    """Subtract the values of a list from the values of another list."""
    if len(lista) != len(listb):
      raise ValueError("Arguments have to be of same length.")
    return map(operator.sub, lista, listb)
  suba = staticmethod(suba)
  def adda(lista, listb):
    """Add the values of a list to the values of another list."""
    if len(lista) != len(listb):
      raise ValueError("Arguments have to be of same length.")
    return map(operator.add, lista, listb)
  adda = staticmethod(adda)
  def diva(lista, listb):
    """Return the element-wise division of 'lista' by 'listb'."""
    if len(lista) != len(listb):
      raise ValueError("Arguments have to be of same length.")
    return map(operator.div, lista, listb)
  diva = staticmethod(diva)
  def mula(lista, listb):
    """Return the element-wise multiplication of 'lista' with 'listb'."""
    if len(lista) != len(listb):
      raise ValueError("Arguments have to be of same length.")
    return map(operator.mul, lista, listb)
  mula = staticmethod(mula)
  def powa(values, exponent):
    """Returns the element-wise exp('exponent') of the given array."""
    if isinstance(exponent, (list, tuple)):
      return map(lambda x, y : x ** y, values, exponent)
    return map(lambda x : x ** exponent, values)
  powa = staticmethod(powa)
  def treemap(fn, tree):
    """Applies `fn' to every value in `tree' which isn't a list and
    returns a list with the same shape as tree and the value of `fn'
    applied to the values in their place.

    """
    result = []
    for node in tree:
      if isinstance(node, (list, tuple)):
        result += [VLAB.treemap(fn, node)]
      else:
        result += [fn(node)]
    return result
  treemap = staticmethod(treemap)
  def makemaskeq(array, value):
    return VLAB.treemap(lambda x : x == value, array)
  makemaskeq = staticmethod(makemaskeq)
  def awhere(mask):
    """Returns the coordinates of the cells which evaluate true."""
    result = []
    if isinstance(mask, (list, tuple)):
      for i, cell in enumerate(mask):
        result += [[i] + sub for sub in VLAB.awhere(cell)
                   if isinstance(sub, (list, tuple))]
      return result
    else:
      if mask:
        return [[]]
      else:
        return []
  awhere = staticmethod(awhere)
  def aclone(tree):
    """Make a deep copy of `tree'."""
    if isinstance(tree, (list, tuple)):
      return list(map(VLAB.aclone, tree))
    return tree
  aclone = staticmethod(aclone)
  def make_chart(title, x_label, y_label, dataset):
    """Helper for creating Charts"""
    if sys.platform.startswith('java'):
      from org.jfree.chart import ChartFactory, ChartFrame, ChartUtilities
      from org.jfree.chart.axis import NumberTickUnit
      from org.jfree.chart.plot import PlotOrientation
      from org.jfree.data.xy import XYSeries, XYSeriesCollection
      chart = ChartFactory.createScatterPlot(title, x_label, y_label, dataset,
                                           PlotOrientation.VERTICAL, True, 
                                           True, False)
      plot = chart.getPlot()
      domain_axis = plot.getDomainAxis()
      domain_axis.setRange(-70, 70)
      range_axis = plot.getRangeAxis()
      range_axis.setRange(0.0, 0.4)
      range_axis.setTickUnit(NumberTickUnit(0.05))
      return chart
    else:
      # TODO: re-merge original python implementation
      raise Exception("original make_chart()")
      return None
  make_chart = staticmethod(make_chart)
  def make_dataset():
    """Dataset helper for supporting chart creation"""
    if sys.platform.startswith('java'):
      from org.jfree.data.xy import XYSeriesCollection
      return XYSeriesCollection()
    else:
      # TODO: re-merge original python implementation
      raise Exception("original make_dataset()")
      return None
  make_dataset = staticmethod(make_dataset)
  def plot(dataset, x, y, label):
    """plot dataset with x and y values and the given label"""
    if sys.platform.startswith('java'):
      from org.jfree.data.xy import XYSeries
      series = XYSeries(label)
      for next_x, next_y in zip(x, y):
        series.add(float(next_x), float(next_y))
      dataset.addSeries(series)
    else:
      # TODO: re-merge original python implementation
      raise Exception("original plot()")
      return None
  plot = staticmethod(plot)
  def save_chart(chart, filename):
    """save the generated chart in the given filename"""
    if sys.platform.startswith('java'):
      from java.io import File
      from org.jfree.chart import ChartUtilities
      ChartUtilities.saveChartAsPNG(File(filename), chart, 800, 600)
    else:
      # TODO: re-merge original python implementation
      raise Exception("save_chart")
      pass
  save_chart = staticmethod(save_chart)
  def maskand(array, mask):
    return map(lambda a, b : a & b, array, mask)
  maskand = staticmethod(maskand)
  def unique(array):
    """return unique values in the given input array"""
    sortedarray = list(array)
    sortedarray.sort()
    result = []
    def addifunique(x, y):
      if x != y:
        result.append(x)
      return y
    result.append(reduce(addifunique, sortedarray))
    return result
  unique = staticmethod(unique)
  def ravel(a):
    def add(a, i):
      if not(isinstance(a, (list, tuple))):
        a = [a]
      if isinstance(i, (list, tuple)):
        return a + ravel(i)
      else:
        return a + [i]
    return reduce(add, a)
  ravel = staticmethod(ravel)
  def min_l_bfgs_b(f, initial_guess, args=(), bounds=None, epsilon=1e-8):
    """The `approx_grad' is not yet implemented, because it's not yet
used."""
    if sys.platform.startswith('java'):
      from java.lang import ClassLoader
      from java.io import File
      for path in sys.path:
        if path.endswith('.jar'):
          libdir = File(path).getParent()
          break
      VLAB.logger.info(libdir)
      # hack taken from: http://blog.cedarsoft.com/2010/11/setting-java-library-path-programmatically
      System.setProperty("java.library.path", libdir)
      syspathfield = ClassLoader.getDeclaredField("sys_paths")
      syspathfield.setAccessible(True);
      syspathfield.set(None, None)
      from lbfgsb import DifferentiableFunction, FunctionValues, Bound, Minimizer
      class function(DifferentiableFunction):
        def __init__(self, function, args=()):
          self.function = function
          self.args = args
        def getValues(self, x):
          f = self.function(x, *self.args)
          g = VLAB.approx_fprime(x, self.function, epsilon, *self.args)
          return FunctionValues(f, g)
      min = Minimizer()
      min.setBounds([Bound(bound[0], bound[1]) for bound in bounds])
      result = min.run(function(f, args), VLAB.ravel(initial_guess))
      return result.point
  min_l_bfgs_b = staticmethod(min_l_bfgs_b)
  def approx_fprime(xk, f, epsilon, *args):
    f0 = f(xk, *args)
    grad = [0. for i in xrange(len(xk))]
    ei = [0. for i in xrange(len(xk))]
    for k in xrange(len(xk)):
      ei[k] = 1.0
      d = map(lambda i : i * epsilon, ei)
      grad[k] = (f(VLAB.adda(xk, d), *args) - f0) / d[k]
      ei[k] = 0.0
    return grad
  approx_fprime = staticmethod(approx_fprime)
  def doLibradtran(args):

    q = {
      'rpv_file':  args['rpv_file'],
      'outfile' :  args['outfile'],
      'infile'  :  args['infile'],
    }

    for k in args:
      q[k] = args[k]

    if args['scene'] == VLAB.K_LAEGERN:
      q['latitude']  = '47.481667'
      q['longitude'] = '-8.394722,17'
    elif args['scene'] == VLAB.K_THARANDT:
      q['latitude']  = '50.9676498'
      q['longitude'] = '-13.520354'

    if args['aerosol'] == VLAB.K_RURAL:
      q['aerosol'] = '1'
    elif args['aerosol'] == VLAB.K_MARITIME:
      q['aerosol'] = '2'
    elif args['aerosol'] == VLAB.K_URBAN:
      q['aerosol'] = '5'
    elif args['aerosol'] == VLAB.K_TROPOSPHERIC:
      q['aerosol'] = '6'

    # TODO: verify these
    if args['sensor'] ==   VLAB.K_SENTINEL2:
      q['wavelength'] = '443 2190'
    elif args['sensor'] == VLAB.K_SENTINEL3:
      q['wavelength'] = '400 2400'
    elif args['sensor'] == VLAB.K_MODIS:
      q['wavelength'] = '405 2155'
    elif args['sensor'] == VLAB.K_MERIS:
      q['wavelength'] = '412 900'
    elif args['sensor'] == VLAB.K_LANDSAT_OLI:
      q['wavelength'] = '433 2300'
    elif args['sensor'] == VLAB.K_LANDSAT_ETM:
      q['wavelength'] = '450 2300'

    if VLAB.osName().startswith('Windows'):
      q['solar_file'] = VLAB.expandEnv('%HOMEDRIVE%%HOMEPATH%\\.beam\\beam-vlab\\auxdata\\libRadtran_win32\\data\\solar_flux\\NewGuey2003.day')
    else:
      q['solar_file'] = VLAB.expandEnv('$HOME/.beam/beam-vlab/auxdata/libRadtran_lin64/data/solar_flux/NewGuey2003.dat')

    sdata = """
solar_file         %s
correlated_k       LOWTRAN
rte_solver         cdisort
rpv_file           %s
deltam             on
nstr               6
zout               TOA
output_user        lambda uu
quiet
""" % (q['solar_file'],
      q['rpv_file'])

    if 'aerosol' in q:
      sdata += "aerosol_default\n"
      sdata += "aerosol_haze       %s\n" % (q['aerosol']) 
    if 'O3' in q:
      sdata += "dens_column        O3 %s\n" % (q['O3']) 
    if 'CO2' in q:
      sdata += "co2_mixing_ratio   %s\n" % (q['CO2']) 
    if 'H2O' in q:
      sdata += "h2o_mixing_ratio   %s\n" % (q['H2O']) 
    if 'umu' in q:
      sdata += "umu                %s\n" % (q['umu']) 
    if 'phi' in q:
      sdata += "phi                %s\n" % (q['phi']) 
    if 'latitude' in q:
      sdata += "latitude           %s\n" % (q['latitude']) 
    if 'longitude' in q:
      sdata += "longitude          %s\n" % (q['longitude']) 
    if 'time' in q:
      sdata += "time               %s\n" % (q['time'])
    if 'sza' in q:
      sdata += "sza                %s\n" % (q['sza'])
    if 'phi0' in q:
      sdata += "phi0               %s\n" % (q['phi0'])
    if 'wavelength' in q:
      sdata += "wavelength         %s\n" % (q['wavelength'])

    fp = open(q['infile'], 'w')
    fp.write(sdata)
    fp.close()

    cmd = {
     'linux' : {
       'cwd'     : '$HOME/.beam/beam-vlab/auxdata/libRadtran_lin64/examples',
       'exe'     : '$HOME/.beam/beam-vlab/auxdata/libRadtran_lin64/bin/uvspec',
       'cmdline' : [],
       'stdin'   : q['infile'],
       'stdout'  : q['outfile'],
       'stderr'  : None,
       'env'     : None
     },
     'windows' : {
       'cwd'     : '%HOMEDRIVE%%HOMEPATH%\\.beam\\beam-vlab\\auxdata\\libRadtran_win32\\examples',
       'exe'     : '%HOMEDRIVE%%HOMEPATH%\\.beam\\beam-vlab\\auxdata\\libRadtran_win32\\uvspec.exe',
       'cmdline' : [],
       'stdin'   : q['infile'],
       'stdout'  : q['outfile'],
       'stderr'  : None,
       'env'     : None
    }
    }
    VLAB.logger.info('%s: spawning libradtran...' % VLAB.me())
    VLAB.doExec(cmd)
  doLibradtran = staticmethod(doLibradtran)


  ###########################################################################
  # 
  # Minimize_* functions...
  #
  # Nelder-Mead simplex minimization of a nonlinear (multivariate) function.
  # 
  # The programming interface is via the minimize() function; see below.
  # 
  # This code has been adapted from the C-coded nelmin.c which was
  # adapted from the Fortran-coded nelmin.f which was, in turn, adapted
  # from the papers
  # 
  #   J.A. Nelder and R. Mead (1965)
  #   A simplex method for function minimization.
  #   Computer Journal, Volume 7, pp 308-313.
  # 
  #   R. O'Neill (1971)
  #   Algorithm AS47. Function minimization using a simplex algorithm.
  #   Applied Statistics, Volume 20, pp 338-345.
  # 
  # and some examples are in
  # 
  #   D.M. Olsson and L.S. Nelson (1975)
  #   The Nelder-Mead Simplex procedure for function minimization.
  #   Technometrics, Volume 17 No. 1, pp 45-51.
  #    
  # For a fairly recent and popular incarnation of this minimizer,
  # see the amoeba function in the famous "Numerical Recipes" text.
  # 
  # P. Jacobs
  # School of Engineering, The University of Queensland
  # 07-Jan-04
  # 
  # Modifications by C. Schenkel
  # Netcetera
  # 31-Oct-13
  # 
  ###########################################################################
  def Minimize_create_new_point(c1, p1, c2, p2):
    """
    Create a new N-dimensional point as a weighting of points p1 and p2.
    """
    p_new = []
    for j in range(len(p1)):
      p_new.append(c1 * p1[j] + c2 * p2[j])
    return p_new
  Minimize_create_new_point = staticmethod(Minimize_create_new_point)

  def Minimize_take_a_step(smplx, Kreflect, Kextend, Kcontract):
    """
    Try to move away from the worst point in the simplex.
  
    The new point will be inserted into the simplex (in place).
    """
    i_low = smplx.lowest()
    i_high = smplx.highest()
    x_high = smplx.vertex_list[i_high]
    f_high = smplx.f_list[i_high]
    # Centroid of simplex excluding worst point.
    x_mid = smplx.centroid(i_high)
    f_mid = smplx.f(x_mid)
    smplx.nfe += 1
  
    # First, try moving away from worst point by
    # reflection through centroid
    x_refl = VLAB.Minimize_create_new_point(1.0+Kreflect, x_mid, -Kreflect, x_high)
    f_refl = smplx.f(x_refl)
    smplx.nfe += 1
    if f_refl < f_mid:
      # The reflection through the centroid is good,
      # try to extend in the same direction.
      x_ext = VLAB.Minimize_create_new_point(Kextend, x_refl, 1.0-Kextend, x_mid)
      f_ext = smplx.f(x_ext)
      smplx.nfe += 1
      if f_ext < f_refl:
        # Keep the extension because it's best.
        smplx.replace_vertex(i_high, x_ext, f_ext)
      else:
        # S…