RealVULS.java

/*
 * Created on 2012/10/09
 * Copyright (C) 2012 Koga Laboratory. All rights reserved.
 *
 */
package org.mklab.cga.vuls;

import org.mklab.cga.interval.matrix.IntervalComplexNumericalMatrix;
import org.mklab.cga.interval.matrix.IntervalRealNumericalMatrix;
import org.mklab.cga.interval.scalar.IntervalComplexNumericalScalar;
import org.mklab.cga.interval.scalar.IntervalRealNumericalScalar;
import org.mklab.cga.linear.LinearEquationVerifier;
import org.mklab.cga.linear.RealIntlabMethod;
import org.mklab.cga.round.FPURoundModeSelector;
import org.mklab.nfc.leq.LUDecomposition;
import org.mklab.nfc.matrix.ComplexNumericalMatrix;
import org.mklab.nfc.matrix.IntMatrix;
import org.mklab.nfc.matrix.RealNumericalMatrix;
import org.mklab.nfc.scalar.ComplexNumericalScalar;
import org.mklab.nfc.scalar.RealNumericalScalar;
import org.mklab.nfc.util.RoundModeManager;


/**
 * VULS means the Verification for Underdetermined Linear Systems. For details see C. Jansson, Rigorous Lower and Upper Bounds in Linear Programming, SIAM J. OPTIM. Vol.14, No.3, pp. 914-935.
 * 
 * @author motoyama
 * @version $Revision$, 2012/10/09
 * @param <RIS> 実区間スカラーの型
 * @param <RIM> 実区間行列の型
 * @param <CIS> 複素区間スカラーの型
 * @param <CIM> 複素区間行列の型
 * @param <RS> 実スカラーの型
 * @param <RM> 実行列の型
 * @param <CS> 複素スカラーの型
 * @param <CM> 複素行列の型
 */
public class RealVULS<RIS extends IntervalRealNumericalScalar<RIS,RIM,CIS,CIM,RS,RM,CS,CM>, RIM extends IntervalRealNumericalMatrix<RIS,RIM,CIS,CIM,RS,RM,CS,CM>, CIS extends IntervalComplexNumericalScalar<RIS,RIM,CIS,CIM,RS,RM,CS,CM>, CIM extends IntervalComplexNumericalMatrix<RIS,RIM,CIS,CIM,RS,RM,CS,CM>, 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 number of matrix
   */
  private int m;
  /**
   * The constant <code>unit</code> based on a general data type
   */
  private RS unit;
  
  private RIS is;

  /**
   * 新しく生成された<code>VULS</code>オブジェクトを初期化します。
   */
  public RealVULS() {
    //do nothing
  }

  /**
   * 新しく生成された<code>VULS</code>オブジェクトを初期化します。
   * 
   * @param unit common data type
   * @param is interval scalar
   * @param m number of matrix
   */
  public RealVULS(RS unit, RIS is, int m) {
    this.unit = unit;
    this.is = is;
    this.m = m;
  }

  /**
   * get a zero vector.
   * 
   * @param dimension the dimension of vector
   * @return a vector
   * 
   */
  private RS[] getVectorZero(int dimension) {
    RS[] vector = this.unit.createArray(dimension);
    for (int i = 0; i < dimension; i++) {
      vector[i] = this.unit.createZero();
    }
    return vector;
  }

  /**
   * get a zero matrix
   * 
   * @param nRow index of row
   * @param nCol index of column
   * @return a matrix
   * 
   */
  public RS[][] getMatrixZero(int nRow, int nCol) {
    RS[][] matrix = this.unit.createArray(nRow, nCol);
    for (int i = 0; i < nRow; i++) {
      for (int j = 0; j < nCol; j++) {
        matrix[i][j] = this.unit.createZero();
      }
    }
    return matrix;
  }
  
  /**
   * Creates interval matrix
   * 
   * @param matrix matrix
   * @return interval matrix
   */
  private RIM createInterval(RM matrix) {
    return this.is.createGrid(new int[] {0}).create(matrix);
  }
  
  /**
   * get a interval matrix from sparse vector-array based on index array.
   * 
   * @param sparseVectorFi sparse vector
   * @param I index array
   * @return interval matrix
   */
  public RIM getIntervalMatrix(RS[][] sparseVectorFi, int[] I) {
    RS[][] B = this.getMatrixZero(this.m, I.length);
    for (int nRow = 0; nRow < this.m; nRow++) {
      for (int nCol = 0; nCol < I.length; nCol++) {
        B[nRow][nCol] = sparseVectorFi[nRow][I[nCol] - 1];
      }
    }
    RM nB = B[0][0].createGrid(B);
    RIM BI = createInterval(nB);
    return BI;
  }

  /**
   * get a index array from P matrix.
   * 
   * @param n dimension of P matrix
   * @param p dimension of index array
   * @param P P matrix of LU decomposition
   * @return index array
   */
  private int[] getIndexArray(final int n, final int p, IntMatrix P) {
    int[] I = new int[p];
    for (int nRow = 0; nRow < p; nRow++) {
      int pRow = 0, pCol = n;
      for (int iRow = 1; iRow <= n; iRow++) {
        for (int iCol = 1; iCol <= n; iCol++) {
          if (P.getIntElement(iRow, iCol) == 1) {
            if (iCol < pCol) {
              pCol = iCol;
              pRow = iRow;
            }
          }
        }
      }

      P.setElement(pRow, pCol, 0);
      if (pRow - 1 == 0) {
        pRow = n + 1;
      }
      I[nRow] = pRow - 1;
    }
    return I;
  }

  /**
   * get a non-index array from index array based on P matrix.
   * 
   * @param n dimension of P matrix
   * @param p p
   * @param I index array
   * @return non-index array from index array
   */
  private int[] getNonIndexArray(final int n, final int p, int[] I) {
    int[] N = new int[n - p];
    int countIndex = 0;

    for (int nRow = 1; nRow <= n; nRow++) {
      if (this.isI(I, nRow) == false) {
        N[countIndex++] = nRow;
      }
    }

    return N;
  }

  /**
   * check data whether is the value of array or not.
   * 
   * @param i array data
   * @param target target value
   * @return if it is the value of array <code>i</code>, true; else, false.
   */
  private boolean isI(int[] i, int target) {
    for (int nRow = 0; nRow < i.length; nRow++) {
      if (i[nRow] == target) {
        return true;
      }
    }

    return false;
  }

  /**
   * @param sparseVectorFi information of sparse block-matrices <code>F</code>
   * @param c2 vector <code>c</code>
   * @param sparseVectorY information of sparse matrix <code>Y</code>
   * @return a verified sparse matrix <code>Y</code>
   */
  public RS[] getVULS(RS[][] sparseVectorFi, RS[] c2, RS[] sparseVectorY) {
    RoundModeManager manager = RoundModeManager.getManager();
    manager.add(new FPURoundModeSelector());

    RM numericalMatrixSparseVectorFi = sparseVectorY[0].createGrid(sparseVectorFi);
    RIM intervalMatrixSparseVectorFi = createInterval(numericalMatrixSparseVectorFi);
    RM numericalMatrix_c2 = c2[0].createGrid(c2);
    RIM intervalMatrix_c2 = createInterval(numericalMatrix_c2.transpose());
    RM numericalMatrixSparseVectorY = sparseVectorY[0] .createGrid(sparseVectorY);
    RIM intervalMatrixSparseVectorY = createInterval(numericalMatrixSparseVectorY);
    final int n = sparseVectorY.length;
    final int p = c2.length;
    IntMatrix P = (numericalMatrixSparseVectorFi.transpose()).luDecomposeWithPermutation(false).getP();

    if (n == p) {
      LinearEquationVerifier<RIS,RIM,RS,RM> verifier = new RealIntlabMethod<>(intervalMatrixSparseVectorFi, intervalMatrix_c2);
      verifier.solve();
      RIM verified_Y = verifier.getSolution();

      int length = verified_Y.length();
      RS[] copy_Y = this.getVectorZero(n);
      for (int index = 0; index < length; index++) {
        copy_Y[index] = verified_Y.getElement(index + 1).getMiddle();
      }
      return copy_Y;
    }

    int[] I = new int[p];
    I = this.getIndexArray(n, p, P);
    
    int[] N = new int[n - p];
    N = this.getNonIndexArray(n, p, I);

    RIM BI = getIntervalMatrix(sparseVectorFi, I);
    RIM intervalMatrixBN = getIntervalMatrix(sparseVectorFi, N);

    RS[] YN = this.getVectorZero(N.length);
    for (int nRow = 0; nRow < N.length; nRow++) {
      YN[nRow] = sparseVectorY[N[nRow] - 1];
    }
    RM numericalMatrixYN =YN[0].createGrid(YN);
    RIM intervalMatrixYN = createInterval(numericalMatrixYN.transpose());

    RIM cI = null;
    if (N.length == 0) {
      cI = intervalMatrix_c2;
    } else {
      cI = intervalMatrix_c2.subtract(intervalMatrixBN.multiply(intervalMatrixYN));
    }

    //    Guarantor<E> guarantor = new Guarantor<E>(new LinearEquationVerifierFactory<E>(BI, cI));
    //    guarantor.solveWithEffectiveDigits(6);
    //    IM verified_Y = ((LinearEquationVerifiedSolution<E>)guarantor.getSolution()).getX();

    LinearEquationVerifier<RIS,RIM,RS,RM> verifier = new RealIntlabMethod<>(BI, cI);
    verifier.solve();
    RIM verified_Y = verifier.getSolution();

    final int length = verified_Y.length();
    RS[] copy_Y = this.getVectorZero(n);
    for (int index = 0; index < n; index++) {
      copy_Y[index] = sparseVectorY[index];
    }

    for (int index = 0; index < length; index++) {
      copy_Y[I[index] - 1] = verified_Y.getElement(index + 1).getMiddle();
    }

    return copy_Y;
  }
  
  /**
   * @param sparseVectorFi information of sparse block-matrices <code>F</code>
   * @param c2 vector <code>c</code>
   * @param sparseVectorY information of sparse matrix <code>Y</code>
   * @param I index vector
   * @param N index vector
   * @param BI submatrix(nonsingular) of B
   * @param intervalMatrixBN submatrix of B 
   * @return a verified sparse matrix <code>Y</code>
   */
  public RIM getEnclosureY(RS[][] sparseVectorFi, RS[] c2, RS[] sparseVectorY, int[] I, int[] N, RIM BI, RIM intervalMatrixBN) {
    RoundModeManager manager = RoundModeManager.getManager();
    manager.add(new FPURoundModeSelector());
//    if (sparseVectorFi[0][0].createUnit() instanceof MPFloat){
//      manager.add(new MPFloatRoundModeSelector());
//    }
    RM numericalMatrixSparseVectorFi = sparseVectorFi[0][0].createGrid(sparseVectorFi);
    RIM intervalMatrixSparseVectorFi = BI.create(numericalMatrixSparseVectorFi);
    RM numericalMatrix_c2 = c2[0].createGrid(c2);
    RIM intervalMatrix_c2 = BI.create(numericalMatrix_c2.transpose());
    RM numericalMatrixSparseVectorY =  sparseVectorY[0].createGrid(sparseVectorY);
    RIM intervalMatrixSparseVectorY = BI.create(numericalMatrixSparseVectorY);
    final int n = sparseVectorY.length;
    final int p = c2.length;
    
    if (n == p) {
      LinearEquationVerifier<RIS,RIM,RS,RM> verifier = new RealIntlabMethod<>(intervalMatrixSparseVectorFi, intervalMatrix_c2);
      verifier.solve();
      RIM xI = verifier.getSolution();

      final int length = xI.length();
      for (int index = 0; index < length; index++) {
        intervalMatrixSparseVectorY.setElement(index, xI.getElement(index + 1));
      }
      return intervalMatrixSparseVectorY;
    }

    RS[] YN = this.getVectorZero(N.length);
    for (int nRow = 0; nRow < N.length; nRow++) {
      YN[nRow] = sparseVectorY[N[nRow] - 1];
    }
    RM numericalMatrixYN = YN[0].createGrid(YN);
    RIM intervalMatrixYN = BI.create(numericalMatrixYN.transpose());

    RIM cI = null;
    if (N.length == 0) {
      cI = intervalMatrix_c2;
    } else {
      cI = intervalMatrix_c2.subtract(intervalMatrixBN.multiply(intervalMatrixYN));
    }
    
    LinearEquationVerifier<RIS,RIM,RS,RM> verifier = new RealIntlabMethod<>(BI, cI);
    verifier.solve();
    RIM xI = verifier.getSolution();

    final int length = xI.length();
    for (int index = 0; index < length; index++) {
      intervalMatrixSparseVectorY.setElement(I[index], xI.getElement(index + 1));
    }
    return intervalMatrixSparseVectorY;
  }

  /**
   * @param sparseVectorFi information of sparse block-matrices <code>F</code>
   * @param n dimension
   * @param p dimension
   * @return I index vector
   */
  public int[] getI(RS[][] sparseVectorFi, final int n, final int p) {
    IntMatrix P;
    int[] I = new int[p];
    I = getSortI(sparseVectorFi,I);
    if (I[p-1]==0) {
      RM numericalMatrixSparseVectorFi = sparseVectorFi[0][0].createGrid(sparseVectorFi);
      LUDecomposition<RS,RM> lud = (numericalMatrixSparseVectorFi.transpose()).luDecomposeWithPermutation(false);
      P = lud.getP();
      I = this.getIndexArray(n, p, P);
    }
    return I;
  }
  
  /**
   * get a non-index array from index array based on P matrix.
   * 
   * @param N dimension of P matrix
   * @param I index array
   * @return non-index array from index array
   */
  public int[] getN(int[] N, int[] I) {
    int countIndex = 0;
    int svl = N.length + I.length;
    for (int nRow = 1; nRow <= svl; nRow++) {
      if (this.isI(I, nRow) == false) {
        N[countIndex++] = nRow;
      }
    }
    return N;
  }
  
//  private int[] getIwithSortY(S[] sparseVectorY) {
//    int[] I = new int[sparseVectorY.length];
//    S[] sYwP = sparseVectorY.clone();
//    sYwP = quickSort(sYwP,0,sYwP.length-1);
//    for (int i = sYwP.length-1; i >=0 ; i--) {
//      System.out.println(sYwP[i]);
//      int k = I.length-1 - i;
//      for (int j = 0; j < sYwP.length; j++) {
//        if (sYwP[i].equals(sparseVectorY[j])) {
//          I[k] = j;
//          break;
//        }
//      }
//    }    
//    return I;
//  }
//  
//  private S[] quickSort(S[] array, int left, int right) {
//    int curleft = left;
//    int curright = right;
//    S pivot = array[(curleft + curright) / 2];
//    
//    do {
//        while (array[curleft].abs().isLessThan(pivot.abs())) {
//            curleft++;
//        }
//        
//        while (array[curright].abs().isGreaterThan(pivot.abs())) {
//            curright--;
//        }
//        if (curleft <= curright) {
//            swap (array, curleft++, curright--);
//        }
//    } while (curleft <= curright);
//    
//    if (left < curright) {
//        quickSort(array, left, curright);
//    }
//    
//    if (curleft < right) {
//      quickSort(array, curleft, right);
//    }
//    return array;
//  }
//
//  private void swap(S[] array, int idx1, int idx2) {
//    S tmp = array[idx1];
//    array[idx1] = array[idx2];
//    array[idx2] = tmp;
//  }
  
  private int[] getSortI(RS[][] sparseVectorFi, int[] I) {
    RS unit1 = sparseVectorFi[0][0].createUnit();
    int p = 0;
    for(int i = 0; i < sparseVectorFi.length; i++) {
      for (int j = 0; j < sparseVectorFi[0].length; j++) {
        if (sparseVectorFi[i][j].abs().equals(unit1) || sparseVectorFi[i][j].abs().equals(unit1.multiply(2))) {
          I[p] = j+1;
          p++;
          break;
        }
      }
    }
    return I;
  }
}