Skip to content

Commit

Permalink
lint check & change to python3 & pre-commit check
Browse files Browse the repository at this point in the history
  • Loading branch information
@stevenhua0320
stevenhua0320 committed Jul 29, 2024
1 parent 3705663 commit ee77d44
Show file tree
Hide file tree
Showing 4 changed files with 181 additions and 149 deletions.
119 changes: 67 additions & 52 deletions diffpy/srmise/baselines/nanospherical.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,8 @@

logger = logging.getLogger("diffpy.srmise")

class NanoSpherical (BaselineFunction):

class NanoSpherical(BaselineFunction):
"""Methods for evaluation of baseline of spherical nanoparticle of uniform density.
Allowed formats are
Expand All @@ -47,29 +48,31 @@ def __init__(self, Cache=None):
evaluations.
"""
# Define parameterdict
parameterdict = {'scale':0, 'radius':1}
formats = ['internal']
default_formats = {'default_input':'internal', 'default_output':'internal'}
parameterdict = {"scale": 0, "radius": 1}
formats = ["internal"]
default_formats = {"default_input": "internal", "default_output": "internal"}
metadict = {}
BaselineFunction.__init__(self, parameterdict, formats, default_formats, metadict, None, Cache)
BaselineFunction.__init__(
self, parameterdict, formats, default_formats, metadict, None, Cache
)

#### Methods required by BaselineFunction ####

# def estimate_parameters(self, r, y):
# """Estimate parameters for spherical baseline. (Not implemented!)
#
# Parameters
# r - array along r from which to estimate
# y - array along y from which to estimate
#
# Returns Numpy array of parameters in the default internal format.
# Raises NotImplementedError if estimation is not implemented for this
# degree, or SrMiseEstimationError if parameters cannot be estimated for
# any other reason.
# """
# if len(r) != len(y):
# emsg = "Arrays r, y must have equal length."
# raise ValueError(emsg)
# def estimate_parameters(self, r, y):
# """Estimate parameters for spherical baseline. (Not implemented!)
#
# Parameters
# r - array along r from which to estimate
# y - array along y from which to estimate
#
# Returns Numpy array of parameters in the default internal format.
# Raises NotImplementedError if estimation is not implemented for this
# degree, or SrMiseEstimationError if parameters cannot be estimated for
# any other reason.
# """
# if len(r) != len(y):
# emsg = "Arrays r, y must have equal length."
# raise ValueError(emsg)

def _jacobianraw(self, pars, r, free):
"""Return the Jacobian of the spherical baseline.
Expand All @@ -83,22 +86,26 @@ def _jacobianraw(self, pars, r, free):
needed. True for evaluation, False for no evaluation.
"""
if len(pars) != self.npars:
emsg = "Argument pars must have "+str(self.npars)+" elements."
emsg = "Argument pars must have " + str(self.npars) + " elements."
raise ValueError(emsg)
if len(free) != self.npars:
emsg = "Argument free must have "+str(self.npars)+" elements."
emsg = "Argument free must have " + str(self.npars) + " elements."
raise ValueError(emsg)
jacobian = [None for p in range(self.npars)]
if (free == False).sum() == self.npars:
return jacobian

if np.isscalar(r):
if r <= 0. or r >= 2.*pars[1]:
if free[0]: jacobian[0] = 0.
if free[1]: jacobian[1] = 0.
if r <= 0.0 or r >= 2.0 * pars[1]:
if free[0]:
jacobian[0] = 0.0
if free[1]:
jacobian[1] = 0.0
else:
if free[0]: jacobian[0] = self._jacobianrawscale(pars, r)
if free[1]: jacobian[1] = self._jacobianrawradius(pars, r)
if free[0]:
jacobian[0] = self._jacobianrawscale(pars, r)
if free[1]:
jacobian[1] = self._jacobianrawradius(pars, r)
else:
s = self._getdomain(pars, r)
if free[0]:
Expand All @@ -120,12 +127,12 @@ def _jacobianrawscale(self, pars, r):
"""
s = np.abs(pars[0])
R = np.abs(pars[1])
rdivR = r/R
rdivR = r / R
# From abs'(s) in derivative, which is equivalent to sign(s) except at 0 where it
# is undefined. Since s=0 is equivalent to the absence of a nanoparticle, sign will
# be fine.
sign = np.sign(pars[1])
return -sign*r*(1-(3./4.)*rdivR+(1./16.)*rdivR**3)
return -sign * r * (1 - (3.0 / 4.0) * rdivR + (1.0 / 16.0) * rdivR**3)

def _jacobianrawradius(self, pars, r):
"""Return partial Jacobian wrt radius without bounds checking.
Expand All @@ -141,7 +148,7 @@ def _jacobianrawradius(self, pars, r):
# From abs'(R) in derivative, which is equivalent to sign(R) except at 0 where it
# is undefined. Since R=0 is a singularity anyway, sign will be fine.
sign = np.sign(pars[1])
return sign*s*(3*r**2*(r**2-4*R**2))/(16*R**4)
return sign * s * (3 * r**2 * (r**2 - 4 * R**2)) / (16 * R**4)

def _transform_parametersraw(self, pars, in_format, out_format):
"""Convert parameter values from in_format to out_format.
Expand All @@ -162,15 +169,17 @@ def _transform_parametersraw(self, pars, in_format, out_format):
temp[0] = np.abs(temp[0])
temp[1] = np.abs(temp[1])
else:
raise ValueError("Argument 'in_format' must be one of %s." \
% self.parformats)
raise ValueError(
"Argument 'in_format' must be one of %s." % self.parformats
)

# Convert to specified output format from "internal" format.
if out_format == "internal":
pass
else:
raise ValueError("Argument 'out_format' must be one of %s." \
% self.parformats)
raise ValueError(
"Argument 'out_format' must be one of %s." % self.parformats
)
return temp

def _valueraw(self, pars, r):
Expand All @@ -185,11 +194,11 @@ def _valueraw(self, pars, r):
r - sequence or scalar over which pars is evaluated.
"""
if len(pars) != self.npars:
emsg = "Argument pars must have "+str(self.npars)+" elements."
emsg = "Argument pars must have " + str(self.npars) + " elements."
raise ValueError(emsg)
if np.isscalar(r):
if r <= 0. or r >= 2.*pars[1]:
return 0.
if r <= 0.0 or r >= 2.0 * pars[1]:
return 0.0
else:
return self._valueraw2(pars, r)
else:
Expand All @@ -209,38 +218,44 @@ def _valueraw2(self, pars, r):
"""
s = np.abs(pars[0])
R = np.abs(pars[1])
rdivR = r/R
return -s*r*(1-(3./4.)*rdivR+(1./16.)*rdivR**3)
rdivR = r / R
return -s * r * (1 - (3.0 / 4.0) * rdivR + (1.0 / 16.0) * rdivR**3)

def _getdomain(self, pars, r):
"""Return slice object for which r > 0 and r < twice the radius"""
low = r.searchsorted(0., side='right')
high = r.searchsorted(2.*pars[1], side='left')
low = r.searchsorted(0.0, side="right")
high = r.searchsorted(2.0 * pars[1], side="left")
return slice(low, high)

def getmodule(self):
return __name__

#end of class NanoSpherical

# end of class NanoSpherical

# simple test code
if __name__ == '__main__':
if __name__ == "__main__":

f = NanoSpherical()
r = np.arange(-5, 10)
pars = np.array([-1., 7.])
pars = np.array([-1.0, 7.0])
free = np.array([False, True])
print "Testing nanoparticle spherical baseline"
print "Scale: %f, Radius: %f" %(pars[0], pars[1])
print "-----------------------------------------"
print("Testing nanoparticle spherical baseline")
print("Scale: %f, Radius: %f" % (pars[0], pars[1]))
print("-----------------------------------------")
val = f._valueraw(pars, r)
jac = f._jacobianraw(pars, r, free)
outjac = [j if j is not None else [None]*len(r) for j in jac]
print "r".center(10), "value".center(10), "jac(scale)".center(10), "jac(radius)".center(10)
jac = f._jacobianraw(pars, r, free)
outjac = [j if j is not None else [None] * len(r) for j in jac]
print(
"r".center(10),
"value".center(10),
"jac(scale)".center(10),
"jac(radius)".center(10),
)
for tup in zip(r, val, *outjac):
for t in tup:
if t is None:
print ("%s" %None).ljust(10),
print("%s" % None).ljust(10),
else:
print ("% .3g" %t).ljust(10),
print("% .3g" % t).ljust(10),
print
93 changes: 51 additions & 42 deletions diffpy/srmise/baselines/polynomial.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,8 @@

logger = logging.getLogger("diffpy.srmise")

class Polynomial (BaselineFunction):

class Polynomial(BaselineFunction):
"""Methods for evaluation and parameter estimation of a polynomial baseline."""

def __init__(self, degree, Cache=None):
Expand All @@ -41,17 +42,19 @@ def __init__(self, degree, Cache=None):
emsg = "Argument degree must be an integer."
raise ValueError(emsg)
if self.degree < 0:
self.degree = -1 # interpreted as negative infinity
self.degree = -1 # interpreted as negative infinity
# Define parameterdict
# e.g. {"a_0":3, "a_1":2, "a_2":1, "a_3":0} if degree is 3.
parameterdict = {}
for d in range(self.degree+1):
parameterdict["a_"+str(d)] = self.degree - d
formats = ['internal']
default_formats = {'default_input':'internal', 'default_output':'internal'}
for d in range(self.degree + 1):
parameterdict["a_" + str(d)] = self.degree - d
formats = ["internal"]
default_formats = {"default_input": "internal", "default_output": "internal"}
metadict = {}
metadict["degree"] = (degree, repr)
BaselineFunction.__init__(self, parameterdict, formats, default_formats, metadict, None, Cache)
BaselineFunction.__init__(
self, parameterdict, formats, default_formats, metadict, None, Cache
)

#### Methods required by BaselineFunction ####

Expand Down Expand Up @@ -81,7 +84,7 @@ def estimate_parameters(self, r, y):
return np.array([])

if self.degree == 0:
return np.array([0.])
return np.array([0.0])

if self.degree == 1:
# Estimate degree=1 baseline.
Expand All @@ -90,15 +93,16 @@ def estimate_parameters(self, r, y):
# lies above the baseline.
# TODO: Make this more sophisticated.
try:
cut = np.max([len(y)/10, 1])
cut = np.max([len(y) / 10, 1])
cut_idx = y.argsort()[:cut]

import numpy.linalg as la

A = np.array([r[cut_idx]]).T
slope = la.lstsq(A, y[cut_idx])[0][0]
return np.array([slope, 0.])
except Exception, e:
emsg = "Error during estimation -- "+str(e)
return np.array([slope, 0.0])
except Exception as e:
emsg = "Error during estimation -- " + str(e)
raise
raise SrMiseEstimationError(emsg)

Expand All @@ -116,10 +120,10 @@ def _jacobianraw(self, pars, r, free):
needed. True for evaluation, False for no evaluation.
"""
if len(pars) != self.npars:
emsg = "Argument pars must have "+str(self.npars)+" elements."
emsg = "Argument pars must have " + str(self.npars) + " elements."
raise ValueError(emsg)
if len(free) != self.npars:
emsg = "Argument free must have "+str(self.npars)+" elements."
emsg = "Argument free must have " + str(self.npars) + " elements."
raise ValueError(emsg)
jacobian = [None for p in range(self.npars)]
if (free == False).sum() == self.npars:
Expand Down Expand Up @@ -148,15 +152,17 @@ def _transform_parametersraw(self, pars, in_format, out_format):
if in_format == "internal":
pass
else:
raise ValueError("Argument 'in_format' must be one of %s." \
% self.parformats)
raise ValueError(
"Argument 'in_format' must be one of %s." % self.parformats
)

# Convert to specified output format from "internal" format.
if out_format == "internal":
pass
else:
raise ValueError("Argument 'out_format' must be one of %s." \
% self.parformats)
raise ValueError(
"Argument 'out_format' must be one of %s." % self.parformats
)
return temp

def _valueraw(self, pars, r):
Expand All @@ -171,51 +177,54 @@ def _valueraw(self, pars, r):
If degree is negative infinity, pars is an empty sequence.
r: sequence or scalar over which pars is evaluated"""
if len(pars) != self.npars:
emsg = "Argument pars must have "+str(self.npars)+" elements."
emsg = "Argument pars must have " + str(self.npars) + " elements."
raise ValueError(emsg)
return np.polyval(pars, r)

def getmodule(self):
return __name__

#end of class Polynomial

# end of class Polynomial

# simple test code
if __name__ == '__main__':
if __name__ == "__main__":

# Test polynomial of degree 3
print "Testing degree 3 polynomial"
print "---------------------------"
f = Polynomial(degree = 3)
print("Testing degree 3 polynomial")
print("---------------------------")
f = Polynomial(degree=3)
r = np.arange(5)
pars = np.array([3, 0, 1, 2])
free = np.array([True, False, True, True])
val = f._valueraw(pars, r)
jac = f._jacobianraw(pars, r, free)
print "Value:\n", val
print "Jacobian: "
for j in jac: print " %s" %j
jac = f._jacobianraw(pars, r, free)
print("Value:\n", val)
print("Jacobian: ")
for j in jac:
print(" %s" % j)

# Test polynomial of degree -oo
print "\nTesting degree -oo polynomial (== 0)"
print "------------------------------------"
f = Polynomial(degree = -1)
print("\nTesting degree -oo polynomial (== 0)")
print("------------------------------------")
f = Polynomial(degree=-1)
r = np.arange(5)
pars = np.array([])
free = np.array([])
val = f._valueraw(pars, r)
jac = f._jacobianraw(pars, r, free)
print "Value:\n", val
print "Jacobian: "
for j in jac: print " %s" %j
jac = f._jacobianraw(pars, r, free)
print("Value:\n", val)
print("Jacobian: ")
for j in jac:
print(" %s" % j)

# Test linear estimation
print "\nTesting linear baseline estimation"
print "------------------------------------"
f = Polynomial(degree = 1)
print("\nTesting linear baseline estimation")
print("------------------------------------")
f = Polynomial(degree=1)
pars = np.array([1, 0])
r = np.arange(0, 10, .1)
y = -r + 10*np.exp(-(r-5)**2) + np.random.rand(len(r))
r = np.arange(0, 10, 0.1)
y = -r + 10 * np.exp(-((r - 5) ** 2)) + np.random.rand(len(r))
est = f.estimate_parameters(r, y)
print "Actual baseline: ", np.array([-1, 0.])
print "Estimated baseline: ", est
print("Actual baseline: ", np.array([-1, 0.0]))
print("Estimated baseline: ", est)
Loading

0 comments on commit ee77d44

Please sign in to comment.