Skip to content

test #39

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

Closed
wants to merge 20 commits into from
Closed

test #39

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
20 commits
Select commit Hold shift + click to select a range
5684c18
Lint check & fix to python3 format (#18)
stevenhua0320 Jul 29, 2024
a4221ff
lint check & change to python3 & pre-commit check (#19)
stevenhua0320 Jul 29, 2024
866e714
lint check and fix print and exception python2 issue (#20)
stevenhua0320 Jul 29, 2024
c28b97b
lint check and fix python2 print and exception issues (#21)
stevenhua0320 Jul 29, 2024
ff0fc58
finish parenthesizing print statements (#24)
sbillinge Jul 29, 2024
7c77b11
fix too many leading #, import modules, and unused var (#29)
stevenhua0320 Jul 30, 2024
21d1a53
requirements (#30)
sbillinge Jul 30, 2024
181a54a
fix import module not used & string check (#25)
stevenhua0320 Jul 30, 2024
8772a95
fix too many leading "#" in string block (#26)
stevenhua0320 Jul 30, 2024
ffe5fc7
lint check, remove unused import modules & remove too many "#". (#27)
stevenhua0320 Jul 30, 2024
edaa35d
remove unused modules, ambiguous variable name (#28)
stevenhua0320 Jul 30, 2024
1f67292
cleaning (#31)
sbillinge Jul 30, 2024
30d0909
Copyright (#32)
sbillinge Jul 30, 2024
7dd69bc
lint check, fix break import modules, remove unused import modules, r…
stevenhua0320 Jul 30, 2024
81b261b
fix break import modules, remove unused import modules, fix docstring…
stevenhua0320 Jul 30, 2024
978dca1
fix formatting issue and typo in copyright (#35)
stevenhua0320 Jul 30, 2024
b61a877
clean out inits
sbillinge Jul 31, 2024
ab27182
[pre-commit.ci] auto fixes from pre-commit hooks
pre-commit-ci[bot] Jul 31, 2024
cf2e82f
dataclusters.py, modelevaluators/aicc and modelparts.py
sbillinge Jul 31, 2024
caeb1fc
Merge branch 'test' of github.com:sbillinge/diffpy.srmise into test
sbillinge Jul 31, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
lint check & change to python3 & pre-commit check (#19)
* lint check & change to python3 & pre-commit check

* [pre-commit.ci] auto fixes from pre-commit hooks

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
  • Loading branch information
stevenhua0320 and pre-commit-ci[bot] authored Jul 29, 2024
commit a4221ff9adeb66d46744f2b59d02d654098b3fea
119 changes: 67 additions & 52 deletions diffpy/srmise/baselines/nanospherical.py
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
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
Loading