#!/usr/bin/env python3
"""
This file is part of the PAFFrontendSim.
The PAFFrontendSim is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License v3 as published by the Free Software Foundation.
The PAFFrontendSim is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with the PAFFrontendSim. If not, see https://www.gnu.org/licenses/.
"""
#############################################################################
# import libs
#############################################################################
import numpy as np
from numpy import random
from copy import copy, deepcopy
from matplotlib import pyplot as plt
try:
from scipy import optimize
_SciPy = True
except:
_SciPy = False
import LogClass
#############################################################################
# constants
#############################################################################
_CLASSNAME = "DistributionData"
# Keywords for the distribution type
_TypeNames = {
"SQUARED" : 'Squared distribution' ,
"HEXAGONAL" : 'Hexagonal distribution' ,
"RANDOM" : 'Random distribution with constrained minumum' ,
"VOGEL" : "Pattern of florets produced by Vogel's model"
}
_DEFAULT_NO_POINTS = 100 # default number of points
_MAX_ITER = 10 # maximum number of iterations to get the right number of points
#############################################################################
# classes
#############################################################################
[docs]
class DistributionData :
"""
This module is used to create a distributions of points on a plane with boundary conditions.
It can be used to generate a list of antenna positions on the focal plane of a PAF.
Input parameters are the size of the plane (XY), the type of distribution (square, hexagonal, random, etc), spacing between antennas and more.
"""
#####################################
# internal data and init method
# init method
def __init__ ( self, log = None ):
"""
Initialize class.
"""
if log == None:
self.log = LogClass.LogClass(_CLASSNAME)
else:
self.log = deepcopy(log) # logging class used for all in and output messages
self.log.addLabel(self)
self.Type = '' # Distribution type (determination algorithm)
self.kwargs = '' # arguments for the distribution type
self.minimumSpacing = 0. # minimum spacing between points
self.xWidth = 10. # x-dimension of the distribution
self.yWidth = 10. # y-dimension of the distribution
self.Circular = False # Circular outlined distribution (xWidth taken as radius
self.NoPoints = _DEFAULT_NO_POINTS # number of points to arrange
self.spacing = -1 # element spacing of the distribution (if negative, spacing is adapted to get the given number of elements)
self.Data = [] # data set (x,y coordinates of the distributed points)
self.TypeCall = {
'SQUARED' : self._squared,
'RANDOM' : self._random,
'HEXAGONAL' : self._hexagonal,
'VOGEL' : self._vogels_pattern
}
#####################################
# internal methods
[docs]
def printOut ( self, line, message_type, noSpace = True, lineFeed = True ) :
"""
Wrapper for the printOut method of the actual log class.
"""
self.log.printOut(line, message_type, noSpace, lineFeed)
#####################################
# methods to set boundary conditions
[docs]
def setMinimumSpacing ( self, minSpace) :
self.printOut("Setting minimum spacing to %.2fmm" % minSpace, self.log.DATA_TYPE)
self.minimumSpacing = minSpace
[docs]
def setGeometry ( self, xw, yw = -1, circ = False) :
if yw < 0 : yw = xw
self.xWidth = xw
self.yWidth = yw
self.Circular = circ
if circ == True :
self.printOut("Using a circular geometry with diameter of %.2fmm" % xw, self.log.DATA_TYPE)
else :
self.printOut("Using a rectangular geometry of %.2f mm x %.2f mm" % (xw, yw), self.log.DATA_TYPE)
[docs]
def setPoints ( self, np = -1) :
if np <= 0 : np = _DEFAULT_NO_POINTS
self.NoPoints = int(np)
self.printOut("Distributing %d points" % np, self.log.DEBUG_TYPE)
[docs]
def setSpacing ( self, spacing = -1) :
self.spacing = spacing
if spacing > 0 :
self.printOut("Distributing points with a spacing of %.2fmm" % spacing, self.log.DATA_TYPE)
[docs]
def setType ( self, type = 'SQUARED', **kwargs ) :
self.printOut("Setting distribution type to '%s'" % type, self.log.DATA_TYPE)
if type.upper() in _TypeNames :
self.Type = type.upper()
self.kwargs = kwargs
else :
self.printOut("type '%s' does not exist" % type, self.log.ERROR_TYPE, False)
#####################################
# Plotting
[docs]
def plotDist ( self, figname = "", show = False ) :
"""
Plot the element usage of the beam-forming over the FoV.
Parameters:
figname (string): file-name of the figure (if given figure is only shown on screen if show == True
show (boolean): flag to show plot on screen even when being written to disk (True)
"""
self.printOut("Plotting the distribution", self.log.FLOW_TYPE)
fig, ax = plt.subplots()
ax.plot(self.Data[:,0], self.Data[:,1], "k+")
# establish labling
ax.set_xlabel("X-position [mm]")
ax.set_ylabel("Y-Position [mm]")
ax.set_title ("Distribution of points" )
# safe of show figure?
if figname != "" :
fig.savefig(figname)
if show :
plt.show()
else :
plt.show()
plt.close()
#####################################
# analytic functions
[docs]
def getDisorder ( self ) :
"""
Calculates the number variance var^2_N(R) of the point-counts in various size circles.
var ~ R^(2-alpha)
Algorithm calculates number variance in a sample of 500 circles of one size. Radius of the circle is then changed. The curve variance versus circle radius is fitted to f(x) = a*x^b and 2-b is returned
This only works reliable for very large sets of distributions!
Returns:
alpha (float): ?
"""
self.printOut("calculating the hyperuniformity of the distribution", self.log.FLOW_TYPE)
if not _SciPy :
self.printOut("SciPy-Lib is not available, can't create the distribution!", self.log.WARNING_TYPE, False)
return
if len(self.Data) < 2 : return 0
# get maximum possible size of the distribution
mx = np.amax(self.Data[:, 0])
my = np.amax(self.Data[:, 1])
r = mx
if my < mx : r = my
self.printOut("maximum radius is %.2fmm and the distribution has %d points" % (r, len(self.Data)), self.log.DATA_TYPE, False)
# loop radii
var = []
rlist = []
number = 20
for ri in range (1, number):
R = 0.5*r/number*ri
self.printOut("using radius %.2f mm" % R, self.log.DEBUG_TYPE, False, False)
rlist.append(R)
n = np.zeros(500, dtype = int)
# loop 1000 random circles
for i in range (0, 500):
rx = (random.rand() - 0.5) * (r - R)
ry = (random.rand() - 0.5) * (r - R)
# calculate number of points within the circle
for dp in (self.Data - np.asarray([rx, ry])) :
rd = np.dot(dp, dp)
if rd < R*R :
n[i] += 1
# calculate variance of the number of points within the circles
n = n-np.mean(n)
var.append(np.dot(n,n)/len(n))
self.printOut("", self.log.DEBUG_TYPE, False)
self.printOut("fitting curve", self.log.DEBUG_TYPE, False)
# fit the exponential function
def f(x, a, b):
return a*x**b
def residual(p, x, y):
return y - f(x, *p)
p0 = [1., 1.]
popt, pcov = optimize.leastsq(residual, p0, args=(rlist, var))
plt.plot(rlist, var, 'k+')
plt.plot(rlist, f(rlist, *popt), 'r-')
plt.show()
plt.close()
self.printOut("got a hyperuniformity of : %.2f " % (2.-popt[1]), self.log.DATA_TYPE, False)
return 2.-popt[1]
[docs]
def getSpacing ( self ) :
"""
Returns the area spacing sqrt(area / noElements).
"""
if self.Circular :
area = np.pi*self.xWidth**2./4.
else :
area = self.xWidth*self.yWidth
points = self.NoPoints
if points < 1 :
points = 1
return np.sqrt(area/points)
#####################################
# internal distribution creation methods
def _squared ( self, step, **kwargs ) :
"""
Create a squared distribution.
Parameters:
step (int): step size between points.
Returns:
tuple (tuple): numpy 2d array of element positions, step-size (int)
"""
dist = []
for x in np.arange(0., self.xWidth/2., step):
for y in np.arange(0., self.yWidth/2., step):
if (not self.Circular) or (np.sqrt(x*x+y*y) < self.xWidth/2.) :
dist.append([x, y])
if x != 0 :
dist.append([-x, y])
if y != 0 :
dist.append([ x, -y])
if x != 0 and y!= 0 :
dist.append([-x, -y])
return np.asarray(dist), step
def _hexagonal ( self, step, **kwargs ) :
"""
Create a list of element positions in a hexagonal arrangement with spacing 1.
Parameters:
rings (int): number of hexagonal rings of the pattern
rot (?): rotation angle between the rings (pure hexagonal: 0). Incorrect!
Returns:
tuple (tuple): numpy 2d array of element positions, step-size (int)
"""
rotangle = kwargs.get("rot", 0.)
self.printOut("using a ring-rotation of %.2f deg\n" % rotangle, self.log.DEBUG_TYPE, False)
dist = [] # empty position list
rawPosList = []
allRings = int(self.xWidth / step) + 1
if self.yWidth > self.xWidth :
allRings = int(self.yWidth / step) + 1
allRings = allRings * 2 # account for corners
rot = -rotangle # ring rotation is zero at start
for rings in range(1, allRings+1) :
for r in range(-rings+1, rings) : # create hex-pattern
y = np.cos(30/180*np.pi)*r*step
for p in range(0, 2*rings-1-np.abs(r)) :
x = (-(2*rings-1-np.abs(r) - 1) /2. + p)*step
if ( ( (not self.Circular) or (np.sqrt(x*x+y*y) < self.xWidth/2.) ) and
(np.absolute(x) <= self.xWidth/2.) and
(np.absolute(y) <= self.yWidth/2.) ) :
xn = np.cos(np.pi*rot/180.)*x - np.sin(np.pi*rot/180.)*y
yn = np.sin(np.pi*rot/180.)*x + np.cos(np.pi*rot/180.)*y
take = True
for p in rawPosList :
if (np.abs(p[0]-x) < 0.001) and (np.abs(p[1]-y) < 0.001) :
take = False
break
if take :
dist.append([xn, yn])
rawPosList.append([x, y])
rot += rotangle # increase rotation angle with every ring
return np.asarray(dist), step
def _random ( self, step, **kwargs ) :
"""
Create a squared distribution.
Parameters:
step (int): step size between points
Returns:
tuple (tuple): numpy 2d array of element positions, step-size (int)
"""
points = int(self.xWidth*self.yWidth/step/step)
dist = []
tries = 0
while (len(dist) < points) and (tries < points**2.) :
x = (random.rand()-0.5)*self.xWidth
y = (random.rand()-0.5)*self.yWidth
if (not self.Circular) or (np.sqrt(x*x+y*y) < self.xWidth/2.) :
# check distance
ok = True
if len(dist) > 0 :
for d in dist:
if (d[0]-x)**2. + (d[1]-y)**2. < self.minimumSpacing**2. :
ok = False
break
if ok :
dist.append([x, y])
tries +=1
self.NoPoints = len(dist)
return np.asarray(dist), self.getSpacing()
def _vogels_pattern ( self, startStep, **kwargs ) :
"""
Pattern of florets produced by Vogel's model (exactly, Vogel's model needs the angle to be set to 135.508deg).
Parameters:
startStep (float?): scaling factor of Vogels model (step-size or distance between points)
angle (float?): angle of the florets. Set at default to Vogel's model
voffset (float?): move center of the florets out of the distribution center [mm] in x-direction
Returns:
tuple (tuple): numpy 2d array of element positions, step-size (int)
"""
angle = kwargs.get("angle", 137.508)
self.printOut("set angle of the florets to %.2fdeg\n" % angle, self.log.DEBUG_TYPE, False)
offset = kwargs.get("voffset", 0.)
self.printOut("shift center of the florets by %.2fmm in x-direction\n" % offset, self.log.DEBUG_TYPE, False)
angle = np.pi*angle/180.
step = startStep
calcStep = 0
i = 0
while (np.abs(calcStep - startStep) > 1.) and (i < _MAX_ITER):
dist = []
r = step
a = 0.
n = 1.
while r*r < ((0.5*self.xWidth+offset)*(0.5*self.xWidth+offset) + 0.25*self.yWidth*self.yWidth) :
x = np.cos(a)*r - offset
y = np.sin(a)*r
if ( (not self.Circular) or (np.sqrt(x*x+y*y) < self.xWidth/2.) ) :
dist.append([x, y])
n+=1
a+=angle
r = np.sqrt(n)*step
self.NoPoints = len(dist)
calcStep = self.getSpacing()
self.printOut("Pattern of florets iteration %d ==> step-size = %.2fmm" % (i, calcStep), self.log.DEBUG_TYPE, False, False)
step = step*(startStep/calcStep)
i += 1
self.printOut("", self.log.DEBUG_TYPE, False)
self.NoPoints = len(dist)
return np.asarray(dist), self.getSpacing()
[docs]
def createDistribution( self ) :
"""
Creates the distribution as set.
Returns:
step (int?): Spacing between the elements NOTE?
"""
self.printOut("Creating a %s" % _TypeNames.get(self.Type), self.log.FLOW_TYPE)
if self.spacing <= 0 :
step = np.sqrt(self.xWidth*self.yWidth/self.NoPoints)
else :
step = self.spacing
self.NoPoints = self.xWidth*self.yWidth/step/step
self.printOut("using a fixed spacing of %.2f mm" % step, self.log.DATA_TYPE, False)
if step < self.minimumSpacing :
self.printOut("Distribution is not possible, spacing needs to be smaller then minimum spacing", self.log.ERROR_TYPE, False)
return False
if self.spacing > 0 :
calcStep = 0
iteration = 0
while (np.abs(calcStep - self.spacing) > 1.) and ( iteration < _MAX_ITER ):
dist, calcStep = self.TypeCall.get(self.Type, self._squared)( step, **self.kwargs )
self.printOut("Iteration %d : getting %d points at step-size %.2fmm\n" % (iteration, len(dist), calcStep), self.log.DEBUG_TYPE, False, False)
step = step*(self.spacing/calcStep)
iteration += 1
self.printOut("", self.log.DEBUG_TYPE, False)
step = calcStep
else :
n = 0
iteration = 0
while ( np.absolute(n - self.NoPoints) > 2) and ( iteration < _MAX_ITER ) :
dist = self.TypeCall.get(self.Type, self._squared)( step, **self.kwargs )
n = len(dist)
self.printOut("Iteration %d : getting %d points at step-size %.2fmm\n" % (iteration,n, step), self.log.DEBUG_TYPE, False, False)
step = step * np.sqrt(n/self.NoPoints)
iteration += 1
self.printOut("", self.log.DEBUG_TYPE, False)
self.NoPoints = len(dist)
self.spacing = step
self.Data = np.asarray(dist)
return step
[docs]
def shuffle ( self, radius ) :
"""
Randomizes an ordered distribution.
Parameters:
radius (float?): Radius in side which the element can be moved around.
"""
for d_idx, d in enumerate(self.Data) :
OK = False
rx = 0.
ry = 0.
while not OK :
rx = (random.rand()-0.5)*radius
ry = (random.rand()-0.5)*radius
OK = True
for e_idx, e in enumerate(self.Data) :
if e_idx != d_idx :
dx = e[0] - d[0] + rx
dy = e[1] - d[1] + ry
if dx**2 + dy**2 < self.minimumSpacing**2. :
OK = False
break
self.Data[d_idx, 0] += rx
self.Data[d_idx, 1] += ry