Skip to content

Commit

Permalink
cleaning up
Browse files Browse the repository at this point in the history
  • Loading branch information
mdoucet committed Jan 19, 2024
1 parent 301bae9 commit 535165b
Show file tree
Hide file tree
Showing 3 changed files with 21 additions and 24 deletions.
2 changes: 1 addition & 1 deletion example/sphere_pytorch_prototype.py
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ def make_kernel(model, q_vectors,device):
print("q",q[6])
kernel = model.make_kernel([q])

pars = {'radius': 200, 'radius_pd': 0.1, 'radius_pd_n':10000, 'scale': 2}
pars = {'radius': 200, 'radius_pd': 0.1, 'radius_pd_n':1000, 'sld':2, 'sld_pd': 0.1, 'sld_pd_n':100, 'scale': 2}

t_before = time.time()
Iq = call_kernel(kernel, pars)
Expand Down
37 changes: 17 additions & 20 deletions sasmodels/kerneltorch.py
Original file line number Diff line number Diff line change
Expand Up @@ -139,37 +139,28 @@ def __init__(self, model_info, q_input, device):
# will use views into this vector, relying on the fact that a
# an array of no dimensions acts like a scalar.
parameter_vector = np.empty(len(partable.call_parameters)-2, np.double)

# Create views into the array to hold the arguments.
offset = 0
kernel_args, volume_args = [], []
kernel_idx = []
for p in partable.kernel_parameters:
if p.length == 1:
# Scalar values are length 1 vectors with no dimensions.
v = parameter_vector[offset:offset+1].reshape(())
else:
# Vector values are simple views.
v = parameter_vector[offset:offset+p.length]
offset += p.length
if p in kernel_parameters:
kernel_args.append(v)
kernel_idx.append(offset)
if p in volume_parameters:
volume_args.append(v)
offset += p.length

# Hold on to the parameter vector so we can use it to call kernel later.
# This may also be required to preserve the views into the vector.
self._parameter_vector = parameter_vector

# Generate a closure which calls the kernel with the views into the
# parameter array.
if q_input.is_2d:
form = model_info.Iqxy
qx, qy = q_input.q[:, 0], q_input.q[:, 1]
self._form = lambda: form(qx, qy, *kernel_args)
else:
form = model_info.Iq
q = q_input.q
self._form = lambda: form(q, *kernel_args)
self._kernel_idx = kernel_idx

# Generate a closure which calls the form_volume if it exists.
self._volume_args = volume_args
Expand All @@ -191,9 +182,11 @@ def _call_kernel(self, call_details, values, cutoff, magnetic, radius_effective_
#call_details.show(values)
radius = ((lambda: 0.0) if radius_effective_mode == 0
else (lambda: self._radius(radius_effective_mode)))

_form = self.info.Iqxy if self.q_input.is_2d else self.info.Iq
self.result = _loops(
self._parameter_vector, self._form, self._volume, radius,
self.q_input.nq, call_details, values, cutoff,self.device)
self._parameter_vector, self._kernel_idx, _form, self._volume, radius,
self.q_input, call_details, values, cutoff, self.device)

def release(self):
# type: () -> None
Expand All @@ -204,7 +197,7 @@ def release(self):
self.q_input = None


def _loops(parameters, form, form_volume, form_radius, nq, call_details,
def _loops(parameters, kernel_idx, form, form_volume, form_radius, q_input, call_details,
values, cutoff,device = 'cpu'):
# type: (np.ndarray, Callable[[], np.ndarray], Callable[[], float], Callable[[], float], int, CallDetails, np.ndarray, float, str) -> None
################################################################
Expand All @@ -225,13 +218,14 @@ def _loops(parameters, form, form_volume, form_radius, nq, call_details,
# calling the respective functions.
n_pars = len(parameters)
parameters = torch.DoubleTensor(parameters)
print(parameters)
q_values = torch.DoubleTensor(q_input.q)
#parameters[:] = values[2:n_pars+2]
parameters[:] = torch.DoubleTensor(values[2:n_pars+2]).to(device)

print("parameters",parameters)
if call_details.num_active == 0:
total = form()
kernel_args = [parameters[i] for i in kernel_idx]
total = form(q_values, *kernel_args)
weight_norm = 1.0
weighted_shell, weighted_form = form_volume()
weighted_radius = form_radius()
Expand Down Expand Up @@ -259,7 +253,7 @@ def _loops(parameters, form, form_volume, form_radius, nq, call_details,
pd_length = call_details.pd_length[:call_details.num_active]

#total = np.zeros(nq, np.float64)
total = torch.zeros(nq, dtype= torch.double).to(device)
total = torch.zeros(q_input.nq, dtype= torch.double).to(device)

#print("ll", range(call_details.num_eval))
#parallel for loop
Expand Down Expand Up @@ -287,7 +281,10 @@ def _loops(parameters, form, form_volume, form_radius, nq, call_details,
# exclude all q for that NaN. Even better would be to have an
# INVALID expression like the C models, but that is expensive.
#Iq = np.asarray(form(), 'd')
Iq = torch.asarray(form()).to(device)
#print(parameters)
#print("kernel_idx:", kernel_idx)
kernel_args = [parameters[i] for i in kernel_idx]
Iq = torch.asarray(form(q_values, *kernel_args)).to(device)
if torch.isnan(Iq).any():
continue

Expand Down
6 changes: 3 additions & 3 deletions sasmodels/models/_spherepy.py
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@
category = "shape:sphere"

# ["name", "units", default, [lower, upper], "type","description"],
parameters = [["sld", "1e-6/Ang^2", 1, [-inf, inf], "",
parameters = [["sld", "1e-6/Ang^2", 1, [-inf, inf], "volume",
"Layer scattering length density"],
["sld_solvent", "1e-6/Ang^2", 6, [-inf, inf], "",
"Solvent scattering length density"],
Expand All @@ -71,7 +71,7 @@
]


def form_volume(radius):
def form_volume(sld, radius):
"""Calculate volume for sphere"""
return 1.333333333333333 * pi * radius ** 3

Expand All @@ -93,7 +93,7 @@ def Iq(q, sld, sld_solvent, radius):
## set numpy to ignore the 0/0 error before we do though...
bes = 3 * (sn - qr * cn) / qr ** 3 # may be 0/0 but we fix that next line
bes[qr == 0] = 1
fq = bes * (sld - sld_solvent) * form_volume(radius)
fq = bes * (sld - sld_solvent) * form_volume(sld, radius)
return 1.0e-4 * fq ** 2
Iq.vectorized = True # Iq accepts an array of q values

Expand Down

0 comments on commit 535165b

Please sign in to comment.