ProblemCreater.java

package org.mklab.sdpj.problem;

import java.io.BufferedReader;
import java.io.FileNotFoundException;
import java.io.FileReader;
import java.io.IOException;

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.Parameter;
import org.mklab.sdpj.algorithm.Solution;
import org.mklab.sdpj.datastructure.BlockSparseMatrix;
import org.mklab.sdpj.datastructure.SparseDenseDiagonal;
import org.mklab.sdpj.datastructure.SparseMatrix;
import org.mklab.sdpj.datastructure.Vector;
import org.mklab.sdpj.iocenter.IO;


/**
 * class <code>ProblemCreater</code> is defined for build a SDP problem, 
 * which can be computed correctly by solver <code>SDPJ</code>. With class 
 * <code>ProblemCreaterOnMemory</code>, it can concern to design how to deal 
 * with inputting data.
 * 
 * @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 ProblemCreater<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>> {
  /** データの入力形式がスパースかどうか */
  private boolean isDataSparse = false;
  /**a setting parameter of computing*/
  private Parameter<RS,RM,CS,CM> parameter;
  /**dimension of vector <code>x</code>*/
  private int m;
  /**the number of block of block-matrices array*/
  private int nBlock;
  /**a integer array save the block information of block-matrices array*/
  private int[] blockStruct;
  /**vector <code>b</code> := <code>c</code>*/
  private Vector<RS,RM,CS,CM> b;
  /**sparse block-matrices array <code>C</code> := <code>F0</code>*/
  private BlockSparseMatrix<RS,RM,CS,CM> C;
  /**sparse block-matrices array <code>A</code> := <code>Fi</code>*/
  private BlockSparseMatrix<RS,RM,CS,CM>[] A;
  /**the number of block matrix <code>C</code> := <code>F0</code>*/
  private int[] CNonZeroCount;
  /**the number of block-matrices array <code>A</code> := <code>Fi</code>*/
  private int[] ANonZeroCount;
  /**initial priaml-dual solution*/
  private Solution<RS,RM,CS,CM> initPt;
  /**current primal-dual solution*/
  private Solution<RS,RM,CS,CM> currentPt;
  /**block-array matrices <code>Fi</code>*/
  private String[][] F;
  /**vector <code>x</code>*/
  private String[] vec;
  
  /** */
  private RS unit;

//  /**
//   * Initialize the generated object of {@link ProblemCreater}.
//   */
//  public ProblemCreater() {
//    //
//  }

  /**
   * create a SDP problem with inputting a dense data file.
   * <br>入力データがスパースの形式ではない場合(***.dat)
   * 
   * @param F block-matrices array <code>Fi</code>
   * @param vec vector <code>x</code>
   * @param m dimension of vector <code>x</code>
   * @param nBlock the number of block of block-matrices array
   * @param blockStruct a integer array save the block information of block-matrices array
   * @param parameter a setting parameter of computing
   * @param unit unit
   */
  private ProblemCreater(String[][] F, String[] vec, int m, int nBlock, int[] blockStruct, Parameter<RS,RM,CS,CM> parameter, RS unit) {
    this.unit = unit;
    this.setDataSparse(false);
    this.F = F;
    this.setVec(vec);
    this.m = m;
    this.nBlock = nBlock;
    this.blockStruct = blockStruct;
    this.parameter = parameter;
    this.CNonZeroCount = new int[nBlock];
    this.ANonZeroCount = new int[nBlock * m];

    this.A = new BlockSparseMatrix[m];
    for (int i = 0; i < m; i++) {
      this.A[i] = new BlockSparseMatrix<>();
    }
    this.b = new Vector<>(m, parameter.getUnit());
    this.C = new BlockSparseMatrix<>(this.nBlock, this.blockStruct);
  }

  /**
   * count the number of nonzero elements of matrix.
   */
  private void countNonzeroNumber() {
    for (int i = 0; i < this.nBlock; i++) {
      this.CNonZeroCount[i] = 0;
      for (int j = 0; j < this.m; j++) {
        this.ANonZeroCount[j * this.nBlock + i] = 0;
      }
    }

    //final RS unit = Tools.getUnitNumber();

    int counter1 = 0;
    for (int k = 0; k < this.nBlock; ++k) {
      int size = this.blockStruct[k];
      if (size > 0) {
        for (int i = 0; i < size; ++i) {
          for (int j = 0; j < size; ++j) {
            RS tmp = this.unit.valueOf(this.F[0][counter1]);
            counter1++;
            if (i <= j && !tmp.isZero()) {
              this.CNonZeroCount[k]++;
            }
          }
        }
      } else {
        for (int j = 0; j < -size; ++j) {
          //E tmp = unit.valueOf(this.F[0][counter1]);
          counter1++;
          this.CNonZeroCount[k]++;
        }
      }
    }

    for (int k = 0; k < this.m; ++k) {
      int counter2 = 0;
      for (int l = 0; l < this.nBlock; ++l) {
        int size = this.blockStruct[l];
        if (size > 0) {
          for (int i = 0; i < size; ++i) {
            for (int j = 0; j < size; ++j) {
              RS tmp = this.unit.valueOf(this.F[k + 1][counter2]);
              counter2++;
              if (i <= j && !tmp.isZero()) {
                this.ANonZeroCount[k * this.nBlock + l]++;
              }
            }
          }
        } else {
          for (int j = 0; j < -size; ++j) {
            //E tmp = unit.valueOf(this.F[k + 1][counter2]);
            counter2++;
            this.ANonZeroCount[k * this.nBlock + l]++;
          }
        }
      }
    }
  }

  /**
   * initialize matrix <code>C</code> and <code>A</code>
   * <br>CとAを初期化します
   * <blockquote><tt>
   * Note:<br>
   *  matrix <code>C</code>  := matrix <code>F0</code><br>
   *  matrix <code>A</code>  := matrix <code>Fi</code><br>
   * </tt></blockquote>
   * 
   */
  private void initializeCandA() {
    for (int i = 0; i < this.nBlock; ++i) {
      int size = this.blockStruct[i];
      SparseMatrix<RS,RM,CS,CM> block;
      if (size > 0) {
        block = new SparseMatrix<>(size, size, SparseDenseDiagonal.SPARSE, this.CNonZeroCount[i], this.unit);
      } else {
        block = new SparseMatrix<>(-size, -size, SparseDenseDiagonal.DIAGONAL, -size, this.unit);
      }
      this.C.setBlock(i, block);
    }
    for (int k = 0; k < this.m; ++k) {
      this.A[k] = new BlockSparseMatrix<>(this.nBlock, this.blockStruct);
      for (int i = 0; i < this.nBlock; ++i) {
        int size = this.blockStruct[i];
        SparseMatrix<RS,RM,CS,CM> block;
        if (size > 0) {
          block = new SparseMatrix<>(size, size, SparseDenseDiagonal.SPARSE, this.ANonZeroCount[k * this.nBlock + i], this.unit);
        } else {
          block = new SparseMatrix<>(-size, -size, SparseDenseDiagonal.DIAGONAL, -size, this.unit);
        }
        this.A[k].setBlock(i, block);
      }
    }
  }

  /**
   * input dense data.
   * <br>スパース形式ではない場合のデータ入力
   * 
   * @param F block-matrices array <code>Fi</code>
   */
  private void setDataCandAWithNotSparse(String[][] F) {
    //final RS unit = Tools.getUnitNumber();

    for (int i = 0; i < this.nBlock; ++i) {
      SparseMatrix<RS,RM,CS,CM> block = this.C.getBlock(i);
      if (block.isSparce()) {
        block.nonZeroCount = 0;
        block.nonZeroEffect = 0;
      }
    }
    for (int k = 0; k < this.m; ++k) {
      for (int i = 0; i < this.nBlock; ++i) {
        SparseMatrix<RS,RM,CS,CM> block = this.A[k].getBlock(i);
        if (block.isSparce()) {
          block.nonZeroCount = 0;
          block.nonZeroEffect = 0;
        }
      }
    }

    int counter = 0;
    for (int k = 0; k < this.nBlock; ++k) {
      int size = this.blockStruct[k];
      if (size > 0) {
        for (int i = 0; i < size; ++i) {
          for (int j = 0; j < size; ++j) {
            RS tmp = this.unit.valueOf(F[0][counter]);
            counter++;
            if (i <= j && !tmp.isZero()) {
              SparseMatrix<RS,RM,CS,CM> target = this.C.getBlock(k);
              int count = target.nonZeroCount;
              target.rowIndex[count] = i;
              target.columnIndex[count] = j;
              // C must be opposite sign.
              target.sparseElements[count] = tmp.unaryMinus();
              target.nonZeroCount++;
              if (i == j) {
                target.nonZeroEffect++;
              } else {
                target.nonZeroEffect += 2;
              }
            }
          }
        }
      } else {
        for (int j = 0; j < -size; ++j) {
          RS tmp = this.unit.valueOf(F[0][counter]);
          counter++;
          // C must be opposite sign.
          this.C.getBlock(k).diagonalElements[j] = tmp.unaryMinus();
        }
      }
    }

    for (int k = 0; k < this.m; ++k) {
      counter = 0;
      for (int p = 0; p < this.nBlock; ++p) {
        int size = this.blockStruct[p];
        if (size > 0) {
          for (int i = 0; i < size; ++i) {
            for (int j = 0; j < size; ++j) {
              RS tmp = this.unit.valueOf(F[k + 1][counter]);
              counter++;
              if (i <= j && !tmp.isZero()) {
                SparseMatrix<RS,RM,CS,CM> target = this.A[k].getBlock(p);
                int count = target.nonZeroCount;
                target.rowIndex[count] = i;
                target.columnIndex[count] = j;
                target.sparseElements[count] = tmp;
                target.nonZeroCount++;
                if (i == j) {
                  target.nonZeroEffect++;
                } else {
                  target.nonZeroEffect += 2;
                }
              }
            }
          }
        } else {
          for (int j = 0; j < -size; ++j) {
            RS tmp = this.unit.valueOf(F[k + 1][counter]);
            counter++;
            this.A[k].getBlock(p).diagonalElements[j] = tmp;
          }
        }
      }
    }
  }

  /**
   * change sparse matrix into dense matrix.
   */
  private void changeToDense() {
    // if possible , change C and A to Dense
    this.C.changeToDense();
    for (int i = 0; i < this.m; ++i) {
      this.A[i].changeToDense();
    }
  }

  /**
   * get a SDP problem that can be computed correctly by solver.
   * 
   * @return a SDP problem
   */
  private Problem<RS,RM,CS,CM> getProblem() {
    return new Problem<>(this.m, this.nBlock, this.blockStruct, this.A, this.b, this.C, this.parameter, this.currentPt, this.initPt);
  }

  /**
   * create a SDP problem with class <code>ProblemCreater</code>.
   * <br>SDPAに渡す問題を生成します。
   * 
 * @param <RS> type of real scalar
 * @param <RM> type of real matrix
 * @param <CS> type of complex scalar
 * @param <CM> type of complex Matrix
   * @param F block-matrices array <code>Fi</code>
   * @param vec vector <code>x</code>
   * @param m dimension of vector <code>x</code>
   * @param nBlock the number of block of block-matrices array
   * @param blockStruct a integer array save the block information of block-matrices array
   * @param parameter パラメータ, a setting parameter of computing
   * @return 問題, a SDP problem can be computed by solver
   */
  public static <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>> Problem<RS,RM,CS,CM> create(String[][] F, String[] vec, int m, int nBlock, int[] blockStruct, Parameter<RS,RM,CS,CM> parameter) {
    final RS unit = parameter.getUnit();
    //Tools.setUnitNumber(unit);

    ProblemCreater<RS,RM,CS,CM> creater = new ProblemCreater<>(F, vec, m, nBlock, blockStruct, parameter, unit);

    for (int i = 0; i < creater.m; i++) {
      creater.b.setElement(i, unit.valueOf(vec[i]));
    }

    creater.countNonzeroNumber();
    creater.initializeCandA();
    creater.setDataCandAWithNotSparse(F);
    creater.changeToDense();
    return creater.getProblem();
  }

  /**
   * make a SDP problem from inputting a sparse data file.
   * <br>Sparseデータ形式から読み込み、問題を作成します
   * 
 * @param <RS> type of real scalar
 * @param <RM> type of real matrix
 * @param <CS> type of complex scalar
 * @param <CM> type of complex Matrix
   * @param dataFileName データファイル名, name of inputting data file
   * @param isSparse if sparse, true; else, false
   * @param unit a unit based on generic data type
   * @return 問題, a SDP problem can be computed by solver
   */
  public static <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>> Problem<RS,RM,CS,CM> createFromFile(String dataFileName, boolean isSparse, RS unit) {
    try (BufferedReader reader = new BufferedReader(new FileReader(dataFileName))) {

      int m = IO.readInteger(reader, System.out);
      int nBlock = IO.readInteger(reader.readLine());
      int[] blockStruct = IO.readBlockStruct(reader, nBlock);
      String[] vec = IO.readVector(reader, m);

      String[][] F;
      if (isSparse) {
        F = IO.readSparseMatrix(reader, m, blockStruct);
      } else {
        F = IO.readDenseMatrix(reader, m, blockStruct);
      }
      
      Parameter<RS,RM,CS,CM> param = new Parameter<>(unit);
      param.setDefaultParameter();

      Problem<RS,RM,CS,CM> problem = create(F, vec, m, nBlock, blockStruct, param);
      return problem;
    } catch (FileNotFoundException e) {
      throw new RuntimeException(e);
    } catch (IOException e) {
      throw new RuntimeException(e);
    }
  }

  /**
   * @param vec vecを設定します。
   */
  public void setVec(String[] vec) {
    this.vec = vec;
  }

  /**
   * @return vecを返します。
   */
  public String[] getVec() {
    return this.vec;
  }

  /**
   * @param isDataSparse isDataSparseを設定します。
   */
  public void setDataSparse(boolean isDataSparse) {
    this.isDataSparse = isDataSparse;
  }

  /**
   * @return isDataSparseを返します。
   */
  public boolean isDataSparse() {
    return this.isDataSparse;
  }
}