Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Allow pure python Fq for beta approximation. #570

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Conversation

pkienzle
Copy link
Contributor

Allow pure python F(q) functions such as the following:

def Fq(q, sld, sld_solvent, radius):
    """Calculate F(q), F^2(q) for sphere"""
    qr = q * radius
    sn, cn = sin(qr), cos(qr)
    with np.errstate(all='ignore'):
        bes = 3 * (sn - qr * cn) / qr ** 3 # may be 0/0 but we fix that next line
        bes[qr == 0] = 1
    fq = bes * (1e-2 * (sld - sld_solvent) * form_volume(radius))
    return fq, fq**2
Fq.vectorized = vectorized  # Fq accepts an array of q value

Test using:

python -m sasmodels.compare -pars -midq sphere@hardsphere,_spherepy@hardsphere radius_pd=0.1 structure_factor_mode=1

@butlerpd
Copy link
Member

@yunliu01 and @smalex-z to speak with @pkienzle about this

@butlerpd butlerpd requested review from lucas-wilkins and butlerpd and removed request for smalex-z August 15, 2023 13:49
@butlerpd
Copy link
Member

Looks like it needs some unit tests.

@wpotrzebowski
Copy link

Not a priority for 6.0 but will be good to have

Copy link
Contributor

@lucas-wilkins lucas-wilkins left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks OK and it looks like it will do what it is supposed to do, but there does seem to be a lot of mordernisation to be done on these files. If assessed it more consistency with existing code, rather than something approximating best practice.

# Note that this is doing slightly too much work since we are composing
# two vector results which are then interlaced via _pack_fq above.
Fq = model_info.Fq
if callable(Fq) and not getattr(Fq, 'vectorized', False):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I take it you're using callable here because model_info.compiled isn't exactly what you want. If it is, consider being consistent with the other changes.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

info.compiled is for info.Fq and info.Iq jointly. Only one of them will be defined, and it may be defined using C source or using a callable python function.

It occurs to me that we could use some numba transforms here so that the code runs as fast as a C model on the CPU (or faster since loop parallelization can run over multiple cores). Python models get much easier to write when we don't have to vectorize across if statements.

Vectorized 1D kernel returning Fq, Fq**2
"""
fq, fqsq = zip(*(Fq(qi, *args) for qi in q))
return np.array(fq), np.array(fqsq)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not certain, but slicing a numpy array might be considerably faster than zip(*)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Surprisingly, the zip version is faster:

In [15]: %timeit x,y = [np.array(v) for v in zip(*(((i, i) for i in np.arange(1000))))]
242 μs ± 7.15 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [16]: %timeit x,y = (a:=np.array(list((i, i) for i in np.arange(1000))))[:,0],a[:,1]
389 μs ± 4.75 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

But in any case, this cost is small compared to evaluating the function at each q.

@@ -1092,13 +1094,12 @@ class ModelInfo(object):
#: the model cannot be run in opencl (e.g., because the model passes
#: functions by reference), then set this to false.
opencl = None # type: bool
#: True if the model is compiled with C or OpenCL
compiled = None # type: bool
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It might make sense to make this a dataclass and avoid having variables initialised to values outside their type

def Iq(q, sld, sld_solvent, radius):
"""Calculate I(q) for sphere"""
# Variants for testing purposes: toggle True/False
vectorized = True # Whether to call a q vector or loop over q scalars.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If they were for prototyping, consider removing them, if they're for testing, they and their dependents should be settable after loading the file.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

prototyping. It could be moved to the explore directory.

@lucas-wilkins
Copy link
Contributor

I think @butlerpd said he would try writing some python models to check with.

@butlerpd
Copy link
Member

yes @lucas-wilkins that bit is now on me.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants