KrawczykMethod.java

/*
 * Created on 2004/10/18
 *
 * $Id: KrawMethod.java,v 1.15 2008/02/02 03:07:13 koga Exp $
 * 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.BooleanMatrix;
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;


/**
 * 線形方程式の精度保証付き解をKrawczyk's method(クラウジック法)で求めるクラスです。
 * 
 * @author hiroki
 * @version $Revision: 1.15 $. 2004/10/18
 * @param <IS> 区間スカラーの型
 * @param <IM> 区間行列の型
 * @param <S> 成分の型
 * @param <M> 行列の型
 */
public class KrawczykMethod<IS extends IntervalNumericalScalar<IS, IM, S, M>, IM extends IntervalNumericalMatrix<IS, IM, S, M>, S extends NumericalScalar<S, M>, M extends NumericalMatrix<S, M>>
    implements LinearEquationVerifier<IS, IM, S, M> {

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

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

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

  //  /**
  //   * コンストラクタ
  //   * 
  //   */
  //  public KrawczykMethod() {
  //    //
  //  }
  //
  //  /**
  //   * 新しく生成された<code>KrawMethod</code>オブジェクトを初期化します。
  //   * 
  //   * @param A 係数行列
  //   * @param b 右辺係数ベクトル
  //   */
  //  public KrawczykMethod(final M A, final M b) {
  //    this.intA = IntervalMatrixFactory.toInterval(A);
  //    this.intB = IntervalMatrixFactory.toInterval(b);
  //  }
  //
  //  /**
  //   * 新しく生成された<code>KrawMethod</code>オブジェクトを初期化します。
  //   * 
  //   * @param A 係数行列
  //   * @param b 右辺係数区間ベクトル
  //   */
  //  public KrawczykMethod(final M A, final IM b) {
  //    this.intA = IntervalMatrixFactory.toInterval(A);
  //    this.intB = b.clone();
  //  }
  //
  //  /**
  //   * 新しく生成された<code>KrawMethod</code>オブジェクトを初期化します。
  //   * 
  //   * @param A 係数区間行列
  //   * @param b 右辺係数ベクトル
  //   */
  //  public KrawczykMethod(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>KrawMethod</code>オブジェクトを初期化します。
   * 
   * @param A 係数行列
   * @param b 右辺係数ベクトル
   */
  public KrawczykMethod(final IM A, final IM b) {
    this.intA = A;
    this.intB = b;
  }

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

    M midA = this.intA.getMiddle();// 区間行列Aの中心
    IM C = createInterval(midA.inverse());

    IM bb = C.multiply(this.intB);// 近似解
    IM CA = C.multiply(this.intA);// inv(midA)*A

    int n = this.intA.getColumnSize();
    IM AA = CA.createUnit(n).subtract(CA);// I-CA
    S beta = AA.norm(NormType.INFINITY).getSupremum();

    IM x = this.intA.createZero(n, 1);

    if (beta.isGreaterThanOrEquals(1)) {
      throw new RuntimeException("Algorithm not suitable for this A"); //$NON-NLS-1$
    }

    IS type = this.intA.getElement(1, 1);
    S alpha = bb.norm(NormType.INFINITY).getSupremum().divide(beta.unaryMinus().add(1)); // (1 - beta);
    IS intervalAlpha = type.createInfSup(alpha.unaryMinus(), alpha);

    for (int i = 0; i < n; i++) {
      x.setElement(i + 1, 1, intervalAlpha);
    }

    final S unit = alpha.createUnit();

    S sOld = unit.getInfinity();
    M rad = x.getSupremum().subtract(x.getInfimum()).divide(2);

    S s = unit.createZero();
    for (int j = 0; j < n; j++) {
      s = s.add(rad.getElement(j + 1, 1));
    }

    S mult = beta.add(1).divide(2); // (1 + beta) / 2;

    S ss = unit.createZero();
    while (s.isLessThan(mult.multiply(sOld))) {
      x = intersect(bb.add(AA.multiply(x)), x);
      sOld = s;

      M radius = x.getSupremum().subtract(x.getInfimum()).divide(2);
      for (int i = 0; i < x.getRowSize(); i++) {
        ss = ss.add(radius.getElement(i + 1, 1));
      }

      s = ss;
    }

    this.solution = x;
  }

  /**
   * Creates interval matrix.
   * 
   * @param matrix matrix
   * @return interval marix
   */
  private IM createInterval(M matrix) {
    return this.intA.create(matrix);
  }

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

  /**
   * AとBの共通の要素を返します。共通でない要素のところにはNaNを入れて返します。
   * 
   * @param A 区間行列
   * @param B 区間行列
   * @return 共通要素を持った区間行列
   */
  public IM intersect(final IM A, final IM B) {
    // TODO 抽象化したCPUの丸めモード指定と行列を使う
    RoundModeManager manager = RoundModeManager.getManager();
    RoundMode oldRoundMode = manager.getRoundMode();

    M inf = (A.getInfimum()).maxElementWise(B.getInfimum());
    M sup = (A.getSupremum()).minElementWise(B.getSupremum());

    BooleanMatrix index = inf.compareElementWise(".>", sup).orElementWise(A.getInfimum().isNanElementWise()).orElementWise(A.getSupremum().isNanElementWise()).orElementWise( //$NON-NLS-1$
        B.getInfimum().isNanElementWise()).orElementWise(B.getSupremum().isNanElementWise());
    if (index.anyTrue()) {
      S unit = A.getMiddle().getElement(1);

      for (int i = 0; i < index.getRowSize(); i++) {
        for (int j = 0; j < index.getColumnSize(); j++) {
          if (index.getElement(i + 1, j + 1)) {
            (inf).setElement(i + 1, j + 1, unit.getNaN());
            (sup).setElement(i + 1, j + 1, unit.getNaN());
          }
        }
      }
    }

    manager.setRoundMode(oldRoundMode);
    return A.createInfSup(inf, sup);
  }

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