RealIntlabMethod.java

/*
 * Created on 2008/01/31
 * Copyright (C) 2008 Koga Laboratory. All rights reserved.
 *
 */
package org.mklab.cga.linear;

import java.util.ArrayList;
import java.util.List;

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.matrix.ComplexNumericalMatrix;
import org.mklab.nfc.matrix.NormType;
import org.mklab.nfc.matrix.RealNumericalMatrix;
import org.mklab.nfc.scalar.ComplexNumericalScalar;
import org.mklab.nfc.scalar.NumericalScalar;
import org.mklab.nfc.scalar.RealNumericalScalar;
import org.mklab.nfc.util.RoundMode;
import org.mklab.nfc.util.RoundModeManager;


/**
 * 線形方程式の精度保証付き解をINTLABに実装されている方法で求めるクラスです。
 * 
 * {@link NumericalScalar}を用いているので多倍長精度に対応しています。
 * 
 * @author yano
 * @version $Revision: 1.7 $, 2008/01/31
 * @param <RIS> 実区間スカラーの型
 * @param <RIM> 実区間行列の型
 * @param <CIS> 複素区間スカラーの型
 * @param <CIM> 複素区間行列の型
 * @param <RS> 実スカラーの型
 * @param <RM> 実行列の型
 * @param <CS> 複素スカラーの型
 * @param <CM> 複素行列の型
 */
public class RealIntlabMethod<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 LinearEquationVerifier<RIS, RIM, RS, RM> {

  /** 区間係数行列 */
  private RIM intA;

  /** 区間右辺ベクトル */
  private RIM intB;

  /** 精度保証付き解 */
  private RIM solution;

  /** 計算途中に使うパラメータ。何なのかよくわからない。 */
  private static int INTERVAL_RESIDUAL = 1;

  /** 計算途中に使うパラメータ。何なのかよくわからない。 */
  private static boolean INTERVAL_FIRST_SECOND_STAGE = true;

  //  /**
  //   * 新しく生成された<code>NumericalMatrixElementIntlabMethod</code>オブジェクトを初期化します。
  //   */
  //  public IntlabMethod() {
  //    //
  //  }
  //
  //  /**
  //   * 新しく生成された<code>NumericalMatrixElementIntlabMethod</code>オブジェクトを初期化します。
  //   * 
  //   * @param A 係数行列
  //   * @param b 右辺係数ベクトル
  //   */
  //  public IntlabMethod(final M A, final M b) {
  //    this.intA = IntervalMatrixFactory.toInterval(A);
  //    this.intB = IntervalMatrixFactory.toInterval(b);
  //  }
  //
  //  /**
  //   * 新しく生成された<code>NumericalMatrixElementIntlabMethod</code>オブジェクトを初期化します。
  //   * 
  //   * @param A 係数行列
  //   * @param b 右辺係数区間ベクトル
  //   */
  //  public IntlabMethod(final M A, final IM b) {
  //    this.intA = IntervalMatrixFactory.toInterval(A);
  //    this.intB = b.clone();
  //  }
  //
  //  /**
  //   * 新しく生成された<code>NumericalMatrixElementIntlabMethod</code>オブジェクトを初期化します。
  //   * 
  //   * @param A 係数区間行列
  //   * @param b 右辺係数ベクトル
  //   */
  //  public IntlabMethod(final IM A, final M b) {
  //    this.intA = A.clone();
  //    this.intB = IntervalMatrixFactory.toInterval(b);
  //  }
  //

  //  /**
  //   * {@inheritDoc}
  //   */
  //  public void solve(final M A, final M b) {
  //    this.intA = IntervalMatrixFactory.toInterval(A);
  //    this.intB = IntervalMatrixFactory.toInterval(b);
  //    solve();
  //  }
  //
  //  /**
  //   * {@inheritDoc}
  //   */
  //  public void solve(final M A, final IM b) {
  //    this.intA = IntervalMatrixFactory.toInterval(A);
  //    this.intB = b.clone();
  //    solve();
  //  }
  //
  //  /**
  //   * {@inheritDoc}
  //   */
  //  public void solve(final IM A, final M b) {
  //    this.intA = A.clone();
  //    this.intB = IntervalMatrixFactory.toInterval(b);
  //    solve();
  //  }

  /**
   * 新しく生成された<code>NumericalMatrixElementIntlabMethod</code>オブジェクトを初期化します。
   * 
   * @param A 係数区間行列
   * @param b 右辺係数区間ベクトル
   */
  public RealIntlabMethod(final RIM A, final RIM b) {
    this.intA = A;
    this.intB = b;
  }

  /**
   * {@inheritDoc}
   */
  public void solve() {
//    this.intA = A;
//    this.intB = b;

    RM midA = this.intA.getMiddle();
    RM midb = this.intB.getMiddle();

    //    String message;
    //    if (((message = Abcdchk.abcdchk(midA, midb)).length()) > 0) {
    //      throw new IllegalArgumentException(message);
    //    }

    int k = midA.getColumnSize();
    int n = midA.getColumnSize();

    RM R = midA.inverse();
    RM xs = R.multiply(midb);// 近似解の計算。

    // Improve residual calculation
    // Interval iteration
    List<RIM> ires_Z_RA = intervalIteration(this.intA, this.intB, R, xs);
    RIM ires = ires_Z_RA.get(0);
    RIM Z = ires_Z_RA.get(1);
    RIM RA = ires_Z_RA.get(2);

    this.solution = stage12(ires, Z, RA, R, xs, n, k);
  }

  //  /**
  //   * @see org.mklab.cga.linear.LinearEquationVerifier#getSolution()
  //   */
  //  public LinearEquationVerifiedSolution<E> getSolution() {
  //    return new LinearEquationVerifiedSolution(this.solution);
  //  }

  /**
   * Creates interval matrix
   * 
   * @param matrix matrix
   * @return interval matrix
   */
  private RIM createInterval(RM matrix) {
    return this.intA.create(matrix);
  }

  /**
   * Creates interval with middle and radius.
   * 
   * @param mid middle
   * @param rad radius
   * @return interval matrxi
   */
  private RIM createMidRad(RM mid, RM rad) {
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    manager.setRoundMode(RoundMode.ROUND_DOWN);
    RM inf = mid.subtract(rad);

    manager.setRoundMode(RoundMode.ROUND_UP);
    RM sup = mid.add(rad);

    RIM ans = this.intA.createInfSup(inf, sup);

    manager.setRoundMode(oldRoundMode);

    return ans;
  }

  /**
   * @param A A
   * @param b b
   * @param R R
   * @param xs xs
   * @return ??
   */
  private List<RIM> intervalIteration(RIM A, RIM b, RM R, RM xs) {
    RIM ires = b.subtract(A.multiply(createInterval(xs)));
    RIM Z = createInterval(R).multiply(ires);
    RIM RA = createInterval(R).multiply(A);

    List<RIM> ans = new ArrayList<>();
    ans.add(ires);
    ans.add(Z);
    ans.add(RA);
    return ans;
  }

  /**
   * @param midA Aの中心
   * @param midb Bの中心
   * @param R R
   * @param xs xs
   * @return ??
   */
  private List<Object> improveResidual(RM midA, RM midb, RM R, RM xs) {
    boolean improvedresidual = (INTERVAL_RESIDUAL == 1) | ((INTERVAL_RESIDUAL == 2) & (midb.getColumnSize() > 1));

    RM localxs = xs;

    // IMPROVE RESIDUAL CALCULATION
    if (improvedresidual) {
      // TODO 表現可能な最大数を返すメソッドがいる?(これよりは小さい値のはず？？ 精度が落ちるから駄目?)
      RS resnorm = midA.getElement(1).getInfinity(); //transformFrom(new MPFloat(1).divide(new MPFloat(0)));

      int i = 0;

      while (i < 15) {
        i = i + 1;
        RS resnormold = resnorm.clone();
        RM res = R.multiply(ResidualImprover.lssresidual(midA, xs, midb));

        // TODO ここで計算失敗（０に値が近いから？？）
        resnorm = (res).norm(NormType.TWO);// 最大特異値
        if (resnorm.isLessThan(resnormold)) {
          localxs = localxs.add(res);
        }

        // TODO 精度が落ちそう
        if (resnorm.isGreaterThanOrEquals(resnormold.divide(10))) {// 0の剰余に注意
          break;
        }
      }

    } else if (INTERVAL_RESIDUAL == 2) {// quadruple precision residual
      // calculation??
      // TODO 最大値を返すメソッドがいる？(この値より小さいはず？ 精度に影響がある？？)
      //S resnorm = (midA).getElement(1).transformFrom(Double.POSITIVE_INFINITY);
      RS resnorm = (midA).getElement(1);
      int i = 0;

      while (i < 15) {
        i = i + 1;

        RS resnormold = resnorm.clone();
        RM res = R.multiply(midb.unaryMinus().add(midA.multiply(localxs)));
        resnorm = (res).norm(NormType.TWO);

        if (resnorm.isLessThan(resnormold)) {
          localxs = localxs.subtract(res);
        }

        // TODO 精度に影響しない？
        if (resnorm.isLessThanOrEquals(resnormold.multiply(1e-1))) {// beware of zero residual
          break;
        }
      }

    } else {
      localxs = localxs.add(R.multiply(midb.subtract(midA.multiply(localxs))));
    }

    List<Object> list = new ArrayList<>();
    list.add(localxs);
    list.add(Boolean.valueOf(improvedresidual));
    return list;
  }

  /**
   * intbalMethod内で使用する計算。
   * 
   * @param ires ires
   * @param Z Z
   * @param RA RA
   * @param R R
   * @param xs xs
   * @param n n
   * @param k k
   * @return result
   */
  private RIM stage12(RIM ires, RIM Z, RIM RA, RM R, RM xs, int n, int k) {
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    manager.setRoundMode(RoundMode.ROUND_NEAR);

    // FIRST STAGE
    RIM X = null;
    if (INTERVAL_FIRST_SECOND_STAGE) {
      RS unit = R.getElement(1, 1).abs();

      // TODO 精度にいぞんする？これより小さい値もありうるよね？
      RS realmin = unit.create(Double.MIN_VALUE); // Z.getMiddle().getElement(1).transformFrom(Double.MIN_VALUE);
      RIM C = createInterval(R.createUnit(n)).subtract(RA);

      RIS type = this.intA.getElement(1, 1);
      RIM Y = Z;
      // TODO 具体的なMPFloatを使わない(精度を落とさずにできるの？)
      RIS minusPlusOneInterval = type.createInfSup(unit.createUnit().unaryMinus(), unit.createUnit());
      RIS disk = type.createInfSup(realmin.multiply(10).unaryMinus(), realmin.multiply(10));
      RM radius = Y.getSupremum().subtract(Y.getInfimum()).divide(2);
      RIM yRadius = createInterval(radius.divide(10));
      RIM E = yRadius.multiply(minusPlusOneInterval).addElementWise(disk);

      int localk = 0;
      final int kmax = 7;
      boolean ready = false;

      while ((!ready) & (localk < kmax) & (!Y.isNanElementWise().anyTrue())) {
        localk = localk + 1;
        X = Y.add(E);
        Y = Z.add(C.multiply(X));
        ready = IntervalUtil.in0(Y, X).allTrue();
      }

      if (ready) {// SUCCESS FIRST STAGE
        X = Y.add(createInterval(xs));
        manager.setRoundMode(oldRoundMode);
        return X;
      }
    }

    // SECOND STAGE STARTS
    RIM b = createInterval(R).multiply(ires);
    RIM dRA = RA.diagonalToVector();
    RM A = RA.comparisonMatrix();
    RM B = A.inverse();
    RM v = (B.multiply(A.createOnes(n, 1))).absElementWise();

    manager.setRoundMode(RoundMode.ROUND_DOWN);

    RM u = A.multiply(v);

    if ((u).minColumnWise().compareElementWise(".>", A.createOnes(1, u.getColumnSize())).allTrue()) { //$NON-NLS-1$
      RM dAc = A.diagonalToVector();
      A = A.multiply(B).subtract(A.createUnit(n));

      manager.setRoundMode(RoundMode.ROUND_UP);

      RM w = (A.unaryMinus().divideElementWise(u.multiply(A.createOnes(1, n)))).maxColumnWise().transpose();
      RM dlow = v.multiplyElementWise(w.conjugateTranspose()).subtract(B.diagonalToVector());
      dlow = (dlow.unaryMinus()).maxElementWise(A.createZero(dlow.getRowSize(), dlow.getColumnSize()));

      B = B.add(v.multiply(w));
      u = B.multiply(b.abssElementWise());
      RM d = B.diagonalToVector();
      RM alpha = dAc.add(A.createOnes(d.getRowSize(), d.getColumnSize()).unaryMinus().divideElementWise(d));

      int localk = b.getColumnSize();

      RIM Ixs = createInterval(xs);

      if (localk == 1) {
        RM beta = u.divideElementWise(dlow).subtract(b.abssElementWise());
        RM mid1 = A.createZero(beta.getRowSize(), beta.getColumnSize());
        RM mid2 = A.createZero(alpha.getRowSize(), alpha.getColumnSize());
        X = Ixs.add((b.add(createMidRad(mid1, beta))).divideElementWise((dRA.add(createMidRad(mid2, alpha)))));
      } else {
        // d and dRA adapted for multiple r.h.s.
        v = A.createOnes(1, localk);
        RM beta = u.divideElementWise(d.multiply(v)).subtract(b.abssElementWise());
        RM mid1 = A.createZero(beta.getRowSize(), beta.getColumnSize());
        RM mid2 = A.createZero(alpha.getRowSize(), alpha.getColumnSize());
        X = Ixs.add(b.add(createMidRad(mid1, beta)).divideElementWise((dRA.add(createMidRad(mid2, alpha)).multiply(createInterval(v)))));

      }
    }

    manager.setRoundMode(oldRoundMode);
    return X;
  }

  /**
   * {@inheritDoc}
   */
  public RIM getSolution() {
    return this.solution;
  }

  /**
   * Returns A.
   * 
   * @return A
   */
  protected RIM getA() {
    return this.intA;
  }
  
  /**
   * Returns B.
   * 
   * @return B
   */
  protected RIM getB() {
    return this.intB;
  }
}