DoubleIntlabMethod.java

/*
 * Created on 2004/10/18
 *
 * $Id: IntlabMethod.java,v 1.25 2008/01/26 04:14:30 koga Exp $
 * Copyringht (C) 2004 Koga Laboratory. All rights reserved.
 */
package org.mklab.cga.linear;

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

import org.mklab.cga.interval.matrix.DoubleIntervalMatrix;
import org.mklab.cga.interval.scalar.DoubleIntervalNumber;
import org.mklab.cga.util.IntervalUtil;
import org.mklab.nfc.matrix.DoubleMatrix;
import org.mklab.nfc.matrix.NormType;
import org.mklab.nfc.scalar.DoubleNumber;
import org.mklab.nfc.util.RoundMode;
import org.mklab.nfc.util.RoundModeManager;


/**
 * 線形方程式の精度保証付き解をINTLABに実装されている方法で求めるクラスです。
 * 
 * @author hiroki
 * @version $Revision: 1.25 $. 2004/10/18
 */
public class DoubleIntlabMethod implements LinearEquationVerifier<DoubleIntervalNumber,DoubleIntervalMatrix,DoubleNumber,DoubleMatrix> {

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

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

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

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

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

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



//  /**
//   * {@inheritDoc}
//   */
//  public void solve(final DoubleMatrix A, final DoubleMatrix b) {
//    this.intA = IntervalMatrixFactory.toInterval(A);
//    this.intB = IntervalMatrixFactory.toInterval(b);
//    solve();
//  }
//
//  /**
//   * @see org.mklab.cga.linear.LinearEquationVerifier#solve(org.mklab.nfc.matrix.NumericalMatrix, org.mklab.cga.interval.matrix.IntervalMatrix)
//   */
//  public void solve(final DoubleMatrix A, final DoubleIntervalMatrix b) {
//    this.intA = IntervalMatrixFactory.toInterval(A);
//    this.intB = b;
//    solve();
//  }
//
//  /**
//   * @see org.mklab.cga.linear.LinearEquationVerifier#solve(org.mklab.cga.interval.matrix.IntervalMatrix, org.mklab.nfc.matrix.NumericalMatrix)
//   */
//  public void solve(final DoubleIntervalMatrix A, final DoubleMatrix b) {
//    this.intA = A;
//    this.intB = IntervalMatrixFactory.toInterval(b);
//    solve();
//  }
  
/**
 * 新しく生成された<code>IntlabMethod</code>オブジェクトを初期化します。
 * 
 * @param A 係数区間行列
 * @param b 右辺係数区間ベクトル
 */
public DoubleIntlabMethod(final DoubleIntervalMatrix A, final DoubleIntervalMatrix b) {
  this.intA = A;
  this.intB = b;
}

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

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

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

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

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

    // Improve residual calculation
    List<Object> improve = improveResidual(midA, midb, R, xs);
    xs = (DoubleMatrix)improve.get(0);

    // Interval iteration
    DoubleIntervalMatrix[] ires_Z_RA = intervalIteration(this.intA, this.intB, R, xs);
    DoubleIntervalMatrix ires = ires_Z_RA[0];
    DoubleIntervalMatrix Z = ires_Z_RA[1];
    DoubleIntervalMatrix RA = ires_Z_RA[2];

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

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

  /**
   * @param A A
   * @param b b
   * @param R R
   * @param xs xs
   * @return ??
   */
  private DoubleIntervalMatrix[] intervalIteration(DoubleIntervalMatrix A, DoubleIntervalMatrix b, DoubleMatrix R, DoubleMatrix xs) {
    DoubleIntervalMatrix ires = b.subtract(A.multiply(new DoubleIntervalMatrix(xs)));
    DoubleIntervalMatrix Z = new DoubleIntervalMatrix(R).multiply(ires);
    DoubleIntervalMatrix RA = new DoubleIntervalMatrix(R).multiply(A);
    return new DoubleIntervalMatrix[] {ires, Z, RA};
  }

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

    // Improve residual calculation
    if (improvedResidual) {
      double resNorm = Double.POSITIVE_INFINITY;

      int i = 0;
      while (i < 15) {
        i = i + 1;
        double resNormOld = resNorm;
        DoubleMatrix res = R.multiply(ResidualImprover.lssresidual(midA, xs, midb));

        resNorm = (res.norm(NormType.TWO)).doubleValue();//最大特異値
        if (resNorm < resNormOld) {
          localxs = localxs.add(res);
        }

        if (resNorm >= (1e-1) * resNormOld) {//0の剰余に注意
          break;
        }
      }

    } else if (INTERVAL_RESIDUAL == 2) {//quadruple precision residual
      double resNorm = Double.POSITIVE_INFINITY;

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

        double resNormOld = resNorm;
        DoubleMatrix res = R.multiply(midb.unaryMinus().add(midA.multiply(localxs)));
        resNorm = (res.norm(NormType.TWO)).doubleValue();

        if (resNorm < resNormOld) {
          localxs = localxs.subtract(res);
        }

        if (resNorm >= 1e-1 * resNormOld) {//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 ??
   */
  private DoubleIntervalMatrix stage12(DoubleIntervalMatrix ires, DoubleIntervalMatrix Z, DoubleIntervalMatrix RA, DoubleMatrix R, DoubleMatrix xs, int n, int k) {
    final RoundModeManager manager = RoundModeManager.getManager();
    final RoundMode oldRoundMode = manager.getRoundMode();

    manager.setRoundMode(RoundMode.ROUND_NEAR);

    // FIRST STAGE
    DoubleIntervalMatrix X = null;
    if (INTERVAL_FIRST_SECOND_STAGE) {
      final double realmin = 2.225073858507201e-308;
      final DoubleIntervalMatrix C = new DoubleIntervalMatrix(R.createUnit(n)).subtract(RA);
      DoubleIntervalMatrix Y = Z;
      final DoubleIntervalMatrix E = new DoubleIntervalMatrix(Y.getRadius().multiply(0.1)).multiply(new DoubleIntervalNumber(-1, 1))
          .addElementWise(DoubleIntervalNumber.createWithMiddleRadius(0, 10 * realmin));

      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(xs);
        manager.setRoundMode(oldRoundMode);
        return X;
      }
    }

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

    manager.setRoundMode(RoundMode.ROUND_DOWN);

    DoubleMatrix u = A.multiply(v);

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

      manager.setRoundMode(RoundMode.ROUND_UP);

      DoubleMatrix w = (A.unaryMinus().divideElementWise(u.multiply(A.createOnes(1, n)))).maxColumnWise().transpose();
      DoubleMatrix 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());
      DoubleMatrix d = B.diagonalToVector();
      DoubleMatrix alpha = dAc.add(A.createOnes(d.getRowSize(), d.getColumnSize()).unaryMinus().divideElementWise(d));

      int localk = b.getColumnSize();
      DoubleIntervalMatrix Ixs = new DoubleIntervalMatrix(xs);

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

      }
    }

    manager.setRoundMode(oldRoundMode);
    return X;
  }

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