[yt-svn] commit/yt: 3 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed Jul 23 04:47:28 PDT 2014
3 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/05e219be80cf/
Changeset: 05e219be80cf
Branch: stable
User: MatthewTurk
Date: 2014-07-22 23:28:48
Summary: Merging from the development branch into stable.
Affected #: 83 files
diff -r c7bae05dd4ef1d7870e5d3e042f5f6ef2cb5c9bf -r 05e219be80cf3afd394c77db58ea8864d06ae410 yt/analysis_modules/absorption_spectrum/absorption_spectrum_fit.py
--- a/yt/analysis_modules/absorption_spectrum/absorption_spectrum_fit.py
+++ b/yt/analysis_modules/absorption_spectrum/absorption_spectrum_fit.py
@@ -4,10 +4,9 @@
from yt.analysis_modules.absorption_spectrum.absorption_line \
import voigt
-
def generate_total_fit(x, fluxData, orderFits, speciesDicts,
- minError=1E-5, complexLim=.999,
- fitLim=.99, minLength=3,
+ minError=1E-4, complexLim=.995,
+ fitLim=.97, minLength=3,
maxLength=1000, splitLim=.99,
output_file=None):
@@ -90,6 +89,7 @@
fluxData[0]=1
fluxData[-1]=1
+
#Find all regions where lines/groups of lines are present
cBounds = _find_complexes(x, fluxData, fitLim=fitLim,
complexLim=complexLim, minLength=minLength,
@@ -111,6 +111,7 @@
yDatBounded=fluxData[b[1]:b[2]]
yFitBounded=yFit[b[1]:b[2]]
+
#Find init redshift
z=(xBounded[yDatBounded.argmin()]-initWl)/initWl
@@ -121,24 +122,33 @@
#Fit Using complex tools
newLinesP,flag=_complex_fit(xBounded,yDatBounded,yFitBounded,
- z,fitLim,minError*(b[2]-b[1]),speciesDict)
+ z,fitLim,minError,speciesDict)
+
+ #If flagged as a bad fit, species is lyman alpha,
+ # and it may be a saturated line, use special tools
+ if flag and species=='lya' and min(yDatBounded)<.1:
+ newLinesP=_large_flag_fit(xBounded,yDatBounded,
+ yFitBounded,z,speciesDict,
+ minSize,minError)
+
+ if na.size(newLinesP)> 0:
+
+ #Check for EXPLOOOOSIIONNNSSS
+ newLinesP = _check_numerical_instability(x, newLinesP, speciesDict,b)
+
#Check existence of partner lines if applicable
if len(speciesDict['wavelength']) != 1:
newLinesP = _remove_unaccepted_partners(newLinesP, x, fluxData,
- b, minError*(b[2]-b[1]),
- x0, xRes, speciesDict)
+ b, minError, x0, xRes, speciesDict)
- #If flagged as a bad fit, species is lyman alpha,
- # and it may be a saturated line, use special tools
- if flag and species=='lya' and min(yDatBounded)<.1:
- newLinesP=_large_flag_fit(xBounded,yDatBounded,
- yFitBounded,z,speciesDict,
- minSize,minError*(b[2]-b[1]))
+
+
#Adjust total current fit
yFit=yFit*_gen_flux_lines(x,newLinesP,speciesDict)
+
#Add new group to all fitted lines
if na.size(newLinesP)>0:
speciesLines['N']=na.append(speciesLines['N'],newLinesP[:,0])
@@ -149,6 +159,7 @@
allSpeciesLines[species]=speciesLines
+
if output_file:
_output_fit(allSpeciesLines, output_file)
@@ -205,10 +216,12 @@
#Setup initial line guesses
if initP==None: #Regular fit
initP = [0,0,0]
- if min(yDat)<.5: #Large lines get larger initial guess
- initP[0] = 10**16
+ if min(yDat)<.01: #Large lines get larger initial guess
+ initP[0] = speciesDict['init_N']*10**2
+ elif min(yDat)<.5:
+ initP[0] = speciesDict['init_N']*10**1
elif min(yDat)>.9: #Small lines get smaller initial guess
- initP[0] = 10**12.5
+ initP[0] = speciesDict['init_N']*10**-1
else:
initP[0] = speciesDict['init_N']
initP[1] = speciesDict['init_b']
@@ -225,9 +238,16 @@
return [],False
#Values to proceed through first run
- errSq,prevErrSq=1,1000
+ errSq,prevErrSq,prevLinesP=1,10*len(x),[]
+ if errBound == None:
+ errBound = len(yDat)*(max(1-yDat)*1E-2)**2
+ else:
+ errBound = errBound*len(yDat)
+
+ flag = False
while True:
+
#Initial parameter guess from joining parameters from all lines
# in lines into a single array
initP = linesP.flatten()
@@ -237,6 +257,7 @@
args=(x,yDat,yFit,speciesDict),
epsfcn=1E-10,maxfev=1000)
+
#Set results of optimization
linesP = na.reshape(fitP,(-1,3))
@@ -247,17 +268,23 @@
#Sum to get idea of goodness of fit
errSq=sum(dif**2)
+ if any(linesP[:,1]==speciesDict['init_b']):
+ # linesP = prevLinesP
+
+ flag = True
+ break
+
#If good enough, break
- if errSq < errBound:
+ if errSq < errBound:
break
#If last fit was worse, reject the last line and revert to last fit
- if errSq > prevErrSq*10:
+ if errSq > prevErrSq*10 :
#If its still pretty damn bad, cut losses and try flag fit tools
if prevErrSq >1E2*errBound and speciesDict['name']=='HI lya':
return [],True
else:
- yNewFit=_gen_flux_lines(x,prevLinesP,speciesDict)
+ linesP = prevLinesP
break
#If too many lines
@@ -266,21 +293,26 @@
if errSq >1E2*errBound and speciesDict['name']=='HI lya':
return [],True
else:
- break
+ flag = True
+ break
#Store previous data in case reject next fit
prevErrSq = errSq
prevLinesP = linesP
-
#Set up initial condition for new line
newP = [0,0,0]
- if min(dif)<.1:
- newP[0]=10**12
- elif min(dif)>.9:
- newP[0]=10**16
+
+ yAdjusted = 1+yFit*yNewFit-yDat
+
+ if min(yAdjusted)<.01: #Large lines get larger initial guess
+ newP[0] = speciesDict['init_N']*10**2
+ elif min(yAdjusted)<.5:
+ newP[0] = speciesDict['init_N']*10**1
+ elif min(yAdjusted)>.9: #Small lines get smaller initial guess
+ newP[0] = speciesDict['init_N']*10**-1
else:
- newP[0]=10**14
+ newP[0] = speciesDict['init_N']
newP[1] = speciesDict['init_b']
newP[2]=(x[dif.argmax()]-wl0)/wl0
linesP=na.append(linesP,[newP],axis=0)
@@ -290,12 +322,12 @@
# acceptable range, as given in dict ref
remove=[]
for i,p in enumerate(linesP):
- check=_check_params(na.array([p]),speciesDict)
+ check=_check_params(na.array([p]),speciesDict,x)
if check:
remove.append(i)
linesP = na.delete(linesP,remove,axis=0)
- return linesP,False
+ return linesP,flag
def _large_flag_fit(x, yDat, yFit, initz, speciesDict, minSize, errBound):
"""
@@ -489,6 +521,9 @@
#List of lines to remove
removeLines=[]
+ #Set error
+
+
#Iterate through all sets of line parameters
for i,p in enumerate(linesP):
@@ -501,16 +536,23 @@
lb = _get_bounds(p[2],b,wl,x0,xRes)
xb,yb=x[lb[0]:lb[1]],y[lb[0]:lb[1]]
+ if errBound == None:
+ errBound = 10*len(yb)*(max(1-yb)*1E-2)**2
+ else:
+ errBound = 10*errBound*len(yb)
+
#Generate a fit and find the difference to data
yFitb=_gen_flux_lines(xb,na.array([p]),speciesDict)
dif =yb-yFitb
+
+
#Only counts as an error if line is too big ---------------<
dif = [k for k in dif if k>0]
err = sum(dif)
#If the fit is too bad then add the line to list of removed lines
- if err > errBound*1E2:
+ if err > errBound:
removeLines.append(i)
break
@@ -640,21 +682,13 @@
#Check if the region needs to be divided
if b[2]-b[1]>maxLength:
- #Find the minimum absorption in the middle two quartiles of
- # the large complex
- q=(b[2]-b[1])/4
- cut = yDat[b[1]+q:b[2]-q].argmax()+b[1]+q
+ split = _split_region(yDat,b,splitLim)
- #Only break it up if the minimum absorption is actually low enough
- if yDat[cut]>splitLim:
-
- #Get the new two peaks
- b1Peak = yDat[b[1]:cut].argmin()+b[1]
- b2Peak = yDat[cut:b[2]].argmin()+cut
+ if split:
#add the two regions separately
- cBounds.insert(i+1,[b1Peak,b[1],cut])
- cBounds.insert(i+2,[b2Peak,cut,b[2]])
+ cBounds.insert(i+1,split[0])
+ cBounds.insert(i+2,split[1])
#Remove the original region
cBounds.pop(i)
@@ -663,7 +697,33 @@
return cBounds
-def _gen_flux_lines(x, linesP, speciesDict):
+
+def _split_region(yDat,b,splitLim):
+ #Find the minimum absorption in the middle two quartiles of
+ # the large complex
+
+ q=(b[2]-b[1])/4
+ cut = yDat[b[1]+q:b[2]-q].argmax()+b[1]+q
+
+ #Only break it up if the minimum absorption is actually low enough
+ if yDat[cut]>splitLim:
+
+ #Get the new two peaks
+ b1Peak = yDat[b[1]:cut].argmin()+b[1]
+ b2Peak = yDat[cut:b[2]].argmin()+cut
+
+ region_1 = [b1Peak,b[1],cut]
+ region_2 = [b2Peak,cut,b[2]]
+
+ return [region_1,region_2]
+
+ else:
+
+ return []
+
+
+
+def _gen_flux_lines(x, linesP, speciesDict,firstLine=False):
"""
Calculates the normalized flux for a region of wavelength space
generated by a set of absorption lines.
@@ -692,6 +752,9 @@
g=speciesDict['Gamma'][i]
wl=speciesDict['wavelength'][i]
y = y+ _gen_tau(x,p,f,g,wl)
+ if firstLine:
+ break
+
flux = na.exp(-y)
return flux
@@ -744,21 +807,25 @@
the difference between the fit generated by the parameters
given in pTotal multiplied by the previous fit and the desired
flux profile, w/ first index modified appropriately for bad
- parameter choices
+ parameter choices and additional penalty for fitting with a lower
+ flux than observed.
"""
pTotal.shape = (-1,3)
yNewFit = _gen_flux_lines(x,pTotal,speciesDict)
error = yDat-yFit*yNewFit
- error[0] = _check_params(pTotal,speciesDict)
+ error_plus = (yDat-yFit*yNewFit).clip(min=0)
+
+ error = error+error_plus
+ error[0] = _check_params(pTotal,speciesDict,x)
return error
-def _check_params(p, speciesDict):
+def _check_params(p, speciesDict,xb):
"""
Check to see if any of the parameters in p fall outside the range
- given in speciesDict.
+ given in speciesDict or on the boundaries
Parameters
----------
@@ -767,6 +834,8 @@
speciesDict : dictionary
dictionary with properties giving the max and min
values appropriate for each parameter N,b, and z.
+ xb : (N) ndarray
+ wavelength array [nm]
Returns
-------
@@ -774,16 +843,137 @@
0 if all values are fine
999 if any values fall outside acceptable range
"""
+
+ minz = (xb[0])/speciesDict['wavelength'][0]-1
+ maxz = (xb[-1])/speciesDict['wavelength'][0]-1
+
check = 0
- if any(p[:,0] > speciesDict['maxN']) or\
- any(p[:,0] < speciesDict['minN']) or\
- any(p[:,1] > speciesDict['maxb']) or\
- any(p[:,1] < speciesDict['minb']) or\
- any(p[:,2] > speciesDict['maxz']) or\
- any(p[:,2] < speciesDict['minz']):
+ if any(p[:,0] >= speciesDict['maxN']) or\
+ any(p[:,0] <= speciesDict['minN']) or\
+ any(p[:,1] >= speciesDict['maxb']) or\
+ any(p[:,1] <= speciesDict['minb']) or\
+ any(p[:,2] >= maxz) or\
+ any(p[:,2] <= minz):
check = 999
+
return check
+def _check_optimization_init(p,speciesDict,initz,xb,yDat,yFit,minSize,errorBound):
+
+ """
+ Check to see if any of the parameters in p are the
+ same as initial paramters and if so, attempt to
+ split the region and refit it.
+
+ Parameters
+ ----------
+ p : (3,) ndarray
+ array with form [[N1, b1, z1], ...]
+ speciesDict : dictionary
+ dictionary with properties giving the max and min
+ values appropriate for each parameter N,b, and z.
+ x : (N) ndarray
+ wavelength array [nm]
+ """
+
+ # Check if anything is a default parameter
+ if any(p[:,0] == speciesDict['init_N']) or\
+ any(p[:,0] == speciesDict['init_N']*10) or\
+ any(p[:,0] == speciesDict['init_N']*100) or\
+ any(p[:,0] == speciesDict['init_N']*.1) or\
+ any(p[:,1] == speciesDict['init_b']) or\
+ any(p[:,1] == speciesDict['maxb']):
+
+ # These are the initial bounds
+ init_bounds = [yDat.argmin(),0,len(xb)-1]
+
+ # Gratitutous limit for splitting region
+ newSplitLim = 1 - (1-min(yDat))*.5
+
+ # Attempt to split region
+ split = _split_region(yDat,init_bounds,newSplitLim)
+
+ # If we can't split it, just reject it. Its unphysical
+ # to just keep the default parameters and we're out of
+ # options at this point
+ if not split:
+ return []
+
+ # Else set up the bounds for each region and fit separately
+ b1,b2 = split[0][2], split[1][1]
+
+ p1,flag = _complex_fit(xb[:b1], yDat[:b1], yFit[:b1],
+ initz, minSize, errorBound, speciesDict)
+
+ p2,flag = _complex_fit(xb[b2:], yDat[b2:], yFit[b2:],
+ initz, minSize, errorBound, speciesDict)
+
+ # Make the final line parameters. Its annoying because
+ # one or both regions may have fit to nothing
+ if na.size(p1)> 0 and na.size(p2)>0:
+ p = na.r_[p1,p2]
+ elif na.size(p1) > 0:
+ p = p1
+ else:
+ p = p2
+
+ return p
+
+
+def _check_numerical_instability(x, p, speciesDict,b):
+
+ """
+ Check to see if any of the parameters in p are causing
+ unstable numerical effects outside the region of fit
+
+ Parameters
+ ----------
+ p : (3,) ndarray
+ array with form [[N1, b1, z1], ...]
+ speciesDict : dictionary
+ dictionary with properties giving the max and min
+ values appropriate for each parameter N,b, and z.
+ x : (N) ndarray
+ wavelength array [nm]
+ b : (3) list
+ list of integers indicating bounds of region fit in x
+ """
+
+ remove_lines = []
+
+
+ for i,line in enumerate(p):
+
+ # First to check if the line is at risk for instability
+ if line[1]<5 or line[0] < 1E12:
+
+
+ # get all flux that isn't part of fit plus a little wiggle room
+ # max and min to prevent boundary errors
+
+ flux = _gen_flux_lines(x,[line],speciesDict,firstLine=True)
+ flux = na.r_[flux[:max(b[1]-10,0)], flux[min(b[2]+10,len(x)):]]
+
+ #Find regions that are absorbing outside the region we fit
+ flux_dif = 1 - flux
+ absorbing_coefficient = max(abs(flux_dif))
+
+
+ #Really there shouldn't be any absorption outside
+ #the region we fit, but we'll give some leeway.
+ #for high resolution spectra the tiny bits on the edges
+ #can give a non negligible amount of flux. Plus the errors
+ #we are looking for are HUGE.
+ if absorbing_coefficient > .1:
+
+ # we just set it to no fit because we've tried
+ # everything else at this point. this region just sucks :(
+ remove_lines.append(i)
+
+ if remove_lines:
+ p = na.delete(p, remove_lines, axis=0)
+
+ return p
def _output_fit(lineDic, file_name = 'spectrum_fit.h5'):
"""
@@ -815,4 +1005,5 @@
f.create_dataset("{0}/z".format(ion),data=params['z'])
f.create_dataset("{0}/complex".format(ion),data=params['group#'])
print 'Writing spectrum fit to {0}'.format(file_name)
+ f.close()
diff -r c7bae05dd4ef1d7870e5d3e042f5f6ef2cb5c9bf -r 05e219be80cf3afd394c77db58ea8864d06ae410 yt/analysis_modules/cosmological_observation/light_ray/light_ray.py
--- a/yt/analysis_modules/cosmological_observation/light_ray/light_ray.py
+++ b/yt/analysis_modules/cosmological_observation/light_ray/light_ray.py
@@ -392,12 +392,9 @@
pf = load(my_segment['filename'])
if self.near_redshift == self.far_redshift:
- h_vel = cm_per_km * pf.units['mpc'] * \
- vector_length(my_segment['start'], my_segment['end']) * \
- self.cosmology.HubbleConstantNow * \
- self.cosmology.ExpansionFactor(my_segment['redshift'])
- next_redshift = np.sqrt((1. + h_vel / speed_of_light_cgs) /
- (1. - h_vel / speed_of_light_cgs)) - 1.
+ next_redshift = my_segment["redshift"] - \
+ self._deltaz_forward(my_segment["redshift"],
+ pf.units["mpc"] * my_segment["traversal_box_fraction"])
elif my_segment['next'] is None:
next_redshift = self.near_redshift
else:
diff -r c7bae05dd4ef1d7870e5d3e042f5f6ef2cb5c9bf -r 05e219be80cf3afd394c77db58ea8864d06ae410 yt/analysis_modules/halo_finding/halo_objects.py
--- a/yt/analysis_modules/halo_finding/halo_objects.py
+++ b/yt/analysis_modules/halo_finding/halo_objects.py
@@ -540,22 +540,23 @@
temp_e2[:,dim] = e2_vector[dim]
length = np.abs(np.sum(rr * temp_e2, axis = 1) * (1 - \
np.sum(rr * temp_e0, axis = 1)**2. * mag_A**-2. - \
- np.sum(rr * temp_e1, axis = 1)**2. * mag_B**-2)**(-0.5))
+ np.sum(rr * temp_e1, axis = 1)**2. * mag_B**-2.)**(-0.5))
length[length == np.inf] = 0.
tC_index = np.nanargmax(length)
mag_C = length[tC_index]
# tilt is calculated from the rotation about x axis
# needed to align e1 vector with the y axis
# after e0 is aligned with x axis
- # find the t1 angle needed to rotate about z axis to align e0 to x
- t1 = np.arctan(e0_vector[1] / e0_vector[0])
- RZ = get_rotation_matrix(-t1, (0, 0, 1)).transpose()
- r1 = (e0_vector * RZ).sum(axis = 1)
+ # find the t1 angle needed to rotate about z axis to align e0 onto x-z plane
+ t1 = np.arctan(-e0_vector[1] / e0_vector[0])
+ RZ = get_rotation_matrix(t1, (0, 0, 1))
+ r1 = np.dot(RZ, e0_vector)
# find the t2 angle needed to rotate about y axis to align e0 to x
- t2 = np.arctan(-r1[2] / r1[0])
- RY = get_rotation_matrix(-t2, (0, 1, 0)).transpose()
+ t2 = np.arctan(r1[2] / r1[0])
+ RY = get_rotation_matrix(t2, (0, 1, 0))
r2 = np.dot(RY, np.dot(RZ, e1_vector))
- tilt = np.arctan(r2[2]/r2[1])
+ # find the tilt angle needed to rotate about x axis to align e1 to y and e2 to z
+ tilt = np.arctan(-r2[2] / r2[1])
return (mag_A, mag_B, mag_C, e0_vector[0], e0_vector[1],
e0_vector[2], tilt)
@@ -771,13 +772,13 @@
Returns
-------
- tuple : (cm, mag_A, mag_B, mag_C, e1_vector, tilt)
+ tuple : (cm, mag_A, mag_B, mag_C, e0_vector, tilt)
The 6-tuple has in order:
#. The center of mass as an array.
#. mag_A as a float.
#. mag_B as a float.
#. mag_C as a float.
- #. e1_vector as an array.
+ #. e0_vector as an array.
#. tilt as a float.
Examples
@@ -808,7 +809,7 @@
def __init__(self, pf, id, size=None, CoM=None,
max_dens_point=None, group_total_mass=None, max_radius=None, bulk_vel=None,
rms_vel=None, fnames=None, mag_A=None, mag_B=None, mag_C=None,
- e1_vec=None, tilt=None, supp=None):
+ e0_vec=None, tilt=None, supp=None):
self.pf = pf
self.gridsize = (self.pf.domain_right_edge - \
@@ -824,7 +825,7 @@
self.mag_A = mag_A
self.mag_B = mag_B
self.mag_C = mag_C
- self.e1_vec = e1_vec
+ self.e0_vec = e0_vec
self.tilt = tilt
# locs=the names of the h5 files that have particle data for this halo
self.fnames = fnames
@@ -902,8 +903,8 @@
def _get_ellipsoid_parameters_basic_loadedhalo(self):
if self.mag_A is not None:
- return (self.mag_A, self.mag_B, self.mag_C, self.e1_vec[0],
- self.e1_vec[1], self.e1_vec[2], self.tilt)
+ return (self.mag_A, self.mag_B, self.mag_C, self.e0_vec[0],
+ self.e0_vec[1], self.e0_vec[2], self.tilt)
else:
return self._get_ellipsoid_parameters_basic()
@@ -917,13 +918,13 @@
Returns
-------
- tuple : (cm, mag_A, mag_B, mag_C, e1_vector, tilt)
+ tuple : (cm, mag_A, mag_B, mag_C, e0_vector, tilt)
The 6-tuple has in order:
#. The center of mass as an array.
#. mag_A as a float.
#. mag_B as a float.
#. mag_C as a float.
- #. e1_vector as an array.
+ #. e0_vector as an array.
#. tilt as a float.
Examples
@@ -996,7 +997,7 @@
max_dens_point=None, group_total_mass=None, max_radius=None, bulk_vel=None,
rms_vel=None, fnames=None, mag_A=None, mag_B=None, mag_C=None,
- e1_vec=None, tilt=None, supp=None):
+ e0_vec=None, tilt=None, supp=None):
self.pf = pf
self.gridsize = (self.pf.domain_right_edge - \
@@ -1012,7 +1013,7 @@
self.mag_A = mag_A
self.mag_B = mag_B
self.mag_C = mag_C
- self.e1_vec = e1_vec
+ self.e0_vec = e0_vec
self.tilt = tilt
self.bin_count = None
self.overdensity = None
@@ -1256,8 +1257,8 @@
"x","y","z", "center-of-mass",
"x","y","z",
"vx","vy","vz","max_r","rms_v",
- "mag_A", "mag_B", "mag_C", "e1_vec0",
- "e1_vec1", "e1_vec2", "tilt", "\n"]))
+ "mag_A", "mag_B", "mag_C", "e0_vec0",
+ "e0_vec1", "e0_vec2", "tilt", "\n"]))
for group in self:
f.write("%10i\t" % group.id)
@@ -1569,17 +1570,17 @@
mag_A = float(line[15])
mag_B = float(line[16])
mag_C = float(line[17])
- e1_vec0 = float(line[18])
- e1_vec1 = float(line[19])
- e1_vec2 = float(line[20])
- e1_vec = np.array([e1_vec0, e1_vec1, e1_vec2])
+ e0_vec0 = float(line[18])
+ e0_vec1 = float(line[19])
+ e0_vec2 = float(line[20])
+ e0_vec = np.array([e0_vec0, e0_vec1, e0_vec2])
tilt = float(line[21])
self._groups.append(LoadedHalo(self.pf, halo, size = size,
CoM = CoM,
max_dens_point = max_dens_point,
group_total_mass = group_total_mass, max_radius = max_radius,
bulk_vel = bulk_vel, rms_vel = rms_vel, fnames = fnames,
- mag_A = mag_A, mag_B = mag_B, mag_C = mag_C, e1_vec = e1_vec,
+ mag_A = mag_A, mag_B = mag_B, mag_C = mag_C, e0_vec = e0_vec,
tilt = tilt))
else:
mylog.error("I am unable to parse this line. Too many or too few items. %s" % orig)
diff -r c7bae05dd4ef1d7870e5d3e042f5f6ef2cb5c9bf -r 05e219be80cf3afd394c77db58ea8864d06ae410 yt/analysis_modules/star_analysis/sfr_spectrum.py
--- a/yt/analysis_modules/star_analysis/sfr_spectrum.py
+++ b/yt/analysis_modules/star_analysis/sfr_spectrum.py
@@ -41,7 +41,7 @@
star_creation_time : Ordered array or list of floats
The creation time for the stars in code units.
volume : Float
- The volume of the region for the specified list of stars.
+ The comoving volume of the region for the specified list of stars.
bins : Integer
The number of time bins used for binning the stars. Default = 300.
@@ -69,7 +69,7 @@
If data_source is not provided, all of these paramters need to be set:
star_mass (array, Msun),
star_creation_time (array, code units),
- volume (float, Mpc**3).
+ volume (float, cMpc**3).
""")
return None
self.mode = 'provided'
@@ -125,14 +125,14 @@
"""
if self.mode == 'data_source':
try:
- vol = self._data_source.volume('mpc')
+ vol = self._data_source.volume('mpccm')
except AttributeError:
# If we're here, this is probably a HOPHalo object, and we
# can get the volume this way.
ds = self._data_source.get_sphere()
- vol = ds.volume('mpc')
+ vol = ds.volume('mpccm')
elif self.mode == 'provided':
- vol = self.volume
+ vol = self.volume('mpccm')
tc = self._pf["Time"]
self.time = []
self.lookback_time = []
@@ -148,6 +148,7 @@
self.redshift.append(self.cosm.ComputeRedshiftFromTime(time * tc))
self.Msol_yr.append(self.mass_bins[i] / \
(self.time_bins_dt[i] * tc / YEAR))
+ # changed vol from mpc to mpccm used in literature
self.Msol_yr_vol.append(self.mass_bins[i] / \
(self.time_bins_dt[i] * tc / YEAR) / vol)
self.Msol.append(self.mass_bins[i])
diff -r c7bae05dd4ef1d7870e5d3e042f5f6ef2cb5c9bf -r 05e219be80cf3afd394c77db58ea8864d06ae410 yt/data_objects/data_containers.py
--- a/yt/data_objects/data_containers.py
+++ b/yt/data_objects/data_containers.py
@@ -215,8 +215,12 @@
if fields == None: fields = []
self.fields = ensure_list(fields)[:]
self.field_data = YTFieldData()
+ self._default_field_parameters = {}
+ self._default_field_parameters["center"] = np.zeros(3, dtype='float64')
+ self._default_field_parameters["bulk_velocity"] = np.zeros(3, dtype='float64')
+ self._default_field_parameters["normal"] = np.array([0,0,1], dtype='float64')
self.field_parameters = {}
- self.__set_default_field_parameters()
+ self._set_default_field_parameters()
self._cut_masks = {}
self._point_indices = {}
self._vc_data = {}
@@ -224,10 +228,13 @@
mylog.debug("Setting %s to %s", key, val)
self.set_field_parameter(key, val)
- def __set_default_field_parameters(self):
- self.set_field_parameter("center",np.zeros(3,dtype='float64'))
- self.set_field_parameter("bulk_velocity",np.zeros(3,dtype='float64'))
- self.set_field_parameter("normal",np.array([0,0,1],dtype='float64'))
+ def _set_default_field_parameters(self):
+ for k,v in self._default_field_parameters.items():
+ self.set_field_parameter(k,v)
+
+ def _is_default_field_parameter(self, parameter):
+ if parameter not in self._default_field_parameters: return False
+ return self._default_field_parameters[parameter] is self.field_parameters[parameter]
def _set_center(self, center):
if center is None:
@@ -783,7 +790,6 @@
@cache_mask
def _get_cut_mask(self, grid):
- #pdb.set_trace()
points_in_grid = np.all(self.positions > grid.LeftEdge, axis=1) & \
np.all(self.positions <= grid.RightEdge, axis=1)
pids = np.where(points_in_grid)[0]
@@ -1772,8 +1778,10 @@
self._distributed = False
self._okay_to_serialize = False
self._check_region = True
+ # Use the data_source's field parameters if they don't exist in the
+ # object or if they are the default values
for k, v in source.field_parameters.items():
- if k not in self.field_parameters:
+ if k not in self.field_parameters or self._is_default_field_parameter(k):
self.set_field_parameter(k,v)
self.source = source
if self._field_cuts is not None:
@@ -2115,7 +2123,7 @@
self._okay_to_serialize = False
self._check_region = True
for k, v in source.field_parameters.items():
- if k not in self.field_parameters:
+ if k not in self.field_parameters or self._is_default_field_parameter(k):
self.set_field_parameter(k,v)
self.source = source
if self._field_cuts is not None:
@@ -2661,14 +2669,14 @@
i = 0
for grid in self._grids:
pointI = self._get_point_indices(grid)
- np = pointI[0].ravel().size
+ npoints = pointI[0].ravel().size
if grid.has_key(field):
new_field = grid[field]
else:
new_field = np.ones(grid.ActiveDimensions, dtype=dtype) * default_val
- new_field[pointI] = self[field][i:i+np]
+ new_field[pointI] = self[field][i:i+npoints]
grid[field] = new_field
- i += np
+ i += npoints
def _is_fully_enclosed(self, grid):
return np.all(self._get_cut_mask)
@@ -3579,23 +3587,23 @@
self._tilt = tilt
# find the t1 angle needed to rotate about z axis to align e0 to x
- t1 = np.arctan(e0[1] / e0[0])
+ t1 = np.arctan(-e0[1] / e0[0])
# rotate e0 by -t1
- RZ = get_rotation_matrix(t1, (0,0,1)).transpose()
- r1 = (e0 * RZ).sum(axis = 1)
+ RZ = get_rotation_matrix(t1, (0,0,1))
+ r1 = np.dot(RZ, e0)
# find the t2 angle needed to rotate about y axis to align e0 to x
- t2 = np.arctan(-r1[2] / r1[0])
+ t2 = np.arctan(r1[2] / r1[0])
"""
calculate the original e1
given the tilt about the x axis when e0 was aligned
to x after t1, t2 rotations about z, y
"""
- RX = get_rotation_matrix(-tilt, (1, 0, 0)).transpose()
- RY = get_rotation_matrix(-t2, (0, 1, 0)).transpose()
- RZ = get_rotation_matrix(-t1, (0, 0, 1)).transpose()
- e1 = ((0, 1, 0) * RX).sum(axis=1)
- e1 = (e1 * RY).sum(axis=1)
- e1 = (e1 * RZ).sum(axis=1)
+ RX = get_rotation_matrix(-tilt, (1, 0, 0))
+ RY = get_rotation_matrix(-t2, (0, 1, 0))
+ RZ = get_rotation_matrix(-t1, (0, 0, 1))
+ e1 = np.dot(RX, (0,1,0))
+ e1 = np.dot(RY, e1)
+ e1 = np.dot(RZ, e1)
e2 = np.cross(e0, e1)
self._e1 = e1
diff -r c7bae05dd4ef1d7870e5d3e042f5f6ef2cb5c9bf -r 05e219be80cf3afd394c77db58ea8864d06ae410 yt/data_objects/profiles.py
--- a/yt/data_objects/profiles.py
+++ b/yt/data_objects/profiles.py
@@ -811,6 +811,7 @@
self.pf = data_source.pf
self.field_data = YTFieldData()
self.weight_field = weight_field
+ ParallelAnalysisInterface.__init__(self, comm=data_source.comm)
def add_fields(self, fields):
fields = ensure_list(fields)
@@ -822,7 +823,9 @@
def _finalize_storage(self, fields, temp_storage):
# We use our main comm here
# This also will fill _field_data
- # FIXME: Add parallelism and combining std stuff
+ temp_storage.values = self.comm.mpi_allreduce(temp_storage.values, op="sum", dtype="float64")
+ temp_storage.weight_values = self.comm.mpi_allreduce(temp_storage.weight_values, op="sum", dtype="float64")
+ temp_storage.used = self.comm.mpi_allreduce(temp_storage.used, op="sum", dtype="bool")
if self.weight_field is not None:
temp_storage.values /= temp_storage.weight_values[...,None]
blank = ~temp_storage.used
diff -r c7bae05dd4ef1d7870e5d3e042f5f6ef2cb5c9bf -r 05e219be80cf3afd394c77db58ea8864d06ae410 yt/data_objects/tests/test_boolean_regions.py
--- a/yt/data_objects/tests/test_boolean_regions.py
+++ b/yt/data_objects/tests/test_boolean_regions.py
@@ -246,10 +246,8 @@
for n in [1, 2, 4, 8]:
pf = fake_random_pf(64, nprocs=n)
pf.h
- ell1 = pf.h.ellipsoid([0.25]*3, 0.05, 0.05, 0.05, np.array([0.1]*3),
- np.array([0.1]*3))
- ell2 = pf.h.ellipsoid([0.75]*3, 0.05, 0.05, 0.05, np.array([0.1]*3),
- np.array([0.1]*3))
+ ell1 = pf.h.ellipsoid([0.25]*3, 0.05, 0.05, 0.05, np.array([0.1]*3), 0.1)
+ ell2 = pf.h.ellipsoid([0.75]*3, 0.05, 0.05, 0.05, np.array([0.1]*3), 0.1)
# Store the original indices
i1 = ell1['ID']
i1.sort()
@@ -287,10 +285,8 @@
for n in [1, 2, 4, 8]:
pf = fake_random_pf(64, nprocs=n)
pf.h
- ell1 = pf.h.ellipsoid([0.45]*3, 0.05, 0.05, 0.05, np.array([0.1]*3),
- np.array([0.1]*3))
- ell2 = pf.h.ellipsoid([0.55]*3, 0.05, 0.05, 0.05, np.array([0.1]*3),
- np.array([0.1]*3))
+ ell1 = pf.h.ellipsoid([0.45]*3, 0.05, 0.05, 0.05, np.array([0.1]*3), 0.1)
+ ell2 = pf.h.ellipsoid([0.55]*3, 0.05, 0.05, 0.05, np.array([0.1]*3), 0.1)
# Get indices of both.
i1 = ell1['ID']
i2 = ell2['ID']
diff -r c7bae05dd4ef1d7870e5d3e042f5f6ef2cb5c9bf -r 05e219be80cf3afd394c77db58ea8864d06ae410 yt/data_objects/tests/test_projection.py
--- a/yt/data_objects/tests/test_projection.py
+++ b/yt/data_objects/tests/test_projection.py
@@ -70,5 +70,8 @@
v1 = proj["Density"].sum()
v2 = (dd["Density"] * dd["d%s" % an]).sum()
yield assert_rel_equal, v1, v2, 10
-
-
+ # test if projections inherit the field parameters of their data sources
+ dd.set_field_parameter("bulk_velocity", np.array([0,1,2]))
+ proj = pf.h.proj(0, "Density", source=dd)
+ yield assert_equal, dd.field_parameters["bulk_velocity"], \
+ proj.field_parameters["bulk_velocity"]
diff -r c7bae05dd4ef1d7870e5d3e042f5f6ef2cb5c9bf -r 05e219be80cf3afd394c77db58ea8864d06ae410 yt/data_objects/time_series.py
--- a/yt/data_objects/time_series.py
+++ b/yt/data_objects/time_series.py
@@ -114,7 +114,12 @@
return self.get_range(key.start, key.stop)
# This will return a sliced up object!
return TimeSeriesData(self._pre_outputs[key], self.parallel)
- o = self._pre_outputs[key]
+ try:
+ o = self._pre_outputs[key]
+ except IndexError:
+ raise InvalidSimulationTimeSeries("Your TimeSeries is empty. \n" +
+ "Confirm you are running from the simulation source directory.")
+
if isinstance(o, types.StringTypes):
o = load(o,**self.kwargs)
return o
diff -r c7bae05dd4ef1d7870e5d3e042f5f6ef2cb5c9bf -r 05e219be80cf3afd394c77db58ea8864d06ae410 yt/utilities/definitions.py
--- a/yt/utilities/definitions.py
+++ b/yt/utilities/definitions.py
@@ -15,7 +15,7 @@
from .physical_constants import \
mpc_per_mpc, kpc_per_mpc, pc_per_mpc, au_per_mpc, rsun_per_mpc, \
- miles_per_mpc, km_per_mpc, cm_per_mpc, sec_per_Gyr, sec_per_Myr, \
+ miles_per_mpc, km_per_mpc, m_per_mpc, cm_per_mpc, sec_per_Gyr, sec_per_Myr, \
sec_per_year, sec_per_day
# The number of levels we expect to have at most
@@ -44,6 +44,7 @@
'rsun' : rsun_per_mpc,
'miles' : miles_per_mpc,
'km' : km_per_mpc,
+ 'm' : m_per_mpc,
'cm' : cm_per_mpc}
# Nicely formatted versions of common length units
diff -r c7bae05dd4ef1d7870e5d3e042f5f6ef2cb5c9bf -r 05e219be80cf3afd394c77db58ea8864d06ae410 yt/utilities/lib/field_interpolation_tables.pxd
--- a/yt/utilities/lib/field_interpolation_tables.pxd
+++ b/yt/utilities/lib/field_interpolation_tables.pxd
@@ -33,6 +33,7 @@
cdef extern from "math.h":
double expf(double x) nogil
+ int isnormal(double x) nogil
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -58,6 +59,7 @@
cdef np.float64_t bv, dy, dd, tf, rv
cdef int bin_id
if dvs[fit.field_id] >= fit.bounds[1] or dvs[fit.field_id] <= fit.bounds[0]: return 0.0
+ if not isnormal(dvs[fit.field_id]): return 0.0
bin_id = <int> ((dvs[fit.field_id] - fit.bounds[0]) * fit.idbin)
bin_id = iclip(bin_id, 0, fit.nbins-2)
dd = dvs[fit.field_id] - (fit.bounds[0] + bin_id * fit.dbin) # x - x0
diff -r c7bae05dd4ef1d7870e5d3e042f5f6ef2cb5c9bf -r 05e219be80cf3afd394c77db58ea8864d06ae410 yt/utilities/physical_constants.py
--- a/yt/utilities/physical_constants.py
+++ b/yt/utilities/physical_constants.py
@@ -48,6 +48,7 @@
mpc_per_rsun = 2.253962e-14
mpc_per_miles = 5.21552871e-20
mpc_per_km = 3.24077929e-20
+mpc_per_m = 3.24077929e-23
mpc_per_cm = 3.24077929e-25
kpc_per_cm = mpc_per_cm / mpc_per_kpc
km_per_pc = 3.08567758e13
@@ -63,6 +64,7 @@
rsun_per_mpc = 1.0 / mpc_per_rsun
miles_per_mpc = 1.0 / mpc_per_miles
km_per_mpc = 1.0 / mpc_per_km
+m_per_mpc = 1.0 / mpc_per_m
cm_per_mpc = 1.0 / mpc_per_cm
cm_per_kpc = 1.0 / kpc_per_cm
cm_per_km = 1.0 / km_per_cm
diff -r c7bae05dd4ef1d7870e5d3e042f5f6ef2cb5c9bf -r 05e219be80cf3afd394c77db58ea8864d06ae410 yt/visualization/eps_writer.py
--- a/yt/visualization/eps_writer.py
+++ b/yt/visualization/eps_writer.py
@@ -15,6 +15,7 @@
import pyx
import numpy as np
from matplotlib import cm
+import matplotlib.pyplot as plt
from _mpl_imports import FigureCanvasAgg
from yt.utilities.logger import ytLogger as mylog
@@ -288,8 +289,6 @@
"""
if isinstance(plot, (PlotWindow, PhasePlot)):
plot.refresh()
- else:
- plot._redraw_image()
if isinstance(plot, (VMPlot, PlotWindow)):
if isinstance(plot, PlotWindow):
data = plot._frb
@@ -344,6 +343,26 @@
_ylabel = plot[k].axes.get_ylabel()
if tickcolor == None:
_tickcolor = None
+ elif isinstance(plot, np.ndarray):
+ ax = plt.gca()
+ _xrange = ax.get_xlim()
+ _yrange = ax.get_ylim()
+ _xlog=False
+ _ylog=False
+ if bare_axes:
+ _xlabel = ""
+ _ylabel = ""
+ else:
+ if xlabel != None:
+ _xlabel = xlabel
+ else:
+ _xlabel = ax.get_xlabel()
+ if ylabel != None:
+ _ylabel = ylabel
+ else:
+ _ylabel = ax.get_ylabel()
+ if tickcolor == None:
+ _tickcolor = None
else:
_xrange = plot._axes.get_xlim()
_yrange = plot._axes.get_ylim()
@@ -461,6 +480,13 @@
# Remove colorbar
_p1 = plot._figure
_p1.delaxes(_p1.axes[1])
+ elif isinstance(plot, np.ndarray):
+ fig = plt.figure()
+ iplot = plt.figimage(plot)
+ _p1 = iplot.figure
+ _p1.set_size_inches(self.figsize[0], self.figsize[1])
+ ax = plt.gca();
+ _p1.add_axes(ax)
else:
raise RuntimeError("Unknown plot type")
@@ -855,7 +881,7 @@
#=============================================================================
- def save_fig(self, filename="test", format="eps"):
+ def save_fig(self, filename="test", format="eps", resolution=250):
r"""Saves current figure to a file.
Parameters
@@ -875,6 +901,10 @@
self.canvas.writeEPSfile(filename)
elif format == "pdf":
self.canvas.writePDFfile(filename)
+ elif format == "png":
+ self.canvas.writeGSfile(filename+".png", "png16m", resolution=resolution)
+ elif format == "jpg":
+ self.canvas.writeGSfile(filename+".jpeg", "jpeg", resolution=resolution)
else:
raise RuntimeError("format %s unknown." % (format))
diff -r c7bae05dd4ef1d7870e5d3e042f5f6ef2cb5c9bf -r 05e219be80cf3afd394c77db58ea8864d06ae410 yt/visualization/plot_container.py
--- a/yt/visualization/plot_container.py
+++ b/yt/visualization/plot_container.py
@@ -7,6 +7,7 @@
from functools import wraps
from matplotlib.font_manager import FontProperties
+from ._mpl_imports import FigureCanvasAgg
from .tick_locators import LogLocator, LinearLocator
from .color_maps import yt_colormaps, is_colormap
from .plot_modifications import \
@@ -108,6 +109,9 @@
font_path = matplotlib.get_data_path() + '/fonts/ttf/STIXGeneral.ttf'
self._font_properties = FontProperties(size=fontsize, fname=font_path)
self._font_color = None
+ self._xlabel = None
+ self._ylabel = None
+ self._colorbarlabel = None
@invalidate_plot
def set_log(self, field, log):
@@ -474,3 +478,67 @@
img = base64.b64encode(self.plots[field]._repr_png_())
ret += '<img src="data:image/png;base64,%s"><br>' % img
return ret
+
+ def set_xlabel(self, x_title, fontsize=18):
+ r"""
+ Allow the user to modify the X-axis title
+ Defaults to the global value. Fontsize defaults
+ to 18.
+
+ Parameters
+ ----------
+ x_title: str
+ The new string for the x-axis. This is a required argument.
+
+ fontsize: float
+ Fontsize for the x-axis title
+
+ >>> plot.set_xtitle("H2I Number Density (cm$^{-3}$)")
+
+ """
+ for f in self.plots:
+ self.plots[f].axes.xaxis.set_label_text(x_title, fontsize=fontsize)
+ self._xlabel = x_title
+
+ def set_ylabel(self, y_title, fontsize=18):
+ r"""
+ Allow the user to modify the Y-axis title
+ Defaults to the global value. Fontsize defaults
+ to 18.
+
+ Parameters
+ ----------
+ y_title: str
+ The new string for the y-axis. This is a required argument.
+ fontsize: float
+ Fontsize for the y-axis title
+
+ >>> plot.set_ytitle("Temperature (K)")
+
+ """
+ for f in self.plots:
+ self.plots[f].axes.yaxis.set_label_text(y_title, fontsize=fontsize)
+ self._ylabel = y_title
+
+ def set_colorbar_label(self, z_title, fontsize=18):
+ r"""
+ Allow the user to modify the Z-axis title
+ Defaults to the global value. Fontsize defaults
+ to 18.
+
+ Parameters
+ ----------
+ z_title: str
+ The new string for the colorbar. This is a required argument.
+ fontsize: float
+ Fontsize for the z-axis title
+
+ >>> plot.set_ztitle("Enclosed Gas Mass ($M_{\odot}$)")
+
+ """
+ for f in self.plots:
+ self.plots[f].cax.yaxis.set_label_text(z_title, fontsize=fontsize)
+ self._colorbarlabel = z_title
+
+ def _get_axes_labels(self):
+ return(self._xlabel, self._ylabel, self._colorbarlabel)
diff -r c7bae05dd4ef1d7870e5d3e042f5f6ef2cb5c9bf -r 05e219be80cf3afd394c77db58ea8864d06ae410 yt/visualization/plot_modifications.py
--- a/yt/visualization/plot_modifications.py
+++ b/yt/visualization/plot_modifications.py
@@ -28,6 +28,8 @@
sec_per_day, sec_per_hr
from yt.visualization.image_writer import apply_colormap
+from matplotlib.colors import colorConverter
+
import _MPL
callback_registry = {}
@@ -337,20 +339,25 @@
class GridBoundaryCallback(PlotCallback):
"""
annotate_grids(alpha=0.7, min_pix=1, min_pix_ids=20, draw_ids=False, periodic=True,
- min_level=None, max_level=None, cmap='B-W LINEAR_r'):
+ min_level=None, max_level=None, cmap='B-W LINEAR_r', edgecolors=None,
+ linewidth=1.0):
Draws grids on an existing PlotWindow object.
Adds grid boundaries to a plot, optionally with alpha-blending. By default,
colors different levels of grids with different colors going from white to
- black, but you can change to any arbitrary colormap with cmap keyword
- (or all black cells for all levels with cmap=None). Cuttoff for display is at
- min_pix wide. draw_ids puts the grid id in the corner of the grid.
+ black, but you can change to any arbitrary colormap with cmap keyword, to all black
+ grid edges for all levels with cmap=None and edgecolors=None, or to an arbitrary single
+ color for grid edges with edgecolors='YourChosenColor' defined in any of the standard ways
+ (e.g., edgecolors='white', edgecolors='r', edgecolors='#00FFFF', or edgecolor='0.3', where
+ the last is a float in 0-1 scale indicating gray).
+ Note that setting edgecolors overrides cmap if you have both set to non-None values.
+ Cutoff for display is at min_pix wide. draw_ids puts the grid id in the corner of the grid.
(Not so great in projections...). One can set min and maximum level of
- grids to display.
+ grids to display, and can change the linewidth of the displayed grids.
"""
_type_name = "grids"
def __init__(self, alpha=0.7, min_pix=1, min_pix_ids=20, draw_ids=False, periodic=True,
- min_level=None, max_level=None, cmap='B-W LINEAR_r'):
+ min_level=None, max_level=None, cmap='B-W LINEAR_r', edgecolors = None, linewidth=1.0):
PlotCallback.__init__(self)
self.alpha = alpha
self.min_pix = min_pix
@@ -359,7 +366,9 @@
self.periodic = periodic
self.min_level = min_level
self.max_level = max_level
+ self.linewidth = linewidth
self.cmap = cmap
+ self.edgecolors = edgecolors
def __call__(self, plot):
x0, x1 = plot.xlim
@@ -399,13 +408,17 @@
( levels >= min_level) & \
( levels <= max_level)
- if self.cmap is not None:
- edgecolors = apply_colormap(levels[(levels <= max_level) & (levels >= min_level)]*1.0,
- color_bounds=[0,plot.data.pf.h.max_level],
- cmap_name=self.cmap)[0,:,:]*1.0/255.
- edgecolors[:,3] = self.alpha
- else:
- edgecolors = (0.0,0.0,0.0,self.alpha)
+ # Grids can either be set by edgecolors OR a colormap.
+ if self.edgecolors is not None:
+ edgecolors = colorConverter.to_rgba(self.edgecolors, alpha=self.alpha)
+ else: # use colormap if not explicity overridden by edgecolors
+ if self.cmap is not None:
+ edgecolors = apply_colormap(levels[(levels <= max_level) & (levels >= min_level)]*1.0,
+ color_bounds=[0,plot.data.pf.h.max_level],
+ cmap_name=self.cmap)[0,:,:]*1.0/255.
+ edgecolors[:,3] = self.alpha
+ else:
+ edgecolors = (0.0,0.0,0.0,self.alpha)
if visible.nonzero()[0].size == 0: continue
verts = np.array(
@@ -414,7 +427,7 @@
verts=verts.transpose()[visible,:,:]
grid_collection = matplotlib.collections.PolyCollection(
verts, facecolors="none",
- edgecolors=edgecolors)
+ edgecolors=edgecolors, linewidth=self.linewidth)
plot._axes.hold(True)
plot._axes.add_collection(grid_collection)
@@ -499,7 +512,7 @@
def get_smallest_appropriate_unit(v, pf):
max_nu = 1e30
good_u = None
- for unit in ['mpc', 'kpc', 'pc', 'au', 'rsun', 'km', 'cm']:
+ for unit in ['mpc', 'kpc', 'pc', 'au', 'rsun', 'km', 'm', 'cm']:
vv = v*pf[unit]
if vv < max_nu and vv > 1.0:
good_u = unit
@@ -1056,7 +1069,8 @@
"""
annotate_particles(width, p_size=1.0, col='k', marker='o', stride=1.0,
ptype=None, stars_only=False, dm_only=False,
- minimum_mass=None, alpha=1.0)
+ minimum_mass=None, alpha=1.0, min_star_age=None,
+ max_star_age=None)
Adds particle positions, based on a thick slab along *axis* with a
*width* along the line of sight. *p_size* controls the number of
@@ -1064,14 +1078,16 @@
restrict plotted particles to only those that are of a given type.
*minimum_mass* will require that the particles be of a given mass,
calculated via ParticleMassMsun, to be plotted. *alpha* determines
- each particle's opacity.
+ each particle's opacity. If stars are plotted, min_star_age and
+ max_star_age filter them (in years).
"""
_type_name = "particles"
region = None
_descriptor = None
def __init__(self, width, p_size=1.0, col='k', marker='o', stride=1.0,
ptype=None, stars_only=False, dm_only=False,
- minimum_mass=None, alpha=1.0):
+ minimum_mass=None, alpha=1.0, min_star_age=None,
+ max_star_age=None):
PlotCallback.__init__(self)
self.width = width
self.p_size = p_size
@@ -1083,10 +1099,12 @@
self.dm_only = dm_only
self.minimum_mass = minimum_mass
self.alpha = alpha
+ self.min_star_age = min_star_age
+ self.max_star_age = max_star_age
def __call__(self, plot):
data = plot.data
- # we construct a recantangular prism
+ # we construct a rectangular prism
x0, x1 = plot.xlim
y0, y1 = plot.ylim
xx0, xx1 = plot._axes.get_xlim()
@@ -1108,6 +1126,16 @@
if self.minimum_mass is not None:
gg &= (reg["ParticleMassMsun"] >= self.minimum_mass)
if gg.sum() == 0: return
+ if self.min_star_age is not None:
+ current_time = plot.pf.current_time
+ age_in_years = (current_time - reg["creation_time"]) * (plot.pf['Time'] / 3.155e7)
+ gg &= (age_in_years >= self.min_star_age)
+ if gg.sum() == 0: return
+ if self.max_star_age is not None:
+ current_time = plot.pf.current_time
+ age_in_years = (current_time - reg["creation_time"]) * (plot.pf['Time'] / 3.155e7)
+ gg &= (age_in_years <= self.max_star_age)
+ if gg.sum() == 0: return
plot._axes.hold(True)
px, py = self.convert_to_plot(plot,
[reg[field_x][gg][::self.stride],
diff -r c7bae05dd4ef1d7870e5d3e042f5f6ef2cb5c9bf -r 05e219be80cf3afd394c77db58ea8864d06ae410 yt/visualization/profile_plotter.py
--- a/yt/visualization/profile_plotter.py
+++ b/yt/visualization/profile_plotter.py
@@ -433,7 +433,7 @@
if field_name is None:
field_name = r'$\rm{'+field+r'}$'
elif field_name.find('$') == -1:
- field_name = r'$\rm{'+field+r'}$'
+ field_name = r'$\rm{'+field_name+r'}$'
if units is None or units == '':
label = field_name
else:
@@ -517,11 +517,10 @@
weight_field=None)
>>> plot.save()
- >>> # Change plot properties.
+ >>> # Change plot properties.
>>> plot.set_cmap("CellMassMsun", "jet")
>>> plot.set_zlim("CellMassMsun", 1e8, 1e13)
>>> plot.set_title("CellMassMsun", "This is a phase plot")
-
"""
x_log = None
y_log = None
@@ -532,6 +531,7 @@
_plot_valid = False
_plot_type = 'Phase'
+
def __init__(self, data_source, x_field, y_field, z_fields,
weight_field="CellMassMsun", x_bins=128, y_bins=128,
accumulation=False, fractional=False,
@@ -540,6 +540,14 @@
self.z_log = {}
self.z_title = {}
self._initfinished = False
+ self._xlimits = [0,0]
+ self._ylimits = [0,0]
+ self._setxlims = False
+ self._setylims = False
+ self._plottext = ""
+ self._textxpos = 0.0
+ self._textypos = 0.0
+
if profile is None:
profile = create_profile(data_source,
@@ -562,10 +570,11 @@
xfi = pf.field_info[field_x]
yfi = pf.field_info[field_y]
zfi = pf.field_info[field_z]
+
x_title = self.x_title or self._get_field_label(field_x, xfi)
y_title = self.y_title or self._get_field_label(field_y, yfi)
z_title = self.z_title.get(field_z, None) or \
- self._get_field_label(field_z, zfi)
+ self._get_field_label(field_z, zfi)
return (x_title, y_title, z_title)
def _get_field_label(self, field, field_info):
@@ -574,13 +583,13 @@
if field_name is None:
field_name = r'$\rm{'+field+r'}$'
elif field_name.find('$') == -1:
- field_name = r'$\rm{'+field+r'}$'
+ field_name = r'$\rm{'+field_name+r'}$'
if units is None or units == '':
label = field_name
else:
label = field_name+r'$\/\/('+units+r')$'
return label
-
+
def _get_field_log(self, field_z, profile):
pf = profile.data_source.pf
zfi = pf.field_info[field_z]
@@ -612,7 +621,16 @@
size = (self.figure_size, self.figure_size)
x_scale, y_scale, z_scale = self._get_field_log(f, self.profile)
+ x_label, y_label, z_label = self._get_axes_labels()
x_title, y_title, z_title = self._get_field_title(f, self.profile)
+ #If the labels are set they take precedence
+ if x_label is not None:
+ x_title = x_label
+ if y_label is not None:
+ y_title = y_label
+ if z_label is not None:
+ z_title = z_label
+
if f in self.plots:
zlim = [self.plots[f].zmin, self.plots[f].zmax]
else:
@@ -627,16 +645,24 @@
x_scale, y_scale, z_scale,
self._colormaps[f], zlim, size, fp.get_size(),
fig, axes, cax)
+
self.plots[f].axes.xaxis.set_label_text(x_title)
self.plots[f].axes.yaxis.set_label_text(y_title)
self.plots[f].cax.yaxis.set_label_text(z_title)
+ if(self._setxlims == True):
+ self.plots[f].axes.set_xlim(self._xlimits[0], self._xlimits[1])
+ if(self._setylims == True):
+ self.plots[f].axes.set_ylim(self._ylimits[0], self._ylimits[1])
+
+ self.plots[f].axes.text(self._textxpos, self._textypos, self._plottext,
+ fontproperties=self._font_properties)
if z_scale == "log":
self._field_transform[f] = log_transform
else:
self._field_transform[f] = linear_transform
if f in self.plot_title:
self.plots[f].axes.set_title(self.plot_title[f])
-
+
if self._font_color is not None:
ax = self.plots[f].axes
cbax = self.plots[f].cb.ax
@@ -648,6 +674,90 @@
label.set_color(self._font_color)
self._plot_valid = True
+
+ def set_xlim(self, xmin=None, xmax=None):
+ r"""
+ Sets the x-axis limits on the Phase plot.
+ Defaults to None leaving the axis unchanged
+ Parameters
+ ----------
+ xmin: float
+ The minimum value on the x-axis
+ xmax: float
+ The maximum value on the x-axis
+
+ >>> plot.set_xlim(5e-21, 1e5)
+ """
+
+ for f, data in self.profile.field_data.items():
+ axes = None
+ if f in self.plots:
+ if self.plots[f].figure is not None:
+ axes = self.plots[f].axes
+
+ self.plots[f].axes.set_xlim(xmin, xmax)
+ self._setxlims = True
+ self._xlimits[0] = xmin
+ self._xlimits[1] = xmax
+
+ def set_ylim(self, ymin=None, ymax=None):
+ r"""
+ Sets the y-axis limits on the Phase plot.
+ Defaults to None leaving the axis unchanged
+ Parameters
+ ----------
+ ymin: float
+ The minimum value on the y-axis
+ ymax: float
+ The maximum value on the y-axis
+
+ >>> plot.set_ylim(1e1, 1e5)
+
+ """
+
+ for f, data in self.profile.field_data.items():
+ axes = None
+ if f in self.plots:
+ if self.plots[f].figure is not None:
+ axes = self.plots[f].axes
+
+ self.plots[f].axes.set_ylim(ymin, ymax)
+ self._setylims = True
+ self._ylimits[0] = ymin
+ self._ylimits[1] = ymax
+
+ def add_text(self, text_str, xpos, ypos, fontsize=18, **kwargs):
+ r"""
+ Allow the user to insert text onto the plot
+ The x-position and y-position must be given as well as the text string.
+ Add text_str plot at location x, y, data coordinates (see example below).
+ Fontsize defaults to 18.
+
+ Parameters
+ ----------
+ text_str: str
+ The text to insert onto the plot. Required argument.
+ xpos: float
+ Position on plot in x-coordinates. Required argument.
+ ypos: float
+ Position on plot in y-coordinates. Required argument.
+ fontsize: float
+ Fontsize for the text (defaults to 18)
+
+ >>> plot.text(1e-15, 5e4, "Hello YT")
+
+ """
+ for f, data in self.profile.field_data.items():
+ axes = None
+ if f in self.plots:
+ if self.plots[f].figure is not None:
+ axes = self.plots[f].axes
+ self.plots[f].axes.text(xpos, ypos, text_str,
+ fontproperties=self._font_properties)
+ self._plottext = text_str
+ self._textxpos = xpos
+ self._textypos = ypos
+
def save(self, name=None, mpl_kwargs=None):
r"""
Saves a 2d profile plot.
https://bitbucket.org/yt_analysis/yt/commits/816186f16396/
Changeset: 816186f16396
Branch: stable
User: MatthewTurk
Date: 2014-07-22 23:30:02
Summary: Updating version to 2.6.3.
Affected #: 1 file
diff -r 05e219be80cf3afd394c77db58ea8864d06ae410 -r 816186f16396a16853810ac9ebcde5057d8d5b1a setup.py
--- a/setup.py
+++ b/setup.py
@@ -156,7 +156,7 @@
# End snippet
######
-VERSION = "2.6.2"
+VERSION = "2.6.3"
if os.path.exists('MANIFEST'):
os.remove('MANIFEST')
https://bitbucket.org/yt_analysis/yt/commits/0b29fa48fa3a/
Changeset: 0b29fa48fa3a
Branch: stable
User: MatthewTurk
Date: 2014-07-22 23:30:09
Summary: Added tag yt-2.6.3 for changeset 816186f16396
Affected #: 1 file
diff -r 816186f16396a16853810ac9ebcde5057d8d5b1a -r 0b29fa48fa3a25d0ed2e1cc28d8f946f878d18ac .hgtags
--- a/.hgtags
+++ b/.hgtags
@@ -5172,3 +5172,4 @@
053487f48672b8fd5c43af992e92bc2f2499f31f yt-2.6
d43ff9d8e20f2d2b8f31f4189141d2521deb341b yt-2.6.1
f1e22ef9f3a225f818c43262e6ce9644e05ffa21 yt-2.6.2
+816186f16396a16853810ac9ebcde5057d8d5b1a yt-2.6.3
Repository URL: https://bitbucket.org/yt_analysis/yt/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list