-
Notifications
You must be signed in to change notification settings - Fork 0
/
convstore.py
215 lines (191 loc) · 6.89 KB
/
convstore.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
#!/usr/bin/env python
import numpy as np
from scipy import special
from scipy.optimize import root_scalar
from subprocess import Popen
def zerocrossings(numarray):
'''
Return the indices of the zero crossings in the given numpy array.
Parameters
----------
numarray : array
The array to find the zero crossings of.
Returns
-------
zc : array
The indices of the elelements right before a zero crossing.
'''
return np.where(np.diff(np.sign(numarray)))[0]
def count_shared_chars(s1, s2):
'''
Given two strings, count the number of characters that are the same
at the beginning of both strings.
Parameters
----------
s1 : str
The first string.
s2 : str
The second string.
Returns
-------
shared_count : int
How many characters are the same at the beginning of both strings
until the first character that is different.
'''
shared_count = 0
for c1, c2 in zip(s1, s2):
if c1 == c2:
shared_count += 1
else:
break
return shared_count
def findallroots(fun, xmin, xmax, dx, dtype,
xtol=1e-9,
max_iter=10000, num_sigfigs=6,
razorsharp=0.00001,
method='secant',
verbose=False,
reverse=False):
'''
Find most roots of a function in a given interval. To accomplish this
the function is sampled with the given resolution dx within the
interval (xmin, xmax) and the zero crossings are found. The zero
crossings are then used as initial guesses for the root finding
algorithm.
At the end of the search, the roots are disambiguated by rounding to
the specified number of significant figures. This is done to avoid
returning multiple roots that are very close to each other.
When a method is given which requires a bracket, this bracket is
estimated by slightly increasing the position for the zero crossing.
This might not give an adequate bracket for the root finding
algorithm, so it it expanded until the function effectively changes
sign. If the function does not change after max_fixes iterations, the
root is skipped.
The roots are returned in a sorted array.
Parameters
----------
fun : function
The function to find the roots of, must be univariate.
xmin : float
The lower bound of the interval to search for roots.
xmax : float
The upper bound of the interval to search for roots.
dx : float
The step size to use when searching for roots.
xtol : float, optional
The tolerance to use when searching for roots.
max_iter : int, optional
The maximum number of iterations to use when searching for roots.
num_sigfigs : int, optional
The number of significant figures to use when disambiguating similar roots.
razorsharp : float, optional
The amount to expand the bracket for the root finding algorithm.
method : str, optional
The method to use when searching for roots. Must be one of the following:
'bisect', 'brentq', 'brenth', 'ridder', 'toms748', 'secant'.
verbose : bool, optional
If True, print information about the roots as they are found.
reverse : bool, optional
If True, return the roots in reverse order.
Returns
-------
zeros : array
'''
rightzor = 1+razorsharp
leftzor = 1-razorsharp
numSamples = int(np.ceil((xmax-xmin)/dx))
x = np.linspace(xmin, xmax, numSamples, dtype=dtype)
y = fun(x)
zcs = zerocrossings(y)
zcs = zcs[~np.isnan(y[zcs])]
if verbose:
print(x[zcs])
zeros = []
for zc in zcs:
if verbose:
print("Evaluating root at x = {:.6f}".format(x[zc]))
bracket = (leftzor*x[zc], rightzor*x[zc])
if method in 'bisect brentq brenth ridder toms748'.split(' '):
l, r = bracket
fl, fr = fun(l), fun(r)
max_fixes = 100
if np.isnan(fl) or np.isnan(fr):
continue
nfix = 0
while np.sign(fun(l) * fun(r)) != -1:
l = l*leftzor
r = r*rightzor
nfix += 1
if nfix > max_fixes:
break
if np.isnan(fun(l)) or np.isnan(fun(r)):
break
bracket = (l,r)
if verbose:
print("bracket", bracket, nfix)
zerofind = root_scalar(fun,
bracket=bracket,
method=method,
xtol=xtol,
maxiter=max_iter)
else:
zerofind = root_scalar(fun,
x0=bracket[0],
x1=bracket[1],
method='secant',
xtol=xtol,
maxiter=max_iter)
the_root = zerofind.root
the_check = fun(the_root)
if verbose:
print(the_root, the_check)
if (zerofind.converged
and (the_root <= xmax)
and (the_root >= xmin)):
zeros.append(the_root)
zeros = np.array(zeros, dtype=dtype)
# remove zeros that are equal to the given number of digits
zeros = np.unique(round2sigfigs(zeros, num_sigfigs))
zeros = np.sort(zeros)
if reverse:
zeros = np.flip(zeros)
return zeros
def find_layout_rectangle(ar, N):
'''
Given an aspect ratio and an integer number of elements,
find two integers such that their products is not larger
than N and whose ratio is closest to the given aspect ratio.
Parameters
----------
ar : float
The aspect ratio.
N : int
The number of elements.
Returns
-------
best_pair : array
The two integers that best fit the given aspect ratio.
'''
w_est = np.sqrt(N / ar)
h_est = w_est * ar
w_candidates = np.array([np.floor(w_est), np.ceil(w_est)], dtype=int)
h_candidates = np.array([np.floor(h_est), np.ceil(h_est)], dtype=int)
valid_pairs = np.array([(w, h) for w in w_candidates for h in h_candidates if w * h >= N])
ar_diffs = np.abs(valid_pairs[:, 0] / valid_pairs[:, 1] - ar)
best_pair = valid_pairs[np.argmin(ar_diffs)]
return best_pair
def round2sigfigs(x, p):
'''
Round x to p significant figures.
Parameters
----------
x : float
The number to round.
p : int
The number of significant figures to round to.
REF: https://stackoverflow.com/questions/18915378/rounding-to-significant-figures-in-numpy
'''
x = np.asarray(x)
x_positive = np.where(np.isfinite(x) & (x != 0), np.abs(x), 10**(p-1))
mags = 10 ** (p - 1 - np.floor(np.log10(x_positive)))
return np.round(x * mags) / mags