code: def impose_expectation(m, f, npts, bounds=None, weights=None, **kwds):
      """impose a given expectation value ``E`` on a given function *f*,
      where ``E = m +/- tol`` and ``E = mean(f(x))`` for ``x`` in *bounds*
  
      Args:
          m (float): target expected mean
          f (func): a function that takes a list and returns a number
          npts (tuple(int)): a tuple of dimensions of the target product measure
          bounds (tuple, default=None): tuple is ``(lower_bounds, upper_bounds)``
          weights (list, default=None): a list of sample weights
          tol (float, default=None): maximum allowable deviation from ``m``
          constraints (func, default=None): a function that takes a nested list of
              ``N x 1D`` discrete measure positions and weights, with the intended
              purpose of kernel-transforming ``x,w`` as ``x' = constraints(x, w)``
          npop (int, default=200): size of the trial solution population
          maxiter (int, default=1000): the maximum number of iterations to perform
          maxfun (int, default=1e+6): the maximum number of function evaluations
  
      Returns:
          a list of sample positions, with expectation ``E``
  
      Notes:
          Expectation value ``E`` is calculated by minimizing ``mean(f(x)) - m``,
          over the given *bounds*, and will terminate when ``E`` is found within
          deviation ``tol`` of the target mean ``m``.  If ``tol`` is not provided,
          then a relative deviation of 1% of ``m`` will be used.
  
          This function does not preserve the mean, variance, or range, as there
          is no initial list of samples to draw the mean, variance, and etc from.
  
          *bounds* is tuple with ``length(bounds) == 2``, composed of all the lower
          bounds, then all the upper bounds, for each positional value (i.e. bounds
          for the weights should not be included).
  
      Examples:
          >>> # provide the dimensions and bounds
          >>> nx = 3;  ny = 2;  nz = 1
          >>> x_lb = [10.0];  y_lb = [0.0];  z_lb = [10.0]
          >>> x_ub = [50.0];  y_ub = [9.0];  z_ub = [90.0]
          >>>
          >>> # prepare the bounds
          >>> lb = (nx * x_lb) + (ny * y_lb) + (nz * z_lb)
          >>> ub = (nx * x_ub) + (ny * y_ub) + (nz * z_ub)
          >>>
          >>> # generate a list of samples with mean +/- dev imposed
          >>> mean = 2.0;  dev = 0.01
          >>> samples = impose_expectation(mean, f, (nx,ny,nz), (lb,ub), tol=dev)
          >>>
          >>> # test the results by calculating the expectation value for the samples
          >>> expectation(f, samples)
          >>> 2.000010010122465"""
      # param[0] is the target mean
      # param[1] is the acceptable deviation from the target mean
      tol = kwds["tol"] if "tol" in kwds else None
      param = (m, tol or 0.01 * m)
  
      # FIXME: the following is a HACK to recover from lost 'weights' information
      #        we 'mimic' discrete measures using the product measure weights
      # plug in the 'constraints' function:  samples' = constrain(samples, weights)
      constrain = None  # default is no constraints
      if "constraints" in kwds:
          constrain = kwds["constraints"]
      if not constrain:  # if None (default), there are no constraints
          constraints = lambda x: x
      else:  # XXX: better to use a standard "xk' = constrain(xk)" interface ?
  
          def constraints(rv):
              coords = _pack(_nested(rv, npts))
              coords = list(zip(*coords))  # 'mimic' a nested list
              coords = constrain(
                  coords, [weights for i in range(len(coords))]
              )
              coords = list(
                  zip(*coords)
              )  # revert back to a packed list
              return _flat(_unpack(coords, npts))
  
      # construct cost function to reduce deviation from expectation value
      def cost(rv):
          """compute cost from a 1-d array of model parameters,
          where:  cost = | E[model] - m |**2"""
          # from mystic.math.measures import _pack, _nested, expectation
          samples = _pack(_nested(rv, npts))
          Ex = expectation(f, samples, weights)
          return (Ex - param[0]) ** 2
  
      # if bounds are not set, use the default optimizer bounds
      if not bounds:
          lower_bounds = []
          upper_bounds = []
          for n in npts:
              lower_bounds += [None] * n
              upper_bounds += [None] * n
      else:
          lower_bounds, upper_bounds = bounds
  
      # construct and configure optimizer
      debug = kwds["debug"] if "debug" in kwds else False
      npop = kwds.pop("npop", 200)
      maxiter = kwds.pop("maxiter", 1000)
      maxfun = kwds.pop("maxfun", 1e6)
      crossover = 0.9
      percent_change = 0.9
  
      def optimize(cost, bounds, tolerance, _constraints):
          (lb, ub) = bounds
          from mystic.solvers import DifferentialEvolutionSolver2
          from mystic.termination import VTR
          from mystic.strategy import Best1Exp
          from mystic.monitors import VerboseMonitor, Monitor
          from mystic.tools import random_seed
  
          if debug:
              random_seed(123)
          evalmon = Monitor()
          stepmon = Monitor()
          if debug:
              stepmon = VerboseMonitor(10)
  
          ndim = len(lb)
          solver = DifferentialEvolutionSolver2(ndim, npop)
          solver.SetRandomInitialPoints(min=lb, max=ub)
          solver.SetStrictRanges(min=lb, max=ub)
          solver.SetEvaluationLimits(maxiter, maxfun)
          solver.SetEvaluationMonitor(evalmon)
          solver.SetGenerationMonitor(stepmon)
          solver.Solve(
              cost,
              termination=VTR(tolerance),
              strategy=Best1Exp,
              CrossProbability=crossover,
              ScalingFactor=percent_change,
              constraints=_constraints,
          )  # XXX: parallel map?
  
          solved = solver.Solution()
          diameter_squared = solver.bestEnergy
          func_evals = len(evalmon)
          return solved, diameter_squared, func_evals
  
      # use optimization to get expectation value
      tolerance = (param[1]) ** 2
      results = optimize(
          cost, (lower_bounds, upper_bounds), tolerance, constraints
      )
  
      # repack the results
      samples = _pack(_nested(results[0], npts))
      return samples
  
imports: from mystic.math.stats import *
  from mystic.math.samples import *
  from mystic.math.integrate import *
  from mystic.symbolic import generate_solvers, generate_constraint, solve
  from mystic.symbolic import generate_conditions, generate_penalty
  from mystic.math import almostEqual
  from mystic.math.distance import Lnorm
  from numpy import asarray
label_desc: Use image node_id 28 for calling a global function with example usage: impose_expectation(m, f, npts, bounds, weights) and returns: samples, unknown, solved,  diameter_squared,  func_evals, _flat