Sedumi2.java

/*
 * Created on 2010/02/12
 * Copyright (C) 2010 Koga Laboratory. All rights reserved.
 *
 */
package org.mklab.sdpj.tool;

import java.io.IOException;
import java.util.List;

import org.mklab.matj.jparser.statement.symbol.JRealMatrixVariable;
import org.mklab.matj.parser.statement.Expression;
import org.mklab.mpf.ap.MPFloat;
import org.mklab.mpf.ap.MPFloatComplexMatrix;
import org.mklab.mpf.ap.MPFloatComplexNumber;
import org.mklab.mpf.ap.MPFloatMatrix;
import org.mklab.nfc.matrix.ComplexNumericalMatrix;
import org.mklab.nfc.matrix.DoubleComplexMatrix;
import org.mklab.nfc.matrix.DoubleMatrix;
import org.mklab.nfc.matrix.Matrix;
import org.mklab.nfc.matrix.RealNumericalMatrix;
import org.mklab.nfc.matx.MatxList;
import org.mklab.nfc.scalar.ComplexNumericalScalar;
import org.mklab.nfc.scalar.DoubleComplexNumber;
import org.mklab.nfc.scalar.DoubleNumber;
import org.mklab.nfc.scalar.RealNumericalScalar;
import org.mklab.sdpj.algorithm.Parameter;
import org.mklab.sdpj.algorithm.ParameterType;
import org.mklab.sdpj.algorithm.Solution;
import org.mklab.sdpj.convert.K;
import org.mklab.sdpj.convert.SDPAM;
import org.mklab.sdpj.convert.Sedumi;
import org.mklab.sdpj.datastructure.DenseMatrix;
import org.mklab.sdpj.datastructure.Vector;
import org.mklab.sdpj.problem.Problem;
import org.mklab.sdpj.problem.ProblemCreater;
import org.mklab.sdpj.solver.Solver;


/**
 * Class <code>Sedumi2</code> can change a sudumi-style inputting SDP problem into <code>SDPJ</code>'s inputting style 
 * and solve the SDP problem with highly accurate based on multiple precision arithmetic simultaneously.
 * The purpose of designing class <code>Sedumi2</code> is that SDPJ can directly solve a sedumi-style inputting SDP 
 * problem. If you want to know more details about <code>Sedumi</code>, which is another SDP solver, please refer to 
 * web site:<a href = "http://sedumi.ie.lehigh.edu/" style=text-decoration:none>
 * <tt>http://sedumi.ie.lehigh.edu/</tt></a>.<br>  
 * 
 * @author Hiroaki Matsuo
 * @version $Revision$, 2010/02/15
 */
public class Sedumi2 {

  /**
   * Solve a sdpam-style SDP problem based on <code>Sdpj</code>.
   * <br>sdpam形式からsdpjをときます。bLOCKsTRUCTの値が -の場合には対応していません。
   * 
 * @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 mDIM dimension of vector
   * @param nBLOCK the number of blocks
   * @param bLOCKsTRUCT0 information of block structure
   * @param c0 a constant matrix <code>c0</code>
   * @param sdpaF a constant matrix <code>Fi</code>
   * @return a SDP solution of vector <code>x</code>, matrix <code>Y</code> and matrix <code>X</code> in order 
   * @throws IOException If I/O exception occurs
   */
  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>> MatxList sdpa(int mDIM, int nBLOCK, Matrix bLOCKsTRUCT0, Matrix c0, MatxList sdpaF) throws IOException {
    //Tools.setUnitNumber(new MPFloat(1));
    int m = mDIM;
    int nBlock = nBLOCK;
    DoubleMatrix bLOCKsTRUCT = (DoubleMatrix)bLOCKsTRUCT0;
    DoubleMatrix c = (DoubleMatrix)c0;
    int[] blockStruct = new int[bLOCKsTRUCT.length()];
    for (int i = 0; i < bLOCKsTRUCT.length(); i++) {
      blockStruct[i] = (int)bLOCKsTRUCT.getElement(i + 1).doubleValue();
    }

    String[] vec = new String[c.length()];
    for (int i = 0; i < c.length(); i++) {
      vec[i] = String.valueOf(c.getElement(i + 1));
    }

    int Fsize = 0;
    for (int i = 0; i < blockStruct.length; i++) {
      Fsize += blockStruct[i] * blockStruct[i];
    }

    String[][] F = new String[sdpaF.size()][Fsize];

    final List<Object> sdpaFList = sdpaF.toList();
    for (int i = 0; i < sdpaF.size(); i++) {
      int start = 1;
      MatxList sdpaFSubList = (MatxList)sdpaFList.get(i);

      for (int j = 0; j < blockStruct.length; j++) {
        int theend = start + blockStruct[j] * blockStruct[j] - 1;
        int n = start;
        final List<Object> sdpaFSubsubList = sdpaFSubList.toList();
        Expression expression = (Expression)sdpaFSubsubList.get(j);
        JRealMatrixVariable popValue = (JRealMatrixVariable)expression.pop();
        DoubleMatrix matrix = (DoubleMatrix)popValue.getValue().getValue();
        for (int k = 0; k < blockStruct[j] * blockStruct[j]; k++) {
          F[i][n - 1] = String.valueOf(matrix.getElement(k + 1)).replace("E", "e"); //$NON-NLS-1$ //$NON-NLS-2$
          n++;
        }
        start = theend + 1;
      }
    }

    Parameter<MPFloat,MPFloatMatrix,MPFloatComplexNumber,MPFloatComplexMatrix> parameter = new Parameter<>(new MPFloat(1));
    parameter.setDefaultParameter(ParameterType.PARAMETER_STABLE);
    Problem<MPFloat,MPFloatMatrix,MPFloatComplexNumber,MPFloatComplexMatrix> problem = ProblemCreater.create(F, vec, m, nBlock, blockStruct, parameter);
    Solver<MPFloat,MPFloatMatrix,MPFloatComplexNumber,MPFloatComplexMatrix> solver = new Solver<>(new MPFloat(1));

    solver.setProblem(problem);
    solver.setOutFile("out.txt"); //$NON-NLS-1$
    solver.setPrintStreamForDisplay(System.out);
    solver.run();
    //   solver.solve(problem, "out.txt", System.out, com); //$NON-NLS-1$
    Solution<RS,RM,CS,CM> s = (Solution<RS, RM, CS, CM>)solver.getSolution();

    MatxList ml = new MatxList();
    Vector<RS,RM,CS,CM> y = s.getYvec();
    ml.add(y.getElement(0).createGrid(y.getElements()));
    DenseMatrix<RS,RM,CS,CM> x = s.getXmat().getBlock(0);
    if (x.isDense()) {
      ml.add(x.denseElements[0].createGrid(x.denseElements));
    } else {
      ml.add(x.diagonalElements[0].createGrid(x.diagonalElements));
    }
    DenseMatrix<RS,RM,CS,CM> z = s.getZmat().getBlock(0);
    if (z.isDense()) {
      ml.add(z.denseElements[0].createGrid(z.denseElements));
    } else {
      ml.add(z.diagonalElements[0].createGrid(z.diagonalElements));
    }

    return ml;
  }

  /**
   * Solve a sedumi-style <code>SDP</code> problem based on <code>Sdpj</code>.
   * <br>sedumi形式からSDPJを起動する関数です。
   * 
 * @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 A0 a constant matrix A0.
   * @param b0 a constant matrix b0.
   * @param c0 a constant matrix c0.
   * @param k0 a constant matrix k0.
   * @return a matrix list of matrix X,Y and Z in order.
   * @throws IOException If I/O exception occurs.
   */
  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>> MatxList sedumi(Matrix A0, Matrix b0, Matrix c0, Matrix k0) throws IOException {
    //Tools.setUnitNumber(new MPFloat(1));
    DoubleMatrix A = (DoubleMatrix)A0;
    DoubleMatrix b = (DoubleMatrix)b0;
    DoubleMatrix c = (DoubleMatrix)c0;
    DoubleMatrix k = (DoubleMatrix)k0;

    int[] k1 = new int[k.length()];
    for (int i = 0; i < k.length(); i++) {
      k1[i] = (int)k.getElement(i + 1).doubleValue();
    }

    Sedumi conv = new Sedumi();
    SDPAM sdpam = conv.convertToSedumi(A, b, c, new K(k1));
    int m = sdpam.getmDIM();
    int nBlock = sdpam.getnBLOCK();

    //List<Integer> bs = new ArrayList<Integer>();
    List<Integer> bs = sdpam.getbLOCKsTRUCT();
    int[] blockStruct = new int[bs.size()];
    for (int i = 0; i < bs.size(); i++) {
      blockStruct[i] = bs.get(i).intValue();
    }

    double[] C = sdpam.getC();
    String[] vec = new String[C.length];
    for (int i = 0; i < C.length; i++) {
      vec[i] = String.valueOf(C[i]);
    }

    int Fsize = 0;
    for (int i = 0; i < blockStruct.length; i++) {
      if (blockStruct[i] > 0) {
        Fsize = Fsize + blockStruct[i] * blockStruct[i];
      } else {
        Fsize = Fsize - blockStruct[i];
      }
    }

    DoubleMatrix[][] sdpamF = sdpam.getSdpamF();
    String[][] F = new String[sdpamF[0].length][Fsize];

    for (int i = 0; i < sdpamF[0].length; i++) {
      int count = 0;
      for (int j = 0; j < sdpamF.length; j++) {
        DoubleMatrix subF = sdpamF[j][i];
        for (int row = 1; row <= subF.getRowSize(); row++) {
          for (int column = 1; column <= subF.getColumnSize(); column++) {
            F[i][count] = String.valueOf(subF.getElement(row, column)).replace("E", "e"); //$NON-NLS-1$ //$NON-NLS-2$
            count++;
          }
        }
      }
    }

    Parameter<DoubleNumber,DoubleMatrix,DoubleComplexNumber,DoubleComplexMatrix> parameter = new Parameter<>(new DoubleNumber(1));
    parameter.setDefaultParameter();
    Problem<DoubleNumber,DoubleMatrix,DoubleComplexNumber,DoubleComplexMatrix> problem = ProblemCreater.create(F, vec, m, nBlock, blockStruct, parameter);
    Solver<DoubleNumber,DoubleMatrix,DoubleComplexNumber,DoubleComplexMatrix> solver = new Solver<>(new DoubleNumber(1));

    solver.setProblem(problem);
    solver.setOutFile("out.txt"); //$NON-NLS-1$
    solver.setPrintStreamForDisplay(System.out);
    solver.run();
    //    solver.solve(problem, "out.txt", System.out, com); //$NON-NLS-1$
    Solution<RS,RM,CS,CM> s = (Solution<RS, RM, CS, CM>)solver.getSolution();

    MatxList ml = new MatxList();
    Vector<RS,RM,CS,CM> y = s.getYvec();
    ml.add(y.getElement(0).createGrid(y.getElements()));
    DenseMatrix<RS,RM,CS,CM> x = s.getXmat().getBlock(0);
    if (x.isDense()) {
      ml.add(x.denseElements[0].createGrid(x.denseElements));
    } else {
      ml.add(x.diagonalElements[0].createGrid(x.diagonalElements));
    }
    DenseMatrix<RS,RM,CS,CM> z = s.getZmat().getBlock(0);
    if (z.isDense()) {
      ml.add(z.denseElements[0].createGrid(z.denseElements));
    } else {
      ml.add(z.diagonalElements[0].createGrid(z.diagonalElements));
    }

    return ml;
  }
}