Ccpair.java

/*
 * $Id: Ccpair.java,v 1.36 2008/07/16 15:40:00 koga Exp $
 * 
 * Copyright (C) 2004 Koga Laboratory. All rights reserved.
 */
package org.mklab.tool.matrix;

import org.mklab.nfc.matrix.DoubleComplexMatrix;
import org.mklab.nfc.matrix.DoubleMatrix;
import org.mklab.nfc.matrix.IndexedMatrix;
import org.mklab.nfc.matrix.IntMatrix;
import org.mklab.nfc.scalar.DoubleComplexNumber;
import org.mklab.nfc.scalar.DoubleNumber;


/**
 * 複素共役の対が並ぶようにソートされた成分からなる行列を求めるクラスです。
 * 
 * <p>Make complex conjugate pairs
 * 
 * @author koga
 * @version $Revision: 1.36 $
 */
public class Ccpair {

  /**
   * 複素ベクトル <code>x</code> の成分について、複素共役の対が並ぶようにソートします。 複素共役対は、実部が昇順になるよう並べられます。実数(虚部がゼロ)は、 複素共役対の後ろに並べなれる。
   * 
   * @param x_ データ
   * @return 複素共役の対が並ぶベクトル (complex conjugate pairs)
   */
  public static DoubleComplexMatrix ccpair(DoubleComplexMatrix x_) {
    DoubleNumber unit = x_.getElement(1, 1).abs().getRealPart().createUnit();
    return ccpair(x_, unit.getMachineEpsilon().multiply(100));
  }

  /**
   * 値を比較するための許容誤差に <code>tolerance</code> を用いる。
   * 
   * @param x_ データ
   * @param tolerance 許容誤差
   * @return 複素共役の対が並ぶべくとる (complex conjugate pairs)
   */
  public static DoubleComplexMatrix ccpair(DoubleComplexMatrix x_, double tolerance) {
    return ccpair(x_, new DoubleNumber(tolerance));
  }

  /**
   * 値を比較するための許容誤差に <code>tolerance</code> を用いる。
   * 
   * @param x_ データ
   * @param tolerance 許容誤差
   * @return 複素共役の対が並ぶべくとる (complex conjugate pairs)
   */
  public static DoubleComplexMatrix ccpair(DoubleComplexMatrix x_, DoubleNumber tolerance) {
    DoubleNumber unit = x_.getElement(1, 1).getRealPart().createUnit();
    DoubleMatrix mat = x_.absElementWise().getRealPart();

    DoubleComplexNumber j = new DoubleComplexNumber(0,1);

    DoubleComplexMatrix x;
    DoubleComplexMatrix y = x_.createZero();
    if (x_.getRowSize() < x_.getColumnSize()) {
      x = x_.transpose();
    } else {
      x = x_.createClone();
    }

    // Find real numbers
    IntMatrix idx_re = x.getImaginaryPart().absElementWise().compareElementWise(".<=", x.absElementWise().getRealPart().multiply(tolerance)).find(); //$NON-NLS-1$
    int nr = idx_re.length();

    DoubleMatrix R = null;
    if (nr != 0) {
      IndexedMatrix<DoubleNumber,DoubleMatrix> temp = x.getSubVector(idx_re).getRealPart().sort();
      R = temp.getMatrix();
      x.setSubVector(idx_re, new DoubleComplexMatrix(mat.createZero(0, 0))); // remove real numbers in x
    }
    int nc = x.length();

    if (nc == 0) {
      if (R != null && R.isReal()) {
        y = new DoubleComplexMatrix(R);
      }

      return y;
    }

    if ((nc % 2) != 0) {
      throw new RuntimeException(Messages.getString("Ccpair.0")); //$NON-NLS-1$
    }

    // Upper half plane
    DoubleComplexMatrix c = new DoubleComplexMatrix(x.getRealPart()).add(new DoubleComplexMatrix(x.getImaginaryPart().absElementWise()).multiply(j));

    // Sort by real part
    IndexedMatrix<DoubleNumber,DoubleMatrix> temp = c.getRealPart().sort();
    DoubleMatrix cr = temp.getMatrix();
    IntMatrix idx_c = temp.getIndices();

    c = new DoubleComplexMatrix(cr).add(new DoubleComplexMatrix(c.getSubVector(idx_c).getImaginaryPart()).multiply(j));

    // Real parts are at least in pair ?
    if (cr.getSubVector(IntMatrix.series(1, nc, 2)).subtract(cr.getSubVector(IntMatrix.series(2, nc, 2))).absElementWise().compareElementWise(".>", //$NON-NLS-1$
        c.getSubVector(IntMatrix.series(1, nc, 2)).absElementWise().getRealPart().multiply(tolerance)).anyTrue()) {
      throw new RuntimeException(Messages.getString("Ccpair.2")); //$NON-NLS-1$
    }

    x = x.getSubVector(idx_c);

    int k = 1;
    while (k <= nc / 2) {
      IntMatrix i1 = cr.getSubVector(IntMatrix.series(1, nc, 2)).subtractElementWise(cr.getElement(2 * k)).absElementWise().compareElementWise(".<=", //$NON-NLS-1$
          (c.getElement(2 * k)).abs().getRealPart().multiply(tolerance)).find();

      if (i1.length() == 0) {
        throw new RuntimeException(Messages.getString("Ccpair.4")); //$NON-NLS-1$
      } else if (i1.length() == 1) {
        // One pair
        if ((x.getElement(2 * k).subtract(x.getElement(2 * k - 1).conjugate())).getRealPart().isGreaterThan(tolerance.multiply(x.getElement(2 * k).abs().getRealPart()))) {
          throw new RuntimeException(Messages.getString("Ccpair.5")); //$NON-NLS-1$
        }
        DoubleComplexMatrix tmp = new DoubleComplexNumber(unit).createGrid(new DoubleComplexNumber[] {j.multiply(new DoubleComplexNumber(x.getElement(2 * k).getImaginaryPart().abs())).unaryMinus().add(new DoubleComplexNumber(x.getElement(2 * k).getRealPart()))});
        x.setSubVector(IntMatrix.series(2 * k - 1, 2 * k), tmp.appendDown(tmp.conjugate()));
      } else {
        // Multi-pairs with same real part
        IntMatrix i2 = cr.getSubVector(1, nc).subtractElementWise(cr.getElement(2 * k)).compareElementWise(".<=", tolerance.multiply(c.getElement(2 * k).abs().getRealPart())).find(); //$NON-NLS-1$

        int n = Math.max(1, i2.length());
        DoubleComplexMatrix x2 = x.getSubVector(i2);
        IndexedMatrix<DoubleNumber,DoubleMatrix> temp2 = x2.getImaginaryPart().sort();
        // DoubleMatrix xi = temp[0];
        IntMatrix xind = temp2.getIndices();
        x2 = x2.getSubVector(xind);

        if (x2.subtract(x2.getSubVector(IntMatrix.series(n, 1, -1)).conjugate()).absElementWise().getRealPart().compareElementWise(".>", x2.absElementWise().getRealPart().multiply(tolerance)).anyTrue()) { //$NON-NLS-1$
          throw new RuntimeException(Messages.getString("Ccpair.8")); //$NON-NLS-1$
        }
        x.setSubVector(IntMatrix.series(i2.getIntElement(1), i2.getIntElement(n - 1), 2), x2.getSubVector(IntMatrix.series(1, n / 2)));
        x.setSubVector(IntMatrix.series(i2.getIntElement(2), i2.getIntElement(n), 2), x2.getSubVector(IntMatrix.series(1, n / 2)).conjugate());
      }

      k = k + i1.length();
    }

    if (nc != 0) {
      y.setSubMatrix(1, nc, 1, 1, x.getSubVector(IntMatrix.series(1, nc)));
    }
    if (nr != 0) {
      y.setSubMatrix(nc + 1, nc + nr, 1, 1, new DoubleComplexMatrix(R));
    }
    return y;
  }

}