Phase.java

package org.mklab.sdpj.algorithm;

import java.io.PrintStream;

import org.mklab.mpf.ap.EFFloat;
import org.mklab.mpf.ap.EFFloatRational;
import org.mklab.mpf.ap.MPFloat;
import org.mklab.nfc.matrix.ComplexNumericalMatrix;
import org.mklab.nfc.matrix.RealNumericalMatrix;
import org.mklab.nfc.scalar.ComplexNumericalScalar;
import org.mklab.nfc.scalar.RealNumericalScalar;
import org.mklab.sdpj.tool.Tools;


/**
 * class <code>Phase</code> manage the information of type of phase, which is the message of solution. About more detail information,
 * please refer to web site:
 *  <a href = "http://sdpa.sourceforge.net/" style=text-decoration:none>
 * <tt>http://sdpa.sourceforge.net/</tt></a>.<br>  
 * 
 * @author koga
 * @version $Revision$, 2009/04/24
 * @param <RS> type of real scalar
 * @param <RM> type of real matrix
 * @param <CS> type of complex scalar
 * @param <CM> type of complex Matrix
 */
public class Phase<RS extends RealNumericalScalar<RS, RM, CS, CM>, RM extends RealNumericalMatrix<RS, RM, CS, CM>, CS extends ComplexNumericalScalar<RS, RM, CS, CM>, CM extends ComplexNumericalMatrix<RS, RM, CS, CM>> {

  /**the dimension of array of block-matrices matrix*/
  private int nDim;
  /**type of phase, which is the message of solution */
  private PhaseType value;
  /***/
  private RS gapValue;

  /**
   * constructor of class <code>Phase</code>.
   * <br>新しく生成された<code>Phase</code>オブジェクトを初期化します。
   * 
   * @param initRes initial residual
   * @param solveInfo information of solution
   * @param param setting parameter
   * @param nDim dimension of array of block-matrices matrix
   */
  public Phase(Residuals<RS,RM,CS,CM> initRes, SolveInfo<RS,RM,CS,CM> solveInfo, Parameter<RS,RM,CS,CM> param, int nDim) {
    this.setnDim(nDim);
    updateCheck(initRes, solveInfo, param);
  }

  /**
   * update the information of computing of each iteration. About more detail of implementing algorithm,
   * please refer to codes of <code>SDPA</code> as below.<br>
   * <code>(rsdpa_parts.cpp : SDPA 6.2.0 Line No. : 2319)</code>
   * 
   * @param currentRes current residual
   * @param solveInfo information of primal-dual solution
   * @param param setting parameter
   * @return boolean if computing stop, false; else, true.
   */
  public boolean updateCheck(Residuals<RS,RM,CS,CM> currentRes, SolveInfo<RS,RM,CS,CM> solveInfo, Parameter<RS,RM,CS,CM> param) {
    final boolean guarantee = false;
    final RS accuracy = param.epsilonDash;
    final RS unit = accuracy.create(1);
    final RS NONZERO = unit.getMachineEpsilon().sqrt();
    this.value = PhaseType.noINFO;

    if (currentRes.getNormPrimalVec().isLessThanOrEquals(accuracy)) {
      if (currentRes.getNormDualMat().isLessThanOrEquals(accuracy)) {
        this.value = PhaseType.pdFEAS;
      } else {
        this.value = PhaseType.pFEAS;
      }
    }
    if (this.getValue().equals(PhaseType.noINFO) && currentRes.getNormDualMat().isLessThanOrEquals(accuracy)) {
      this.value = PhaseType.dFEAS;
    }
    if (this.getValue() == PhaseType.pdFEAS) {
//      RS mean =((solveInfo.getObjValPrimal().abs()).add(solveInfo.getObjValDual().abs())).divide(2);
//      RS PDgap =(solveInfo.getObjValPrimal().subtract(solveInfo.getObjValDual())).abs();
      
      RS relativeGap;
      if(solveInfo.getObjValPrimal().createUnit() instanceof MPFloat){
        EFFloat objValPrimal = new EFFloat((MPFloat)solveInfo.getObjValPrimal());
        EFFloat objValDual = new EFFloat((MPFloat)solveInfo.getObjValDual());
        EFFloat effloatMean = objValPrimal.abs().add(objValDual.abs()).divide(2);
        EFFloat effloatPDgap = objValPrimal.subtract(objValDual).abs();
        if(effloatMean.isLessThan(1)){
          relativeGap =(RS)effloatPDgap.toMPFloatWithDefaultPrecision();
        } else {
          EFFloatRational eM = new EFFloatRational(effloatMean);
          EFFloatRational eG = new EFFloatRational(effloatPDgap);
          EFFloatRational rG = eG.divide(eM);
          relativeGap =(RS)rG.toMPFloat();
          //relativeGap =effloatPDgap.divide(effloatMean).toMPFloatWithDefaultPrecision();
        }
      } else {
        RS mean =((solveInfo.getObjValPrimal().abs()).add(solveInfo.getObjValDual().abs())).divide(2);
        RS PDgap =(solveInfo.getObjValPrimal().subtract(solveInfo.getObjValDual())).abs();
        if(mean.isLessThan(1)){
          relativeGap =PDgap;
        } else {
          relativeGap =PDgap.divide(mean);
        }
      }
      this.gapValue = relativeGap;

      if (relativeGap.isLessThanOrEquals(param.epsilonStar)) {
        this.value = PhaseType.pdOPT;
        return false;
      }

//      RS dominator;
//
//      if (mean.equals(0)) {
//        dominator = unit.create(1);
//      } else {
//        dominator = mean;
//      }
//
//      if ((PDgap.divide(dominator)).isLessThanOrEquals(param.epsilonStar)) {
//        this.value = PhaseType.pdOPT;
//        return false;
//      }
    }
    if (!guarantee) {
      if (this.getValue().equals(PhaseType.noINFO) && solveInfo.getRho().isGreaterThan(NONZERO.add(1))) {
        this.value = PhaseType.pdINF;
        return false;
      }
    }
    if (this.getValue().equals(PhaseType.pFEAS)) {
      if (Tools.REVERSE_PRIMAL_DUAL == 1) {
        if (solveInfo.getObjValPrimal().isLessThanOrEquals(param.upperBound.unaryMinus())) {
          this.value = PhaseType.pUNBD;
          return false;
        }
      } else {
        if (solveInfo.getObjValPrimal().isLessThanOrEquals(param.lowerBound)) {
          this.value = PhaseType.pUNBD;
          return false;
        }
      }
      if (solveInfo.getRho().isGreaterThan(NONZERO.add(1))) {
        this.value = PhaseType.pFEAS_dINF;
        return false;
      }
    }

    if (this.getValue().equals(PhaseType.dFEAS)) {
      if (Tools.REVERSE_PRIMAL_DUAL == 1) {
        if (solveInfo.getObjValDual().isGreaterThanOrEquals(param.lowerBound.unaryMinus())) {
          this.value = PhaseType.dUNBD;
          return false;
        }
      } else {
        if (solveInfo.getObjValDual().isGreaterThanOrEquals(param.upperBound)) {
          this.value = PhaseType.dUNBD;
          return false;
        }
      }
      if (solveInfo.getRho().isGreaterThan(NONZERO.add(1))) {
        this.value = PhaseType.pINF_dFEAS;
        return false;
      }
    }
    return true;
  }

  /**
   * reverse the computing data. About more detail information, please the codes of <code>SDPA</code>.<br> 
   * <code>(rsdpa_parts.cpp : SDPA 6.2.0 Line No. 2407)</code>
   */
  public void reverse() {
    if (Tools.REVERSE_PRIMAL_DUAL == 1) {
      switch (this.getValue()) {
        case noINFO:
          break;
        case pFEAS:
          this.value = PhaseType.dFEAS;
          break;
        case dFEAS:
          this.value = PhaseType.pFEAS;
          break;
        case pdFEAS:
          break;
        case pdINF:
          break;
        case pFEAS_dINF:
          this.value = PhaseType.pINF_dFEAS;
          break;
        case pINF_dFEAS:
          this.value = PhaseType.pFEAS_dINF;
          break;
        case pdOPT:
          break;
        case pUNBD:
          this.value = PhaseType.dUNBD;
          break;
        case dUNBD:
          this.value = PhaseType.pUNBD;
          break;
        default:
          break;
      }
    }
  }

  //  /**
  //   * 標準出力に出力します。
  //   */
  //  public void display() {
  //    Tools.message(toString());
  //  }

  /**
   * @param ps 出力先
   */
  public void display(PrintStream ps) {
    Tools.message(toString(), ps);
  }

  /**
   * @see java.lang.Object#toString()
   */
  @Override
  public String toString() {

    String str;
    switch (this.getValue()) {
      case noINFO:
        str = "noINFO    "; //$NON-NLS-1$
        break;
      case pFEAS:
        str = "pFEAS     "; //$NON-NLS-1$
        break;
      case dFEAS:
        str = "dFEAS     "; //$NON-NLS-1$
        break;
      case pdFEAS:
        str = "pdFEAS    "; //$NON-NLS-1$
        break;
      case pdINF:
        str = "pdINF     "; //$NON-NLS-1$
        break;
      case pFEAS_dINF:
        str = "pFEAS_dINF"; //$NON-NLS-1$
        break;
      case pINF_dFEAS:
        str = "pINF_dFEAS"; //$NON-NLS-1$
        break;
      case pdOPT:
        str = "pdOPT     "; //$NON-NLS-1$
        break;
      case pUNBD:
        str = "pUNBD     "; //$NON-NLS-1$
        break;
      case dUNBD:
        str = "dUNBD     "; //$NON-NLS-1$
        break;
      default:
        str = "phase error"; //$NON-NLS-1$
        Tools.message("rPhase:: phase error"); //$NON-NLS-1$
        break;
    }
    return "phase.value = " + str + "\n"; //$NON-NLS-1$ //$NON-NLS-2$
  }

  /**
   * valueを返します。
   * 
   * @return value
   */
  public PhaseType getValue() {
    return this.value;
  }

  /**
   * @param nDim nDimを設定します。
   */
  public void setnDim(int nDim) {
    this.nDim = nDim;
  }

  /**
   * @return nDimを返します。
   */
  public int getnDim() {
    return this.nDim;
  }

  /**
   * @return gapValue value of relativeGap
   */
  public RS getGapValue() {
    return this.gapValue;
  }
}