Solver.java

/*
 * Created on 2011/02/21
 * Copyright (C) 2011 Koga Laboratory. All rights reserved.
 *
 */
package org.mklab.sdpj.solver;

import java.io.File;
import java.io.FileOutputStream;
import java.io.IOException;
import java.io.PrintStream;
import java.util.ArrayList;

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.algorithm.AverageComplementarity;
import org.mklab.sdpj.algorithm.DirectionParameter;
import org.mklab.sdpj.algorithm.Lanczos;
import org.mklab.sdpj.algorithm.Newton;
import org.mklab.sdpj.algorithm.Parameter;
import org.mklab.sdpj.algorithm.ParameterType;
import org.mklab.sdpj.algorithm.Phase;
import org.mklab.sdpj.algorithm.PhaseType;
import org.mklab.sdpj.algorithm.RatioInitResCurrentRes;
import org.mklab.sdpj.algorithm.Residuals;
import org.mklab.sdpj.algorithm.Solution;
import org.mklab.sdpj.algorithm.SolveInfo;
import org.mklab.sdpj.algorithm.StepLength;
import org.mklab.sdpj.algorithm.Switch;
import org.mklab.sdpj.algorithm.SwitchType;
import org.mklab.sdpj.algorithm.WHICH_DIRECTION;
import org.mklab.sdpj.datastructure.BlockSparseMatrix;
import org.mklab.sdpj.datastructure.Vector;
import org.mklab.sdpj.iocenter.SolverIO;
import org.mklab.sdpj.iocenter.Vsolution;
import org.mklab.sdpj.problem.Problem;
import org.mklab.sdpj.tool.Tools;


/**
 * Class <code>Solver</code>, which is the main class of <code>SDPJ</code> package, is implemented for solving the SDP problem based on PDIPM algorithm. The PDIPM algorithm is considered the most
 * power and efficiency algorithm of solving SDP problem in polynomial time. If you want to know more details about PDIPM algorithm, please refer to the book, "interior method" wrote by professor
 * kojima of Tokyo Institute of Technology. And the other SDP solver <code>SDPA</code> is developed by the group around professor kojima. Our group referred codes of <code>SDPA</code> for implementing
 * the main class <code>Solver</code> more efficiently and correctly based on multiple precision arthimetic. About information of <code>SDPA</code>, please refer to web site: <a href =
 * "http://sdpa.sourceforge.net/" style=text-decoration:none><tt>http://sdpa.sourceforge.net/</tt></a>.<br>
 * 
 * @author wu
 * @version $Revision$, 2011/02/21
 * @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 Solver<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>> {

  /** a constant parameter, <i>2.2</i> by default */
  private RS kappa;// KAPPA = 2.2 
  /** information of SDP solution */
  private SolveInfo<RS,RM,CS,CM> solutionInformaton;
  /** variable <code>mu := X ・Y</code> */
  private AverageComplementarity<RS,RM,CS,CM> mu;
  /** current primal-dual problem solution */
  private Solution<RS,RM,CS,CM> currentPt;
  /** variable <code>relativeGap := objValPrimal-objValDual </code> */
  private Phase<RS,RM,CS,CM> relativeGap;
  /** temporary current primal-dual problem solution */
  private Solution<RS,RM,CS,CM> tempCurrentPt;
  /** a unit based on a specified data type */
  private RS unit;
  /** SDP problem */
  private Problem<RS,RM,CS,CM> problem;
  /** データファイル */
  private String dataFile;
  /** 初期値のファイルnullでも可（現状だとむしろnullで。） */
  private String initFile;
  /** 出力ファイル */
  private String outFile;
  /** parameter file data */
  private String paraFile;
  /** あればtrue、なければfalse */
  private boolean isInitFile;
  /** 初期値のファイルのデータ形式 sparse:true */
  private boolean isInitSparse;
  /** 問題の入力データの形式 sparse:true */
  private boolean isDataSparse;
  /** パラメータファイルがあるかどうか */
  private boolean isParameter;
  /** type of parameter, refer to <code>{@link org.mklab.sdpj.algorithm.ParameterType <tt>ParameterType</tt>}</code> */
  private ParameterType parameterType;
  /** 表示先 */
  private PrintStream display;
  /** 出力先のストリーム */
  private PrintStream fpOut;
  /**
   * save outputting and inputting information of SDP solution, which is a object of class <code>{@link org.mklab.sdpj.iocenter.SolverIO <tt>SolverIO</tt>}</code>
   */
  private SolverIO<RS,RM,CS,CM> io;
  /**
   * save matrices data by array list,which is a object of class <code>{@link java.util.ArrayList <tt>ArrayList</tt>}</code>
   */
  private ArrayList<String> matrixData = new ArrayList<>(1);
  /** if calling class <code>Vsolution</code> to verify SDP solution, true; else false. */
  private boolean isVsolution;
  /** if show the process of solving SDP problem, true; else false. default is true. */
  private boolean isDisplay = true;
  /**
   * status of SDP solution, which is a object of class <code>{@link org.mklab.sdpj.algorithm.PhaseType <tt>PhaseType</tt>}</code>
   */
  private PhaseType phaseValue;
  /** if call class <code>Vsdpj</code> to compute the lower and upper bounds of SDP solution, true; else false. */
  private boolean isVsdpj;
  /** if call class <code>Solver</code> directly to solve SDP problem, true; else, false. */
  private boolean isSolver = false;
  /** save information of class <code>Vsolution</code> */
  private Vsolution<RS,RM,CS,CM> vsolution;

  /**
   * set SDP problem
   * 
   * @param problem SDP problem
   */
  public void setProblem(Problem<RS,RM,CS,CM> problem) {
    this.problem = problem;
  }

  /**
   * set local dataFile
   * 
   * @param dataFile local dataFile
   */
  public void setDataFile(String dataFile) {
    this.dataFile = dataFile;
  }

  /**
   * set the initial file
   * 
   * @param initFile initial file
   */
  public void setInitFile(String initFile) {
    this.initFile = initFile;
  }

  /**
   * set the name for outputting a file
   * 
   * @param outFile name for outputting a file
   */
  public void setOutFile(String outFile) {
    this.outFile = outFile;
  }

  /**
   * set the parameter file
   * 
   * @param paraFile parameter file
   */
  public void setParaFile(String paraFile) {
    this.paraFile = paraFile;
  }

  /**
   * @param isInitFile if there is a initial file, true.Otherwise, false.
   */
  public void setIsInitFile(boolean isInitFile) {
    this.isInitFile = isInitFile;
  }

  /**
   * @param isInitSparse if initial file is sparse, true. Otherwise, false.
   */
  public void setIsInitSparse(boolean isInitSparse) {
    this.isInitSparse = isInitSparse;
  }

  /**
   * @param isDataSparse if data is sparse, true. Otherwise, false.
   */
  public void setIsDataSparse(boolean isDataSparse) {
    this.isDataSparse = isDataSparse;
  }

  /**
   * @param isParameter if it is parameter, true. Otherwise, false.
   */
  public void setIsParameter(boolean isParameter) {
    this.isParameter = isParameter;
  }

  /**
   * set the type of parameter.
   * 
   * @param parameterType type of parameter.
   */
  public void setParameterType(ParameterType parameterType) {
    this.parameterType = parameterType;
  }

  /**
   * set the instance of PrintStream
   * 
   * @param display instance of PrintStream
   */
  public void setPrintStreamForDisplay(PrintStream display) {
    this.display = display;
  }

  /**
   * set the instance of PrintStream
   * 
   * @param fpOut instance of PrintStream
   */
  public void setPrintStreamForSave(PrintStream fpOut) {
    this.fpOut = fpOut;
  }

  /**
   * 新しく生成された{@link Solver}オブジェクトを初期化します。
   * 
   * @param unit 単位数
   */
  public Solver(RS unit) {
    this.unit = unit;
    //Tools.setUnitNumber(unit);
    this.kappa = unit.create(2).divide(10).add(2);

    this.problem = null;
    this.dataFile = null;
    this.initFile = null;
    this.outFile = "result/others/result.sdpj"; //$NON-NLS-1$
    this.paraFile = null;
    this.isInitFile = false;
    this.isInitSparse = false;
    this.isDataSparse = false;
    this.isParameter = false;
    this.display = System.out;
  }

  /**
   * Kappaを設定します。
   * 
   * @param kappa Kappa
   */
  public void setKappa(RS kappa) {
    this.kappa = kappa;
  }

  /**
   * set <code>mu</code>
   * 
   * @param mu mu
   */
  private void setMu(AverageComplementarity<RS,RM,CS,CM> mu) {
    this.mu = mu;
  }

  /**
   * get <code>mu</code>
   * 
   * @return mu
   */
  public AverageComplementarity<RS,RM,CS,CM> getMu() {
    return this.mu;
  }

  /**
   * set <code>relativeGap</code>
   * 
   * @param relativeGap relativeGap
   */
  private void setRelativeGap(Phase<RS,RM,CS,CM> relativeGap) {
    this.relativeGap = relativeGap;
  }

  /**
   * get <code>relativeGap</code>
   * 
   * @return relativeGap
   */
  public Phase<RS,RM,CS,CM> getRelativeGap() {
    return this.relativeGap;
  }
  /**
   * set <code>solutionInformation</code>
   * 
   * @param solutionInformation information of solution
   */
  public void setSolutionInformation(SolveInfo<RS,RM,CS,CM> solutionInformation) {
    this.solutionInformaton = solutionInformation;
  }

  /**
   * get <code>solutionInformation</code>
   * 
   * @return information of solution
   */
  public SolveInfo<RS,RM,CS,CM> getSolutionInformation() {
    return this.solutionInformaton;
  }

  /**
   * This function is a solver of SDP problems
   * 
   * @return if successful, return true. Otherwise, false.
   * @throws IOException if reading file is failure, IOException is thrown.
   */
  public boolean run() throws IOException {
    if (this.fpOut == null) {
      this.fpOut = new PrintStream(new FileOutputStream(new File(this.outFile)));
    }

    if (this.io == null) {
      this.io = new SolverIO<>(this.fpOut, this.display, this.unit);
      //Set necessary variance of class io
      this.io.setParameterType(this.parameterType);
      this.io.setIsParameter(this.isParameter);
      this.io.setParaFile(this.paraFile);
      this.io.setIsDataSparse(this.isDataSparse);
      this.io.setIsInitFile(this.isInitFile);
      this.io.setInitFile(this.initFile);
      this.io.setIsInitSparse(this.isInitSparse);
      this.io.setMatrixData(this.matrixData);
    }

    this.io.setVsolution(this.isVsolution);
    this.io.setDisplay(this.isDisplay);
    this.io.setVsdpj(this.isVsdpj);

    // Problem<RS,RM,CS,CM> problem = readFile(dataFile, initFile, outFile, fpOut, paraFile, isInitFile, isInitSparse, isDataSparse, isParameter, parameterType, Display, com);

    if (this.problem == null) {
      this.problem = this.io.readFile(this.dataFile);
    }

    // pinpal(problem, fpOut, Display, com, outFile);
    // Problem<RS,RM,CS,CM> problem = prob;

    int m = this.problem.getM();
    int blockSize = this.problem.getBlockSize();
    int[] blockStruct = this.problem.getBlockStruct();
    Vector<RS,RM,CS,CM> b = this.problem.getB();
    BlockSparseMatrix<RS,RM,CS,CM> C = this.problem.getC();
    BlockSparseMatrix<RS,RM,CS,CM>[] A = this.problem.getA();
    //  if(this.isSolver == false) { 
    if (this.currentPt == null) {
      this.currentPt = this.problem.getCurrentPt();
      this.tempCurrentPt = this.currentPt.createClone();
    }

    Solution<RS,RM,CS,CM> initPt = this.problem.getInitPt();
    Parameter<RS,RM,CS,CM> param = this.problem.getParam();
    int nDim = this.problem.getNDim();

    // the end of file read
    // 初期値の計算

    Residuals<RS,RM,CS,CM> initRes = new Residuals<>(m, blockSize, blockStruct, b, C, A, this.currentPt);
    Residuals<RS,RM,CS,CM> currentRes = initRes.createClone();

    Newton<RS,RM,CS,CM> newton = new Newton<>(m, blockSize, blockStruct, this.unit);
    newton.computeFormula(m, A, this.kappa);

    StepLength<RS,RM,CS,CM> alpha = new StepLength<>(this.unit.create(1), this.unit.create(1), blockSize, blockStruct);
    DirectionParameter<RS,RM,CS,CM> beta = new DirectionParameter<>(param.betaStar);
    Switch<RS,RM,CS,CM> reduction = new Switch<>(SwitchType.ON);
    AverageComplementarity<RS,RM,CS,CM> mu2 = new AverageComplementarity<>(param.lambdaStar);
    Lanczos<RS,RM,CS,CM> lanczos = new Lanczos<>();

    if (this.problem.isInitializeInitPt()) {
      mu2 = new AverageComplementarity<>(nDim, initPt);
    }
    RatioInitResCurrentRes<RS,RM,CS,CM> theta = new RatioInitResCurrentRes<>(param, initRes);

    SolveInfo<RS,RM,CS,CM> solveInfo;
    //  if(this.isSolver ==false) { 
    solveInfo = new SolveInfo<>(nDim, b, C, A, initPt, mu2.getInitialValue(), param.omegaStar);
    //   }else {
    //    solveInfo = this.solutionInformaton; 
    //   }

    Phase<RS,RM,CS,CM> phase = new Phase<>(initRes, solveInfo, param, nDim);
    int pIteration = 0;
    //  IO.printHeader(this.fpOut, this.Display);
    this.io.printHeader(this.fpOut, this.display);

    //Starting main loop

    final boolean modifiedAlgorithm = false;

    while (phase.updateCheck(currentRes, solveInfo, param) && pIteration < param.maxIteration) {

      // Mehrotra's Predictor, set variable of Mehrotra
      reduction.MehrotraPredictor(phase);
      beta.MehrotraPredictor(phase, reduction, param);
      boolean isSuccessCholesky;
      isSuccessCholesky = newton.Mehrotra(WHICH_DIRECTION.PREDICTOR, m, A, mu2, beta, phase, this.currentPt, currentRes);
      if (isSuccessCholesky == false) {
        break;
      }
      // the end of Predictor

      if (!modifiedAlgorithm) {
        alpha.MehrotraPredictor(b, C, phase, newton);
        //    alpha.MehrotraCorrector(nDim, b, C, this.currentPt, phase, reduction, newton, mu, theta, lanczos, param, this.com);
      }
      // Mehrotra's Corrector
      beta.MehrotraCorrector(nDim, phase, alpha, this.currentPt, newton, mu2, param);
      newton.Mehrotra(WHICH_DIRECTION.CORRECTOR, m, A, mu2, beta, phase, this.currentPt, currentRes);

      alpha.MehrotraCorrector(nDim, b, C, this.currentPt, phase, reduction, newton, mu2, theta, lanczos, param);
      // the end of Corrector

      // IO.printOneIteration(pIteration, mu, theta, solveInfo, alpha, beta, currentRes, this.fpOut, this.Display);

      this.io.setPIteration(pIteration);
      this.io.setMu(mu2);
      this.io.setTheta(theta);
      this.io.setSolveInfo(solveInfo);
      this.io.setAlpha(alpha);
      this.io.setBeta(beta);
      this.io.setCurrentRes(currentRes);
      this.io.printOneIteration();

      // Compute condition number, singular value of matrix B  
      // newton.getBMat().getConditionNumber(condBMat, minBMat, maxBMat);

      // if step length is too short, we finish algorithm
      if (this.currentPt.update(alpha, newton) == false) {
        Tools.message("cannot move"); //$NON-NLS-1$
        pIteration++;
        break;
      }
      //   theta.update(reduction, alpha);
      mu2.update(nDim, this.currentPt);
      currentRes.update(m, b, C, A, this.currentPt, mu2);

      theta.update_exact(initRes, currentRes);
      solveInfo.update(nDim, b, C, initPt, this.currentPt, currentRes, mu2, theta, param);
      pIteration++;
    } // end of MAIN_LOOP

    this.currentPt.update_last();
    currentRes.compute(m, blockSize, blockStruct, b, C, A, this.currentPt);
    if (Tools.REVERSE_PRIMAL_DUAL == 1) {
      phase.reverse();
    }

    this.io.setPIteration(pIteration);
    //if(this.isSolver == false) { 
    this.io.setMu(mu2);
    //   }
    this.io.setTheta(theta);
    this.io.setSolveInfo(solveInfo);
    this.io.setAlpha(alpha);
    this.io.setBeta(beta);
    this.io.setCurrentRes(currentRes);
    this.io.printOneIteration();

    this.io.setBlockStruct(blockStruct);
    this.io.setPhase(phase);
    this.io.setCurrentPt(this.currentPt);
    this.io.setNDim(nDim);
    this.io.setM(this.problem.getM());
    this.io.setB(b);
    this.io.setC(C);
    this.io.setA(A);
    this.io.setParam(param);
    this.io.setOutFile(this.outFile);
    this.io.printLastInfo();

    //ready for the running of Vsdpj
    if (this.isVsdpj == true) {
      this.phaseValue = this.io.getPhase().getValue();
      this.vsolution = this.io.getVsolution();
    }
    
    Phase<RS,RM,CS,CM> rGap = new Phase<>(currentRes, solveInfo, param, nDim);
    this.setRelativeGap(rGap);
    
    // Compute condition number, singular value of matrix Y,X
    //    this.currentPt.zMat.getCondtionNumber(condzMat, minzMat, maxzMat);
    //    this.currentPt.xMat.getCondtionNumber(condxMat, minxMat, maxxMat);

    this.setSolutionInformation(solveInfo);
    this.setMu(mu2);

    //    Tools.message("  main loop time = " + this.com.getMainLoopTime(), this.Display); //$NON-NLS-1$
    //    Tools.message("    main loop time = " + this.com.getMainLoopTime(), this.fpOut); //$NON-NLS-1$
    //    Tools.message("      total time = " + this.com.getTotalTime(), this.Display); //$NON-NLS-1$
    //    Tools.message("        total time = " + this.com.getTotalTime(), this.fpOut); //$NON-NLS-1$
    //    Tools.message("file   read time = " + this.com.getFileReadTime() + "\n", this.Display); //$NON-NLS-1$ //$NON-NLS-2$
    //    Tools.message("  file   read time = " + this.com.getFileReadTime(), this.fpOut); //$NON-NLS-1$
    if (Sdpj.getElapsedTime() == 0.0 && this.isDisplay == true) {
      Tools.message("\n"); //$NON-NLS-1$
    }
    this.fpOut.close();
    //  SDPJ.unit = null;

    return true;
  }

  /**
   * 解を返します。
   * 
   * @return 解
   */
  public Solution<RS,RM,CS,CM> getSolution() {
    return this.currentPt;
  }

  /**
   * Return io.
   * 
   * @return io
   */
  public SolverIO<RS,RM,CS,CM> getIo() {
    return this.io;
  }

  /**
   * Set io.
   * 
   * @param io io
   */
  public void setIo(SolverIO<RS,RM,CS,CM> io) {
    this.io = io;
  }

  /**
   * set <code>currentPt</code>
   * 
   * @param currentPt reset current solution.
   */
  public void setCurrentPt(Solution<RS,RM,CS,CM> currentPt) {
    this.currentPt = currentPt;
  }

  /**
   * set <code>vsolution</code>
   * 
   * @param vsolution vsolutionを設定します。
   */
  public void setVsolution(boolean vsolution) {
    this.isVsolution = vsolution;
  }

  /**
   * get <code>isVsolution</code>
   * 
   * @return vsolutionを返します。
   */
  public boolean isVsolution() {
    return this.isVsolution;
  }

  /**
   * set <code>isDisplay</code>
   * 
   * @param isDisplay isDisplayを設定します。
   */
  public void setDisplay(boolean isDisplay) {
    this.isDisplay = isDisplay;
  }

  /**
   * get <code>isDisplay</code>
   * 
   * @return isDisplayを返します。
   */
  public boolean isDisplay() {
    return this.isDisplay;
  }

  /**
   * set <code>phaseValue</code>
   * 
   * @param phaseValue phaseValueを設定します。
   */
  public void setPhaseValue(PhaseType phaseValue) {
    this.phaseValue = phaseValue;
  }

  /**
   * get <code>phaseValue</code>
   * 
   * @return phaseValueを返します。
   */
  public PhaseType getPhaseValue() {
    return this.phaseValue;
  }

  /**
   * set <code>isVsdpj</code>
   * 
   * @param isVsdpj isVsdpjを設定します。
   */
  public void setVsdpj(boolean isVsdpj) {
    this.isVsdpj = isVsdpj;
  }

  /**
   * get <code>isVsdpj</code>
   * 
   * @return isVsdpjを返します。
   */
  public boolean isVsdpj() {
    return this.isVsdpj;
  }

  /**
   * set <code>vsolution</code>
   * 
   * @param vsolution vsolutionを設定します。
   */
  public void setVsolution(Vsolution<RS,RM,CS,CM> vsolution) {
    this.vsolution = vsolution;
  }

  /**
   * get <code>vsolution</code>
   * 
   * @return vsolutionを返します。
   */
  public Vsolution<RS,RM,CS,CM> getVsolution() {
    return this.vsolution;
  }

  /**
   * get <code>problem</code>
   * 
   * @return information of SDP problem
   */
  public Problem<RS,RM,CS,CM> getProblem() {
    return this.problem;
  }

  /**
   * set <code>isSolver</code>
   * 
   * @param isSolver isSolverを設定します。
   */
  public void setSolver(boolean isSolver) {
    this.isSolver = isSolver;
  }

  /**
   * get <code>isSolver</code>
   * 
   * @return isSolverを返します。
   */
  public boolean isSolver() {
    return this.isSolver;
  }

  /**
   * set <code>tempCurrentPt</code>
   * 
   * @param tempCurrentPt tempCurrentPtを設定します。
   */
  public void setTempCurrentPt(Solution<RS,RM,CS,CM> tempCurrentPt) {
    this.tempCurrentPt = tempCurrentPt;
  }

  /**
   * get <code>tempCurrentPt</code>
   * 
   * @return tempCurrentPtを返します。
   */
  public Solution<RS,RM,CS,CM> getTempCurrentPt() {
    return this.tempCurrentPt;
  }
  
  
  /**
   * @return fpOutを返します。
   */
  public PrintStream getFpOut(){
    return this.fpOut;
  }

}