Tanh-Sinh Quadrature via F2PY – Part 3

For the last of the current series on combining Fortran, Python and Excel with F2PY and ExcelPython I have updated the xlSciPy spreadsheet to include two other variants of the Tanh-Sinh function:

  • Quad_TSi integrates a non-periodic function over a semi-infinite interval (a to infinity)
  • Quad_TSo integrates a periodic function over a semi-infinite interval

The new version can be downloaded from:

Download xlScipy.zip

and now includes all the Fortran source code (tsc.f90), as well as the compiled Fortran file, and Python and VBA code.

The Quad_TSi function input is the same as for Quad_TS, except that the limits range is replaced by a single value, a, the start of the integration interval.

F2PY2-1

The Quad_TSo function has an additional optional “omega” argument. For most functions use of the default value of 1 will give satisfactory results, but in some cases (such as Function 3 below) a different value will give greater precision with far fewer evaluations.  I have not yet found any clear statement of how omega should be optimised, other than trial and error.

F2PY2-2

In modifying the original Fortran code to work with F2PY I found the following links very useful:

and for more information on the background to Tanh-Sinh Quadrature see:

Also, don’t forget the Tanh-Sinh Quadrature spreadsheets from Graeme Dennes posted here, which include VBA versions of all the functions presented here, and many more, together with examples of use in real applications:   

Posted in Excel, Fortran, Link to Python, Maths, Newton, Numerical integration, NumPy and SciPy, UDFs, VBA | Tagged , , , , , , , | Leave a comment

Tanh-Sinh Quadrature via F2PY – Part 2

  • Lets look at some details of the coding.

The original code (available here) includes routines “intdeini” to generate an array of points and weights of the quadrature formula, and “intde” to carry out the integration of any supplied function “f”.

Modifications to intdeini (renamed getaw) were straightforward:

  • an “implicit none” statement was added
  • !f2py lines were added for the input and output
  • The eps and tinyln constants were hard coded

 

      subroutine getaw(lenaw, aw)
        implicit none
        integer lenaw
!f2py intent(in) lenaw
        real*8 aw(0:lenaw)
!f2py intent(out)  aw
        real*8 eps
        real*8 efs, hoff
        integer noff, nk, k, j
        real*8 pi2, tinyln, epsln, h0, ehp, ehm, h, t, ep, em, xw, wg
    ! ---- adjustable parameter ----
          efs = 0.1d0
          hoff = 8.45d0
    ! ------------------------------
          eps =   2.22044604925031E-16
          pi2 = 2 * atan(1.0d0)
          tinyln =  708.396418532264
...

The intde routine (renamed intde3) needed a little more work. The code below shows the f2py comment lines required to allow a python function to be called from the Fortran code:

      subroutine intde3(a, b, aw, x, lenaw)
      implicit none
!f2py intent(callback) f
      external f
      real*8 f
!     The following lines define the signature of func for F2PY:
!f2py real*8 y
!f2py y = fun(y)

      real*8 x(3)
!f2py intent(out) x 

      integer lenaw
!f2py intent(in) lenaw
      real*8 aw(0:lenaw)
!f2py intent(in)  aw
      real*8 a, b
!f2py intent(in) a, b

...

The other major change was:

  •  Where the calls to the target function included a calculation, this was replaced with a variable.

So lines in the original code such as:

...
        xa = ba * aw(j)
          fa = f(a + xa)
          fb = f(b - xa)
...

were replaced with:

...
          xa = ba * aw(j)
          c = a + xa
          fa = f(c)
          c = b - xa
          fb = f(c)
...

Finally 11 short functions were added to the Fortran code for testing purposes:

  •  ...       real*8 function f1(x)
    real*8 x
    !f2py intent(in)  x
    f1 = 1 / sqrt(x)
    end function f1
    ... 

    The Fortran code in tsc2.f90 was then compiled in a form that could be accessed from Python using F2PY:
    Open a Command Prompt window (right-click the Windows Icon and select Command Prompt)

  • If necessary, go to the directory containing your Fortran code
  • Enter f2py -c -m tsc2 tsc2.f90

F2PY1-1

F2PY generates a large volume of text during the compilation process, but will terminate either with an error message, or if compilation was successful, with the message shown below:

F2PY1-2

The resulting compiled file (tsc2.pyd) can be called from Python, in the same way as any other Python module:

F2PY1-3

To call from Excel(or from another Python function) it is convenient to create an interface function, to first call the getaw function, then pass the function to be integrated in alternative ways:

def Quad_TS(Func, limits, ftype, lenaw, var):
    res = np.zeros(6)
    stime = time.clock()
    x = np.zeros(3)
    a = limits[0]
    b = limits[1]
  
    try:
        aw[0] != 0.
    except:
        aw = tsc2.getaw(lenaw)
        
    stime2 = time.clock()
    if ftype == 0:
        c = tsc2.intde3(a, b, aw, tsc2.f1)
    elif ftype > 0:
        Func = getattr(tsc2, Func)
        c = tsc2.intde3(a, b, aw, Func)
    else:
        if ftype == -2: Func = 'lambda ' + var + ': ' + Func
        c = tsc2.intde3(a, b, aw, eval(Func))
    etime = time.clock()
    res[0] = c[0]
    res[1] = c[1]
    res[2] = c[2]
    res[3] = etime - stime2
    res[4] = etime - stime    
    return res

The ftype argument currently has four options for the function to be integrated:

  • ftype = 0 uses the hard coded function tsc2.f1
  • ftype > 0 converts the string Func argument into the named function from the tsc2 module.  In future versions this will be extended to allow the module name to be passed, as well as the function name.
  • ftype = -2 converts the passed text string into a Python lambda function, with the variable symbol passed in the var argument.
  • ftype = any other negative number converts the string Func to a Python Function object, using Eval.

Finally the Python interface function (Quad-TS) is called from Excel, using Excel-Python:

Function Quad_TS(Func As String, Limits As Variant, Optional FType As Long = 1, Optional LenAW As Long = 11150, Optional VarSym As String = "x")
Dim Methods As Variant, Result As Variant, STime As Double, Rtn As Variant, IFunc As String, RunAW As Boolean
    On Error GoTo rtnerr
    STime = MicroTimer
    
    GetArray Limits, 1
        
    Set Methods = Py.Module(modname1)
    Set Result = Py.Call(Methods, "Quad_TS", Py.Tuple(Func, Limits, FType, LenAW, VarSym))
    
    Rtn = Py.Var(Result)
    Rtn(4) = MicroTimer - STime
    Quad_TS = Rtn
    Exit Function

rtnerr:
    Quad_TS = Err.Description

End Function

 

 

Posted in Excel, Link to dll, Link to Python, Maths, Newton, Numerical integration, NumPy and SciPy, UDFs, VBA | Tagged , , , , , , , , | Leave a comment

Tanh-Sinh Quadrature via F2PY

Fast and accurate numerical integration routines based on the Tanh-Sinh method are available as Excel VBA based functions in the Tanh_Sinh Quadrature spreadsheet by Graeme Dennes.  Graeme’s VBA code is based on open source Fortran code by Takuya Ooura, which I thought would form the basis of a useful exercise in the use of F2PY to link the compiled Fortran code to Python, and hence to Excel, via Excel-Python.

The procedure turned out not to be as straightforward as I had hoped, but is now complete, and the  Quad_TS function has been added to the xlScipy spreadsheet, including full open source code:

Download xlScipy.zip

The spreadsheet requires Excel plus an installed copy of Python, including Numpy and Scipy.  All other required files are included in the download.  Unzip all the contents of the download zip file to any convenient folder, and it should work.

This post will look at use of the new function, and the next will look at the detailed changes required to the Fortran code to compile successfully with F2PY.

The screenshots below shows use of the function with four different “function type” options (click on any image for a full-size view):

Quad_TS

Option -1 allows the integration of a function entered as text on the spreadsheet, either a Python lambda function, or a named Python function, including the module name.  Option 0 evaluates a hard-coded function (in this case tsc2.f1), and Option 1 (default) evaluates any named function in the tsc2 module.

Quad_TS2

Option -2 evaluates a text string entered in a worksheet cell, without the Python lambda notation.   The default variable symbol is “x”, but any other symbol may be used as an optional argument, as illustrated above.  The formula may include the Python math functions: log, exp, sin, cos, tan, or sqrt.  Any other function must be preceded by the Python module name, as in Python code.

Quad_TS3

The integration process involves the generation of an array of weights, which defaults to a length of 11150 items.  This is sufficient for the accurate integration of most functions, but for the function shown above, which oscillates rapidly about the x axis, a longer table is required.  This may be specified as the optional fourth argument, allowing the accurate evaluation of the integral.

The function output shows the number of function evaluations, the execution time of the compiled code, and the total execution time, including overhead associated with the call of the function from Excel.  For a single evaluation of a single integral the overhead swamps any performance gain, and the compiled code provides little or no benefit over the VBA function.  With repeated calls from a Python routine however the saving would become significant, and even a single evaluation of an oscillating function can also show a huge performance improvement.

The screen shot below shows the evaluation of the F11 function ( exp(x) * sin( 2*x**3)) using the SciPy integration function.  The compiled Quad_TS function has given a 58 times improvement in the evaluation time and a 20 times improvement in the overall execution time for a single evaluation.

Quad_TS4

Posted in Excel, Fortran, Link to Python, Maths, Newton, Numerical integration, NumPy and SciPy, UDFs, VBA | Tagged , , , , , , , , | Leave a comment

The Big Short and When the Levee Breaks

The new film the Big Short features Led Zeppelin’s “When the Levee Breaks” in the closing credits.

The story of how the complete song (with 1 minute 30 seconds of intro) came to be included in the film is also interesting:

The crazy story of how ‘The Big Short’ got Led Zeppelin to approve song rights

 

Posted in Bach | Tagged , , | Leave a comment

Running Fortran from Python with F2PY

The last post looking at using the Python F2Py package said:

It is also necessary to install a Fortran compiler.  After some trial and error I found that the MinGW (Minimalist GNU for Windows) package worked well.

Trying some simple examples on a new computer however yielded a variety of “unable to find” errors, and a search on each returned a range of conflicting suggestions on how they might be fixed.  Eventually the following (relatively straightforward) process worked for me:

  • Uninstall Python and install the latest version of Anaconda Python, appropriate to the operating system.  I actually had 2 versions of 32 bit Python in different folders, on 64 bit Windows, which may or may not have contributed to the problems.
  • Install the version of MinGW provided with the Anaconda download, from the command line using:
    conda install mingw
    from the top level Anaconda folder.  There is an installation package available on the MinGW site, but I couldn’t get F2Py to work with that installation.
  • That’s all

I had previously installed Visual Studio Express 2015.  It is not clear to me whether this is required or not.

After this process I found that:

  • The simple examples from the F2Py tutorial worked without a problem.
  • Any example with a PRINT* statement in the Fortran code still does not compile.  I haven’t found a solution to this one yet, other than avoiding PRINT statements.
  • The callback example from the tutorial worked after adding a line to declare the FUN variable as REAL*8.  More on this one later.

I hope to have some more detailed and useful examples to present in the near future!

 

Posted in Excel, Link to Python, NumPy and SciPy | Tagged , , , , | Leave a comment

xlSciPy – Python SciPy for Excel; Update with new functions

The xlSciPy spreadsheet, previously described here, has been updated with new functions for integration, finding equation roots and maxima and minima, solving systems of non-linear equations, and evaluation of equations entered as text.  The new version (including full open-source VBA and Python code) may be downloaded from:

xlSciPy.zip

As before, the spreadsheet requires Python, SciPy, and Numpy to be installed, but all other files (including the required ExcelPython files) are included in the download zip file.  Unzip the download file to its own folder, and enable macros when Excel is started, and it should work.

All the available functions (about 350) are listed on the Index tab.  The new functions introduced in this versions are illustrated in the screenshots below:

The integration functions allow the numerical integration of a function entered as text, a named Python function, or a Python function hard coded in the “xl_IntFunc1” function:

SciPy3-01

The first function below oscillates rapidly about the x axis, and requires an increased number of iterations to converge to an accurate result.  The second has converged to machine precision in a much shorter time.

Scipy3-02

The xl_Brent function finds a root of any function of one variable (between specified limits), using Brent’s method.

Scipy3-03

The examples below show the evaluation of the integral of a function entered as text, the same function with additional fixed parameters, and Python versions of the two functions.

Scipy3-04

The xl_MinimizeFS function calls the SciPy minimize_scalar function, which will find the minimum of a Python function of one variable, or an equivalent lambda function entered as text.

Scipy3-05

The xl_MinimizeF function calls the Python minimize function to find the unconstrained or constrained minimum of a scalar function of one or more variables, using one of 11 alternative methods.

Scipy3-06

The examples below show the same root finding problem as in previous examples, followed by an example from the tutorial in the SciPy Manual, finding the minimum of the Rosenbrock Function, using different methods.

Scipy3-07

For constrained minimization, as shown below, the constraints may either be specified with Python functions or text Lambda functions. Change the “Constraint Type” to 2 on the spreadsheet to see the lambda function constraint input.

Scipy3-08

The xl_SolveF function solves a system of non-linear equations with an equal number of unknowns and objective values, using one of 10 alternative methods.

Scipy3-09

The example below shows a system of non-linear equations, from the SciPy Manual reference section.

Scipy3-10

The example below analyses a reinforced concrete section to find depth of neutral axis (DNA), compression face strain (epsc), and tension bar diameter for specified axial force, bending moment, and tensile steel stress values.  A VBA routine to solve the same problem was presented in: https://newtonexcelbach.wordpress.com/2015/10/20/solving-non-linear-equations-with-two-or-more-unknowns-5/

Scipy3-11

The example below, from the SciPy Manual tutorial, illustrates the use of Krylov’s Method (Method 9), with and without the use of pre-conditioners.

Scipy3-12

The results shown in the SciPy manual show faster solution time with the use of a preconditioner, but my results for this example were slightly slower with the preconditioner.

Scipy3-13

Posted in Excel, Link to Python, Maths, Newton, Numerical integration, NumPy and SciPy, UDFs, VBA | Tagged , , , , , , , , , , | Leave a comment

Worksheetfunction.Transpose – changed behaviour in Excel 2013 and 2016

It is (fairly) well known that when using Worksheetfunction.Transpose on an array in VBA the array must be no more than 65,536 (2^16) rows long, or it will return an error.

Recently when using Excel 2016 (32 bit) I thought this limitation had been fixed; I was able to transpose arrays with up to the maximum number of rows allowed by Excel (1,048,576 0r 16 * 2^16), without generating any error.  Then I checked the end of the returned arrays, and discovered that they had in fact been truncated at 65,536 rows.  A closer look revealed that the truncation is in fact done in blocks of 2^16 rows, so if you transpose an array of 65,537 rows the operation will delete all but 1 of them, with no indication that anything unusual has happened.

To investigate what is happening I wrote 7 functions to split an array into single columns, then transpose them, the intent being to generate 1D arrays, which are much faster when transferring large data sets to other programs.  The variations were:

  • Splitarray: Split a 3 column array into 3 single columns, with no transpose.
  • Splitarray-a: Split and transpose in a single operation.
  • Splitarray-b: Redim the output arrays first, then split and transpose in a single operation.
  • Splitarray-c: Redim the output arrays first, then split into three variant arrays, then transpose.
  • Splitarray-d: Transpose without splitting, creating a 3 row array.
  • Splitarray-e: Split into single columns, then copy to 1D arrays with a loop.
  • Splitarray-f:  Convert the input range to a variant array, then split and transpose with a loop.

The functions were modified to return the following information:

  • Row 1:  Run time
  • Rows 2-4:  The number of values in the returned arrays, or in the case of Splitarray-d the  number of rows, then the number of columns.

In the case of Splitarray-c and Splitarray-f, which had a two stage operation, two columns were returned, with times and array sizes at the end of each step.

The results for Excel 2016 are shown below for four different array sizes:

65,536 rows; all functions return all the input data

Transpose1-1

1,048,576 rows: the first function (which did not transpose) and the last two (which used loops) return all the data, but all the others truncate the output to 65,536 values.

Transpose1-2

65,537 rows: The second to fifth functions have all truncated the data to just one value.

Transpose1-3

1,048577 rows: Now all the functions return #VALUE! errors (as in previous versions of Excel), because the specified range size exceeds the maximum number of rows.  Note that VBA arrays can have a much greater number of rows (limited by available memory), but if you are reading data from the spreadsheet the range limit on number of rows applies.

Transpose1-4

Comparing with Excel 2010; results for 65536 rows are the same:

Transpose1-5

With 1 additional row up to 1,048, 576 rows, all the functions that use Transpose return #VALUE! errors:

Transpose1-6

Transpose1-7

And with 1 or more rows past the range limit all the functions return errors:

Transpose1-8

It was my recollection that Excel 2013 worked the same as 2010, but checking with 65,537 rows I found that the arrays were truncated to one value, when Transpose was used, as for Excel 2016.

Transpose1-9

So the conclusions are:

  1. In Excel 2013 and 2016 (32 bit) using Worksheetfunction.Transpose in any VBA routine will truncate the array without warning if it has more than 65536 rows, so if you must use this function be sure to check the size of the array first.
  2. But really, there is no good reason for using it at all.  It can be replaced by a short VBA routine using loops, which will return the full array up to 1,048,576 rows (or larger if working entirely in VBA), and is actually faster than the Transpose version (at least for arrays of the size used in these tests).
Posted in Arrays, Excel, UDFs, VBA | Tagged , , , , , , , | Leave a comment