code: def flux_randomization(model, threshold, tfba, solver):
      
  
      optimize = {}
      for reaction_id in model.reactions:
          if model.is_reversible(reaction_id):
              optimize[reaction_id] = 2*random.random() - 1.0
          else:
              optimize[reaction_id] = random.random()
  
      fba = _get_fba_problem(model, tfba, solver)
      for reaction_id, value in iteritems(threshold):
          fba.prob.add_linear_constraints(fba.get_flux_var(reaction_id) >= value)
  
      fba.maximize(optimize)
      for reaction_id in model.reactions:
          yield reaction_id, fba.get_flux(reaction_id)
docs: Find a random flux solution on the boundary of the solution space.
  
      The reactions in the threshold dictionary are constrained with the
      associated lower bound.
  
      Args:
          model: MetabolicModel to solve.
          threshold: dict of additional lower bounds on reaction fluxes.
          tfba: If True enable thermodynamic constraints.
          solver: LP solver instance to use.
  
      Returns:
          An iterator of reaction ID and reaction flux pairs.