RealRumpMethod.java

/*
 * Created on 2004/08/03
 * 
 * Copyright 2004. Created by hiroki
 */
package org.mklab.cga.eigen;

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.util.IntervalUtil;
import org.mklab.nfc.eig.EigenSolution;
import org.mklab.nfc.matrix.BooleanMatrix;
import org.mklab.nfc.matrix.ComplexNumericalMatrix;
import org.mklab.nfc.matrix.ElementHolder;
import org.mklab.nfc.matrix.IndexedMatrix;
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.RoundMode;
import org.mklab.nfc.util.RoundModeManager;


/**
 * 多重および近接固有値の精度保証を行うクラスです。 Rumpの方法を実装しています。
 * 
 * @author hiroki
 * @param <RIS> 実区間スカラーの型
 * @param <RIM> 実区間行列の型
 * @param <CIS> 複素区間スカラーの型
 * @param <CIM> 複素区間行列の型
 * @param <RS> 実スカラーの型
 * @param <RM> 実行列の型
 * @param <CS> 複素スカラーの型
 * @param <CM> 複素行列の型
 */
public class RealRumpMethod<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>>
    implements RealEigenVerifier<RIS, RIM, CIS, CIM, RS, RM, CS, CM> {

  /** 対象となる行列 */
  private RM A;

  /** 対象となる区間行列 */
  private RIM IA;

  /** 精度保証付き固有値 */
  private CIM eigenValue;

  /** 精度保証付き固有ベクトル(不変部分空間の基底ベクトル) */
  private CIM eigenVector;

  /**
   * 新しく生成された<code>RumpMethod</code>オブジェクトを初期化します。
   * 
   * @param IA interval matrix
   */
  public RealRumpMethod(RIM IA) {
    this.IA = IA;
    this.A = IA.getMiddle();
  }
  
  /**
   * Creates interval matrix
   * 
   * @param matrix real matrix
   * @return interval real matrix
   */
  private RIM createInterval(RM matrix) {
    return this.IA.create(matrix);
  }

  /**
   * Creates interval matrix
   * 
   * @param mid middle
   * @param rad radius
   * @return rim
   */
  private RIM createMidRad(RM mid, RM rad) {
    return this.IA.createMidRad(mid, rad);
  }

  /**
   * Creates interval matrix
   * 
   * @param mid middle
   * @param rad radius
   * @return rad radius
   */
  private CIM createMidRad(CM mid, RM rad) {
    return this.IA.toComplex().createMidRad(mid, rad);
  }

  /**
   * Creates interval matrix
   * 
   * @param matrix complex matrix
   * @return interval complex matrix
   */
  private CIM createInterval(CM matrix) {
    return this.IA.toComplex().create(matrix);
  }

  /**
   * Creates interval matrix
   * 
   * @param matrix complex matrix
   * @return interval complex matrix
   */
  private CIS createInterval(CS matrix) {
    return this.IA.getElement(1,1).toComplex().create(matrix);
  }

  /**
   * Creates interval scalar with infimum and supremum.
   * 
   * @param inf infimum
   * @param sup supremum
   * @return interval scalara
   */
  private RIS createIntervalInfSup(double inf, double sup) {
    RS type = this.A.getElement(1,1);
    RIS itype = this.IA.getElement(1,1);
    return itype.createInfSup(type.create(inf), type.create(sup));
  }

  /**
   * Creates interval scalar with middle and radius.
   * 
   * @param inf infimum
   * @param sup supremum
   * @return interval scalara
   */
  private RIS createMidRad(double inf, double sup) {
    RS type = this.A.getElement(1,1);
    RIS itype = this.IA.getElement(1,1);
    return itype.createMidRad(type.create(inf), type.create(sup));
  }

  //  /**
  //   * @see org.mklab.cga.eigen.EigenVerifier#getSolution()
  //   */
  //  public EigenVerifiedSolution<RIS,RIM,CIS,CIM,RS,RM,CS,CM> getSolution() {
  //    return new EigenVerifiedSolution<>(this.eigenValue, this.eigenVector);
  //  }

  /**
   * {@inheritDoc}
   */
  public CIM getEigenValue() {
    return this.eigenValue;
  }

  /**
   * {@inheritDoc}
   */
  public CIM getEigenVector() {
    return this.eigenVector;
  }

//  /**
//   * {@inheritDoc}
//   */
//  public void solve() {
//  }

//  /**
//   * {@inheritDoc}
//   */
//  public void solveAsSparse() {
//    if (this.A != null) {
//      solveAsSparse(this.A);
//    } else if (this.IA != null) {
//      solveAsSparse(this.IA);
//    } else {
//      throw new RuntimeException("Target matrix is not given."); //$NON-NLS-1$
//    }
//  }

//  /**
//   * {@inheritDoc}
//   */
//  public void solve(final RM matrix) {
//    final EigenSolution<RS, RM, CS, CM> eigen = (matrix).eigenDecompose();
//    CM value = eigen.getValues();
//    CM vector = eigen.getVectors();
//    this.eigenValue = this.IA.toComplex().create(value.createZero(value.getRowSize(), 1));
//    this.eigenVector = this.IA.toComplex().create(vector.createZero(vector.getRowSize(), vector.getColumnSize()));
//
//    for (int i = 0; i < value.getRowSize(); i++) {
//      final CS rambda = value.getElement(i + 1, i + 1);
//      final CM xs = vector.getColumnVector(i + 1);
//      final CIM[] ans = solve(matrix, rambda, xs);
//      this.eigenValue.setRowVector(i + 1, ans[0]);// 固有値
//      this.eigenVector.setColumnVector(i + 1, ans[1]);// 不変部分空間の基底
//    }
//  }

  /**
   * {@inheritDoc}
   */
  public void solve() {
    final EigenSolution<RS, RM, CS, CM> eigen = this.A.eigenDecompose();
    CM value = eigen.getValue();
    CM vector = eigen.getVector();
    this.eigenValue = this.IA.toComplex().create(value.createZero(value.getRowSize(), 1));
    this.eigenVector = this.IA.toComplex().create(vector.createZero(vector.getRowSize(), vector.getColumnSize()));

    for (int i = 0; i < value.getRowSize(); i++) {
      final CS rambda = value.getElement(i + 1, i + 1);
      final CM xs = vector.getColumnVector(i + 1);
      final CIM[] ans = solve(this.A, rambda, xs);
      this.eigenValue.setRowVector(i + 1, ans[0]);// 固有値
      this.eigenVector.setColumnVector(i + 1, ans[1]);// 不変部分空間の基底
    }
  }

  /**
   * 多重および近接固有値の精度保証を行います。
   * 
   * 一つの固有値に対し、精度保証を行います。
   * 
   * 近接する固有値の境界点を求めたい場合は、 <code>xs</code> には 近接する固有値の固有ベクトル列を設定してください。
   * 
   * @param matrix 固有値を求めたい行列
   * @param lambda Aの固有値のひとつ
   * @param xs lambdaに対応する固有ベクトル
   * @return {精度保証付き固有値,精度保証付き固有ベクトル(不変部分空間の基底)}
   */
  private CIM[] solve(final RM matrix, CS lambda, final CM xs) {
    return solve(this.IA.create(matrix), lambda, xs);
  }

  /**
   * 精度保証付き固有値と精度保証付き固有ベクトル求めます。
   * 
   * @param IA 固有値を求めたい区間行列
   * @param lambda Aの固有値のひとつ
   * @param xs lambdaに対応する固有ベクトル
   * @return {精度保証付き固有値,精度保証付き固有ベクトル(不変部分空間の基底)}
   */
  @SuppressWarnings("unchecked")
  private CIM[] solve(RIM IA, CS lambda, CM xs) {
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    int n = xs.getRowSize();// 行数
    int k = xs.getColumnSize();// 列数

    manager.setRoundMode(RoundMode.ROUND_NEAR);

    CM xss = (xs).absElementWise().sumRowWise();
    IndexedMatrix<CS, CM> NI = (xss).sortColumnWise();

    IntMatrix I = NI.getIndices();
    IntMatrix u = I.getRowVectors(1, n - k).transpose();
    IntMatrix v = I.getRowVectors(n - k + 1, n).transpose();
    RM midA = IA.getMiddle();

    // TODO Aが実数行列じゃないときもこれでいいの？
    CM R = midA.toComplex().subtract(midA.createUnit(n).toComplex().multiply(lambda));// R=midA-lambda*I
    R.setColumnVectors(v, xs.unaryMinus()); // R(:, v) = -xs
    CM y = R.leftDivide(midA.toComplex().multiply(xs).subtract(xs.multiply(lambda)));// y = R\(midA*xs-lambda*xs)

    xs.setRowVectors(u, (xs.getRowVectors(u).subtract(y.getRowVectors(u))));// xs(u,:)= xs(u,:)-y(u,:)

    // lambda - sum(diag(y(v,:)))/k
    CS locallambda = lambda.subtract((y.getRowVectors(v).diagonalToVector().sumColumnWise()).getElement(1, 1).divide(k));
    R = midA.toComplex().subtract(midA.toComplex().createUnit(n).multiply(locallambda));// R=midA-lambda*I
    R.setColumnVectors(v, xs.unaryMinus());// R(:,v)=-xs
    R = R.inverse();// R=inv(R)

    RIM II =  this.IA.create(midA.createUnit(n));// 単位行列のIntervalMatrixの作成
    CIS I_lambda = createInterval(locallambda);

    CIM IR = this.IA.toComplex().create(R);
    CIM IC = IA.toComplex().subtract(II.toComplex().multiply(I_lambda));// C=A-intval(lambda)*I
    CIM IZ = (createInterval(R.unaryMinus())).multiply(IC.multiply(createInterval(xs)));// Z=-R*(C*xs)
    IC.setColumnVectors(v, createInterval(xs.unaryMinus()));// C(:, v)=-xs
    IC = II.toComplex().subtract(IR.multiply(IC));// C=I-RC
    CIM IY = IZ.createClone();// Y=Z

    CIM Eps = createInterval(IY.abssElementWise().multiply(0.1)).multiply(createIntervalInfSup(-1.0, 1.0).toComplex()).addElementWise(createMidRad(0, Double.MIN_VALUE).toComplex());
    // TODO こうしたらいいのかな？(NumericalMatrixElementIntervalの実装が必要みたい)
    CIM Z_row = IZ.getRowVectors(v);
    RM Z_row_abs = Z_row.abssElementWise().getRealPart();
    BooleanMatrix Z_flag = Z_row_abs.compareElementWise(".>", 0.1); //$NON-NLS-1$
    int z_sum = Z_flag.find().length();

    double m = 0;
    double mmax = Math.min(15 * (z_sum + 1), 20);

    boolean ready = false;

    CIM X = null;
    CIM XX = createInterval(R.createZero(IY.getRowSize(), IY.getColumnSize()));

    while ((!ready) && (m < mmax) && (!IY.isNanElementWise().anyTrue())) {
      m = m + 1;
      X = IY.add(Eps);
      XX.copy(X);
      XX.setRowVectors(v,createInterval(R.createZero(v.getColumnSize(), XX.getColumnSize())));
      IY = IZ.add(IC.multiply(X)).add(IR.multiply(XX.multiply(X.getRowVectors(v))));
      ready = IntervalUtil.in0forRump(IY, X).allTrue();
    }

    CIM L = null;

    if (ready) {
      RM M = (IY.getRowVectors(v)).abssElementWise().getRealPart();

      if (v.length() == 1) {
        CS[] value = locallambda.createArray(1);
        value[0] = locallambda;
        CM mid = locallambda.createGrid(value);
        L =  createMidRad(mid, M);
      } else {
        EigenSolution<RS, RM, CS, CM> EV = (M).eigenDecompose();
        ElementHolder<?> list = ((EV.getValue().diagonalToVector()).absElementWise()).maximum();
        int index = list.getRow();

        CM Perronx = EV.getVector().getColumnVector(index);

        manager.setRoundMode(RoundMode.ROUND_UP);

        /* upper bound for Perron root */
        RM rad = ((M.toComplex().multiply(Perronx).divideElementWise(Perronx))).maxColumnWise().getRealPart();

        manager.setRoundMode(RoundMode.ROUND_NEAR);
        CS[] value = locallambda.createArray(1);
        value[0] = locallambda;
        CM mid = locallambda.createGrid(value);
        L = createMidRad(mid, rad);
      }

      IY.setRowVectors(v, IY.createZero(v.getColumnSize(), IY.getColumnSize()));
      X = createInterval(xs).add(IY);
    } else {
      RS nanElement = IA.getMiddle().getElement(1).getNaN();
      RS[] value = nanElement.createArray(1);
      value[0] = nanElement;
      CM nan = IA.getMiddle().getElement(1).createGrid(value).toComplex();
      CM nanL = xs.createZero(1, 1);
      nanL.setRowVector(1, nan);

      CM nanX = xs.createZero(IA.getRowSize(), 1);
      for (int i = 0; i < IA.getColumnSize(); i++) {
        nanX.setRowVector(i + 1, nanL);
      }

      L = createMidRad(nanL, nanL.getRealPart());
      X = createMidRad(nanX, nanX.getRealPart());

      System.err.println("no inclusion achieved"); //$NON-NLS-1$
      //throw new RuntimeException("no inclusion achieved"); //$NON-NLS-1$
    }

    manager.setRoundMode(oldRoundMode);
    CIM[] ans = (CIM[])new IntervalComplexNumericalMatrix[2];
    ans[0] = L;
    ans[1] = X;
    return ans;
  }

  /**
   * @param IA 固有値を求めたい区間行列
   * @param lambda Aの固有値のひとつ
   * @param xs lambdaに対応する固有ベクトル
   * @return {精度保証付き固有値,精度保証付き固有ベクトル(不変部分空間の基底)}
   */
  private CIM[] solveSparse(RIM IA, CS lambda, CM xs) {
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    int n = xs.getRowSize();// 行数
    int k = xs.getColumnSize();// 列数

    manager.setRoundMode(RoundMode.ROUND_NEAR);

    CM xss = (xs).absElementWise().sumRowWise();
    IndexedMatrix<CS, CM> NI = (xss).sortColumnWise();
    IntMatrix I = NI.getIndices();
    IntMatrix v = I.getRowVectors(n - k + 1, n).transpose();
    RM midA = IA.getMiddle();

    CS locallambda = lambda.clone();
    CM R = midA.toComplex().subtract(midA.createUnit(n).toComplex().multiply(locallambda));// R=midA-lambda*I
    R.setColumnVectors(v, xs.unaryMinus());// R(:,v)=-xs
    CIM G =  IA.toComplex().create(R).inverse();

    CIM IZ = IA.toComplex().unaryMinus().multiply(createInterval(xs)).add(createInterval(xs.multiply(lambda)));
    CIM IY = IZ.createClone();// Y=Z

    // TODO ここ精度に依存しそう
    RIM Eps = createInterval(IY.abssElementWise().getRealPart().multiply(0.1)).multiply(createIntervalInfSup(-1.0, 1.0)).addElementWise(createMidRad(0, Double.MIN_VALUE));
    CIM Z_row = IZ.getRowVectors(v);
    CM Z_row_abs = Z_row.abssElementWise();
    BooleanMatrix Z_flag = Z_row_abs.compareElementWise(".>", 0.1); //$NON-NLS-1$
    int z_sum = Z_flag.find().length();

    double m = 0;
    double mmax = Math.min(15 * (z_sum + 1), 20);

    boolean ready = false;

    CIM X = null;
    CIM XX = createInterval(R.createZero(IY.getRowSize(), IY.getColumnSize()));

    while ((!ready) && (m < mmax) && (!IY.isNanElementWise().anyTrue())) {
      m = m + 1;
      X = IY.add(Eps.toComplex());
      XX.copy(X);
      XX.setRowVectors(v, XX.createZero(v.getColumnSize(), XX.getColumnSize()));
      CIM b = IZ.add(XX.multiply(X.getRowVectors(v)));
      IY = G.inverse().multiply(b);
      ready = IntervalUtil.in0forRump(IY, X).allTrue();
    }

    CIM L = null;

    if (ready) {
      RM M = (IY.getRowVectors(v)).abssElementWise().getRealPart();

      if (v.length() == 1) {
        CS[] value = locallambda.createArray(1);
        value[0] = locallambda;
        CM mid = locallambda.createGrid(value);
        L = createMidRad(mid, M);
      } else {
        EigenSolution<RS, RM, CS, CM> EV = (M).eigenDecompose();
        ElementHolder<?> list = (EV.getValue().diagonalToVector()).absElementWise().maximum();
        int index = list.getRow();

        CM Perronx = EV.getVector().getColumnVector(index);

        manager.setRoundMode(RoundMode.ROUND_UP);

        /* upper bound for Perron root */
        RM rad = ((M.toComplex().multiply(Perronx).divideElementWise(Perronx))).maxColumnWise().getRealPart();

        manager.setRoundMode(RoundMode.ROUND_NEAR);
        CS[] value = locallambda.createArray(1);
        value[0] = locallambda;
        CM mid = locallambda.createGrid(value);
        L = createMidRad(mid, rad);
      }

      IY.setRowVectors(v, createInterval(R.createZero(v.getColumnSize(), IY.getColumnSize())));
      X = createInterval(xs).add(IY);
    } else {
      CS nanElement = IA.getMiddle().getElement(1).getNaN().toComplex();
      CS[] value = nanElement.createArray(1);
      value[0] = nanElement;
      CM nan = IA.getMiddle().getElement(1).toComplex().createGrid(value);
      CM nanL = xs.createZero(1, 1);
      nanL.setRowVector(1, nan);

      CM nanX = xs.createZero(IA.getRowSize(), 1);
      for (int i = 0; i < IA.getColumnSize(); i++) {
        nanX.setRowVector(i + 1, nanL);
      }

      L = createMidRad(nanL, nanL.getRealPart());
      X = createMidRad(nanX, nanX.getRealPart());

      System.err.println("no inclusion achieved"); //$NON-NLS-1$
      //throw new RuntimeException("no inclusion achieved"); //$NON-NLS-1$
    }

    manager.setRoundMode(oldRoundMode);
    CIM[] ans = (CIM[])new IntervalComplexNumericalMatrix[2];
    ans[0] = L;
    ans[1] = X;
    return ans;
  }

//  /**
//   * @see org.mklab.cga.eigen.EigenVerifier#solveAsSparse(org.mklab.nfc.matrix.NumericalMatrix)
//   */
//  public void solveAsSparse(RM A) {
//    throw new UnsupportedOperationException();
//  }

  /**
   * {@inheritDoc}
   */
  public void solveAsSparse() {
    final EigenSolution<RS, RM, CS, CM> eigen = (this.IA.getMiddle()).eigenDecompose();
    CM value = eigen.getValue();
    CM vector = eigen.getVector();
    this.eigenValue = createInterval(value.createZero(value.getRowSize(), 1));
    this.eigenVector = createInterval(vector.createZero(vector.getRowSize(), vector.getColumnSize()));

    for (int i = 0; i < value.getRowSize(); i++) {
      final CS rambda = (value).getElement(i + 1, i + 1);
      final CM xs = vector.getColumnVector(i + 1);
      final CIM[] ans = solveSparse(this.IA, rambda, xs);
      this.eigenValue.setRowVector(i + 1, ans[0]);// 固有値
      this.eigenVector.setColumnVector(i + 1, ans[1]);// 不変部分空間の基底
    }
  }

//  /**
//   * 不変部分空間の基底を求めます。
//   * 
//   * @param a 不変部分空間の基底を求める行列
//   * @return 不変部分空間の基底
//   */
//  CIM getBasis(RM a) {
//    solve(a);
//    return this.eigenVector;
//  }

//  /**
//   * 不変部分空間の基底を求めます。
//   * 
//   * @return 不変部分空間の基底
//   */
//  CIM getBasis() {
//    solve();
//    return this.eigenVector;
//  }
  
  /**
   * Returns A.
   * @return A
   */
  protected RIM getA() {
    return this.IA;
  }

}