Number alpha_dual_aff = IpCq().dual_frac_to_the_bound(1.0,
                            *step->z_L(),
                            *step->z_U(),
                            *step->v_L(),
                            *step->v_U());

    Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
                   "  The affine maximal step sizes are\n"
                   "   alpha_primal_aff = %23.16e\n"
                   "   alpha_dual_aff = %23.16e\n",
                   alpha_primal_aff,
                   alpha_dual_aff);

    // now compute the average complementarity at the affine step
    // ToDo shoot for mu_min instead of 0?
    Number mu_aff = CalculateAffineMu(alpha_primal_aff, alpha_dual_aff, *step);
    Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
                   "  The average complementariy at the affine step is %23.16e\n",
                   mu_aff);

    // get the current average complementarity
    Number mu_curr = IpCq().curr_avrg_compl();
    Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
                   "  The average complementariy at the current point is %23.16e\n",
                   mu_curr);
    DBG_ASSERT(mu_curr>0.);

    // Apply Mehrotra's rule
    Number sigma = pow((mu_aff/mu_curr),3);
    // Make sure, sigma is not too large
    sigma = Min(sigma, sigma_max_);

    Number mu = sigma*mu_curr;

    // Store the affine search direction (in case it is needed in the
    // line search for a corrector step)
    IpData().set_delta_aff(step);
    IpData().SetHaveAffineDeltas(true);

    char ssigma[40];
    sprintf(ssigma, " sigma=%8.2e", sigma);
    IpData().Append_info_string(ssigma);
    //sprintf(ssigma, " xi=%8.2e ", IpCq().curr_centrality_measure());
    //IpData().Append_info_string(ssigma);

    new_mu = Max(Min(mu, mu_max), mu_min);
    return true;
  }

  Number ProbingMuOracle::CalculateAffineMu
  (
    Number alpha_primal,
    Number alpha_dual,
    const IteratesVector& step)
  {
    // Get the current values of the slack variables and bound multipliers
    SmartPtr<const Vector> slack_x_L = IpCq().curr_slack_x_L();
    SmartPtr<const Vector> slack_x_U = IpCq().curr_slack_x_U();
    SmartPtr<const Vector> slack_s_L = IpCq().curr_slack_s_L();
    SmartPtr<const Vector> slack_s_U = IpCq().curr_slack_s_U();

    SmartPtr<const Vector> z_L = IpData().curr()->z_L();
    SmartPtr<const Vector> z_U = IpData().curr()->z_U();
    SmartPtr<const Vector> v_L = IpData().curr()->v_L();
    SmartPtr<const Vector> v_U = IpData().curr()->v_U();

    SmartPtr<Vector> tmp_slack;
    SmartPtr<Vector> tmp_mult;
    SmartPtr<const Matrix> P;
    Index ncomp = 0;
    Number sum =0.;

    // For each combination of slack and multiplier, compute the new
    // values and their dot products.

    // slack_x_L
    if (slack_x_L->Dim()>0) {
      ncomp += slack_x_L->Dim();

      P = IpNLP().Px_L();
      tmp_slack = slack_x_L->MakeNew();
      tmp_slack->Copy(*slack_x_L);
      P->TransMultVector(alpha_primal, *step.x(), 1.0, *tmp_slack);

      tmp_mult = z_L->MakeNew();
      tmp_mult->Copy(*z_L);