Java Code Examples for org.apache.commons.math3.complex.Complex#multiply()

public Complex response(final double normalizedFrequency) {
    final double w = 2 * Math.PI * normalizedFrequency;
    final Complex czn1 = ComplexUtils.polar2Complex(1., -w);
    final Complex czn2 = ComplexUtils.polar2Complex(1., -2 * w);
    Complex ch = new Complex(1);
    Complex cbot = new Complex(1);

    for (int i = 0; i < mNumBiquads; i++) {
        final Biquad stage = mBiquads[i];

        Complex ct = new Complex(stage.getB0() / stage.getA0()); // NOPMD
        ct = ct.add(czn1.multiply(stage.getB1() / stage.getA0()));
        ct = ct.add(czn2.multiply(stage.getB2() / stage.getA0()));

        Complex cb = new Complex(1); // NOPMD
        cb = cb.add(czn1.multiply(stage.getA1() / stage.getA0()));
        cb = cb.add(czn2.multiply(stage.getA2() / stage.getA0()));

        ch = ch.multiply(ct);
        cbot = cbot.multiply(cb);
    }

    return ch.divide(cbot);
}
 

/**
 * Compute Zernike moments at specified order.
 *
 * @param w Width of the bounding box of the shape.
 * @param h Height of the bounding box of the shape.
 * @param n 1st order of the moment.
 * @param m 2nd order of the moment.
 *
 * @return Zernike moment of data in f[][].
 */
public static Complex calculateZernikeMoment(double[][] f, int w, int h, int n, int m){
    int diff = n-Math.abs(m);
    if ((n<0) || (Math.abs(m) > n) || (diff%2!=0)){
        throw new IllegalArgumentException("zer_mom: n="+n+", m="+m+", n-|m|="+diff);
    }

    final double c = -1;
    final double d = 1;


    ZernikeBasisFunction zernike = new ZernikeBasisFunction(n,m);
    Complex res = new Complex(0.0, 0.0);
    for (int i=0;i<w;i++){
        for (int j=0;j<h;j++) {
            Complex v = new Complex(c+(i*(d-c))/(w-1), d-(j*(d-c))/(h-1));
            res = res.add(zernike.value(v).conjugate().multiply(f[i][j]));
        }
    }

    return res.multiply((n+1)/Math.PI);
}
 

private ComplexPair transform(Complex c) {
	if (c.isInfinite()) {
		return new ComplexPair(new Complex(-1), new Complex(1));
	}

	c = ((new Complex(1)).add(c)).divide((new Complex(1)).subtract(c)); // bilinear

	Complex v = new Complex(0);
	v = MathSupplement.addmul(v, 4 * (b2 * (a2 - 1) + 1), c);
	v = v.add(8 * (b2 * (a2 - 1) - 1));
	v = v.multiply(c);
	v = v.add(4 * (b2 * (a2 - 1) + 1));
	v = v.sqrt();

	Complex u = v.multiply(-1);
	u = MathSupplement.addmul(u, ab_2, c);
	u = u.add(ab_2);

	v = MathSupplement.addmul(v, ab_2, c);
	v = v.add(ab_2);

	Complex d = new Complex(0);
	d = MathSupplement.addmul(d, 2 * (b - 1), c).add(2 * (1 + b));

	return new ComplexPair(u.divide(d), v.divide(d));
}
 

@Test
public void testButterworthResponse() {
    final int filterOrder = 4;
    final Butterworth iirLowPass = new Butterworth();
    iirLowPass.lowPass(filterOrder, 1.0, F_CUT_LOW);
    final DataSet magLowPass = filterAndGetMagnitudeSpectrum(iirLowPass, demoDataSet);
    assertThat("low-pass cut-off", magLowPass.getValue(DIM_X, F_CUT_LOW), lessThan(-3.0 + EPSILON_DB));
    assertThat("low-pass rejection", magLowPass.getValue(DIM_X, 10 * F_CUT_LOW), lessThan(-3.0 + EPSILON_DB - 20 * filterOrder));
    assertThat("low-pass pass-band ripple", getRange(magLowPass, 0, magLowPass.getIndex(DIM_X, F_CUT_LOW * 0.1)), lessThan(EPSILON_DB));

    // response
    assertEquals(filterOrder / 2, iirLowPass.getNumBiquads());
    Complex cornerFrequency = new Complex(1, 0);
    for (int i = 0; i < iirLowPass.getNumBiquads(); i++) {
        cornerFrequency = cornerFrequency.multiply(iirLowPass.getBiquad(i).response(F_CUT_LOW));
    }
    final double valueAtCutOff = 20 * TMathConstants.Log10(cornerFrequency.abs());
    assertThat("low-pass cut-off (response calculation)", valueAtCutOff, lessThan(-3.0 + EPSILON_DB));
}
 

@Override
public CharacteristicFunction apply(final double time) {
	final double logDiscountFactorForForward = this.getLogDiscountFactorForForward(time);
	final double logDiscountFactorForDiscounting = this.getLogDiscountFactorForDiscounting(time);

	return new CharacteristicFunction() {
		@Override
		public Complex apply(final Complex argument) {
			final Complex iargument = argument.multiply(Complex.I);
			final Complex denominator = ((Complex.ONE).subtract(iargument.multiply(theta*nu))).add(argument.multiply(argument).multiply(0.5*sigma*sigma*nu));
			final Complex firstLevyExponent = (((Complex.ONE).divide(denominator)).log()).multiply(time/nu);
			final Complex compensator =  iargument.multiply(time/nu * Math.log(1/(1.0-theta*nu-0.5*sigma*sigma*nu)));

			return (firstLevyExponent.subtract(compensator)
					.add(iargument.multiply(Math.log(initialValue)-logDiscountFactorForForward))
					.add(logDiscountFactorForDiscounting))
					.exp();
		}
	};
}
 

@Override
public CharacteristicFunction apply(final double time) {
	final double logDiscountFactorForForward = this.getLogDiscountFactorForForward(time);
	final double logDiscountFactorForDiscounting = this.getLogDiscountFactorForDiscounting(time);

	return new CharacteristicFunction() {
		@Override
		public Complex apply(final Complex argument) {
			final Complex iargument = argument.multiply(Complex.I);
			final Complex denominator = ((Complex.ONE).subtract(iargument.multiply(theta*nu))).add(argument.multiply(argument).multiply(0.5*sigma*sigma*nu));
			final Complex firstLevyExponent = (((Complex.ONE).divide(denominator)).log()).multiply(time/nu);
			final Complex compensator =  iargument.multiply(time/nu * Math.log(1/(1.0-theta*nu-0.5*sigma*sigma*nu)));

			return (firstLevyExponent.subtract(compensator)
					.add(iargument.multiply(Math.log(initialValue)-logDiscountFactorForForward))
					.add(logDiscountFactorForDiscounting))
					.exp();
		}
	};
}
 

@Override
public CharacteristicFunction apply(final double time) {
	final double logDiscountFactorForForward		= this.getLogDiscountFactorForForward(time);
	final double logDiscountFactorForDiscounting	= this.getLogDiscountFactorForDiscounting(time);

	return new CharacteristicFunction() {
		@Override
		public Complex apply(final Complex argument) {
			final Complex iargument = argument.multiply(Complex.I);
			return	iargument
					.multiply(
							iargument
							.multiply(0.5*volatility*volatility*time)
							.add(Math.log(initialValue)-0.5*volatility*volatility*time-logDiscountFactorForForward))
					.add(logDiscountFactorForDiscounting)
					.exp();
		}
	};
}
 

public SV otherSide(double r, double x, double g, double b, double ratio) {
    Complex z = new Complex(r, x); // z=r+jx
    Complex y = new Complex(g, b); // y=g+jb
    Complex s1 = new Complex(p, q); // s1=p1+jq1
    Complex u1 = ComplexUtils.polar2Complex(u, Math.toRadians(a));
    Complex v1 = u1.divide(Math.sqrt(3f)); // v1=u1/sqrt(3)

    Complex v1p = v1.multiply(ratio); // v1p=v1*rho
    Complex i1 = s1.divide(v1.multiply(3)).conjugate(); // i1=conj(s1/(3*v1))
    Complex i1p = i1.divide(ratio); // i1p=i1/rho
    Complex i2 = i1p.subtract(y.multiply(v1p)); // i2=i1p-y*v1p
    Complex v2 = v1p.subtract(z.multiply(i2)); // v2=v1p-z*i2
    Complex s2 = v2.multiply(3).multiply(i2.conjugate()); // s2=3*v2*conj(i2)

    Complex u2 = v2.multiply(Math.sqrt(3f));
    return new SV(-s2.getReal(), -s2.getImaginary(), u2.abs(), Math.toDegrees(u2.getArgument()));
}
 

private static Complex transform(final Complex in, final double f) {
    if (in.isInfinite()) {
        return new Complex(1, 0);
    }
    Complex c;
    // frequency transform
    c = in.multiply(f);

    // bilinear high pass transform
    return new Complex(-1).multiply(new Complex(1).add(c)).divide(new Complex(1).subtract(c));
}
 

@Override
public CharacteristicFunction apply(final double time) {
	final double logDiscountFactorForForward		= this.getLogDiscountFactorForForward(time);
	final double logDiscountFactorForDiscounting	= this.getLogDiscountFactorForDiscounting(time);
	final double transformedMean =  jumpSizeMean - 0.5 * jumpSizeStdDev*jumpSizeStdDev;
	return new CharacteristicFunction() {
		@Override
		public Complex apply(final Complex argument) {
			final Complex iargument = argument.multiply(Complex.I);

			final Complex exponent = (iargument.multiply(transformedMean))
					.add(iargument.multiply(iargument.multiply(jumpSizeStdDev*jumpSizeStdDev/2.0)));

			final Complex jumpTransform = ((exponent.exp()).subtract(1.0)).multiply(jumpIntensity*time);

			final double jumpTransformCompensator = jumpIntensity*time*(Math.exp(transformedMean+jumpSizeStdDev*jumpSizeStdDev/2.0)-1.0);

			return	iargument
					.multiply(
							iargument
							.multiply(0.5*volatility*volatility*time)
							.add(Math.log(initialValue)-0.5*volatility*volatility*time-logDiscountFactorForForward))
					.add(logDiscountFactorForDiscounting).add(jumpTransform.subtract(jumpTransformCompensator))
					.exp();
		}
	};

}
 

private Complex transform(Complex c) {
	if (c.isInfinite())
		return new Complex(-1, 0);

	// frequency transform
	c = c.multiply(f);

	Complex one = new Complex(1, 0);

	// bilinear low pass transform
	return (one.add(c)).divide(one.subtract(c));
}
 

private Complex[] multiply(double[] t, Complex c) {
    Complex[] ret = new Complex[t.length];
    for (int i = 0; i < ret.length; i++) {
        ret[i] = c.multiply(t[i]);
    }
    return ret;
}
 

@Override
public CharacteristicFunction apply(final double time) {

	final double logDiscountFactorForForward		= this.getLogDiscountFactorForForward(time);
	final double logDiscountFactorForDiscounting	= this.getLogDiscountFactorForDiscounting(time);

	return new CharacteristicFunction() {
		@Override
		public Complex apply(final Complex argument) {

			final Complex iargument = argument.multiply(Complex.I);

			final Complex gamma = iargument.multiply(rho * xi).subtract(kappa).pow(2)
					.subtract(
							iargument.multiply(iargument)
							.add(iargument.multiply(-1)).multiply(0.5)
							.multiply(2 * xi * xi))
					.sqrt();

			final Complex a = iargument
					.multiply(rho * xi)
					.subtract(kappa)
					.subtract(gamma).multiply((-theta*kappa * time) / (xi * xi))
					.subtract(iargument.multiply(rho * xi).subtract(kappa).subtract(gamma)
							.multiply(new Complex(1).divide(gamma.multiply(time).exp()).subtract(1).divide(gamma))
							.multiply(0.5).add(new Complex(1).divide(gamma.multiply(time).exp())).log()
							.add(gamma.multiply(time)).multiply((2 * theta*kappa) / (xi * xi)));

			final Complex b = iargument.multiply(iargument).add(iargument.multiply(-1)).multiply(-1)
					.divide(iargument.multiply(rho * xi).subtract(kappa)
							.add(gamma.multiply(new Complex(1).divide(gamma.multiply(time).exp()).add(1)
									.divide(new Complex(1).divide(gamma.multiply(time).exp()).subtract(1)))));

			return a.add(b.multiply(volatility*volatility)).add(iargument.multiply(Math.log(initialValue) - logDiscountFactorForForward)).add(logDiscountFactorForDiscounting).exp();
		}
	};
}
 

@Override
public Complex apply(final Complex argument) {
	final Complex iargument = argument.multiply(Complex.I);
	final Complex exponent = (iargument).add(1);
	final Complex numerator = (new Complex(strike)).pow(exponent);
	final Complex denominator = (argument.multiply(argument)).subtract(iargument);

	return numerator.divide(denominator).negate();
}
 

public Complex response(final double normalizedFrequency) {
    final double a0 = getA0();
    final double a1 = getA1();
    final double a2 = getA2();
    final double b0 = getB0();
    final double b1 = getB1();
    final double b2 = getB2();

    final double w = 2 * Math.PI * normalizedFrequency;
    final Complex czn1 = ComplexUtils.polar2Complex(1., -w);
    final Complex czn2 = ComplexUtils.polar2Complex(1., -2 * w);
    Complex ch = new Complex(1);
    Complex cbot = new Complex(1);

    Complex ct = new Complex(b0 / a0);

    ct = ct.add(czn1.multiply(b1 / a0));
    ct = ct.add(czn2.multiply(b2 / a0));

    Complex cb = new Complex(1);
    cb = cb.add(czn1.multiply(a1 / a0));
    cb = cb.add(czn2.multiply(a2 / a0));

    ch = ch.multiply(ct);
    cbot = cbot.multiply(cb);

    return ch.divide(cbot);
}
 

private static Complex transform(final Complex in, final double f) {
    if (in.isInfinite()) {
        return new Complex(-1, 0);
    }
    Complex c;
    // frequency transform
    c = in.multiply(f);

    final Complex one = new Complex(1, 0);

    // bilinear low pass transform
    return one.add(c).divide(one.subtract(c));
}
 

@Override
public Complex apply(final Complex argument) {
	final Complex iargument = argument.multiply(Complex.I);
	final Complex exponent = iargument.add(1.0);
	final Complex numerator = (new Complex(strike)).pow(exponent.subtract(1.0)).multiply(exponent);
	final Complex denominator = (argument.multiply(argument)).subtract(iargument);

	return numerator.divide(denominator);
}
 

/**
 * Complex Paul wavelet function
 *
 * @param x input
 * @param m input
 * @return complex Paul wavelet function
 */
public Complex Paul(final double x, final int m) {
    final double val = Math.pow(2, m) * TMath.Factorial(m)
                       / Math.sqrt(TMathConstants.Pi() * TMath.Factorial(2 * m));
    Complex c1 = new Complex(1, 0);
    Complex c2 = new Complex(1, 0);
    for (int i = 0; i < m + 1; i++) {
        c1 = c1.multiply(new Complex(1, -x));
        if (i < m) {
            c2 = c2.multiply(new Complex(0, 1));
        }
    }
    return c1.multiply(c2).multiply(val);
}
 

private double gotoLimit2(AclfNetwork net) throws InterpssException {
		
		int busi = 15;
		System.out.println("============ Go to Limit ============");
		LoadflowAlgorithm algo = CoreObjectFactory.createLoadflowAlgorithm(net);
		algo.loadflow();
//		LFOut.showlf(net);
		
		
		AclfBus bus = net.getBusList().get(busi);
		Complex base = bus.getLoadPQ();
		oBase = bus.getLoadPQ();
		Complex delta = base.multiply(100);
		Complex result = base;

//		System.out.println("base = "+base+", delta = "+delta);
		bus.getContributeLoadList().get(0).setLoadCP(base.add(delta));
		System.out.println("\tbus 16 load pq = "+bus.getLoadPQ());
		
		int iter = 0;
		while(delta.abs() > 1e-10) {
			iter += 1;
			//��ǰ��������
//			System.out.println("base = "+base+", delta = "+delta);
			
			//��̽�е�
			delta = delta.multiply(0.5);
			bus.getContributeLoadList().get(0).setLoadCP(base.add(delta));
			//�����
//			System.out.println("bus 9 load pq = "+bus.getLoadPQ());
			
			//����̽�ɹ���...
			if (algo.loadflow()) {
				base = base.add(delta);
				result = base;
//				System.out.println("a success bus 9 load PQ = "+bus.getLoadPQ());
			}
//			System.out.println(AclfOutFunc.loadFlowSummary(net));
		}
		
//		System.out.println("iter = "+iter);
		bus.getContributeLoadList().get(0).setLoadCP(result);
		System.out.println("after searching, bus 16 load PQ = "+result+", lf result = "+algo.loadflow());
//		System.out.println(AclfOutFunc.loadFlowSummary(net));
		
		return result.divide(oBase).abs();
	}
 

private void gotoLimit(AclfNetwork net) throws InterpssException {
		
		int busi = 9;
		System.out.println("============ Go to Limit ============");
		LoadflowAlgorithm algo = CoreObjectFactory.createLoadflowAlgorithm(net);
		
		AclfBus bus = net.getBusList().get(9);
		Complex base = bus.getLoadPQ();
		Complex load0 = bus.getLoadPQ();
		Complex delta = base.multiply(10);
		Complex result = base;

		System.out.println("base = "+base+", delta = "+delta);
		bus.getContributeLoadList().get(0).setLoadCP(base.add(delta));
		System.out.println("\tbus load pq = "+bus.getLoadPQ());
		
		int iter = 0;
		while(delta.abs() > 1e-10) {
			iter += 1;
			//��ǰ��������
			System.out.println("base = "+base+", delta = "+delta);
			
			//��̽�е�
			delta = delta.multiply(0.5);
			bus.getContributeLoadList().get(0).setLoadCP(base.add(delta));
			//�����
			System.out.println("bus load pq = "+bus.getLoadPQ());
			
			//����̽�ɹ���...
			if (algo.loadflow()) {
				base = base.add(delta);
				result = base;
				System.out.println("a success bus 9 load PQ = "+bus.getLoadPQ());
			}
//			System.out.println(AclfOutFunc.loadFlowSummary(net));
		}
		
		System.out.println("iter = "+iter);
		bus.getContributeLoadList().get(0).setLoadCP(result);
		System.out.println("after searching, bus 9 load PQ = "+result+", lf result = "+algo.loadflow()+", result/load0 = "+result.divide(load0));
		System.out.println(AclfOutFunc.loadFlowSummary(net));
		System.out.println(AclfOutFunc.loadLossAllocation(net));
		LFOut.showlf(net);
	}