Remove scipy dependency by replacing eigvals() function

[smith120bh]

[smith120bh]

It appears that the nonlinear solver is getting random, uncaught failures after this switch... Random failures are always concerning, but I'm also concerned that anaStruct is returning a buckling factor of 0.0 rather than throwing an actual error. This'll require some more digging and investigation before this can be safely merged!

[smith120bh]

Okay, looks like this is actually a known but long-running bug in Numpy and/or Conda and/or OpenBlas: https://github.com/numpy/numpy/issues/20233

Unfortunately, this goes beyond my comprehension of the intricacies of Python. I would love to get any help that others could provide here for working around this bug!

I'd like to drop the scipy dependency, but I'm not willing to do so if the alternative numpy implementation of inv() / eigvals() is unreliable.

[someparsa]

Hi Brooks, I suggest that instead of clearing dependencies at this stage, just please comment the lines. This way it will be easier to see the changes on the go. Just a suggestion.

[someparsa]

An easy way to add the eigen function from the SciPy to our package. The link is here (+) and this is a copy. Mapping all the dependencies here may be tough but will get us rid of all failures and broken dependencies.

def eigvals(a, b=None, overwrite_a=False, check_finite=True,
            homogeneous_eigvals=False):

    return eig(a, b=b, left=0, right=0, overwrite_a=overwrite_a,
               check_finite=check_finite,
               homogeneous_eigvals=homogeneous_eigvals)

And this is the eig function as an objective reference (+):

def eig(a, b=None, left=False, right=True, overwrite_a=False,
        overwrite_b=False, check_finite=True, homogeneous_eigvals=False):

    a1 = _asarray_validated(a, check_finite=check_finite)
    if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
        raise ValueError('expected square matrix')
    overwrite_a = overwrite_a or (_datacopied(a1, a))
    if b is not None:
        b1 = _asarray_validated(b, check_finite=check_finite)
        overwrite_b = overwrite_b or _datacopied(b1, b)
        if len(b1.shape) != 2 or b1.shape[0] != b1.shape[1]:
            raise ValueError('expected square matrix')
        if b1.shape != a1.shape:
            raise ValueError('a and b must have the same shape')
        return _geneig(a1, b1, left, right, overwrite_a, overwrite_b,
                       homogeneous_eigvals)

    geev, geev_lwork = get_lapack_funcs(('geev', 'geev_lwork'), (a1,))
    compute_vl, compute_vr = left, right

    lwork = _compute_lwork(geev_lwork, a1.shape[0],
                           compute_vl=compute_vl,
                           compute_vr=compute_vr)

    if geev.typecode in 'cz':
        w, vl, vr, info = geev(a1, lwork=lwork,
                               compute_vl=compute_vl,
                               compute_vr=compute_vr,
                               overwrite_a=overwrite_a)
        w = _make_eigvals(w, None, homogeneous_eigvals)
    else:
        wr, wi, vl, vr, info = geev(a1, lwork=lwork,
                                    compute_vl=compute_vl,
                                    compute_vr=compute_vr,
                                    overwrite_a=overwrite_a)
        t = {'f': 'F', 'd': 'D'}[wr.dtype.char]
        w = wr + _I * wi
        w = _make_eigvals(w, None, homogeneous_eigvals)

    _check_info(info, 'eig algorithm (geev)',
                positive='did not converge (only eigenvalues '
                         'with order >= %d have converged)')

    only_real = numpy.all(w.imag == 0.0)
    if not (geev.typecode in 'cz' or only_real):
        t = w.dtype.char
        if left:
            vl = _make_complex_eigvecs(w, vl, t)
        if right:
            vr = _make_complex_eigvecs(w, vr, t)
    if not (left or right):
        return w
    if left:
        if right:
            return w, vl, vr
        return w, vl
    return w, vr

We will need some polish on the codes but another concern will be the copy-right of the SciPy package. Although I think it is under open-source free-source licenses and we won't have any specific problem. Other solution will be using the GUN's pieces of codes. BTW if we only demand this for eigen values, doing it from scratch won't be much of stress. Random thoughts.

[smith120bh]

Thanks, @someparsa ! I did look at the option of extracting the eigvals() function from scipy - but if you keep tracing back dependent functions, you'll find that the eigenvalue calculation in scipy is not a pure Python function. The hard maths of it is actually run through a series of low-level C and Fortran functions... Copying those into anaStruct is beyond my Python expertise, but AFAIK would mean bringing in conda and slew of other things as as dependencies. I'd suggest that that'd be worse than just keeping scipy as a dependency!

Another option could be to make scipy into an extras_require dependency, since the scipy eigvals() function is only needed for non-linear buckling calculations. We could just fail the non-linear buckling calc if scipy isn't available.

All that said, removing scipy as a dependency would be nice, but there is also only so much work I'm willing to put into it at this point. I'm more interested in building some better validation tests, and in adding modal analysis functionality to anaStruct!