Solution.java

package org.mklab.sdpj.algorithm;

import java.io.PrintStream;

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.algebra.Algebra;
import org.mklab.sdpj.datastructure.BlockDenseMatrix;
import org.mklab.sdpj.datastructure.BlockVector;
import org.mklab.sdpj.datastructure.Vector;
import org.mklab.sdpj.tool.Tools;


/**
 * class <code>Solution</code> is defined to save the solution information of <code>SDP</code> problem solved by solver.
 * 
 * @author takafumi
 * @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 Solution<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>> {

  /**unit value based on a specified data type */
  private RS unit;
  /**dense block-matrices matrix <code>X</code> := <code>Y</code> */
  private BlockDenseMatrix<RS,RM,CS,CM> xMat;
  /**dense block-matrices matrix <code>Z</code> := <code>X</code>*/
  private BlockDenseMatrix<RS,RM,CS,CM> zMat;
  /**vector <code>y</code> := <code>x</code> */
  private Vector<RS,RM,CS,CM> yVec;
  /**<code>mu</code> := <code> X・Z</code> := <code> Y・X</code> */
  BlockDenseMatrix<RS,RM,CS,CM> xMatzMat;
  /**cholesky decomposition of block-matrices matrix <code>X</code> := <code>Y</code> */
  BlockDenseMatrix<RS,RM,CS,CM> choleskyX;
  /**inverse matrix of cholesky decomposition of block-matrices matrix <code>X</code> := <code>Y</code> */
  BlockDenseMatrix<RS,RM,CS,CM> invCholeskyX;
  /**cholesky decomposition of block-matrices matrix <code>Z</code> := <code>X</code>*/
  BlockDenseMatrix<RS,RM,CS,CM> choleskyZ;
  /**inverse matrix of cholesky decomposition of block-matrices matrix <code>Z</code> := <code>X</code>*/
  BlockDenseMatrix<RS,RM,CS,CM> invCholeskyZ;
  /**eigenvalue of <code>mu</code> := <code> X・Z</code> := <code> Y・X</code> */
  BlockVector<RS,RM,CS,CM> xzEigenValues;
  /**minimum eigenvalue of <code>mu</code> := <code> X・Z</code> := <code> Y・X</code> */
  RS xzMinEigenValue;
  /**working dense block-matrices matrix */
  private BlockDenseMatrix<RS,RM,CS,CM> workMat;
  /**working dense block-matrices matrix */
  private BlockDenseMatrix<RS,RM,CS,CM> workMat2;
  /**working vector */
  private BlockVector<RS,RM,CS,CM> workVec;

  /**
   * 新しく生成された<code>Solution</code>オブジェクトを初期化します。
   * 
   * @param unit unit value based on a specified data type
   */
  public Solution(RS unit) {
    this.unit = unit;
    this.xMat = new BlockDenseMatrix<>(unit);
    this.zMat = new BlockDenseMatrix<>(unit);
    this.yVec = new Vector<>(0, unit.create(0));
    this.xMatzMat = new BlockDenseMatrix<>(unit);
    this.choleskyX = new BlockDenseMatrix<>(unit);
    this.choleskyZ = new BlockDenseMatrix<>(unit);
    this.invCholeskyX = new BlockDenseMatrix<>(unit);
    this.invCholeskyZ = new BlockDenseMatrix<>(unit);
    this.xzEigenValues = new BlockVector<>(unit);
    this.workMat = new BlockDenseMatrix<>(unit);
    this.workMat2 = new BlockDenseMatrix<>(unit);
    this.workVec = new BlockVector<>(unit);
  }

  /**
   * constructor of class <code>Solution</code>.<br>
   * 新しく生成された<code>Solution</code>オブジェクトを初期化します。
   * 
   * @param m dimension of vector <code>y</code> := <code>x</code>
   * @param blockSize dimension of integer array, which is <code>blockStruct</code>
   * @param blockStruct structure information of block-matrices matrix
   * @param lambda setting parameter <code>lambda</code>
   */
  public Solution(int m, int blockSize, int[] blockStruct, RS lambda) {
    initialize(m, blockSize, blockStruct, lambda);
  }

  /**
   * initialize matrices and arrays. <br>
   * The more detail of implementing of algorithm, 
   * please refer to the codes of <code>SDPA</code>.<br>  
   * <code>(rsdpa_parts.cpp: SDPA 6.2.0 Line No.: 282)</code>
   * 
   * @param m dimension of vector <code>y</code> := <code>x</code>
   * @param blockSize dimension of integer array, which is <code>blockStruct</code>
   * @param blockStruct structure information of block-matrices matrix
   * @param lambda setting parameter <code>lambda</code>
   */
  public void initialize(int m, int blockSize, int[] blockStruct, RS lambda) {
    this.unit = lambda.createUnit();
    this.xMat = new BlockDenseMatrix<>(blockSize, blockStruct, this.unit);
    this.xMat.setIdentity(lambda);
    this.zMat = new BlockDenseMatrix<>(blockSize, blockStruct, this.unit);
    this.zMat.setIdentity(lambda);
    this.yVec = new Vector<>(m, this.unit.create(0));
    this.xMatzMat = new BlockDenseMatrix<>(blockSize, blockStruct, this.unit);
    this.xMatzMat.setIdentity(lambda.multiply(lambda));
    this.choleskyX = new BlockDenseMatrix<>(blockSize, blockStruct, this.unit);
    this.choleskyX.setIdentity(lambda.sqrt());
    this.choleskyZ = new BlockDenseMatrix<>(blockSize, blockStruct, this.unit);
    this.choleskyZ.setIdentity(lambda.sqrt());
    this.invCholeskyX = new BlockDenseMatrix<>(blockSize, blockStruct, this.unit);
    this.invCholeskyX.setIdentity(this.unit.create(1).divide(lambda.sqrt()));
    this.invCholeskyZ = new BlockDenseMatrix<>(blockSize, blockStruct, this.unit);
    this.invCholeskyZ.setIdentity(this.unit.create(1).divide(lambda.sqrt()));
    this.xzEigenValues = new BlockVector<>(blockSize, blockStruct, this.unit);
    this.workMat = new BlockDenseMatrix<>(blockSize, blockStruct, this.unit);
    this.workMat2 = new BlockDenseMatrix<>(blockSize, blockStruct, this.unit);

    int[] workStruct = new int[blockSize];
    for (int i = 0; i < blockSize; ++i) {
      if (blockStruct[i] > 0) {
        workStruct[i] = 3 * blockStruct[i] - 1;
      } else {
        workStruct[i] = -3 * blockStruct[i] - 1;
      }
    }
    this.workVec = new BlockVector<>(blockSize, workStruct, this.unit);
  }

  /**
   * initialize matrices and arrays with zero by default.
   * 
   * @param m dimension of vector <code>y</code> := <code>x</code>
   * @param blockSize dimension of integer array, which is <code>blockStruct</code>
   * @param blockStruct structure information of block-matrices matrix
   */
  public void initializeZero(int m, int blockSize, int[] blockStruct) {
    this.xMat.initialize(blockSize, blockStruct);
    this.xMat.setZero();
    this.zMat.initialize(blockSize, blockStruct);
    this.zMat.setZero();
    this.yVec.initialize(m);
    this.xMatzMat.initialize(blockSize, blockStruct);
    this.choleskyX.initialize(blockSize, blockStruct);
    this.choleskyZ.initialize(blockSize, blockStruct);
    this.invCholeskyX.initialize(blockSize, blockStruct);
    this.invCholeskyZ.initialize(blockSize, blockStruct);
    this.xzEigenValues.initialize(blockSize, blockStruct);
    this.workMat.initialize(blockSize, blockStruct);
    this.workMat2.initialize(blockSize, blockStruct);
  }

  /**
   * initialize the matrices (<code>X,Z &amp; X・Z</code> := <code>Y,X &amp; Y・X</code>) again.<br>
   * The more detail of implementing of algorithm, 
   * please refer to the codes of <code>SDPA</code>.<br>  
   * <code>(rsdpa_parts.cpp: SDPA 6.2.0 Line No.: 373)</code>
   * 
   * @param blockSize dimension of integer array, which is <code>blockStruct</code>
   * @param blockStruct structure information of block-matrices matrix
   */
  public void initializeResetup(int blockSize, int[] blockStruct) {
    final int[] workStruct = new int[blockSize];
    for (int i = 0; i < blockSize; ++i) {
      if (blockStruct[i] > 0) {
        workStruct[i] = 3 * blockStruct[i] - 1;
      } else {
        workStruct[i] = -3 * blockStruct[i] - 1;
      }
    }
    this.workVec.initialize(blockSize, workStruct);

    BlockDenseMatrix<RS,RM,CS,CM>[] choleskyInverseX = Algebra.getCholeskyAndInv(this.xMat);
    this.choleskyX = choleskyInverseX[0];
    this.invCholeskyX = choleskyInverseX[1];

    BlockDenseMatrix<RS,RM,CS,CM>[] choleskyInverseZ = Algebra.getCholeskyAndInv(this.zMat);
    this.choleskyZ = choleskyInverseZ[0];
    this.invCholeskyZ = choleskyInverseZ[1];

    this.xMatzMat = Algebra.let(this.xMat, '*', this.zMat);
  }

  /**
   * get a judgment of next computing iteration whether can be done or not. 
   * If step length is not too small, true; else, false. The more detail of implementing of algorithm, 
   * please refer to the codes of <code>SDPA</code>.<br>  
   * <code>(rsdpa_parts.cpp: SDPA 6.2.0 Line No.: 464)</code>
   * 
   * @param alpha length of step
   * @param newton information of newton method
   * @return boolean if step length is not too small, true; else, false. 
   */
  public boolean update(StepLength<RS,RM,CS,CM> alpha, Newton<RS,RM,CS,CM> newton) {
    final boolean guarantee = false;
    boolean total_judge = true;

    this.xMat = Algebra.let(this.xMat, '+', newton.getDxMat(), alpha.getPrimal());

    this.yVec = Algebra.let(this.yVec, '+', newton.getDyVec(), alpha.getDual());

    this.zMat = Algebra.let(this.zMat, '+', newton.getDzMat(), alpha.getDual());

    // TODO パラメータとして抽出すべきです。
    final RS minimumStepLength = this.unit.create(10).power(-4);

    if (alpha.getPrimal().isLessThan(minimumStepLength) && alpha.getDual().isLessThan(minimumStepLength)) {
      Tools.message("Step length is too small. "); //$NON-NLS-1$
      return false;
    }

    if (!guarantee) {
      BlockDenseMatrix<RS,RM,CS,CM>[] choleskyInverseX = Algebra.getCholeskyAndInv(this.xMat);
      this.choleskyX = choleskyInverseX[0];
      this.invCholeskyX = choleskyInverseX[1];

      BlockDenseMatrix<RS,RM,CS,CM>[] choleskyInverseZ = Algebra.getCholeskyAndInv(this.zMat);
      this.choleskyZ = choleskyInverseZ[0];
      this.invCholeskyZ = choleskyInverseZ[1];
    }
    this.xMatzMat = Algebra.let(this.xMat, '*', this.zMat);

    this.xzMinEigenValue = this.unit.create(1);

    return total_judge;
  }

  /**
   * update last step information. The more detail of implementing of algorithm, 
   * please refer to the codes of <code>SDPA</code>.<br>
   * <code>(rsdpa_parts.cpp: SDPA 6.2.0 Line No.: 527)</code>
   * 
   */
  public void update_last() {
    this.xzMinEigenValue = this.unit.create(0);
  }

  /**
   * ファイルへ出力
   * 
   * @param ps プリントストリーム
   */
  public void display(PrintStream ps) {
    if (ps == null) {
      return;
    }
    Tools.message(this.toString(), ps);
  }

  /**
   * コンソールへ出力
   */
  public void display() {
    System.out.println(this.toString());
  }

  /**
   * @see java.lang.Object#toString()
   */
  @Override
  public String toString() {
    StringBuffer buff = new StringBuffer();
    buff.append("xMat = \n"); //$NON-NLS-1$
    buff.append(this.xMat.toString());

    buff.append("yVec = \n"); //$NON-NLS-1$
    buff.append(this.yVec.toString());

    buff.append("zMat = \n"); //$NON-NLS-1$
    buff.append(this.zMat.toString());

    return buff.toString();
  }

  /**
   * get a vector <code>y</code> := <code>x</code>.
   * 
   * @return vector <code>y</code> := <code>x</code>
   */
  public Vector<RS,RM,CS,CM> getYvec() {
    return this.yVec;
  }

  /**
   * get a block-matrices matrix <code>X</code>.
   * 
   * @return block-matrices matrix X
   */
  public BlockDenseMatrix<RS,RM,CS,CM> getXmat() {
    return this.xMat;
  }

  /**
   * get a block-matrices matrix <code>Z</code>.
   * 
   * @return block-matrices matrix <code>Z</code> := <code>Y</code>
   */
  public BlockDenseMatrix<RS,RM,CS,CM> getZmat() {
    return this.zMat;
  }

  /**
   * return a copy of solution of <code>SDP</code> problem.
   * 
   * @return a copy
   */
  public Solution<RS,RM,CS,CM> createClone() {
    Solution<RS,RM,CS,CM> inst = new Solution<>(this.unit);
    inst.unit = this.unit == null ? null : this.unit;
    inst.xMat = this.xMat == null ? null : (BlockDenseMatrix<RS,RM,CS,CM>)this.xMat.createClone();
    inst.zMat = this.zMat == null ? null : (BlockDenseMatrix<RS,RM,CS,CM>)this.zMat.createClone();
    inst.yVec = this.yVec == null ? null : (Vector<RS,RM,CS,CM>)this.yVec.createClone();
    inst.xMatzMat = this.xMatzMat == null ? null : (BlockDenseMatrix<RS,RM,CS,CM>)this.xMatzMat.createClone();
    inst.choleskyX = this.choleskyX == null ? null : (BlockDenseMatrix<RS,RM,CS,CM>)this.choleskyX.createClone();
    inst.invCholeskyX = this.invCholeskyX == null ? null : (BlockDenseMatrix<RS,RM,CS,CM>)this.invCholeskyX.createClone();
    inst.choleskyZ = this.choleskyZ == null ? null : (BlockDenseMatrix<RS,RM,CS,CM>)this.choleskyZ.createClone();
    inst.invCholeskyZ = this.invCholeskyZ == null ? null : (BlockDenseMatrix<RS,RM,CS,CM>)this.invCholeskyZ.createClone();
    inst.xzEigenValues = this.xzEigenValues == null ? null : (BlockVector<RS,RM,CS,CM>)this.xzEigenValues.createClone();
    inst.xzMinEigenValue = this.xzMinEigenValue == null ? null : this.xzMinEigenValue;
    inst.workMat = this.workMat == null ? null : (BlockDenseMatrix<RS,RM,CS,CM>)this.workMat.createClone();
    inst.workMat2 = this.workMat2 == null ? null : (BlockDenseMatrix<RS,RM,CS,CM>)this.workMat2.createClone();
    inst.workVec = this.workVec == null ? null : (BlockVector<RS,RM,CS,CM>)this.workVec.createClone();
    return inst;
  }
}