OishiMethod.java

/*
 * 作成日: 2003/08/22
 *
 * Copyright (C) 2004 Koga Laboratory. All rights reserved.
 */
package org.mklab.cga.linear;

import org.mklab.cga.interval.matrix.IntervalNumericalMatrix;
import org.mklab.cga.interval.scalar.IntervalNumericalScalar;
import org.mklab.nfc.matrix.NormType;
import org.mklab.nfc.matrix.NumericalMatrix;
import org.mklab.nfc.scalar.NumericalScalar;
import org.mklab.nfc.util.RoundMode;
import org.mklab.nfc.util.RoundModeManager;


/**
 * 線形方程式の精度保証付き解を(丸めの変更を行わない)Oishiの方法で求めるクラスです。
 * 
 * @author Hiroki
 * @param <IS> 区間スカラーの型
 * @param <IM> 区間行列の型
 * @param <S> 成分の型
 * @param <M> 行列の型
 */
public class OishiMethod<IS extends IntervalNumericalScalar<IS, IM, S, M>, IM extends IntervalNumericalMatrix<IS, IM, S, M>, S extends NumericalScalar<S, M>, M extends NumericalMatrix<S, M>> {

  /** 係数行列 */
  private IM IA;
  /** 係数行列 */
  private M A;
  /** 右辺係数ベクトル */
  private M b;
  /** Aの逆行列 */
  private M R;
  /** 連立一次方程式Ax=bの近似解 */
  private M x;
  /** Aのサイズ * */
  /** error */
  private S error;
  private M merror;
  private int length;

  //  /**
  //   * Ax = bのAとbを設定します。 コンストラクタではAx=bの近似解xが計算します。
  //   * 
  //   * @param A 係数行列
  //   * @param b 右辺係数ベクトル
  //   */
  //  public OishiMethod(final M A, final M b) {
  //    this.A = A;
  //    this.b = b;
  //    this.length = A.length();
  //
  //    if (A.getColumnSize() == A.getRowSize()) {
  //      this.R = A.inverse();
  //    } else {
  //      this.R = (A).pseudoInverse();
  //    }
  //
  //    this.x = this.R.multiply(b);
  //  }

  /**
   * Ax=bのAとbを設定します。 係数行列Aは区間行列です。コンストラクタでは近似解を計算します。
   * 
   * @param A 係数区間行列
   * @param b 右辺係数ベクトル
   */
  public OishiMethod(final IM A, M b) {
    this.IA = A;
    this.b = b;
    this.A = A.getInfimum().add(A.getSupremum()).divide(2);
    this.R = this.A.inverse();
    this.x = this.R.multiply(b);
    this.length = A.getRowSize();
  }
  
  /**
   * Creates interval matrix.
   * 
   * @param matrix matrix
   * @return interval matrix
   */
  private IM createInterval(M matrix) {
    return this.IA.create(matrix);
  }

  /**
   * 連立一次方程式の近似解とその誤差を求めます。
   */
  public void  solveError() {
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    IM i_A = this.IA;
    IM i_R = this.IA.create(this.R);

    IM i_G = ((i_R.multiply(this.IA)).subtract(this.IA.createUnit(this.length))).absElementWise();
    IM i_r = ((i_A.multiply(this.IA.create(this.x))).subtract(this.IA.create(this.b))).absElementWise();

    M ru = i_r.getSupremum();
    M rd = i_r.getInfimum();
    M Gu = i_G.getSupremum();
    M Gd = i_G.getInfimum();

    manager.setRoundMode(RoundMode.ROUND_UP);

    M r_max = (rd).maxElementWise(ru);
    M G_max = (Gd).maxElementWise(Gu);

    S R_norm = ((this.R).norm(NormType.INFINITY));
    S G_norm = ((G_max).norm(NormType.INFINITY));
    S r_norm = ((r_max).norm(NormType.INFINITY));

    manager.setRoundMode(RoundMode.ROUND_DOWN);

    S D = G_norm.unaryMinus().add(1);

    if (D.isLessThanOrEquals(0)) {
      manager.setRoundMode(oldRoundMode);
      throw new RuntimeException("false"); //$NON-NLS-1$
    }

    manager.setRoundMode(RoundMode.ROUND_UP);

    S A_inv = R_norm.divide(D);
    this.error = A_inv.multiply(r_norm);

    manager.setRoundMode(oldRoundMode);
  }

  /**
   * Returns x.
   * 
   * @return x
   */
  public M getSolution() {
    return this.x;
  }

  /**
   * Returns error.
   * 
   * @return error;
   */
  public S getError() {
    return this.error;
  }
  
  /**
   * Returns matrix error.
   * 
   * @return matrix error
   */
  public M getMatrixError() {
    return this.merror;
  }
  
  /**
   * 連立一次方程式の近似解とそのそれぞれに対する誤差を求めます。 まだ問題あり。
   */
  public void solveElementsError() {
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    IM i_R = this.IA.create(this.R);
    IM i_A = this.IA;

    IM i_G = (i_R.multiply(this.IA).subtract(this.IA.createUnit(this.length))).absElementWise();
    IM i_r = (i_A.multiply(this.IA.create(this.x)).subtract(this.IA.create(this.b))).absElementWise();

    M Gd = i_G.getInfimum();
    M Gu = i_G.getSupremum();
    M rd = i_r.getInfimum();
    M ru = i_r.getSupremum();

    M rd_abs = (rd).absElementWise();
    M ru_abs = (ru).absElementWise();

    manager.setRoundMode(RoundMode.ROUND_UP);

    M center = (ru.add(rd)).divide(2);
    M radius = center.subtract(rd);

    M r_max = (rd_abs).maxElementWise(ru_abs);
    M G_max = (Gd).maxElementWise(Gu);

    S G_norm = ((G_max).norm(NormType.INFINITY));
    S R_norm = ((this.R).norm(NormType.INFINITY));
    S r_norm = ((r_max).norm(NormType.INFINITY));

    IM i_Rru = ((i_R.multiply(createInterval(center)).add((createInterval(this.R)).absElementWise().multiply(createInterval(radius))))).absElementWise();

    manager.setRoundMode(RoundMode.ROUND_DOWN);

    S D = G_norm.unaryMinus().add(1);

    // 1 - ||G|| > 0
    if (D.isLessThanOrEquals(0)) {
      manager.setRoundMode(oldRoundMode);
      throw new RuntimeException("false"); //$NON-NLS-1$
    }

    manager.setRoundMode(RoundMode.ROUND_UP);

    M Rr = (i_Rru.getSupremum()).maxElementWise(i_Rru.getInfimum());
    this.error = R_norm.multiply(r_norm).divide(D);
    M kappa = G_max.multiply(G_max.createOnes(this.length, 1));
    M error2 = Rr.add(kappa.multiply(this.error));
    S cwerror = max(error2);
    this.error = cwerror;
    this.merror = error2;

    manager.setRoundMode(oldRoundMode);
  }

  /**
   */
  public void solveElementsError2() {
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    manager.setRoundMode(RoundMode.ROUND_DOWN);

    M Gd = (this.R.multiply(this.A).subtract(this.A.createUnit(this.length))).absElementWise();
    M rd = this.A.multiply(this.x).subtract(this.b);

    manager.setRoundMode(RoundMode.ROUND_UP);

    M Gu = (this.R.multiply(this.A).subtract(this.A.createUnit(this.length))).absElementWise();
    M ru = this.A.multiply(this.x).subtract(this.b);
    M center = (ru.add(rd)).divide(2);
    M radius = center.subtract(rd);

    M G_max = (Gd).maxElementWise(Gu);
    S G_norm = ((G_max).norm(NormType.INFINITY));
    M Rru = (this.R.multiply(center).add((this.R).absElementWise().multiply(radius))).absElementWise();

    manager.setRoundMode(RoundMode.ROUND_DOWN);

    M Rrd = (this.R.multiply(center).add((this.R).absElementWise().multiply(radius))).absElementWise();
    S D = G_norm.unaryMinus().add(1);

    // 1 - ||G|| > 0
    if (D.isLessThanOrEquals(0)) {
      manager.setRoundMode(oldRoundMode);
      throw new RuntimeException("false"); //$NON-NLS-1$
    }

    manager.setRoundMode(RoundMode.ROUND_UP);

    M Rr = (Rru).maxElementWise(Rrd);
    S Rr_norm = ((Rr).norm(NormType.INFINITY));
    this. error = Rr_norm.divide(D);
    M kappa = G_max.multiply(G_max.createOnes(this.length, 1));
    M error2 = Rr.add(kappa.multiply(this.error));
    S cwerror = max(error2);
    this.error = cwerror;
    this.merror = error2;

    manager.setRoundMode(oldRoundMode);
  }

  /**
   * 誤差ベクトルの要素の最大値を返します。
   * 
   * @param errorVector 誤差ベクトル
   * @return 最大値
   */
  private S max(M errorVector) {
    S cwmax = (errorVector).getElement(1, 1);
    for (int i = 1; i < this.A.length(); i++) {
      if (cwmax.isLessThan(errorVector.getElement(i + 1, 1))) {
        cwmax = (errorVector).getElement(i + 1, 1);
      }

    }
    return cwmax;
  }
}