I have added a new py_Fact function to the py_Scipy module, that returns the factorial of an integer as a float if it is up to 15 digits long, or a string for longer numbers.
The revised code, and the example shown below in the py_Special-Dist.xlsb spreadsheet, can be downloaded from:
@xl_func(category = "Scipy", help_topic="https://docs.scipy.org/doc/scipy/reference/special.factorial.html!0")
@xl_arg('n', 'numpy_array<int>', ndim=1)
@xl_arg('exact', 'bool')
@xl_return('numpy_array<var>')
def py_fact(n, exact = False):
res = scipy.special.factorial(n, exact = exact)
if exact:
for i in range(0, res.shape[0]):
if res[i] > 1E15: res[i] = str(res[i])
return res
If the “exact” argument is set to True, or omitted, the Python function returns a long integer of the required length, but if this is returned to Excel as a number it would be converted to a float, so all the results greater than 1E15 are converted to a string so no digits are lost, and the resulting array is returned to Excel as a variant array.
In the example below the last 13 digits of the strings have been compared with the last 13 digits of the values returned as floats. Because the factorial numbers always have trailing zeros the integer and float results are exactly equal up to 22!, after which the difference rapidly starts to increase.
Recently someone told me that they were going to a birthday party of a neighbour who was turning 106! and I can now point out that this is equal to: 114628056373470835453434738414834942870388487424139673389282723476762012382449946252660360871841673476016298287096435143747350528228224302506311680000000000000000000000000.
I have recently needed code to rotate 2D coordinates, either about the origin, or some specified other point, using pyxll to call the code from Excel. My first attempt was a direct translation of some VBA code, that also worked with 3D coordinates, with rotation about any of the 3 axes:
@xl_func
@xl_arg('Rcoords', 'numpy_array')
@xl_arg('Axis', 'int')
def py_Rotate(Rcoords, Rotation, Axis, result = 0):
return Rotate(Rcoords, Rotation, Axis, result)
def Rotate(Rcoords, Rotation, Axis, result = 0):
# # Rotate about the origin, Rotation in radians
Plane = np.zeros(2, dtype = int)
if Rotation != 0:
if Axis == 1:
Plane[0] = 1
Plane[1] = 2
elif Axis == 2:
Plane[0] = 2
Plane[1] = 0
elif Axis == 3:
Plane[0] = 0
Plane[1] = 1
else:
return "Invalid Axis"
NumRows = Rcoords.shape[0]
# Rotate
for i in range(0, NumRows):
try:
if (Rcoords[i, Plane[0]] == 0 and Rcoords[i, Plane[1]] == 0) == False:
r = (Rcoords[i, Plane[0]]** 2 + Rcoords[i, Plane[1]]** 2)** 0.5
Theta = atan2(Rcoords[i, Plane[1]], Rcoords[i, Plane[0]])
Theta = Theta + Rotation
Rcoords[i, Plane[0]] = r * cos(Theta)
Rcoords[i, Plane[1]] = r * sin(Theta)
except:
pass
if result == 0:
return Rcoords
else:
return Rcoords[0, result-1]
This was greatly simplified, and speeded up, by using 2D coordinates for the input, simplifying the rotation formula, and using Numpy procedures to operate on entire columns, rather than looping through row by row:
To rotate about a point other than the origin the 3 functions above were modified by:
Translating the coordinates to move the rotation point to the origin.
Rotate the coordinates about the origin.
Translate the rotated coordinates back by the same amount.
@xl_func()
@xl_arg('Rcoords', 'numpy_array')
@xl_arg('RotnPt', 'numpy_array')
@xl_arg('Axis', 'int')
def RotateC(Rcoords, Rotation, RotnPt, Axis, result = 0):
# Rotate about RotnPt, Rotation in radians
Tol = 0.000000000001
NumRows = Rcoords.shape[0]
if abs(Rcoords[0, 0] - Rcoords[NumRows-1, 0]) + abs(Rcoords[0, 1] - Rcoords[NumRows-1, 1]) < Tol: NumRows = NumRows - 1
NumRCols = Rcoords.shape[1]
if RotnPt.shape[0] > RotnPt.shape[1]:
Rotnpt2 = RotnPt
NumCCols = RotnPt.shape[0]
RotnPt = np.zeros((1, NumCCols))
RotnPt[0, :] = Rotnpt2[:, 0]
NumCCols = RotnPt.shape[1]
if NumRCols < NumCCols:
NumCols = NumRCols
else:
NumCols = NumCCols
OCoords = np.zeros((NumRows, NumRCols))
for i in range(0, NumRows):
for j in range(0, NumCols):
OCoords[i, j] = Rcoords[i, j] - RotnPt[0, j]
OCoords = Rotate(OCoords, Rotation, Axis)
for i in range(0, NumRows):
for j in range(0, NumCols):
OCoords[i, j] = OCoords[i, j] + RotnPt[0, j]
OCoords[0:NumRows, NumCols:NumRCols] = Rcoords[0:NumRows, NumCols:NumRCols]
return OCoords
@xl_func()
@xl_arg('Rcoords', 'numpy_array')
@xl_arg('RotnPt', 'numpy_array')
@xl_arg('Axis', 'int')
def RotateC2(Rcoords, Rotation, RotnPt, Axis, result = 0):
# Rotate about RotnPt, Rotation in radians; with numpy array calculations
Tol = 0.000000000001
NumRows = Rcoords.shape[0]
if abs(Rcoords[0, 0] - Rcoords[NumRows-1, 0]) + abs(Rcoords[0, 1] - Rcoords[NumRows-1, 1]) < Tol: NumRows = NumRows - 1
NumRCols = Rcoords.shape[1]
if RotnPt.shape[0] > RotnPt.shape[1]:
Rotnpt2 = RotnPt
NumCCols = RotnPt.shape[0]
RotnPt = np.zeros((1, NumCCols))
RotnPt[0, :] = Rotnpt2[:, 0]
NumCCols = RotnPt.shape[1]
if NumRCols < NumCCols:
NumCols = NumRCols
else:
NumCols = NumCCols
OCoords = Rcoords - RotnPt
OCoords = Rotate2(OCoords, Rotation)
OCoords = OCoords + RotnPt
OCoords[0:NumRows, NumCols:NumRCols] = Rcoords[0:NumRows, NumCols:NumRCols]
return OCoords
@xl_func()
@xl_arg('Rcoords', 'numpy_array')
@xl_arg('RotnPt', 'numpy_array')
def RotateCM(Rcoords, Rotation, RotnPt):
# Rotate about RotnPt, Rotation in radians; with numpy array calculations
OCoords = Rcoords - RotnPt
OCoordsr = RotateM(OCoords, Rotation)
OCoordsr = OCoordsr + RotnPt
return OCoordsr
Finally for comparison and to check results some code was adapted from https://gist.github.com/LyleScott/ . Note that this code treats clockwise rotations as positive, whereas the other functions use the standard approach of anti-clockwise rotation being positive.
@xl_func()
@xl_arg('xy', 'numpy_array')
@xl_arg('radians', 'float')
@xl_arg('origin', 'numpy_array', ndim=1)
def rotate_around_point_highperf(xy, radians, origin=(0, 0)):
"""Rotate a point around a given point.
From: https://gist.github.com/LyleScott/
I call this the "high performance" version since we're caching some
values that are needed >1 time. It's less readable than the previous
function but it's faster.
"""
x = xy[:,0]
y = xy[:,1]
offset_x, offset_y = origin[:]
adjusted_x = (x - offset_x)
adjusted_y = (y - offset_y)
cos_rad = cos(radians)
sin_rad = sin(radians)
qx = offset_x + cos_rad * adjusted_x + sin_rad * adjusted_y
qy = offset_y + -sin_rad * adjusted_x + cos_rad * adjusted_y
res = np.zeros((qx.shape[0],2))
res[:,0] = qx
res[:,1] = qy
return res
Results and execution times for rotation of a simple shape are shown in the screen-shot below (click on the image for full-size view):
The functions using matrix multiplication are significantly faster, and this improvement is increased with a longer list of coordinates:
The py_Evalu spreadsheet was last presented at Scipy Functions with Excel and pyxll 5 – Evaluate with Units and included the Plot_Math function, which converts a text string to Latex format, then uses Matplotlib to convert this to a graphic image, and optionally evaluates the function. To extend this functionality I have now added the Plot_Quad function to evaluate and display definite integrals.
The updated files are included in the py_SciPy zip file that can be downloaded from:
Details of the required pyxll package (including download, free trial, and full documentation) can be found at:pyxll
For those installing a new copy of pyxll, a 10% discount on the first year’s fees is available using the coupon code “NEWTONEXCELBACH10”.
Examples of the Plot_Quad function are shown on the Latex sheet.
Required inputs are:
The function to be integrated as a text string.
The integration limits.
Any fixed parameters and their values.
When the function is entered the result of the integration is returned in the function cell, and an image of the definite integral is displayed in Latex format immediately below. The image may be selected and dragged to any desired new location, as shown in the second example above. This example returns the area of a circle quadrant, which in this case is accurate to 15 significant figures.
Help on the function can be displayed by clicking on the “insert function” icon, immediately to the left of the edit line, to open the “function wizard”. This shows brief help on each function argument. Detailed help on the SciPy Quad function (used for the evaluation of the integral) can be displayed by clicking on “Help on this function”, at the bottom left of the function wizard display:
By default, the integration variable is “x”. Any other symbol may be specified with the optional “var” argument. Optional arguments for the quad function may also be specified with the “kwargs” argument (see the Scipy help for details of all options). In the example below the epsrel argument has been entered to allow a greater relative error in the result, with in this case only a small reduction in execution time:
In the next example the function has been used to find the force on a circular concrete section generated by a parabolic-rectangular stress distribution, as defined by the Eurocode 2 concrete design code. In this case two separate functions are required, for the rectangular and parabolic part of the stress distribution:
The resulting total force is compared with the results of the py_Circu concrete design function, showing near exact agreement.
The spreadsheet also has an on-sheet check using Simpson’s method with each block divided into 20 layers, showing good (but less exact) agreement:
The py_RC Elastic spreadsheet (last discussed here) has been updated with a new function py_Fatigue3600, providing design for fatigue to the latest AS 3600 concrete structures code. The new spreadsheet and associated Python code can be downloaded from:
The download file includes full open-source Python code. For details of the pyxll package required to link the Python code to Excel see: https://newtonexcelbach.com/python-and-pyxll/
Input for the new py_Fatigue3600 function is as before. Note that the new function is intended for use on bridge structures, and the number of fatigue cycles is calculated to AS 5100.2 requirements, based on the current number of heavy vehicles/lane/day:
The output of both functions is shown below:
The main differences between the calculations to the two codes are:
Different variation in the steel fatigue stress limit related to the number of stress cycles. AS 3600 has a higher stress limit for mid-range cycles, but for low cycle numbers (as in this case), or very high cycles, the AS 5100.5 stress limits are higher.
For the concrete compressive stress limit AS 5100.5 has a conservative limit, not related to the number of cycles. AS 3600 has a procedure for calculating the number of cycles to failure for any specified compressive stress. This has been adapted in the function to calculate a stress limit for any given number of cycles. In most cases the AS 3600 limit is significantly higher than the AS 5100.5 simplified approach.
A check has been added to the spreadsheet to verify that the concrete stress limit calculated by the AS 3600 function does correspond with the required number of stress cycles:
VBA code to import text files to Excel was presented at Importing text files with VBA – 2 and following posts. We now look at Python code for a user defined function (UDF) for the same task.
The code requires Python and Numpy, together with pyxll to link the code to Excel. See Python and pyxll for information on the pyxll package, including a discount coupon for new users.
The code uses the Python readlines function to read the file to a list of text lines, then optionally the Numpy fromstring function to extract the numbers from each line.
import os
import numpy as np
@xl_func
@xl_arg('filename', 'str')
@xl_arg('dir', 'str')
@xl_arg('ExtractVals', 'bool')
@xl_return('numpy_array<var>')
def py_ReadText(filename, dir = '', ExtractVals = True ):
"""
Returns the contents of a text file
:param filename: Text file name.
:param dir: Full directory path.
:param ExtractVals: True = extract numerical values from each line and return as a numpy array.
"""
if dir == '':
dir = os.getcwd() # + r'/'
filename = os.path.join(dir, filename)
try:
with open(filename, 'r') as file:
txta = file.readlines()
except:
return np.array(['File not found at ' + dir])
if ExtractVals:
numlist = []
nrows = len(txta)
numlist = []
maxn = 0
for row in txta:
numline = np.fromstring(row, sep = ' ')
numlist.append(numline)
n = numline.shape[0]
if n > maxn: maxn = n
numa = np.zeros((nrows, maxn))
for i in range(0, nrows):
n = numlist[i].shape[0]
numa[i, 0:n] = numlist[i]
return numa
else:
return np.array(txta, dtype = str)
Sample output is shown in the screenshot below:
Column B imports the file as text, with Column F importing the same file with the numbers split into separate columns.
The text files to be imported may be in any folder, if the full path is specified as the second function argument. Alternatively, the active folder will be used. For the example shown the text files were in the same folder as the spreadsheet, but to make this the active folder the file had to be saved to the same location after opening, using Save-As-Replace.
Newton Excel Bach, not (just) an Excel Blog 