"""
Utility functions for complex analysis of plane geometry.
Uses and generates (numpy) complex numbers.
"""
import numpy as np
import ddg.math.euclidean2d as euclidean2d
###########
# sampling
###########
[docs]def complex_sample(size=None, seed=None):
"""
Get a sample of random complex numbers or numpy.ndarray.
Parameters
----------
size : array_like of ints (default=None)
Shape of numpy.ndarray.
If `size=None` then a single `complex` will be returned.
seed : int (default=None)
Seed for obtaining a deterministic function.
Returns
-------
complex or numpy.ndarray of numpy.complex128
If `size=None` then `complex` will be returned,
else `numpy.ndarray` of shape `size`.
Notes
-----
Real and imaginary part are standard normally distributed
Examples
--------
>>> import numpy as np
>>> from ddg.math.complex import complex_sample
>>> complex_sample(seed=15)
(-1.4308730228590871-0.9365477163197146j)
>>> complex_sample(1, 15)
array([-1.43087302-0.93654772j])
>>> complex_sample([2, 1], 15)
array([[-1.43087302+0.39393836j],
[-0.93654772-0.52408663j]])
>>> type(complex_sample(1)[0])
<class 'numpy.complex128'>
>>> type(complex_sample())
<class 'complex'>
"""
gen = np.random.default_rng(seed=seed)
re = gen.standard_normal(size=size)
im = gen.standard_normal(size=size)
return re + 1j * im
#############
# conversion
#############
[docs]def to_array(z):
"""
Convert a complex number to a numpy.ndarray.
Parameters
----------
z : complex
Returns
-------
numpy.ndarray of shape (2, input_shape)
An array with a shape of (2, input_shape), where input_shape is the shape
of the input numpy array. If the input is `complex`, the shape is `(2,)`.
Examples
--------
>>> from ddg.math.complex import to_array, complex_sample
>>> to_array(complex_sample(seed=15))
array([-1.43087302, -0.93654772])
>>> to_array(complex_sample(1, 15))
array([[-1.43087302],
[-0.93654772]])
"""
return np.array([z.real, z.imag])
[docs]def to_complex(x):
"""
Convert an array_like to a complex number.
Parameters
----------
x : array_like of length shape (2,)
Returns
-------
complex
Examples
--------
>>> from ddg.math.complex import to_complex
>>> to_complex((1, 2))
(1+2j)
>>> to_complex(np.array([1, 0.5]))
(1+0.5j)
>>> to_complex((1, 2, 11, 2, 3, 4, 512, 1412))
(1+2j)
"""
return complex(x[0], x[1])
[docs]def homogeneous(z):
"""
Get complex number in homogeneous coordinates.
z is interpreted as the affine coordinate of an element in CP^1.
Parameters
----------
z : complex
Returns
-------
numpy.ndarray of shape (2,)
Examples
--------
>>> from ddg.math.complex import homogeneous, to_complex
>>> homogeneous(to_complex((1, 2)))
array([1.+2.j, 1.+0.j])
"""
if np.isinf(z):
return np.array([1.0, 0.0])
else:
return np.array([z, 1.0])
#################
# plane geometry
#################
[docs]def scalar_product(z1, z2):
"""
Calculate scalar product of two complex numbers of corresponding elements in R^2.
Also accepts numpy.ndarray as input.
Parameters
----------
z1,z2 : complex or numpy.ndarray of dtype numpy.complex128
Returns
-------
float or numpy.ndarray
Returns a float if inputs are of type complex.
Otherwise a numpy.ndarray is returned.
Examples
--------
>>> from ddg.math.complex import scalar_product, complex_sample
>>> a = complex_sample((2, 2), 15)
>>> a
array([[-1.43087302+0.5256162j , -0.93654772+0.80732362j],
[ 0.39393836-1.44353139j, -0.52408663+1.01706379j]])
>>> b = a[0]
>>> z = a[0][0]
>>> scalar_product(z, z)
2.323670001997294
>>> scalar_product(b, b)
array([2.32367 , 1.52889306])
>>> scalar_product(b, z)
array([2.32367 , 1.76442324])
>>> scalar_product(a, a)
array([[2.32367 , 1.52889306],
[2.2389703 , 1.30908554]])
>>> scalar_product(a, b)
array([[ 2.32367 , 1.52889306],
[-1.32241926, 1.31193176]])
"""
return (z1.conjugate() * z2).real
[docs]def determinant(z1, z2):
"""
For two complex numbers get determinant of matrix of corresponding elements in R^2.
Also accepts numpy.ndarray as input.
Parameters
----------
z1,z2 : complex or numpy.ndarray of dtype numpy.complex128
Returns
-------
float or numpy.ndarray
Returns a float if inputs are of type complex.
Otherwise a numpy.ndarray is returned.
Examples
--------
>>> from ddg.math.complex import determinant, complex_sample
>>> a = complex_sample((2, 2), 15)
>>> a
array([[-1.43087302+0.5256162j , -0.93654772+0.80732362j],
[ 0.39393836-1.44353139j, -0.52408663+1.01706379j]])
>>> b = a[0]
>>> z = b[0]
>>> determinant(z, z)
0.0
>>> determinant(a, a)
array([[0., 0.],
[0., 0.]])
>>> determinant(z, a)
array([[ 0. , -0.66291294],
[ 1.85844974, -1.17982071]])
"""
return (z1.conjugate() * z2).imag
[docs]def rel_angle(z1, z2):
"""Get angle in [0,2*pi) between complex numbers z1,z2.
Parameters
----------
z1,z2 : complex
Returns
-------
float
Examples
--------
>>> from ddg.math.complex import rel_angle, complex_sample
>>> z1 = complex_sample(seed=15)
>>> z2 = complex_sample(seed=16)
>>> rel_angle(z1, z2)
4.888818869307385
"""
alpha = np.angle(z2 / z1)
return alpha if alpha >= 0 else 2 * np.pi + alpha
###############
# cross-ratios
###############
[docs]def cr(z1, z2, z3, z4):
"""Compute complex cross ratio.
cr(z1,z2,z3,z4) = (z1 - z2)/(z2 - z3) * (z3 - z4)/(z4 - z1)
Parameters
----------
z1,z2,z3,z4 : complex
Returns
-------
complex
Examples
--------
>>> from ddg.math.complex import cr, complex_sample
>>> Z = complex_sample(4, 15)
>>> Z
array([-1.43087302+0.5256162j , -0.93654772+0.80732362j,
0.39393836-1.44353139j, -0.52408663+1.01706379j])
>>> cr(Z[0], Z[1], Z[2], Z[3])
(0.5474067528516228-0.08577542856779714j)
"""
z1, z2, z3, z4 = map(homogeneous, [z1, z2, z3, z4])
a = np.linalg.det(np.array([z1, z2]))
b = np.linalg.det(np.array([z2, z3]))
c = np.linalg.det(np.array([z3, z4]))
d = np.linalg.det(np.array([z4, z1]))
if b == 0.0 or d == 0.0:
return complex("infinity")
else:
return (a * c) / (b * d)
[docs]def fourth_point_from_cross_ratio(z, z1, z2, q):
"""
Get z12 such that cr(z, z1, z12, z2) = q
For a quadrilateral of the form ::
z2_____z12
| |
| |
z_____z1
with given z,z1,z2 and q, compute z12 .
Parameters
----------
z,z1,z2,q : complex
Returns
-------
complex
Examples
--------
>>> from ddg.math.complex import fourth_point_from_cross_ratio, complex_sample
>>> Z = complex_sample(4, 15)
>>> Z
array([-1.43087302+0.5256162j , -0.93654772+0.80732362j,
0.39393836-1.44353139j, -0.52408663+1.01706379j])
>>> q = complex(0.5474067528516228 - 0.08577542856779719j)
>>> fourth_point_from_cross_ratio(Z[0], Z[1], Z[2], q)
(-1.580286264659597+0.09784352065927185j)
"""
a = (z - z1) / (z2 - z)
return (a * z2 + q * z1) / (a + q)
################
# circumcircles
################
[docs]def circumcenter(z1, z2, z3):
"""
Get center of circle through three complex numbers.
Parameters
----------
z1,z2,z3 : complex
Returns
-------
numpy.complex128
Examples
--------
>>> from ddg.math.complex import circumcenter, complex_sample
>>> Z = complex_sample(3, 15)
>>> Z
array([-1.43087302-0.52408663j, -0.93654772+0.5256162j ,
0.39393836+0.80732362j])
>>> circumcenter(Z[0], Z[1], Z[2])
(-0.013691724033015891-0.5502195548404201j)
"""
x, y, z = map(lambda x: np.array([x.real, x.imag]), [z1, z2, z3])
D = 2 * (x[0] * (y[1] - z[1]) + y[0] * (z[1] - x[1]) + z[0] * (x[1] - y[1]))
s_x = x[0] ** 2 + x[1] ** 2
s_y = y[0] ** 2 + y[1] ** 2
s_z = z[0] ** 2 + z[1] ** 2
u = (s_x * (y[1] - z[1]) + s_y * (z[1] - x[1]) + s_z * (x[1] - y[1])) / D
v = (s_x * (z[0] - y[0]) + s_y * (x[0] - z[0]) + s_z * (y[0] - x[0])) / D
return u + 1j * v
[docs]def circumradius(z1, z2, z3):
"""
Get radius of circle through three complex numbers.
Parameters
----------
z1,z2,z3 : complex
Returns
-------
float
Examples
--------
>>> from ddg.math.complex import circumradius, complex_sample
>>> Z = complex_sample(3, 15)
>>> Z
array([-1.43087302-0.52408663j, -0.93654772+0.5256162j ,
0.39393836+0.80732362j])
>>> circumradius(Z[0], Z[1], Z[2])
1.4174222248902117
"""
x, y, z = map(lambda x: np.array([x.real, x.imag]), [z1, z2, z3])
l_xy, l_yz, l_zx = map(np.linalg.norm, (x - y, y - z, z - x))
return (l_xy * l_yz * l_zx) / (4 * euclidean2d.triangle_area(l_xy, l_yz, l_zx))
########
# quads
########
[docs]def intersect_diags(z1, z2, z3, z4):
r"""
For a complex quadrilateral (z1,z2,z3,z4)
get the intersection point of the diagonals.
The four complex points are given as an argument in positive cyclic order. ::
z4-----z3
|\ /|
| x |
|/ \|
z1-----z2
Parameters
----------
z1,z2,z3,z4 : complex
Returns
-------
complex
Examples
--------
>>> from ddg.math.complex import intersect_diags, complex_sample
>>> Z = complex_sample(4, 15)
>>> Z
array([-1.43087302+0.5256162j , -0.93654772+0.80732362j,
0.39393836-1.44353139j, -0.52408663+1.01706379j])
>>> intersect_diags(Z[0], Z[1], Z[2], Z[3])
(-1.449982659344593+0.5462373465733343j)
"""
n1 = 1j * (z3 - z1)
n2 = 1j * (z4 - z2)
A = np.array([[n1.real, n1.imag], [n2.real, n2.imag]])
b = np.array([scalar_product(n1, z1), scalar_product(n2, z2)])
x = np.linalg.solve(A, b)
return complex(*x)
[docs]def intersect_edges(z1, z2, z3, z4):
r"""
For a complex quadrilateral (z1,z2,z3,z4)
get the intersection of the two opposite edges (z1,z2) and (z3,z4).
Will return complex infinity if intersection does not exist. ::
z4-----z3--
| | \
| | x
| | /
z1-----z2--
Parameters
----------
z1,z2,z3,z4 : complex
Returns
-------
complex
Examples
--------
>>> from ddg.math.complex import intersect_edges, complex_sample
>>> Z = complex_sample(4, 15)
>>> Z
array([-1.43087302+0.5256162j , -0.93654772+0.80732362j,
0.39393836-1.44353139j, -0.52408663+1.01706379j])
>>> intersect_edges(Z[0], Z[1], Z[2], Z[3])
(-0.5318751583626388+1.037939495771862j)
>>> z = Z[0]
>>> intersect_edges(z, z, z, z)
(inf+0j)
"""
n1 = 1j * (z2 - z1)
n2 = 1j * (z4 - z3)
A = np.array([[n1.real, n1.imag], [n2.real, n2.imag]])
b = np.array([scalar_product(n1, z1), scalar_product(n2, z3)])
if np.linalg.cond(A) < 1 / np.finfo(float).eps:
x = np.linalg.solve(A, b)
else:
return complex("infinity")
return complex(*x)