zmath.java

Java实现的各种数学算法

    if (dmin < dmax) singular = true;
    /****************
    if (dmin < dmax) {
      System.out.println("singular");
      v[1] = 0.f;
      return;
    }
    ****************/

    // diagonal matrix, eigenalues are the diagonal elements
    float babs = (b >= 0.f) ? b : (-b);
    if (babs < dmax) {
      if (a > c) {
	v[0] = a;
	v[1] = c;
	// think this is right, but be careful
	if (e != null) {  e[0] = 1.f;  e[1] = 0.f; }
      }
      else {
	v[0] = c;
	v[1] = a;
	if (e != null) {  e[0] = 0.f;  e[1] = 1.f; }
      }
      return !singular;
    }

    double det = a * c - b * b;
    double apc = a + c;
    double discr = apc * apc - 4. * det;
    if (!(discr >= 0.)) {
      System.out.println("\nproblem.  eigvalues abc = "+a+" "+b+" "+c);
      System.out.println("det = "+det);
      zliberror._assert(discr >= 0., "discr < 0: "+discr);
    }
    double sqrt = Math.sqrt(discr);
    v[0] = (float)((apc + sqrt) / 2.);
    v[1] = (float)((apc - sqrt) / 2.);
    zliberror._assert(v[0] >= v[1], ("eigs wrong order "+v[0]+" "+v[1]));

    if (e != null) {
      e[0] = 2.f * b;
      e[1] = (float)((c - a) + sqrt);
      double len = Math.sqrt(e[0]*e[0] + e[1]*e[1]);
      e[0] = (float)(e[0] / len);
      e[1] = (float)(e[1] / len);
    }

    return !singular;
  } //eigenvalues22



  // matlab: m = [0.9, 0.35; 0.35, 0.4];  eig(m)
  static void testEigenValues22()
  {
    float[][] m = new float[2][2];
    m[0][0] = 0.9f;
    m[0][1] = m[1][0] = 0.35f;
    m[1][1] = 0.4f;

    float[] l = new float[2];
    float[] e = new float[2];
    boolean singular = eigvalues22(m, l, e);

    matrix.print(new PrintWriter(System.out), "m=", m);
    System.out.println("l1="+l[0]+"  l2="+l[1] + " singular="+!singular);
    System.out.println("e1="+e[0]+"  e2="+e[1]);

    // second test: diagonal matrix, eigenvalues are the diagonal elements
    m[0][0] = 0.3f;
    m[0][1] = m[1][0] = 0.0f;
    m[1][1] = 0.4f;
    singular=eigvalues22(m, l, e);
    matrix.print(new PrintWriter(System.out), "m=", m);
    System.out.println("l1="+l[0]+"  l2="+l[1] + " singular="+!singular);
    System.out.println("e1="+e[0]+"  e2="+e[1]);

    // third test
    m[0][0] = 0.99f;
    m[0][1] = m[1][0] = 0.2f; 
    m[1][1] = 0.001f;
    singular=eigvalues22(m, l, e);
    matrix.print(new PrintWriter(System.out), "m=", m);
    System.out.println("l1="+l[0]+"  l2="+l[1] + " singular="+!singular);
    System.out.println("e1="+e[0]+"  e2="+e[1]);

  } //testEigenValues22


  /**
   * Return eigenvalues of a 3x3 SYMMETRIC matrix in v[].
   * Return value is false if the matrix is near singular.
   * Tested, smallest eigenvalue is correct at least.
   */
  public static final boolean eigvalues33(final float[][] M, double[] lambda)
  {
    float trace = M[0][0] + M[1][1] + M[2][2];
    float trace3 = trace / 3.f;
    M[0][0] -= trace3;
    M[1][1] -= trace3;
    M[2][2] -= trace3;

    float a = M[0][0];
    float b = M[1][1];
    float c = M[2][2];
    float d = M[0][1];
    float e = M[0][2];
    float f = M[1][2];

    double p = a*b + a*c + b*c - d*d - e*e - f*f;
    double q = a*f*f + b*e*e + c*d*d - 2.*d*e*f - a*b*c;
    zliberror._assert(p <= 0.);
    //System.out.println("p="+p);
    //System.out.println("q="+q);

    double epsilon = 0.00001;	// Double.MIN_VALUE.  in case p,q are both 0
    double beta = Math.sqrt(-4.*p / 3.);
    double alphaarg = 3.*q / (p*beta - epsilon);
    if (alphaarg > 1.) alphaarg = 1.;
    double alpha = (1./3.)*Math.acos(alphaarg);
    double twopio3 = 2.*Math.PI/3.;
    //System.out.println("alphaarg="+alphaarg);
    //System.out.println("alpha="+alpha);

    lambda[0] = (float)(beta * Math.cos(alpha));
    lambda[1] = (float)(beta * Math.cos(alpha - twopio3));
    lambda[2] = (float)(beta * Math.cos(alpha + twopio3));

    return true;	// todo
  } //eigvalues33


  /**
   * this matrix fails with alphaarg slightly greater than 1,
   * causing alpha to be NaN.
   * Thus check for alphaarg > 1 above
   */
  private static void eigvalues33fails()
  {
    float[][] m = new float[3][3];
    m[0][0] = -59.2151f;
    m[0][1] = -4.6064f;
    m[0][2] = 14.0243f;

    m[1][0] = -4.6064f;
    m[1][1] = -42.9025f;
    m[1][2] = -53.3784f;

    m[2][0] = 14.0243f;
    m[2][1] = -53.3784f;
    m[2][2] = 102.1176f;

    double[] l = new double[3];
    zmath.eigvalues33(m, l);

    if (ArrUtils.hasNaN(l) >= 0) {
      matrix.print("got nan here: m=", m);
      matrix.print("got nan here: l=", l);
      System.exit(1);
    }
  } //eigvalues33fails


  /** W.Clocksin implemented this version.  TODO double check
   public static double gaussEval(int nd, double [] data, double [] means, double [][] covar)
   {
     int i;
     double s = 0, a = 1;
     double [] arg;
     double [] x = new double[nd];
     for (i = 0; i < nd; i++) x[i] = data[i] - means[i];
     try
     {
       a = Math.pow(2.0 * Math.PI, nd / 2.0) * Math.sqrt(DoubleMatrix.determinant(covar));
       arg = DoubleMatrix.multiply(x, DoubleMatrix.inverse(covar));
       s = 0.0;
       for (i = 0; i < nd; i++) s += arg[i] * x[i];
     } catch (MathException e) {};
     return  Math.exp(-0.5 * s) / a;
   }
  */


  /**
   * evaluate multi-dimensional gaussian.
   */
  public static final double gaussEval(double[] data,
				       double[] mean,
				       double[][] invcovariance)
    throws Exception
  {
    zliberror._assert(mean.length == data.length);

    int nd = mean.length;
    double det = DoubleMatrix.determinant(invcovariance);
    double f1 = Math.pow(2.*Math.PI, nd/2.);
    double f2 = 1. / Math.sqrt(det);  // det(inverse) = 1/det(original)
    double scale = 1. / (f1 * f2);

    return gaussEval(data, mean, invcovariance, scale);
  } //gaussEval


  /**
   * Return the scale part of the gaussian - requires sqrt,pow,
   * detminant of covariance.  If the probability of a bunch of points
   * wrt a fixed gaussian is being evaluated it is faster to evaluate
   * the scale outside the loop.
   */
  public static final double gaussEvalScale(double[] mean,
					    double[][] invcovariance)
    throws Exception
  {
    int nd = mean.length;
    double det = DoubleMatrix.determinant(invcovariance);
    double f1 = Math.pow(2.*Math.PI, nd/2.);
    double f2 = 1. / Math.sqrt(det);  // det(inverse) = 1/det(original)
    double scale = 1. / (f1 * f2);

    return scale;
  } //gaussEvalScale


  /**
   * evaluate multi-dimensional gaussian.
   */
  public static final double gaussEval(double[] data,
				       double[] mean,
				       double[][] invcovariance,
				       double scale)
  {
    int nd = data.length;
    double sum = 0.;

    for( int r=0; r < nd; r++ ) {
      double rsum = 0.;
      for( int c=0; c < nd; c++ ) {
	rsum += invcovariance[r][c] * (data[c] - mean[c]);
      }
      sum += (rsum * (data[r] - mean[r]));
    }

    return scale * Math.exp(-0.5 * sum);
  } //gaussEval


  /**
   * quadratic, float version
   * @return number of real roots (0,1,2)
   */
  public static int quadsolve(float a, float b, float c, float[] roots)
  {
    float discrim = b*b - 4.f*a*c;

    if (discrim < 0.f)
      return 0;
    else if (discrim == 0.f) {
      roots[0] = -b/(2.f*a);
      return 1;
    }
    else {
      discrim = (float)Math.sqrt(discrim);
      roots[0] = (-b + discrim)/(2.f*a);
      roots[1] = (-b - discrim)/(2.f*a);
      return 2;
    }
  } //quadsolve


  /**
   * quadratic
   * @return number of real roots (0,1,2)
   */
  public static int quadsolve(double a, double b, double c, double[] roots)
  {
    double discrim = b*b - 4.*a*c;

    if (discrim < 0.)
      return 0;
    else if (discrim == 0.) {
      roots[0] = -b/(2.*a);
      return 1;
    }
    else {
      discrim = Math.sqrt(discrim);
      roots[0] = (-b + discrim)/(2.*a);
      roots[1] = (-b - discrim)/(2.*a);
      return 2;
    }
  } //quadsolve

  //----------------------------------------------------------------

  // some code wants that name pythag also
   public final static double hypot(double a, double b) {
     return Math.sqrt(a*a + b*b);
   }

   /**
    * sqrt(a^2 + b^2) without under/overflow.
    * this approach recommended by e.g. jama
    */
  public final static double Hypot(double a, double b)
  {
     double r;
     if (Math.abs(a) > Math.abs(b)) {
       r = b/a;
       r = Math.abs(a)*Math.sqrt(1.+r*r);
     } else if (b != 0.) {
       r = a/b;
       r = Math.abs(b)*Math.sqrt(1.+r*r);
     } else {
       r = 0.0;
     }
     return r;
   }

  // was called len
  public final static double dist(double ix, double iy, double lx, double ly)
  {
      double dx2 = ix - lx;  dx2 = dx2*dx2;
      double dy2 = iy - ly;  dy2 = dy2*dy2;
      double len2 = dx2 + dy2;
      return Math.sqrt(len2);
  } //dist


  

  //================ TEST ================

  public static void main(String[] args)
  {
    //testBinom();
    //testEigenValues22();
    //testInv2x2();
    testBeta();
  } //main

} //zmath