Here are the examples of the python api numpy.polyfit taken from open source projects. By voting up you can indicate which examples are most useful and appropriate.
50 Examples
3
Example 1
Project: chaco Source File: regression_lasso.py
def _selection_changed_fired(self, event):
indices = self.selection_datasource.metadata["selection"]
if any(indices):
x = compress(indices, self.component.index.get_data())
y = compress(indices, self.component.value.get_data())
if len(x) < 2 or len(y) < 2:
self.fit_params = None
self.centroid = None
else:
self.fit_params = tuple(polyfit(x,y,1))
self.centroid = (sum(x)/len(x)), (sum(y)/len(y))
else:
self.fit_params = None
self.centroid = None
return
3
Example 2
Project: wordgraph Source File: analysers.py
def get_validity(self):
x_values = [point.x for point in self.points]
y_values = [point.y for point in self.points]
par = np.polyfit(x_values, y_values, 1, full=True)
self.gradient = par[0][0]
self.constant = par[0][1]
residuals = np.var([(self.gradient * point.x + self.constant - point.y)
for point in self.points])
# variance = np.var(y_values)
# Rsqr = np.round(1 - residuals / variance, decimals=2)
return 1 - residuals
3
Example 3
Project: spinalcordtoolbox Source File: sct_smooth_spinal_cord_shifting_centerline.py
Function: polynome_centerline
Function: polynome_centerline
def polynome_centerline(x_centerline,y_centerline,z_centerline):
"""Fit polynomial function through centerline"""
# Fit centerline in the Z-X plane using polynomial function
print '\nFit centerline in the Z-X plane using polynomial function...'
coeffsx = np.polyfit(z_centerline, x_centerline, deg=5)
polyx = np.poly1d(coeffsx)
x_centerline_fit = np.polyval(polyx, z_centerline)
#Fit centerline in the Z-Y plane using polynomial function
print '\nFit centerline in the Z-Y plane using polynomial function...'
coeffsy = np.polyfit(z_centerline, y_centerline, deg=5)
polyy = np.poly1d(coeffsy)
y_centerline_fit = np.polyval(polyy, z_centerline)
return x_centerline_fit,y_centerline_fit
3
Example 4
def polynome_centerline(x_centerline,y_centerline,z_centerline):
"""Fit polynomial function through centerline"""
# Fit centerline in the Z-X plane using polynomial function
print '\nFit centerline in the Z-X plane using polynomial function...'
coeffsx = numpy.polyfit(z_centerline, x_centerline, deg=param.deg_poly)
polyx = numpy.poly1d(coeffsx)
x_centerline_fit = numpy.polyval(polyx, z_centerline)
#Fit centerline in the Z-Y plane using polynomial function
print '\nFit centerline in the Z-Y plane using polynomial function...'
coeffsy = numpy.polyfit(z_centerline, y_centerline, deg=param.deg_poly)
polyy = numpy.poly1d(coeffsy)
y_centerline_fit = numpy.polyval(polyy, z_centerline)
return x_centerline_fit,y_centerline_fit,polyx,polyy
3
Example 5
def slope(self):
"""
:return: the esitmated slope for points in the current window
"""
x = range(self.window_size)
y = self.vals
slope, bias = np.polyfit(x, y, 1)
return slope
3
Example 6
Project: waveform-analyzer Source File: common.py
def parabolic_polyfit(f, x, n):
"""
Use the built-in polyfit() function to find the peak of a parabola
f is a vector and x is an index for that vector.
n is the number of samples of the curve used to fit the parabola.
"""
a, b, c = polyfit(arange(x-n//2, x+n//2+1), f[x-n//2:x+n//2+1], 2)
xv = -0.5 * b/a
yv = a * xv**2 + b * xv + c
return (xv, yv)
3
Example 7
Project: spinalcordtoolbox Source File: msct_smooth.py
def polynomial_fit(x,y,degree):
import numpy as np
coeffs = np.polyfit(x, y, degree)
poly = np.poly1d(coeffs)
y_fit = np.polyval(poly, x)
return y_fit, poly
3
Example 8
def get_profile(self):
"""
get location of profile by linear regression
"""
if self.Model_suite.x is None:
print "can't get profile, no x y locations"
return
x = self.Model_suite.x
y = self.Model_suite.y
self.profile = np.polyfit(x,y,1)
3
Example 9
def fit_parabola(xs, ys):
a, b, c = np.polyfit(xs, ys, 2)
if a < 0: raise BadFitShape()
return a, b, c
3
Example 10
def test_residplot(self):
x, y = self.df.x, self.df.y
ax = lm.residplot(x, y)
resid = y - np.polyval(np.polyfit(x, y, 1), x)
x_plot, y_plot = ax.collections[0].get_offsets().T
npt.assert_array_equal(x, x_plot)
npt.assert_array_almost_equal(resid, y_plot)
3
Example 11
Project: quant_at Source File: hurst.py
def hurst(ts):
"""Returns the Hurst Exponent of the time series vector ts"""
# Create the range of lag values
lags = range(2, 100)
# Calculate the array of the variances of the lagged differences
tau = [sqrt(std(subtract(ts[lag:], ts[:-lag]))) for lag in lags]
# Use a linear fit to estimate the Hurst Exponent
poly = polyfit(log(lags), log(tau), 1)
# Return the Hurst exponent from the polyfit output
return poly[0]*2.0
3
Example 12
def polynome_centerline(x_centerline,y_centerline,z_centerline):
"""Fit polynomial function through centerline"""
# Fit centerline in the Z-X plane using polynomial function
print '\nFit centerline in the Z-X plane using polynomial function...'
coeffsx = numpy.polyfit(z_centerline, x_centerline, deg=5)
polyx = numpy.poly1d(coeffsx)
x_centerline_fit = numpy.polyval(polyx, z_centerline)
# Fit centerline in the Z-Y plane using polynomial function
print '\nFit centerline in the Z-Y plane using polynomial function...'
coeffsy = numpy.polyfit(z_centerline, y_centerline, deg=5)
polyy = numpy.poly1d(coeffsy)
y_centerline_fit = numpy.polyval(polyy, z_centerline)
return x_centerline_fit, y_centerline_fit
3
Example 13
Project: battery-status Source File: battery-status-graph.py
def guess_expiry(x, y, zero = 0):
fit = np.polyfit(data['energy_full'], data['timestamp'], 1)
#print "fit: %s" % fit
fit_fn = np.poly1d(fit)
#print "fit_fn: %s" % fit_fn
return datetime.datetime.fromtimestamp(fit_fn(zero))
3
Example 14
Project: seisflows Source File: math.py
def lsq2(x, f):
# parabolic least squares fit
p = np.polyfit(x, f, 2)
if p[0] > 0:
return -p[1]/(2*p[0])
else:
print -1
raise Exception()
3
Example 15
Project: TheCannon Source File: snr_test.py
def quad_fit(x, y, yerr):
""" Fit a qudratic to the SNR to make a lookup error table """
qfit = np.polyfit(x, y, deg=2, w = 1 / yerr)
plt.figure()
plt.scatter(x, y)
plt.errorbar(x, y, yerr=yerr)
xvals = np.linspace(min(x), max(x), 100)
yvals = qfit[2] + qfit[1]*xvals + qfit[0]*xvals**2
plt.plot(xvals, yvals, color='r', lw=2)
plt.show()
print(qfit)
0
Example 16
Project: pymatgen Source File: elastic.py
@classmethod
def from_stress_dict(cls, stress_dict, tol=0.1, vasp=True, symmetry=False):
"""
Constructs the elastic tensor from IndependentStrain-Stress dictionary
corresponding to legacy behavior of elasticity package.
Args:
stress_dict (dict): dictionary of stresses indexed by corresponding
IndependentStrain objects.
tol (float): tolerance for zeroing small values of the tensor
vasp (boolean): flag for whether the stress tensor should be
converted based on vasp units/convention for stress
symmetry (boolean): flag for whether or not the elastic tensor
should fit from data based on symmetry
"""
inds = [(0, 0), (1, 1), (2, 2), (1, 2), (0, 2), (0, 1)]
c_ij = np.array([[np.polyfit([strain[ind1] for strain in list(stress_dict.keys())
if (strain.i, strain.j) == ind1],
[stress_dict[strain][ind2] for strain
in list(stress_dict.keys())
if (strain.i, strain.j) == ind1], 1)[0]
for ind1 in inds] for ind2 in inds])
if vasp:
c_ij *= -0.1 # Convert units/sign convention of vasp stress tensor
c_ij[0:, 3:] = 0.5 * c_ij[0:, 3:] # account for voigt doubling of e4,e5,e6
c = cls.from_voigt(c_ij)
c = c.zeroed()
return c
0
Example 17
Project: distalgo Source File: plot.py
def load_graphline(graphline, ax=plt):
"""Plot a line on the graph, return its handle."""
datafile = graphline.dataset.data_file
results = load_data(datafile)
xset = []
yset = []
for xaxis, datapoint in graphline.dataset:
datas = [(data[graphline.key] if not graphline.avg_over_procs
else data[graphline.key] / data['Total_processes'])
/ graphline.avg_factor
for config, data, ts in results if config == datapoint]
if len(datas) == 0:
print("No data for ", graphline.key, datapoint)
xset.append(xaxis)
yset.append(graphline.aggregate(datas))
if graphline.fit_degree is not None:
pol = np.polyfit(xset, yset, graphline.fit_degree)
sample_xs = np.linspace(xset[0], xset[-1], len(xset) * 10)
fitline_args = graphline.line_properties
datapoint_args = graphline.line_properties
fitline_args.pop('marker')
datapoint_args.pop('linestyle')
datapoint_args.pop('label')
ax.plot(sample_xs, np.polyval(pol, sample_xs), **fitline_args)
ax.plot(xset, yset, label='_nolegend_', linestyle="", **datapoint_args)
return Line2D([0, 1], [0, 1], **graphline.line_properties)
else:
return ax.plot(xset, yset, **graphline.line_properties)
0
Example 18
def deltafactor_polyfit(volumes, energies):
"""
This is the routine used to compute V0, B0, B1 in the deltafactor code.
"""
fitdata = np.polyfit(volumes**(-2./3.), energies, 3, full=True)
deriv0 = np.poly1d(fitdata[0])
deriv1 = np.polyder(deriv0, 1)
deriv2 = np.polyder(deriv1, 1)
deriv3 = np.polyder(deriv2, 1)
for x in np.roots(deriv1):
if x > 0 and deriv2(x) > 0:
v0 = x**(-3./2.)
break
else:
raise EOSError("No minimum could be found")
derivV2 = 4./9. * x**5. * deriv2(x)
derivV3 = (-20./9. * x**(13./2.) * deriv2(x) - 8./27. * x**(15./2.) * deriv3(x))
b0 = derivV2 / x**(3./2.)
b1 = -1 - x**(-3./2.) * derivV3 / derivV2
# e0, b0, b1, v0, eos_params
return np.poly1d(fitdata[0])(v0**(-2./3.)), b0, b1, v0, fitdata[0]
0
Example 19
def fit_parabola(self, deltas, perf, output_prefix=None, create_pdf=True):
"""
for best k, fit parabola to 3 points in delta dimension and determine optimum accordingly,
assume convex function
"""
#self.run_once() - don't call this because we don't need snp data loaded
assert len(perf) == len(deltas)
assert perf[1] <= perf[0]
assert perf[1] <= perf[2]
assert deltas[0] < deltas[1] < deltas[2]
coef = np.polyfit(deltas, perf, 2)
xfit = np.linspace(min(deltas), max(deltas), num=100)
yfit = np.polyval(coef, xfit)
best_idx = np.argmin(yfit)
best_delta = xfit[best_idx]
best_obj = yfit[best_idx]
#if create_pdf and (output_prefix!=None):
# pylab.figure()
# pylab.plot(deltas, perf, '.')
# pylab.plot(xfit, yfit, '-')
# pylab.grid(True)
# plot_fn = output_prefix + "_parabola.pdf"
# util.create_directory_if_necessary(plot_fn)
# pylab.savefig(plot_fn)
return best_delta, best_obj
0
Example 20
def simple(args):
"""
%prog simple anchorfile --qbed=qbedfile --sbed=sbedfile [options]
Write the block ends for each block in the anchorfile.
GeneA1 GeneA2 GeneB1 GeneB2 +/- score
Optional additional columns:
orderA1 orderA2 orderB1 orderB2 sizeA sizeB size block_id
With base coordinates (--coords):
block_id seqidA startA endA bpSpanA GeneA1 GeneA2 geneSpanA
block_id seqidB startB endB bpSpanB GeneB1 GeneB2 geneSpanB
"""
p = OptionParser(simple.__doc__)
p.add_option("--rich", default=False, action="store_true", \
help="Output additional columns [default: %default]")
p.add_option("--coords", default=False, action="store_true",
help="Output columns with base coordinates [default: %default]")
p.add_option("--bed", default=False, action="store_true",
help="Generate BED file for the blocks")
p.add_option("--noheader", default=False, action="store_true",
help="Don't output header [default: %default]")
p.set_beds()
opts, args = p.parse_args(args)
if len(args) != 1:
sys.exit(not p.print_help())
anchorfile, = args
additional = opts.rich
coords = opts.coords
header = not opts.noheader
bed = opts.bed
if bed:
coords = True
bbed = Bed()
ac = AnchorFile(anchorfile)
simplefile = anchorfile.rsplit(".", 1)[0] + ".simple"
qbed, sbed, qorder, sorder, is_self = check_beds(anchorfile, p, opts)
pf = "-".join(anchorfile.split(".", 2)[:2])
blocks = ac.blocks
if coords:
h = "Block|Chr|Start|End|Span|StartGene|EndGene|GeneSpan|Orientation"
else:
h = "StartGeneA|EndGeneA|StartGeneB|EndGeneB|Orientation|Score"
if additional:
h += "|StartOrderA|EndOrderA|StartOrderB|EndOrderB|"\
"SizeA|SizeB|Size|Block"
fws = open(simplefile, "w")
if header:
print >> fws, "\t".join(h.split("|"))
atotalbase = btotalbase = 0
for i, block in enumerate(blocks):
a, b, scores = zip(*block)
a = [qorder[x] for x in a]
b = [sorder[x] for x in b]
ia, oa = zip(*a)
ib, ob = zip(*b)
astarti, aendi = min(ia), max(ia)
bstarti, bendi = min(ib), max(ib)
astart, aend = min(a)[1].accn, max(a)[1].accn
bstart, bend = min(b)[1].accn, max(b)[1].accn
sizeA = len(set(ia))
sizeB = len(set(ib))
size = len(block)
slope, intercept = np.polyfit(ia, ib, 1)
orientation = "+" if slope >= 0 else '-'
aspan = aendi - astarti + 1
bspan = bendi - bstarti + 1
score = int((aspan * bspan) ** .5)
score = str(score)
block_id = pf + "-block-{0}".format(i)
if coords:
aseqid, astartbase, aendbase = \
get_boundary_bases(astart, aend, qorder)
bseqid, bstartbase, bendbase = \
get_boundary_bases(bstart, bend, sorder)
abase = aendbase - astartbase + 1
bbase = bendbase - bstartbase + 1
atotalbase += abase
btotalbase += bbase
# Write dual lines
aargs = [block_id, aseqid, astartbase, aendbase,
abase, astart, aend, aspan, "+"]
bargs = [block_id, bseqid, bstartbase, bendbase,
bbase, bstart, bend, bspan, orientation]
if bed:
bbed.append(BedLine("\t".join(str(x) for x in \
(bseqid, bstartbase - 1, bendbase,
"{}:{}-{}".format(aseqid, astartbase, aendbase),
size, orientation))))
for args in (aargs, bargs):
print >> fws, "\t".join(str(x) for x in args)
continue
args = [astart, aend, bstart, bend, score, orientation]
if additional:
args += [astarti, aendi, bstarti, bendi,
sizeA, sizeB, size, block_id]
print >> fws, "\t".join(str(x) for x in args)
fws.close()
logging.debug("A total of {0} blocks written to `{1}`.".format(i + 1, simplefile))
if coords:
print >> sys.stderr, "Total block span in {0}: {1}".format(qbed.filename, \
human_size(atotalbase, precision=2))
print >> sys.stderr, "Total block span in {0}: {1}".format(sbed.filename, \
human_size(btotalbase, precision=2))
print >> sys.stderr, "Ratio: {0:.1f}x".format(\
max(atotalbase, btotalbase) * 1. / min(atotalbase, btotalbase))
if bed:
bedfile = simplefile + ".bed"
bbed.print_to_file(filename=bedfile, sorted=True)
logging.debug("Bed file written to `{}`".format(bedfile))
0
Example 21
Project: alchemical-analysis Source File: parser_amber.py
def _extrapol(x, y, scheme):
"""Simple extrapolation schemes."""
y0 = None
y1 = None
if scheme == 'linear':
if 0.0 not in x:
y0 = (x[0] * y[1] - x[1] * y[0]) / (x[0] - x[1])
if 1.0 not in x:
y1 = (((x[-2] - 1.0) * y[-1] + ((1.0 - x[-1]) * y[-2])) /
(x[-2] - x[-1]))
elif scheme == 'polyfit':
nval = len(x)
if nval < 6:
deg = nval - 1
else:
deg = 6
coeffs = np.polyfit(x, y, deg)
if 0.0 not in x:
y0 = coeffs[-1]
if 1.0 not in x:
y1 = sum(coeffs)
else:
raise SystemExit('ERROR: Unsupported extrapolation scheme: %s' %
scheme)
return y0, y1
0
Example 22
Project: AWS-Lambda-ML-Microservice-Skeleton Source File: test_regression.py
Function: test_polyfit_build
Function: test_polyfit_build
def test_polyfit_build(self):
# Ticket #628
ref = [-1.06123820e-06, 5.70886914e-04, -1.13822012e-01,
9.95368241e+00, -3.14526520e+02]
x = [90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103,
104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115,
116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 129,
130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141,
146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157,
158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169,
170, 171, 172, 173, 174, 175, 176]
y = [9.0, 3.0, 7.0, 4.0, 4.0, 8.0, 6.0, 11.0, 9.0, 8.0, 11.0, 5.0,
6.0, 5.0, 9.0, 8.0, 6.0, 10.0, 6.0, 10.0, 7.0, 6.0, 6.0, 6.0,
13.0, 4.0, 9.0, 11.0, 4.0, 5.0, 8.0, 5.0, 7.0, 7.0, 6.0, 12.0,
7.0, 7.0, 9.0, 4.0, 12.0, 6.0, 6.0, 4.0, 3.0, 9.0, 8.0, 8.0,
6.0, 7.0, 9.0, 10.0, 6.0, 8.0, 4.0, 7.0, 7.0, 10.0, 8.0, 8.0,
6.0, 3.0, 8.0, 4.0, 5.0, 7.0, 8.0, 6.0, 6.0, 4.0, 12.0, 9.0,
8.0, 8.0, 8.0, 6.0, 7.0, 4.0, 4.0, 5.0, 7.0]
tested = np.polyfit(x, y, 4)
assert_array_almost_equal(ref, tested)
0
Example 23
Project: workload-automation Source File: __init__.py
def fit_polynomial(s, n):
# pylint: disable=no-member
coeffs = np.polyfit(s.index, s.values, n)
poly = np.poly1d(coeffs)
return poly(s.index)
0
Example 24
def detrend(self, method='linear', order=5):
"""
Detrend series data with linear or nonlinear detrending.
Preserve intercept so that subsequent operations can adjust the baseline.
Parameters
----------
method : str, optional, default = 'linear'
Detrending method
order : int, optional, default = 5
Order of polynomial, for non-linear detrending only
"""
check_options(method, ['linear', 'nonlinear'])
if method == 'linear':
order = 1
def func(y):
x = arange(len(y))
p = polyfit(x, y, order)
p[-1] = 0
yy = polyval(p, x)
return y - yy
return self.map(func)
0
Example 25
Project: spinalcordtoolbox Source File: sct_get_centerline_automatic.py
def main():
# Initialization
fname_anat = ''
fname_point = ''
slice_gap = param.gap
remove_tmp_files = param.remove_tmp_files
gaussian_kernel = param.gaussian_kernel
start_time = time.time()
verbose = 1
# get path of the toolbox
status, path_sct = commands.getstatusoutput('echo $SCT_DIR')
path_sct = sct.slash_at_the_end(path_sct, 1)
# Parameters for debug mode
if param.debug == 1:
sct.printv('\n*** WARNING: DEBUG MODE ON ***\n\t\t\tCurrent working directory: '+os.getcwd(), 'warning')
status, path_sct_testing_data = commands.getstatusoutput('echo $SCT_TESTING_DATA_DIR')
fname_anat = path_sct_testing_data+'/t2/t2.nii.gz'
fname_point = path_sct_testing_data+'/t2/t2_centerline_init.nii.gz'
slice_gap = 5
else:
# Check input param
try:
opts, args = getopt.getopt(sys.argv[1:],'hi:p:g:r:k:')
except getopt.GetoptError as err:
print str(err)
usage()
if not opts:
usage()
for opt, arg in opts:
if opt == '-h':
usage()
elif opt in ('-i'):
fname_anat = arg
elif opt in ('-p'):
fname_point = arg
elif opt in ('-g'):
slice_gap = int(arg)
elif opt in ('-r'):
remove_tmp_files = int(arg)
elif opt in ('-k'):
gaussian_kernel = int(arg)
# display usage if a mandatory argument is not provided
if fname_anat == '' or fname_point == '':
usage()
# check existence of input files
sct.check_file_exist(fname_anat)
sct.check_file_exist(fname_point)
# extract path/file/extension
path_anat, file_anat, ext_anat = sct.extract_fname(fname_anat)
path_point, file_point, ext_point = sct.extract_fname(fname_point)
# extract path of schedule file
# TODO: include schedule file in sct
# TODO: check existence of schedule file
file_schedule = path_sct + param.schedule_file
# Get input image orientation
input_image_orientation = get_orientation(fname_anat)
# Display arguments
print '\nCheck input arguments...'
print ' Anatomical image: '+fname_anat
print ' Orientation: '+input_image_orientation
print ' Point in spinal cord: '+fname_point
print ' Slice gap: '+str(slice_gap)
print ' Gaussian kernel: '+str(gaussian_kernel)
print ' Degree of polynomial: '+str(param.deg_poly)
# create temporary folder
print('\nCreate temporary folder...')
path_tmp = 'tmp.'+time.strftime("%y%m%d%H%M%S")
sct.create_folder(path_tmp)
print '\nCopy input data...'
sct.run('cp '+fname_anat+ ' '+path_tmp+'/tmp.anat'+ext_anat)
sct.run('cp '+fname_point+ ' '+path_tmp+'/tmp.point'+ext_point)
# go to temporary folder
os.chdir(path_tmp)
# convert to nii
convert('tmp.anat'+ext_anat, 'tmp.anat.nii')
convert('tmp.point'+ext_point, 'tmp.point.nii')
# Reorient input anatomical volume into RL PA IS orientation
print '\nReorient input volume to RL PA IS orientation...'
#sct.run(sct.fsloutput + 'fslswapdim tmp.anat RL PA IS tmp.anat_orient')
set_orientation('tmp.anat.nii', 'RPI', 'tmp.anat_orient.nii')
# Reorient binary point into RL PA IS orientation
print '\nReorient binary point into RL PA IS orientation...'
# sct.run(sct.fsloutput + 'fslswapdim tmp.point RL PA IS tmp.point_orient')
set_orientation('tmp.point.nii', 'RPI', 'tmp.point_orient.nii')
# Get image dimensions
print '\nGet image dimensions...'
nx, ny, nz, nt, px, py, pz, pt = Image('tmp.anat_orient.nii').dim
print '.. matrix size: '+str(nx)+' x '+str(ny)+' x '+str(nz)
print '.. voxel size: '+str(px)+'mm x '+str(py)+'mm x '+str(pz)+'mm'
# Split input volume
print '\nSplit input volume...'
split_data('tmp.anat_orient.nii', 2, '_z')
file_anat_split = ['tmp.anat_orient_z'+str(z).zfill(4) for z in range(0, nz, 1)]
split_data('tmp.point_orient.nii', 2, '_z')
file_point_split = ['tmp.point_orient_z'+str(z).zfill(4) for z in range(0, nz, 1)]
# Extract coordinates of input point
# sct.printv('\nExtract the slice corresponding to z='+str(z_init)+'...', verbose)
#
data_point = Image('tmp.point_orient.nii').data
x_init, y_init, z_init = unravel_index(data_point.argmax(), data_point.shape)
sct.printv('Coordinates of input point: ('+str(x_init)+', '+str(y_init)+', '+str(z_init)+')', verbose)
# Create 2D gaussian mask
sct.printv('\nCreate gaussian mask from point...', verbose)
xx, yy = mgrid[:nx, :ny]
mask2d = zeros((nx, ny))
radius = round(float(gaussian_kernel+1)/2) # add 1 because the radius includes the center.
sigma = float(radius)
mask2d = exp(-(((xx-x_init)**2)/(2*(sigma**2)) + ((yy-y_init)**2)/(2*(sigma**2))))
# Save mask to 2d file
file_mask_split = ['tmp.mask_orient_z'+str(z).zfill(4) for z in range(0,nz,1)]
nii_mask2d = Image('tmp.anat_orient_z0000.nii')
nii_mask2d.data = mask2d
nii_mask2d.setFileName(file_mask_split[z_init]+'.nii')
nii_mask2d.save()
#
# # Get the coordinates of the input point
# print '\nGet the coordinates of the input point...'
# data_point = Image('tmp.point_orient.nii').data
# x_init, y_init, z_init = unravel_index(data_point.argmax(), data_point.shape)
# print '('+str(x_init)+', '+str(y_init)+', '+str(z_init)+')'
# x_init, y_init, z_init = (data > 0).nonzero()
# x_init = x_init[0]
# y_init = y_init[0]
# z_init = z_init[0]
# print '('+str(x_init)+', '+str(y_init)+', '+str(z_init)+')'
#
# numpy.unravel_index(a.argmax(), a.shape)
#
# file = nibabel.load('tmp.point_orient.nii')
# data = file.get_data()
# x_init, y_init, z_init = (data > 0).nonzero()
# x_init = x_init[0]
# y_init = y_init[0]
# z_init = z_init[0]
# print '('+str(x_init)+', '+str(y_init)+', '+str(z_init)+')'
#
# # Extract the slice corresponding to z=z_init
# print '\nExtract the slice corresponding to z='+str(z_init)+'...'
# file_point_split = ['tmp.point_orient_z'+str(z).zfill(4) for z in range(0,nz,1)]
# nii = Image('tmp.point_orient.nii')
# data_crop = nii.data[:, :, z_init:z_init+1]
# nii.data = data_crop
# nii.setFileName(file_point_split[z_init]+'.nii')
# nii.save()
#
# # Create gaussian mask from point
# print '\nCreate gaussian mask from point...'
# file_mask_split = ['tmp.mask_orient_z'+str(z).zfill(4) for z in range(0,nz,1)]
# sct.run(sct.fsloutput+'fslmaths '+file_point_split[z_init]+' -s '+str(gaussian_kernel)+' '+file_mask_split[z_init])
#
# # Obtain max value from mask
# print '\nFind maximum value from mask...'
# file = nibabel.load(file_mask_split[z_init]+'.nii')
# data = file.get_data()
# max_value_mask = numpy.max(data)
# print '..'+str(max_value_mask)
#
# # Normalize mask beween 0 and 1
# print '\nNormalize mask beween 0 and 1...'
# sct.run(sct.fsloutput+'fslmaths '+file_mask_split[z_init]+' -div '+str(max_value_mask)+' '+file_mask_split[z_init])
## Take the square of the mask
#print '\nCalculate the square of the mask...'
#sct.run(sct.fsloutput+'fslmaths '+file_mask_split[z_init]+' -mul '+file_mask_split[z_init]+' '+file_mask_split[z_init])
# initialize variables
file_mat = ['tmp.mat_z'+str(z).zfill(4) for z in range(0,nz,1)]
file_mat_inv = ['tmp.mat_inv_z'+str(z).zfill(4) for z in range(0,nz,1)]
file_mat_inv_cuemul = ['tmp.mat_inv_cuemul_z'+str(z).zfill(4) for z in range(0,nz,1)]
# create identity matrix for initial transformation matrix
fid = open(file_mat_inv_cuemul[z_init], 'w')
fid.write('%i %i %i %i\n' %(1, 0, 0, 0) )
fid.write('%i %i %i %i\n' %(0, 1, 0, 0) )
fid.write('%i %i %i %i\n' %(0, 0, 1, 0) )
fid.write('%i %i %i %i\n' %(0, 0, 0, 1) )
fid.close()
# initialize centerline: give value corresponding to initial point
x_centerline = [x_init]
y_centerline = [y_init]
z_centerline = [z_init]
warning_count = 0
# go up (1), then down (2) in reference to the binary point
for iUpDown in range(1, 3):
if iUpDown == 1:
# z increases
slice_gap_signed = slice_gap
elif iUpDown == 2:
# z decreases
slice_gap_signed = -slice_gap
# reverse centerline (because values will be appended at the end)
x_centerline.reverse()
y_centerline.reverse()
z_centerline.reverse()
# initialization before looping
z_dest = z_init # point given by user
z_src = z_dest + slice_gap_signed
# continue looping if 0 < z < nz
while 0 <= z_src and z_src <= nz-1:
# print current z:
print 'z='+str(z_src)+':'
# estimate transformation
sct.run(fsloutput+'flirt -in '+file_anat_split[z_src]+' -ref '+file_anat_split[z_dest]+' -schedule '+file_schedule+ ' -verbose 0 -omat '+file_mat[z_src]+' -cost normcorr -forcescaling -inweight '+file_mask_split[z_dest]+' -refweight '+file_mask_split[z_dest])
# display transfo
status, output = sct.run('cat '+file_mat[z_src])
print output
# check if transformation is bigger than 1.5x slice_gap
tx = float(output.split()[3])
ty = float(output.split()[7])
norm_txy = linalg.norm([tx, ty],ord=2)
if norm_txy > 1.5*slice_gap:
print 'WARNING: Transformation is too large --> using previous one.'
warning_count = warning_count + 1
# if previous transformation exists, replace current one with previous one
if os.path.isfile(file_mat[z_dest]):
sct.run('cp '+file_mat[z_dest]+' '+file_mat[z_src])
# estimate inverse transformation matrix
sct.run('convert_xfm -omat '+file_mat_inv[z_src]+' -inverse '+file_mat[z_src])
# compute cuemulative transformation
sct.run('convert_xfm -omat '+file_mat_inv_cuemul[z_src]+' -concat '+file_mat_inv[z_src]+' '+file_mat_inv_cuemul[z_dest])
# apply inverse cuemulative transformation to initial gaussian mask (to put it in src space)
sct.run(fsloutput+'flirt -in '+file_mask_split[z_init]+' -ref '+file_mask_split[z_init]+' -applyxfm -init '+file_mat_inv_cuemul[z_src]+' -out '+file_mask_split[z_src])
# open inverse cuemulative transformation file and generate centerline
fid = open(file_mat_inv_cuemul[z_src])
mat = fid.read().split()
x_centerline.append(x_init + float(mat[3]))
y_centerline.append(y_init + float(mat[7]))
z_centerline.append(z_src)
#z_index = z_index+1
# define new z_dest (target slice) and new z_src (moving slice)
z_dest = z_dest + slice_gap_signed
z_src = z_src + slice_gap_signed
# Reconstruct centerline
# ====================================================================================================
# reverse back centerline (because it's been reversed once, so now all values are in the right order)
x_centerline.reverse()
y_centerline.reverse()
z_centerline.reverse()
# fit centerline in the Z-X plane using polynomial function
print '\nFit centerline in the Z-X plane using polynomial function...'
coeffsx = polyfit(z_centerline, x_centerline, deg=param.deg_poly)
polyx = poly1d(coeffsx)
x_centerline_fit = polyval(polyx, z_centerline)
# calculate RMSE
rmse = linalg.norm(x_centerline_fit-x_centerline)/sqrt( len(x_centerline) )
# calculate max absolute error
max_abs = max( abs(x_centerline_fit-x_centerline) )
print '.. RMSE (in mm): '+str(rmse*px)
print '.. Maximum absolute error (in mm): '+str(max_abs*px)
# fit centerline in the Z-Y plane using polynomial function
print '\nFit centerline in the Z-Y plane using polynomial function...'
coeffsy = polyfit(z_centerline, y_centerline, deg=param.deg_poly)
polyy = poly1d(coeffsy)
y_centerline_fit = polyval(polyy, z_centerline)
# calculate RMSE
rmse = linalg.norm(y_centerline_fit-y_centerline)/sqrt( len(y_centerline) )
# calculate max absolute error
max_abs = max( abs(y_centerline_fit-y_centerline) )
print '.. RMSE (in mm): '+str(rmse*py)
print '.. Maximum absolute error (in mm): '+str(max_abs*py)
# display
if param.debug == 1:
import matplotlib.pyplot as plt
plt.figure()
plt.plot(z_centerline,x_centerline,'.',z_centerline,x_centerline_fit,'r')
plt.legend(['Data','Polynomial Fit'])
plt.title('Z-X plane polynomial interpolation')
plt.show()
plt.figure()
plt.plot(z_centerline,y_centerline,'.',z_centerline,y_centerline_fit,'r')
plt.legend(['Data','Polynomial Fit'])
plt.title('Z-Y plane polynomial interpolation')
plt.show()
# generate full range z-values for centerline
z_centerline_full = [iz for iz in range(0, nz, 1)]
# calculate X and Y values for the full centerline
x_centerline_fit_full = polyval(polyx, z_centerline_full)
y_centerline_fit_full = polyval(polyy, z_centerline_full)
# Generate fitted transformation matrices and write centerline coordinates in text file
print '\nGenerate fitted transformation matrices and write centerline coordinates in text file...'
file_mat_inv_cuemul_fit = ['tmp.mat_inv_cuemul_fit_z'+str(z).zfill(4) for z in range(0,nz,1)]
file_mat_cuemul_fit = ['tmp.mat_cuemul_fit_z'+str(z).zfill(4) for z in range(0,nz,1)]
fid_centerline = open('tmp.centerline_coordinates.txt', 'w')
for iz in range(0, nz, 1):
# compute inverse cuemulative fitted transformation matrix
fid = open(file_mat_inv_cuemul_fit[iz], 'w')
fid.write('%i %i %i %f\n' %(1, 0, 0, x_centerline_fit_full[iz]-x_init) )
fid.write('%i %i %i %f\n' %(0, 1, 0, y_centerline_fit_full[iz]-y_init) )
fid.write('%i %i %i %i\n' %(0, 0, 1, 0) )
fid.write('%i %i %i %i\n' %(0, 0, 0, 1) )
fid.close()
# compute forward cuemulative fitted transformation matrix
sct.run('convert_xfm -omat '+file_mat_cuemul_fit[iz]+' -inverse '+file_mat_inv_cuemul_fit[iz])
# write centerline coordinates in x, y, z format
fid_centerline.write('%f %f %f\n' %(x_centerline_fit_full[iz], y_centerline_fit_full[iz], z_centerline_full[iz]) )
fid_centerline.close()
# Prepare output data
# ====================================================================================================
# write centerline as text file
for iz in range(0, nz, 1):
# compute inverse cuemulative fitted transformation matrix
fid = open(file_mat_inv_cuemul_fit[iz], 'w')
fid.write('%i %i %i %f\n' %(1, 0, 0, x_centerline_fit_full[iz]-x_init) )
fid.write('%i %i %i %f\n' %(0, 1, 0, y_centerline_fit_full[iz]-y_init) )
fid.write('%i %i %i %i\n' %(0, 0, 1, 0) )
fid.write('%i %i %i %i\n' %(0, 0, 0, 1) )
fid.close()
# write polynomial coefficients
savetxt('tmp.centerline_polycoeffs_x.txt',coeffsx)
savetxt('tmp.centerline_polycoeffs_y.txt',coeffsy)
# apply transformations to data
print '\nApply fitted transformation matrices...'
file_anat_split_fit = ['tmp.anat_orient_fit_z'+str(z).zfill(4) for z in range(0,nz,1)]
file_mask_split_fit = ['tmp.mask_orient_fit_z'+str(z).zfill(4) for z in range(0,nz,1)]
file_point_split_fit = ['tmp.point_orient_fit_z'+str(z).zfill(4) for z in range(0,nz,1)]
for iz in range(0, nz, 1):
# forward cuemulative transformation to data
sct.run(fsloutput+'flirt -in '+file_anat_split[iz]+' -ref '+file_anat_split[iz]+' -applyxfm -init '+file_mat_cuemul_fit[iz]+' -out '+file_anat_split_fit[iz])
# inverse cuemulative transformation to mask
sct.run(fsloutput+'flirt -in '+file_mask_split[z_init]+' -ref '+file_mask_split[z_init]+' -applyxfm -init '+file_mat_inv_cuemul_fit[iz]+' -out '+file_mask_split_fit[iz])
# inverse cuemulative transformation to point
sct.run(fsloutput+'flirt -in '+file_point_split[z_init]+' -ref '+file_point_split[z_init]+' -applyxfm -init '+file_mat_inv_cuemul_fit[iz]+' -out '+file_point_split_fit[iz]+' -interp nearestneighbour')
# Merge into 4D volume
print '\nMerge into 4D volume...'
# sct.run(fsloutput+'fslmerge -z tmp.anat_orient_fit tmp.anat_orient_fit_z*')
# sct.run(fsloutput+'fslmerge -z tmp.mask_orient_fit tmp.mask_orient_fit_z*')
# sct.run(fsloutput+'fslmerge -z tmp.point_orient_fit tmp.point_orient_fit_z*')
concat_data(glob.glob('tmp.anat_orient_fit_z*.nii'), 'tmp.anat_orient_fit.nii', dim=2)
concat_data(glob.glob('tmp.mask_orient_fit_z*.nii'), 'tmp.mask_orient_fit.nii', dim=2)
concat_data(glob.glob('tmp.point_orient_fit_z*.nii'), 'tmp.point_orient_fit.nii', dim=2)
# Copy header geometry from input data
print '\nCopy header geometry from input data...'
copy_header('tmp.anat_orient.nii', 'tmp.anat_orient_fit.nii')
copy_header('tmp.anat_orient.nii', 'tmp.mask_orient_fit.nii')
copy_header('tmp.anat_orient.nii', 'tmp.point_orient_fit.nii')
# Reorient outputs into the initial orientation of the input image
print '\nReorient the centerline into the initial orientation of the input image...'
set_orientation('tmp.point_orient_fit.nii', input_image_orientation, 'tmp.point_orient_fit.nii')
set_orientation('tmp.mask_orient_fit.nii', input_image_orientation, 'tmp.mask_orient_fit.nii')
# Generate output file (in current folder)
print '\nGenerate output file (in current folder)...'
os.chdir('..') # come back to parent folder
#sct.generate_output_file('tmp.centerline_polycoeffs_x.txt','./','centerline_polycoeffs_x','.txt')
#sct.generate_output_file('tmp.centerline_polycoeffs_y.txt','./','centerline_polycoeffs_y','.txt')
#sct.generate_output_file('tmp.centerline_coordinates.txt','./','centerline_coordinates','.txt')
#sct.generate_output_file('tmp.anat_orient.nii','./',file_anat+'_rpi',ext_anat)
#sct.generate_output_file('tmp.anat_orient_fit.nii', file_anat+'_rpi_align'+ext_anat)
#sct.generate_output_file('tmp.mask_orient_fit.nii', file_anat+'_mask'+ext_anat)
fname_output_centerline = sct.generate_output_file(path_tmp+'/tmp.point_orient_fit.nii', file_anat+'_centerline'+ext_anat)
# Delete temporary files
if remove_tmp_files == 1:
print '\nRemove temporary files...'
sct.run('rm -rf '+path_tmp)
# print number of warnings
print '\nNumber of warnings: '+str(warning_count)+' (if >10, you should probably reduce the gap and/or increase the kernel size'
# display elapsed time
elapsed_time = time.time() - start_time
print '\nFinished! \n\tGenerated file: '+fname_output_centerline+'\n\tElapsed time: '+str(int(round(elapsed_time)))+'s\n'
0
Example 26
def BdRate(group1, group2):
"""Compute the BD-rate between two score groups.
The returned object also contains the range of PSNR values used
to compute the result.
Bjontegaard's metric allows to compute the average % saving in bitrate
between two rate-distortion curves [1].
rate1,psnr1 - RD points for curve 1
rate2,psnr2 - RD points for curve 2
adapted from code from: (c) 2010 Giuseppe Valenzise
copied from code by [email protected], [email protected]
"""
# pylint: disable=too-many-locals
metric_set1 = group1.dataPoints()
metric_set2 = group2.dataPoints()
# numpy plays games with its exported functions.
# pylint: disable=no-member
# pylint: disable=bad-builtin
psnr1 = [x[1] for x in metric_set1]
psnr2 = [x[1] for x in metric_set2]
log_rate1 = map(math.log, [x[0] for x in metric_set1])
log_rate2 = map(math.log, [x[0] for x in metric_set2])
# Best cubic poly fit for graph represented by log_ratex, psrn_x.
poly1 = numpy.polyfit(psnr1, log_rate1, 3)
poly2 = numpy.polyfit(psnr2, log_rate2, 3)
# Integration interval.
min_int = max([min(psnr1), min(psnr2)])
max_int = min([max(psnr1), max(psnr2)])
# find integral
p_int1 = numpy.polyint(poly1)
p_int2 = numpy.polyint(poly2)
# Calculate the integrated value over the interval we care about.
int1 = numpy.polyval(p_int1, max_int) - numpy.polyval(p_int1, min_int)
int2 = numpy.polyval(p_int2, max_int) - numpy.polyval(p_int2, min_int)
# Calculate the average improvement.
avg_exp_diff = (int2 - int1) / (max_int - min_int)
# In really bad formed data the exponent can grow too large.
# clamp it.
if avg_exp_diff > 200:
avg_exp_diff = 200
# Convert to a percentage.
avg_diff = (math.exp(avg_exp_diff) - 1) * 100
return {'difference': avg_diff, 'psnr':[min_int, max_int]}
0
Example 27
Project: iris Source File: polynomial_fit.py
def main():
# Enable a future option, to ensure that the netcdf load works the same way
# as in future Iris versions.
iris.FUTURE.netcdf_promote = True
# Load some test data.
fname = iris.sample_data_path('A1B_north_america.nc')
cube = iris.load_cube(fname)
# Extract a single time series at a latitude and longitude point.
location = next(cube.slices(['time']))
# Calculate a polynomial fit to the data at this time series.
x_points = location.coord('time').points
y_points = location.data
degree = 2
p = np.polyfit(x_points, y_points, degree)
y_fitted = np.polyval(p, x_points)
# Add the polynomial fit values to the time series to take
# full advantage of Iris plotting functionality.
long_name = 'degree_{}_polynomial_fit_of_{}'.format(degree, cube.name())
fit = iris.coords.AuxCoord(y_fitted, long_name=long_name,
units=location.units)
location.add_aux_coord(fit, 0)
qplt.plot(location.coord('time'), location, label='data')
qplt.plot(location.coord('time'),
location.coord(long_name),
'g-', label='polynomial fit')
plt.legend(loc='best')
plt.title('Trend of US air temperature over time')
qplt.show()
0
Example 28
Project: compare-codecs Source File: visual_metrics.py
def bdsnr(metric_set1, metric_set2):
"""
BJONTEGAARD Bjontegaard metric calculation
Bjontegaard's metric allows to compute the average gain in psnr between two
rate-distortion curves [1].
rate1,psnr1 - RD points for curve 1
rate2,psnr2 - RD points for curve 2
returns the calculated Bjontegaard metric 'dsnr'
code adapted from code written by : (c) 2010 Giuseppe Valenzise
http://www.mathworks.com/matlabcentral/fileexchange/27798-bjontegaard-metric/content/bjontegaard.m
"""
# pylint: disable=too-many-locals
# numpy seems to do tricks with its exports.
# pylint: disable=no-member
# map() is recommended against.
# pylint: disable=bad-builtin
rate1 = [x[0] for x in metric_set1]
psnr1 = [x[1] for x in metric_set1]
rate2 = [x[0] for x in metric_set2]
psnr2 = [x[1] for x in metric_set2]
log_rate1 = map(math.log, rate1)
log_rate2 = map(math.log, rate2)
# Best cubic poly fit for graph represented by log_ratex, psrn_x.
poly1 = numpy.polyfit(log_rate1, psnr1, 3)
poly2 = numpy.polyfit(log_rate2, psnr2, 3)
# Integration interval.
min_int = max([min(log_rate1), min(log_rate2)])
max_int = min([max(log_rate1), max(log_rate2)])
# Integrate poly1, and poly2.
p_int1 = numpy.polyint(poly1)
p_int2 = numpy.polyint(poly2)
# Calculate the integrated value over the interval we care about.
int1 = numpy.polyval(p_int1, max_int) - numpy.polyval(p_int1, min_int)
int2 = numpy.polyval(p_int2, max_int) - numpy.polyval(p_int2, min_int)
# Calculate the average improvement.
if max_int != min_int:
avg_diff = (int2 - int1) / (max_int - min_int)
else:
avg_diff = 0.0
return avg_diff
0
Example 29
Project: spinalcordtoolbox Source File: sct_get_centerline.py
def get_centerline_from_point(input_image, point_file, gap=4, gaussian_kernel=4, remove_tmp_files=1):
# Initialization
fname_anat = input_image
fname_point = point_file
slice_gap = gap
remove_tmp_files = remove_tmp_files
gaussian_kernel = gaussian_kernel
start_time = time()
verbose = 1
# get path of the toolbox
status, path_sct = commands.getstatusoutput('echo $SCT_DIR')
path_sct = sct.slash_at_the_end(path_sct, 1)
# Parameters for debug mode
if param.debug == 1:
sct.printv('\n*** WARNING: DEBUG MODE ON ***\n\t\t\tCurrent working directory: '+os.getcwd(), 'warning')
status, path_sct_testing_data = commands.getstatusoutput('echo $SCT_TESTING_DATA_DIR')
fname_anat = path_sct_testing_data+'/t2/t2.nii.gz'
fname_point = path_sct_testing_data+'/t2/t2_centerline_init.nii.gz'
slice_gap = 5
# check existence of input files
sct.check_file_exist(fname_anat)
sct.check_file_exist(fname_point)
# extract path/file/extension
path_anat, file_anat, ext_anat = sct.extract_fname(fname_anat)
path_point, file_point, ext_point = sct.extract_fname(fname_point)
# extract path of schedule file
# TODO: include schedule file in sct
# TODO: check existence of schedule file
file_schedule = path_sct + param.schedule_file
# Get input image orientation
input_image_orientation = get_orientation_3d(fname_anat, filename=True)
# Display arguments
print '\nCheck input arguments...'
print ' Anatomical image: '+fname_anat
print ' Orientation: '+input_image_orientation
print ' Point in spinal cord: '+fname_point
print ' Slice gap: '+str(slice_gap)
print ' Gaussian kernel: '+str(gaussian_kernel)
print ' Degree of polynomial: '+str(param.deg_poly)
# create temporary folder
print('\nCreate temporary folder...')
path_tmp = 'tmp.'+strftime('%y%m%d%H%M%S')
sct.create_folder(path_tmp)
print '\nCopy input data...'
sct.run('cp '+fname_anat+ ' '+path_tmp+'/tmp.anat'+ext_anat)
sct.run('cp '+fname_point+ ' '+path_tmp+'/tmp.point'+ext_point)
# go to temporary folder
os.chdir(path_tmp)
# convert to nii
im_anat = convert('tmp.anat'+ext_anat, 'tmp.anat.nii')
im_point = convert('tmp.point'+ext_point, 'tmp.point.nii')
# Reorient input anatomical volume into RL PA IS orientation
print '\nReorient input volume to RL PA IS orientation...'
set_orientation(im_anat, 'RPI')
im_anat.setFileName('tmp.anat_orient.nii')
# Reorient binary point into RL PA IS orientation
print '\nReorient binary point into RL PA IS orientation...'
# sct.run(sct.fsloutput + 'fslswapdim tmp.point RL PA IS tmp.point_orient')
set_orientation(im_point, 'RPI')
im_point.setFileName('tmp.point_orient.nii')
# Get image dimensions
print '\nGet image dimensions...'
nx, ny, nz, nt, px, py, pz, pt = Image('tmp.anat_orient.nii').dim
print '.. matrix size: '+str(nx)+' x '+str(ny)+' x '+str(nz)
print '.. voxel size: '+str(px)+'mm x '+str(py)+'mm x '+str(pz)+'mm'
# Split input volume
print '\nSplit input volume...'
im_anat_split_list = split_data(im_anat, 2)
file_anat_split = []
for im in im_anat_split_list:
file_anat_split.append(im.absolutepath)
im.save()
im_point_split_list = split_data(im_point, 2)
file_point_split = []
for im in im_point_split_list:
file_point_split.append(im.absolutepath)
im.save()
# Extract coordinates of input point
data_point = Image('tmp.point_orient.nii').data
x_init, y_init, z_init = unravel_index(data_point.argmax(), data_point.shape)
sct.printv('Coordinates of input point: ('+str(x_init)+', '+str(y_init)+', '+str(z_init)+')', verbose)
# Create 2D gaussian mask
sct.printv('\nCreate gaussian mask from point...', verbose)
xx, yy = mgrid[:nx, :ny]
mask2d = zeros((nx, ny))
radius = round(float(gaussian_kernel+1)/2) # add 1 because the radius includes the center.
sigma = float(radius)
mask2d = exp(-(((xx-x_init)**2)/(2*(sigma**2)) + ((yy-y_init)**2)/(2*(sigma**2))))
# Save mask to 2d file
file_mask_split = ['tmp.mask_orient_Z'+str(z).zfill(4) for z in range(0, nz, 1)]
nii_mask2d = Image('tmp.anat_orient_Z0000.nii')
nii_mask2d.data = mask2d
nii_mask2d.setFileName(file_mask_split[z_init]+'.nii')
nii_mask2d.save()
# initialize variables
file_mat = ['tmp.mat_Z'+str(z).zfill(4) for z in range(0, nz, 1)]
file_mat_inv = ['tmp.mat_inv_Z'+str(z).zfill(4) for z in range(0, nz, 1)]
file_mat_inv_cuemul = ['tmp.mat_inv_cuemul_Z'+str(z).zfill(4) for z in range(0, nz, 1)]
# create identity matrix for initial transformation matrix
fid = open(file_mat_inv_cuemul[z_init], 'w')
fid.write('%i %i %i %i\n' % (1, 0, 0, 0))
fid.write('%i %i %i %i\n' % (0, 1, 0, 0))
fid.write('%i %i %i %i\n' % (0, 0, 1, 0))
fid.write('%i %i %i %i\n' % (0, 0, 0, 1))
fid.close()
# initialize centerline: give value corresponding to initial point
x_centerline = [x_init]
y_centerline = [y_init]
z_centerline = [z_init]
warning_count = 0
# go up (1), then down (2) in reference to the binary point
for iUpDown in range(1, 3):
if iUpDown == 1:
# z increases
slice_gap_signed = slice_gap
elif iUpDown == 2:
# z decreases
slice_gap_signed = -slice_gap
# reverse centerline (because values will be appended at the end)
x_centerline.reverse()
y_centerline.reverse()
z_centerline.reverse()
# initialization before looping
z_dest = z_init # point given by user
z_src = z_dest + slice_gap_signed
# continue looping if 0 <= z < nz
while 0 <= z_src < nz:
# print current z:
print 'z='+str(z_src)+':'
# estimate transformation
sct.run(fsloutput+'flirt -in '+file_anat_split[z_src]+' -ref '+file_anat_split[z_dest]+' -schedule ' +
file_schedule + ' -verbose 0 -omat ' + file_mat[z_src] +
' -cost normcorr -forcescaling -inweight ' + file_mask_split[z_dest] +
' -refweight '+file_mask_split[z_dest])
# display transfo
status, output = sct.run('cat '+file_mat[z_src])
print output
# check if transformation is bigger than 1.5x slice_gap
tx = float(output.split()[3])
ty = float(output.split()[7])
norm_txy = linalg.norm([tx, ty], ord=2)
if norm_txy > 1.5*slice_gap:
print 'WARNING: Transformation is too large --> using previous one.'
warning_count = warning_count + 1
# if previous transformation exists, replace current one with previous one
if os.path.isfile(file_mat[z_dest]):
sct.run('cp '+file_mat[z_dest]+' '+file_mat[z_src])
# estimate inverse transformation matrix
sct.run('convert_xfm -omat '+file_mat_inv[z_src]+' -inverse '+file_mat[z_src])
# compute cuemulative transformation
sct.run('convert_xfm -omat '+file_mat_inv_cuemul[z_src]+' -concat '+file_mat_inv[z_src]+' '+file_mat_inv_cuemul[z_dest])
# apply inverse cuemulative transformation to initial gaussian mask (to put it in src space)
sct.run(fsloutput+'flirt -in '+file_mask_split[z_init]+' -ref '+file_mask_split[z_init]+' -applyxfm -init '+file_mat_inv_cuemul[z_src]+' -out '+file_mask_split[z_src])
# open inverse cuemulative transformation file and generate centerline
fid = open(file_mat_inv_cuemul[z_src])
mat = fid.read().split()
x_centerline.append(x_init + float(mat[3]))
y_centerline.append(y_init + float(mat[7]))
z_centerline.append(z_src)
#z_index = z_index+1
# define new z_dest (target slice) and new z_src (moving slice)
z_dest = z_dest + slice_gap_signed
z_src = z_src + slice_gap_signed
# Reconstruct centerline
# ====================================================================================================
# reverse back centerline (because it's been reversed once, so now all values are in the right order)
x_centerline.reverse()
y_centerline.reverse()
z_centerline.reverse()
# fit centerline in the Z-X plane using polynomial function
print '\nFit centerline in the Z-X plane using polynomial function...'
coeffsx = polyfit(z_centerline, x_centerline, deg=param.deg_poly)
polyx = poly1d(coeffsx)
x_centerline_fit = polyval(polyx, z_centerline)
# calculate RMSE
rmse = linalg.norm(x_centerline_fit-x_centerline)/sqrt( len(x_centerline) )
# calculate max absolute error
max_abs = max(abs(x_centerline_fit-x_centerline))
print '.. RMSE (in mm): '+str(rmse*px)
print '.. Maximum absolute error (in mm): '+str(max_abs*px)
# fit centerline in the Z-Y plane using polynomial function
print '\nFit centerline in the Z-Y plane using polynomial function...'
coeffsy = polyfit(z_centerline, y_centerline, deg=param.deg_poly)
polyy = poly1d(coeffsy)
y_centerline_fit = polyval(polyy, z_centerline)
# calculate RMSE
rmse = linalg.norm(y_centerline_fit-y_centerline)/sqrt( len(y_centerline) )
# calculate max absolute error
max_abs = max( abs(y_centerline_fit-y_centerline) )
print '.. RMSE (in mm): '+str(rmse*py)
print '.. Maximum absolute error (in mm): '+str(max_abs*py)
# display
if param.debug == 1:
import matplotlib.pyplot as plt
plt.figure()
plt.plot(z_centerline,x_centerline,'.',z_centerline,x_centerline_fit,'r')
plt.legend(['Data','Polynomial Fit'])
plt.title('Z-X plane polynomial interpolation')
plt.show()
plt.figure()
plt.plot(z_centerline,y_centerline,'.',z_centerline,y_centerline_fit,'r')
plt.legend(['Data','Polynomial Fit'])
plt.title('Z-Y plane polynomial interpolation')
plt.show()
# generate full range z-values for centerline
z_centerline_full = [iz for iz in range(0, nz, 1)]
# calculate X and Y values for the full centerline
x_centerline_fit_full = polyval(polyx, z_centerline_full)
y_centerline_fit_full = polyval(polyy, z_centerline_full)
# Generate fitted transformation matrices and write centerline coordinates in text file
print '\nGenerate fitted transformation matrices and write centerline coordinates in text file...'
file_mat_inv_cuemul_fit = ['tmp.mat_inv_cuemul_fit_z'+str(z).zfill(4) for z in range(0,nz,1)]
file_mat_cuemul_fit = ['tmp.mat_cuemul_fit_z'+str(z).zfill(4) for z in range(0,nz,1)]
fid_centerline = open('tmp.centerline_coordinates.txt', 'w')
for iz in range(0, nz, 1):
# compute inverse cuemulative fitted transformation matrix
fid = open(file_mat_inv_cuemul_fit[iz], 'w')
fid.write('%i %i %i %f\n' % (1, 0, 0, x_centerline_fit_full[iz]-x_init))
fid.write('%i %i %i %f\n' % (0, 1, 0, y_centerline_fit_full[iz]-y_init))
fid.write('%i %i %i %i\n' % (0, 0, 1, 0))
fid.write('%i %i %i %i\n' % (0, 0, 0, 1))
fid.close()
# compute forward cuemulative fitted transformation matrix
sct.run('convert_xfm -omat '+file_mat_cuemul_fit[iz]+' -inverse '+file_mat_inv_cuemul_fit[iz])
# write centerline coordinates in x, y, z format
fid_centerline.write('%f %f %f\n' %(x_centerline_fit_full[iz], y_centerline_fit_full[iz], z_centerline_full[iz]) )
fid_centerline.close()
# Prepare output data
# ====================================================================================================
# write centerline as text file
for iz in range(0, nz, 1):
# compute inverse cuemulative fitted transformation matrix
fid = open(file_mat_inv_cuemul_fit[iz], 'w')
fid.write('%i %i %i %f\n' % (1, 0, 0, x_centerline_fit_full[iz]-x_init))
fid.write('%i %i %i %f\n' % (0, 1, 0, y_centerline_fit_full[iz]-y_init))
fid.write('%i %i %i %i\n' % (0, 0, 1, 0))
fid.write('%i %i %i %i\n' % (0, 0, 0, 1))
fid.close()
# write polynomial coefficients
savetxt('tmp.centerline_polycoeffs_x.txt',coeffsx)
savetxt('tmp.centerline_polycoeffs_y.txt',coeffsy)
# apply transformations to data
print '\nApply fitted transformation matrices...'
file_anat_split_fit = ['tmp.anat_orient_fit_z'+str(z).zfill(4) for z in range(0, nz, 1)]
file_mask_split_fit = ['tmp.mask_orient_fit_z'+str(z).zfill(4) for z in range(0, nz, 1)]
file_point_split_fit = ['tmp.point_orient_fit_z'+str(z).zfill(4) for z in range(0, nz, 1)]
for iz in range(0, nz, 1):
# forward cuemulative transformation to data
sct.run(fsloutput+'flirt -in '+file_anat_split[iz]+' -ref '+file_anat_split[iz]+' -applyxfm -init '+file_mat_cuemul_fit[iz]+' -out '+file_anat_split_fit[iz])
# inverse cuemulative transformation to mask
sct.run(fsloutput+'flirt -in '+file_mask_split[z_init]+' -ref '+file_mask_split[z_init]+' -applyxfm -init '+file_mat_inv_cuemul_fit[iz]+' -out '+file_mask_split_fit[iz])
# inverse cuemulative transformation to point
sct.run(fsloutput+'flirt -in '+file_point_split[z_init]+' -ref '+file_point_split[z_init]+' -applyxfm -init '+file_mat_inv_cuemul_fit[iz]+' -out '+file_point_split_fit[iz]+' -interp nearestneighbour')
# Merge into 4D volume
print '\nMerge into 4D volume...'
# im_anat_list = [Image(fname) for fname in glob.glob('tmp.anat_orient_fit_z*.nii')]
fname_anat_list = glob.glob('tmp.anat_orient_fit_z*.nii')
im_anat_concat = concat_data(fname_anat_list, 2)
im_anat_concat.setFileName('tmp.anat_orient_fit.nii')
im_anat_concat.save()
# im_mask_list = [Image(fname) for fname in glob.glob('tmp.mask_orient_fit_z*.nii')]
fname_mask_list = glob.glob('tmp.mask_orient_fit_z*.nii')
im_mask_concat = concat_data(fname_mask_list, 2)
im_mask_concat.setFileName('tmp.mask_orient_fit.nii')
im_mask_concat.save()
# im_point_list = [Image(fname) for fname in glob.glob('tmp.point_orient_fit_z*.nii')]
fname_point_list = glob.glob('tmp.point_orient_fit_z*.nii')
im_point_concat = concat_data(fname_point_list, 2)
im_point_concat.setFileName('tmp.point_orient_fit.nii')
im_point_concat.save()
# Copy header geometry from input data
print '\nCopy header geometry from input data...'
im_anat = Image('tmp.anat_orient.nii')
im_anat_orient_fit = Image('tmp.anat_orient_fit.nii')
im_mask_orient_fit = Image('tmp.mask_orient_fit.nii')
im_point_orient_fit = Image('tmp.point_orient_fit.nii')
im_anat_orient_fit = copy_header(im_anat, im_anat_orient_fit)
im_mask_orient_fit = copy_header(im_anat, im_mask_orient_fit)
im_point_orient_fit = copy_header(im_anat, im_point_orient_fit)
for im in [im_anat_orient_fit, im_mask_orient_fit, im_point_orient_fit]:
im.save()
# Reorient outputs into the initial orientation of the input image
print '\nReorient the centerline into the initial orientation of the input image...'
set_orientation('tmp.point_orient_fit.nii', input_image_orientation, 'tmp.point_orient_fit.nii')
set_orientation('tmp.mask_orient_fit.nii', input_image_orientation, 'tmp.mask_orient_fit.nii')
# Generate output file (in current folder)
print '\nGenerate output file (in current folder)...'
os.chdir('..') # come back to parent folder
fname_output_centerline = sct.generate_output_file(path_tmp+'/tmp.point_orient_fit.nii', file_anat+'_centerline'+ext_anat)
# Delete temporary files
if remove_tmp_files == 1:
print '\nRemove temporary files...'
sct.run('rm -rf '+path_tmp, error_exit='warning')
# print number of warnings
print '\nNumber of warnings: '+str(warning_count)+' (if >10, you should probably reduce the gap and/or increase the kernel size'
# display elapsed time
elapsed_time = time() - start_time
print '\nFinished! \n\tGenerated file: '+fname_output_centerline+'\n\tElapsed time: '+str(int(round(elapsed_time)))+'s\n'
0
Example 30
Project: librosa Source File: spectral.py
def poly_features(y=None, sr=22050, S=None, n_fft=2048, hop_length=512,
order=1, freq=None):
'''Get coefficients of fitting an nth-order polynomial to the columns
of a spectrogram.
Parameters
----------
y : np.ndarray [shape=(n,)] or None
audio time series
sr : number > 0 [scalar]
audio sampling rate of `y`
S : np.ndarray [shape=(d, t)] or None
(optional) spectrogram magnitude
n_fft : int > 0 [scalar]
FFT window size
hop_length : int > 0 [scalar]
hop length for STFT. See `librosa.core.stft` for details.
order : int > 0
order of the polynomial to fit
freq : None or np.ndarray [shape=(d,) or shape=(d, t)]
Center frequencies for spectrogram bins.
If `None`, then FFT bin center frequencies are used.
Otherwise, it can be a single array of `d` center frequencies,
or a matrix of center frequencies as constructed by
`librosa.core.ifgram`
Returns
-------
coefficients : np.ndarray [shape=(order+1, t)]
polynomial coefficients for each frame
Examples
--------
>>> y, sr = librosa.load(librosa.util.example_audio_file())
>>> S = np.abs(librosa.stft(y))
Line features
>>> line = librosa.feature.poly_features(S=S, sr=sr)
>>> line
array([[ -2.406e-08, -5.051e-06, ..., -1.103e-08, -5.651e-09],
[ 3.445e-04, 3.834e-02, ..., 2.661e-04, 2.239e-04]])
Quadratic features
>>> quad = librosa.feature.poly_features(S=S, order=2)
>>> quad
array([[ 6.276e-12, 2.010e-09, ..., 1.493e-12, 1.000e-13],
[ -9.325e-08, -2.721e-05, ..., -2.749e-08, -6.754e-09],
[ 4.715e-04, 7.902e-02, ..., 2.963e-04, 2.259e-04]])
Plot the results for comparison
>>> import matplotlib.pyplot as plt
>>> plt.figure()
>>> plt.subplot(3, 1, 1)
>>> librosa.display.specshow(line)
>>> plt.colorbar()
>>> plt.title('Line coefficients')
>>> plt.subplot(3, 1, 2)
>>> librosa.display.specshow(quad)
>>> plt.colorbar()
>>> plt.title('Quadratic coefficients')
>>> plt.subplot(3, 1, 3)
>>> librosa.display.specshow(librosa.logamplitude(S**2, ref_power=np.max),
... y_axis='log', x_axis='time')
>>> plt.title('log Power spectrogram')
>>> plt.colorbar(format='%+2.0f dB')
>>> plt.tight_layout()
'''
S, n_fft = _spectrogram(y=y, S=S, n_fft=n_fft, hop_length=hop_length)
# Compute the center frequencies of each bin
if freq is None:
freq = fft_frequencies(sr=sr, n_fft=n_fft)
# If frequencies are constant over frames, then we only need to fit once
if freq.ndim == 1:
coefficients = np.polyfit(freq, S, order)
else:
# Else, fit each frame independently and stack the results
coefficients = np.concatenate([[np.polyfit(freq[:, i], S[:, i], order)]
for i in range(S.shape[1])], axis=0).T
return coefficients
0
Example 31
def deltafactor_polyfit(volumes, energies):
"""
This is the routine used to compute V0, B0, B1 in the deltafactor code.
Taken from deltafactor/eosfit.py
"""
fitdata = np.polyfit(volumes**(-2./3.), energies, 3, full=True)
ssr = fitdata[1]
sst = np.sum((energies - np.average(energies))**2.)
residuals0 = ssr/sst
deriv0 = np.poly1d(fitdata[0])
deriv1 = np.polyder(deriv0, 1)
deriv2 = np.polyder(deriv1, 1)
deriv3 = np.polyder(deriv2, 1)
v0 = 0
x = 0
for x in np.roots(deriv1):
if x > 0 and deriv2(x) > 0:
v0 = x**(-3./2.)
break
else:
raise EOSError("No minimum could be found")
derivV2 = 4./9. * x**5. * deriv2(x)
derivV3 = (-20./9. * x**(13./2.) * deriv2(x) - 8./27. * x**(15./2.) * deriv3(x))
b0 = derivV2 / x**(3./2.)
b1 = -1 - x**(-3./2.) * derivV3 / derivV2
#print('deltafactor polyfit:')
#print('e0, b0, b1, v0')
#print(fitdata[0], b0, b1, v0)
n = collections.namedtuple("DeltaFitResults", "v0 b0 b1 poly1d")
return n(v0, b0, b1, fitdata[0])
0
Example 32
def get_extent(self, i, order, debug=True):
# Some blocks file, such as ones manually edited, will have garbled
# order, which prompts the hack below
acol = [order[x][0] for x in self.columns[0] if x in order]
bcol = [order[x][0] for x in self.columns[i] if x in order]
elen = min(len(acol), len(bcol))
ia, ib = acol[:elen], bcol[:elen]
slope, intercept = np.polyfit(ia, ib, 1)
orientation = '+' if slope >= 0 else '-'
ocol = [order[x] for x in self.columns[i] if x in order]
#orientation = '+' if ocol[0][0] <= ocol[-1][0] else '-'
si, start = min(ocol)
ei, end = max(ocol)
same_chr = (start.seqid == end.seqid)
chr = start.seqid if same_chr else None
ngenes = ei - si + 1
if debug:
r = "{0}:{1}-{2}".format(chr, start.start, end.end)
print >> sys.stderr, "Column {0}: {1} - {2} ({3})".\
format(i, start.accn, end.accn, r)
print >> sys.stderr, " {0} .. {1} ({2}) features .. {3}".\
format(chr, ngenes, len(ocol), orientation)
span = abs(start.start - end.end)
return start, end, si, ei, chr, orientation, span
0
Example 33
def initial_guess(self):
"""
Quadratic fit to get an initial guess for the parameters.
Returns:
tuple: (minimum energy(e0), buk modulus(b0),
derivative of bulk modulus wrt pressure(b1),
minimum volume(v0))
"""
a, b, c = np.polyfit(self.volumes, self.energies, 2)
self.eos_params = [a, b, c]
v0 = -b / (2 * a)
e0 = a * v0 ** 2 + b * v0 + c
b0 = 2 * a * v0
b1 = 4 # b1 is usually a small number like 4
vmin, vmax = self.volumes.min(), self.volumes.max()
if not vmin < v0 and v0 < vmax:
raise EOSError('The minimum volume of a fitted parabola is '
'not in the input volumes\n.')
return e0, b0, b1, v0
0
Example 34
Project: pandashells Source File: p_regplot.py
def main():
msg = textwrap.dedent(
"""
Create a single variable regression plot of specified order.
-----------------------------------------------------------------------
Examples:
* Fit a line to synthetic data with boostrap errors.
p.linspace 0 10 20 \\
| p.df 'df["y_true"] = .2 * df.x' \\
'df["noise"] = np.random.randn(20)' \\
'df["y"] = df.y_true + df.noise' --names x \\
| p.regplot -x x -y y
* Fit a quadratic to synthetic data with boostrap errors.
p.linspace 0 10 40 \\
| p.df 'df["y_true"] = .5 * df.x + .3 * df.x ** 2'\\
'df["noise"] = np.random.randn(40)' \\
'df["y"] = df.y_true + df.noise' --names x \\
| p.regplot -x x -y y --order 2
* Fit sealevel data with no bootstrap
p.example_data -d sealevel\\
| p.regplot -x year -y sealevel_mm --n_boot 1
-----------------------------------------------------------------------
"""
)
# read command line arguments
parser = argparse.ArgumentParser(
formatter_class=argparse.RawDescriptionHelpFormatter, description=msg)
arg_lib.add_args(parser, 'io_in', 'io_out', 'decorating')
msg = 'Column for dependent variable'
parser.add_argument('-x', nargs=1, type=str, dest='x', metavar='col',
help=msg, required=True)
msg = 'Column for independent variable'
parser.add_argument('-y', nargs=1, type=str, dest='y',
metavar='col', help=msg, required=True)
msg = 'The order of the polynomial to fit (default = 1)'
parser.add_argument('--order', help=msg, nargs=1, default=[1], type=int)
msg = 'Number of bootstrap samples for uncertainty region (default=1000)'
parser.add_argument(
'--n_boot', help=msg, nargs=1, default=[1000], type=int)
parser.add_argument('-a', '--alpha', help='Set opacity',
nargs=1, default=[0.5], type=float)
# parse arguments
args = parser.parse_args()
# get the input dataframe
df = io_lib.df_from_input(args)
# extract command line params
x = df[args.x[0]].values
y = df[args.y[0]].values
# do a polyfit with the specified order
coeffs = np.polyfit(x, y, args.order[0])
label = make_label(coeffs, args.savefig)
sns.regplot(
x, y, order=args.order[0], n_boot=args.n_boot[0],
line_kws={'label': label, 'color': CC[2], 'alpha': .5},
scatter_kws={'alpha': args.alpha[0], 'color': CC[0]})
pl.legend(loc='best')
pl.xlabel(args.x[0])
pl.ylabel(args.y[0])
plot_lib.refine_plot(args)
plot_lib.show(args)
0
Example 35
Project: odl Source File: operator.py
def _adjoint_definition(self):
"""Verify ``<Ax, y> == <x, A^* y>``."""
left_inner_vals = []
right_inner_vals = []
with FailCounter(
test_name='Verifying the identity <Ax, y> = <x, A^T y>',
err_msg='error = |<Ax, y< - <x, A^* y>| / ||A|| ||x|| ||y||',
logger=self.log) as counter:
for [name_x, x], [name_y, y] in samples(self.operator.domain,
self.operator.range):
x_norm = x.norm()
y_norm = y.norm()
l_inner = self.operator(x).inner(y)
r_inner = x.inner(self.operator.adjoint(y))
denom = self.operator_norm * x_norm * y_norm
error = 0 if denom == 0 else abs(l_inner - r_inner) / denom
if error > self.tol:
counter.fail('x={:25s} y={:25s} : error={:6.5f}'
''.format(name_x, name_y, error))
left_inner_vals.append(l_inner)
right_inner_vals.append(r_inner)
scale = np.polyfit(left_inner_vals, right_inner_vals, 1)[0]
self.log('\nThe adjoint seems to be scaled according to:')
self.log('(x, A^T y) / (Ax, y) = {}. Should be 1.0'.format(scale))
0
Example 36
Project: FreeCAD_assembly2 Source File: lineSearches.py
def quadraticLineSearch( f, x1, f1, intialStep, it, debugPrintLevel, printF, tol_stag=3, tol_x=10**-6):
if norm(intialStep) == 0:
printF(' quadraticLineSearch: norm search direction is 0, aborting!')
return x1
def LSEval(lam, fv=None):
return LineSearchEvaluation( f, x1, intialStep, lam, fv )
Y = [ LSEval( 0.0, f1), LSEval( 1 ), LSEval( 2 )]
y_min_prev = min(Y)
count_stagnation = 0
tol_lambda = tol_x / norm(intialStep)
for k in range(it):
Y.sort()
if debugPrintLevel > 0:
printF(' quadratic line search it %i, fmin %1.2e, lam %1.6f %1.6f %1.6f, f(lam) %1.2e %1.2e %1.2e'%( k+1, Y[0].fv, Y[0].lam,Y[1].lam,Y[2].lam,Y[0].fv,Y[1].fv,Y[2].fv ))
#``p[0]*x**(N-1) + p[1]*x**(N-2) + ... + p[N-2]*x + p[N-1]``
quadraticCoefs, residuals, rank, singular_values, rcond = numpy.polyfit( [y.lam for y in Y], [y.fv for y in Y], 2, full=True)
if quadraticCoefs[0] > 0 and rank == 3:
lam_c = -quadraticCoefs[1] / (2*quadraticCoefs[0]) #diff poly a*x**2 + b*x + c -> grad_poly = 2*a*x + b
lam_c = min( max( [y.lam for y in Y])*4, lam_c)
if lam_c < 0:
if debugPrintLevel > 1: printF(' quadratic line search lam_c < 0')
lam_c = 1.0 / (k + 1) ** 2
else:
if debugPrintLevel > 1: printF(' quadratic fit invalid, using interval halving instead')
lam_c = ( Y[0].lam + Y[1].lam )/2
del Y[2] # Y sorted at start of each iteration
Y.append( LSEval( lam_c ))
y_min = min(Y)
if y_min == y_min_prev:
count_stagnation = count_stagnation + 1
if count_stagnation > tol_stag:
if debugPrintLevel > 0: printF(' terminating quadratic line search as count_stagnation > tol_stag')
break
else:
y_min_prev = y_min
count_stagnation = 0
Lam = [y.lam for y in Y]
if max(Lam) - min(Lam) < tol_lambda:
if debugPrintLevel > 0: printF(' terminating quadratic max(Lam)-min(Lam) < tol_lambda (%e < %e)' % (max(Lam) - min(Lam), tol_lambda))
break
return min(Y).xv
0
Example 37
@staticmethod
def fit(position, fwhm, temperature=None, minmax=None):
if minmax and len(minmax) >= 2:
idxs = (fwhm >= minmax[0]) & (fwhm <= minmax[1])
position = position[idxs]
fwhm = fwhm[idxs]
A, B, C = N.polyfit(position, fwhm, 2)
fwhm_fit = N.polyval([A, B, C], position)
err = sqrt(sum((fwhm_fit - fwhm)**2) / len(position))
fit = FocusFit()
fit.position = position
fit.fwhm = fwhm
fit.temperature = temperature
fit.minmax = minmax
fit.A, fit.B, fit.C = A, B, C
fit.err = err
fit.fwhm_fit = fwhm_fit
return fit
0
Example 38
def test_polyfit_build(self,level=rlevel):
"""Ticket #628"""
ref = [-1.06123820e-06, 5.70886914e-04, -1.13822012e-01,
9.95368241e+00, -3.14526520e+02]
x = [90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103,
104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115,
116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 129,
130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141,
146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157,
158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169,
170, 171, 172, 173, 174, 175, 176]
y = [9.0, 3.0, 7.0, 4.0, 4.0, 8.0, 6.0, 11.0, 9.0, 8.0, 11.0, 5.0,
6.0, 5.0, 9.0, 8.0, 6.0, 10.0, 6.0, 10.0, 7.0, 6.0, 6.0, 6.0,
13.0, 4.0, 9.0, 11.0, 4.0, 5.0, 8.0, 5.0, 7.0, 7.0, 6.0, 12.0,
7.0, 7.0, 9.0, 4.0, 12.0, 6.0, 6.0, 4.0, 3.0, 9.0, 8.0, 8.0,
6.0, 7.0, 9.0, 10.0, 6.0, 8.0, 4.0, 7.0, 7.0, 10.0, 8.0, 8.0,
6.0, 3.0, 8.0, 4.0, 5.0, 7.0, 8.0, 6.0, 6.0, 4.0, 12.0, 9.0,
8.0, 8.0, 8.0, 6.0, 7.0, 4.0, 4.0, 5.0, 7.0]
tested = np.polyfit(x, y, 4)
assert_array_almost_equal(ref, tested)
0
Example 39
def guess_starting_values(self, y, x):
sval, oval = np.polyfit(x, y, 1)
self.set_initval(self.params.offset, oval)
self.set_initval(self.params.slope, sval)
0
Example 40
def estimate_parameters(self, signal, x1, x2, only_current=False):
"""Estimate the parameters by the two area method
Parameters
----------
signal : Signal1D instance
x1 : float
Defines the left limit of the spectral range to use for the
estimation.
x2 : float
Defines the right limit of the spectral range to use for the
estimation.
only_current : bool
If False estimates the parameters for the full dataset.
Returns
-------
bool
"""
super(Polynomial, self)._estimate_parameters(signal)
axis = signal.axes_manager.signal_axes[0]
binned = signal.metadata.Signal.binned
i1, i2 = axis.value_range_to_indices(x1, x2)
if only_current is True:
estimation = np.polyfit(axis.axis[i1:i2],
signal()[i1:i2],
self.get_polynomial_order())
if binned is True:
self.coefficients.value = estimation / axis.scale
else:
self.coefficients.value = estimation
return True
else:
if self.coefficients.map is None:
self._create_arrays()
nav_shape = signal.axes_manager._navigation_shape_in_array
with signal.unfolded():
dc = signal.data
# For polyfit the spectrum goes in the first axis
if axis.index_in_array > 0:
dc = dc.T # Unfolded, so simply transpose
cmaps = np.polyfit(axis.axis[i1:i2], dc[i1:i2, :],
self.get_polynomial_order())
if axis.index_in_array > 0:
cmaps = cmaps.T # Transpose back if needed
# Shape needed to fit coefficients.map:
cmap_shape = nav_shape + (self.get_polynomial_order() + 1, )
self.coefficients.map['values'][:] = cmaps.reshape(cmap_shape)
if binned is True:
self.coefficients.map["values"] /= axis.scale
self.coefficients.map['is_set'][:] = True
self.fetch_stored_values()
return True
0
Example 41
Project: emailinsight Source File: kerasPlottingUtils.py
def plot_regression_line(self, plot_ax, data, batches, future_batches,
n_backwards, linestyle, label):
"""Calculates and plots a regression line.
Args:
plot_ax: The ax on which to plot the line.
data: The data used to perform the regression
(e.g. training loss values).
batchs: List of all batchs (0, 1, 2, ...).
future_batchs: List of the future batchs for which values are
ought to be predicted.
n_backwards: How far back to go in time (in batchs) in order
to compute the regression. (E.g. 10 = calculate it on the
last 10 values max.)
linestyle: Linestyle of the regression line.
label: Label of the regression line.
"""
# dont try to draw anything if the data list is empty or it's the
# first batch
if len(data) > 1:
poly = np.poly1d(np.polyfit(batches[-n_backwards:],
data[-n_backwards:], self.poly_degree))
future_values = [poly(i) for i in future_batches]
plot_ax.plot(future_batches, future_values, linestyle, label=label)
0
Example 42
def guess_starting_values(self, y, x):
qval, sval, oval = np.polyfit(x, y, 2)
self.set_initval(self.params.offset, oval)
self.set_initval(self.params.slope, sval)
self.set_initval(self.params.quad, qval)
0
Example 43
Project: landlab Source File: fault_facet_finder.py
def analyse_fault_trace(self, fault_trace_node_ids):
"""
This method takes the grid and an array listing the (contiguous) node
ids of cells that contain a single fault segment trace of interest.
It sets and returns the azimuth of the fault trace, az,
-pi/2 < az <= pi/2.
(i.e., smallest angle between north and the trace).
"""
self.ft_trace_node_ids = fault_trace_node_ids
x = self.grid.node_x[fault_trace_node_ids]
y = self.grid.node_y[fault_trace_node_ids]
(grad, offset) = np.polyfit(x, y, 1)
angle = np.arctan(grad)
if grad >= 0.:
self.az = np.pi / 2. - angle
else:
self.az = -np.pi / 2. - angle
return self.az
0
Example 44
Project: odl Source File: operator.py
def self_adjoint(self):
"""Verify ``<Ax, y> == <x, Ay>``."""
left_inner_vals = []
right_inner_vals = []
with FailCounter(
test_name='Verifying the identity <Ax, y> = <x, Ay>',
err_msg='error = |<Ax, y> - <x, Ay>| / ||A|| ||x|| ||y||',
logger=self.log) as counter:
for [name_x, x], [name_y, y] in samples(self.operator.domain,
self.operator.range):
x_norm = x.norm()
y_norm = y.norm()
l_inner = self.operator(x).inner(y)
r_inner = x.inner(self.operator(y))
denom = self.operator_norm * x_norm * y_norm
error = 0 if denom == 0 else abs(l_inner - r_inner) / denom
if error > self.tol:
counter.fail('x={:25s} y={:25s} : error={:6.5f}'
''.format(name_x, name_y, error))
left_inner_vals.append(l_inner)
right_inner_vals.append(r_inner)
scale = np.polyfit(left_inner_vals, right_inner_vals, 1)[0]
self.log('\nThe adjoint seems to be scaled according to:')
self.log('(x, Ay) / (Ax, y) = {}. Should be 1.0'.format(scale))
0
Example 45
def bdrate(metric_set1, metric_set2):
"""
BJONTEGAARD Bjontegaard metric calculation
Bjontegaard's metric allows to compute the average % saving in bitrate
between two rate-distortion curves [1].
rate1,psnr1 - RD points for curve 1
rate2,psnr2 - RD points for curve 2
adapted from code from: (c) 2010 Giuseppe Valenzise
"""
# numpy plays games with its exported functions.
# pylint: disable=no-member
# pylint: disable=too-many-locals
# pylint: disable=bad-builtin
rate1 = [x[0] for x in metric_set1]
psnr1 = [x[1] for x in metric_set1]
rate2 = [x[0] for x in metric_set2]
psnr2 = [x[1] for x in metric_set2]
log_rate1 = map(math.log, rate1)
log_rate2 = map(math.log, rate2)
# Best cubic poly fit for graph represented by log_ratex, psrn_x.
poly1 = numpy.polyfit(psnr1, log_rate1, 3)
poly2 = numpy.polyfit(psnr2, log_rate2, 3)
# Integration interval.
min_int = max([min(psnr1), min(psnr2)])
max_int = min([max(psnr1), max(psnr2)])
# find integral
p_int1 = numpy.polyint(poly1)
p_int2 = numpy.polyint(poly2)
# Calculate the integrated value over the interval we care about.
int1 = numpy.polyval(p_int1, max_int) - numpy.polyval(p_int1, min_int)
int2 = numpy.polyval(p_int2, max_int) - numpy.polyval(p_int2, min_int)
# Calculate the average improvement.
avg_exp_diff = (int2 - int1) / (max_int - min_int)
# In really bad formed data the exponent can grow too large.
# clamp it.
if avg_exp_diff > 200:
avg_exp_diff = 200
# Convert to a percentage.
avg_diff = (math.exp(avg_exp_diff) - 1) * 100
return avg_diff
0
Example 46
Project: histwords Source File: plothelper.py
def trendline(xd, yd, order=1, c='r', alpha=1, plot_r=False, text_pos=None):
"""Make a line of best fit"""
#Calculate trendline
coeffs = np.polyfit(xd, yd, order)
intercept = coeffs[-1]
slope = coeffs[-2]
if order == 2: power = coeffs[0]
else: power = 0
minxd = np.min(xd)
maxxd = np.max(xd)
xl = np.array([minxd, maxxd])
yl = power * xl ** 2 + slope * xl + intercept
#Plot trendline
plt.plot(xl, yl, color=c, alpha=alpha)
#Calculate R Squared
r = sp.stats.pearsonr(xd, yd)[0]
if plot_r == False:
#Plot R^2 value
if text_pos == None:
text_pos = (0.9 * maxxd + 0.1 * minxd, 0.9 * np.max(yd) + 0.1 * np.min(yd),)
plt.text(text_pos[0], text_pos[1], '$R = %0.2f$' % r)
else:
#Return the R^2 value:
return r
0
Example 47
Project: homer Source File: fingerprint.py
def fingerprint_page(page, grid_size=16):
"Extract scale-invariant fingerprint from the systems on the page"
# Find the average angle of each system to undo rotation.
# This may be a useful step to undo rotation on the actual image.
# Also, fit a quadratic for deskewing?
angles = []
for sys in page.systems:
xs = sys['barlines'][:, :, 0].mean(1)
ys = sys['barlines'][:, :, 1].mean(1)
coefs = np.polyfit(xs, ys, 1)
theta = np.arctan(coefs[0])
angles.append(theta)
theta = np.mean(angles)
# Rotate all x values by -theta
rotation = np.array([[np.cos(-theta), -np.sin(-theta)],
[np.sin(-theta), np.cos(-theta)]])
systems_xs = []
for sys in page.systems:
def barline_rotated_x(barline):
center = barline.mean(0)
rotated_center = rotation.dot(center[:, None])
return rotated_center[0, 0]
systems_xs.append(map(barline_rotated_x, sys['barlines']))
# Make margin-invariant
xmin = min(min(xs) for xs in systems_xs)
xmax = max(max(xs) for xs in systems_xs)
width = xmax - xmin
systems_xs = [[x - xmin for x in xs] for xs in systems_xs]
systems_str = StringIO()
for xs in systems_xs:
grid_pos = ((np.asarray(xs) * grid_size) / width).astype(int)
grid_pos[grid_pos >= grid_size] = grid_size - 1
systems_str.write(','.join(map(str, grid_pos)))
systems_str.write('\n')
systems_str.seek(0)
systems_str = systems_str.read()
sha = hashlib.sha1(systems_str).digest()
fingerprint, = struct.unpack("<L", sha[:4])
return fingerprint
0
Example 48
Project: memopol2 Source File: views.py
def mep(request, mep_id):
mep_ = MEP.get(mep_id)
mep_['achievements']=autoTrophies(mep_)
positions = Position.objects.filter(mep_id=mep_id)
score_list = mep_.scores
for score in score_list:
score['color'] = score_to_color(int(score['value']))
score_list.sort(key = lambda k : datetime.strptime(k['date'], "%d/%m/%Y"))
scores = [s['value'] for s in mep_.scores]
if score_list:
try:
import numpy
import matplotlib
matplotlib.use("Agg")
from matplotlib import pyplot
pyplot.plot(scores, 'bo')
a, b = numpy.polyfit(range(len(scores)), [int(x) for x in scores], 1)
pyplot.plot([a*int(x) + b for x in range(len(scores))])
pyplot.legend(('Scores', 'Mediane'), 'best', shadow=True)
pyplot.plot(scores)
pyplot.axis([0, len(scores) - 1, 0, 102])
pyplot.title("%s - Votes notes evolution over time" % (mep_.infos['name']['full']))
pyplot.xticks(range(len(scores)), [k['date'] for k in score_list])
pyplot.xlabel("Votes dates")
pyplot.ylabel("Scores on votes")
pyplot.savefig(realpath("./memopol2/%simg/trends/meps/%s-scores.png" % (settings.MEDIA_URL, mep_id)), format="png")
pyplot.clf()
except ImportError:
pass
context = {
'mep_id': mep_id,
'mep': mep_,
'positions': positions,
'visible_count': len([x for x in positions if x.visible]),
'average': sum(scores)/len(scores) if len(scores) > 0 else "",
'score_list' : score_list,
}
return direct_to_template(request, 'meps/mep.html', context)
0
Example 49
def fit_background(q,I):
## Working on background calculation
## mkak 2016.09.28
x = q
y = I
pfit = np.polyfit(x,y,4)
yfit = np.polyval(pfit,x)
#panel.plot(xrd_spectra[0], xrd_spectra[1]-yfit, label='no bkg')
#panel.plot(xrd_spectra[0], yfit, color='blue', label='bkg')
### calculation works, but plotting here wipes previous plots - only shows last
import matplotlib as plt
plt.figure()
plt.plot(x,y,label='raw data')
plt.plot(x,yfit,label='background')
plt.plot(x,y-yfit,label='background subtracted')
plt.legend()
plt.show()
0
Example 50
Project: landlab Source File: fault_facet_finder.py
def fit_slopes_to_facet_lines(self, polynomial_degree=4, curvature_threshold=0.0004):
"""
Fits (linear) lines of best fit to extracted profiles, already stored as
class properties.
"""
avg_slopes_linear = []
avg_slopes_poly = []
curv_of_flattest_part_list = []
slope_min_curv = []
rsqd_list = []
big_slope_small_curv = []
elev_at_bssc = []
for i in range(len(self.profile_x_facet_pts)):
x = self.profile_x_facet_pts[i]
z = self.profile_z_facet_pts[i]
(grad, offset) = np.polyfit(x, z, 1)
coeffs, residuals = np.polyfit(
x, z, polynomial_degree, full=True)[:2]
rsqd = 1. - residuals / (z.size * z.var())
# differentiate the coeffs to get slope:
diff_multiplier = np.arange(polynomial_degree + 1)[::-1]
curv_multiplier = np.arange(polynomial_degree)[::-1]
z_equ = np.poly1d(coeffs)
S_equ = np.poly1d((coeffs * diff_multiplier)[:-1])
curv_equ = np.poly1d(
((coeffs * diff_multiplier)[:-1] * curv_multiplier)[:-1])
S_at_each_pt = S_equ(x)
curv_at_each_pt = curv_equ(x)
avg_slopes_linear.append(abs(grad))
avg_slopes_poly.append(np.amax(np.fabs(S_at_each_pt)))
loc_of_flattest_part = np.argmin(
np.fabs(curv_at_each_pt[2:-2])) + 2
curv_of_flattest_part = curv_at_each_pt[loc_of_flattest_part]
S_at_min_curve_untested = abs(S_at_each_pt[loc_of_flattest_part])
small_curves = np.less(
np.fabs(curv_at_each_pt[2:-2]), curvature_threshold)
try:
big_slope_small_curv.append(
np.amax(S_at_each_pt[small_curves]))
elev_at_bssc.append(z[np.argmax(S_at_each_pt[small_curves])])
except ValueError:
big_slope_small_curv.append(np.nan)
elev_at_bssc.append(np.nan)
slope_min_curv.append(S_at_min_curve_untested)
curv_of_flattest_part_list.append(curv_of_flattest_part)
rsqd_list.append(rsqd)
# figure(8)
#synthetic_z = grad*x + offset
synthetic_z = z_equ(x)
plot(x, z, 'x')
plot(x, synthetic_z, '-')
self.avg_slopes_linear = np.array(avg_slopes_linear)
self.avg_slopes_poly = np.array(avg_slopes_poly)
self.curv_of_flattest_part = np.array(curv_of_flattest_part_list)
self.slope_min_curv = np.array(slope_min_curv)
self.big_slope_small_curv = np.array(big_slope_small_curv)
self.elev_at_bssc = np.array(elev_at_bssc)
self.rsqd = np.array(rsqd_list)