[Yt-svn] yt-commit r917 - trunk/yt/extensions/lightcone
britton at wrangler.dreamhost.com
britton at wrangler.dreamhost.com
Tue Nov 11 13:58:37 PST 2008
Author: britton
Date: Tue Nov 11 13:58:37 2008
New Revision: 917
URL: http://yt.spacepope.org/changeset/917
Log:
The light cone generator will now calculate the fractional volume in common
between an original solution and a recycled one made using the original. This
will help to avoid making light cone realizations that sample the same data.
Added:
trunk/yt/extensions/lightcone/Common_nVolume.py
Modified:
trunk/yt/extensions/lightcone/LightCone.py
Added: trunk/yt/extensions/lightcone/Common_nVolume.py
==============================================================================
--- (empty file)
+++ trunk/yt/extensions/lightcone/Common_nVolume.py Tue Nov 11 13:58:37 2008
@@ -0,0 +1,101 @@
+import numpy as na
+
+def commonNVolume(nCube1,nCube2,periodic=None):
+ "Return the n-volume in common between the two n-cubes."
+
+ # Check for proper args.
+ if ((len(na.shape(nCube1)) != 2) or
+ (na.shape(nCube1)[1] != 2) or
+ (na.shape(nCube1) != na.shape(nCube2))):
+ print "Arguments must be 2 (n,2) numpy array."
+ return 0
+
+ if ((periodic is not None) and
+ (na.shape(nCube1) != na.shape(periodic))):
+ print "periodic argument must be (n,2) numpy array."
+ return 0
+
+ nCommon = 1.0
+ for q in range(na.shape(nCube1)[0]):
+ if (periodic is None):
+ nCommon *= commonSegment(nCube1[q],nCube2[q])
+ else:
+ nCommon *= commonSegment(nCube1[q],nCube2[q],periodic=periodic[q])
+
+ return nCommon
+
+def commonSegment(seg1,seg2,periodic=None):
+ "Return the length of the common segment."
+
+ # Check for proper args.
+ if ((len(seg1) != 2) or (len(seg2) != 2)):
+ print "Arguments must be arrays of size 2."
+ return 0
+
+ # If not periodic, then this is very easy.
+ if periodic is None:
+ seg1.sort()
+ len1 = seg1[1] - seg1[0]
+ seg2.sort()
+ len2 = seg2[1] - seg2[0]
+
+ common = 0.0
+
+ add = seg1[1] - seg2[0]
+ if ((add > 0) and (add <= max(len1,len2))):
+ common += add
+ add = seg2[1] - seg1[0]
+ if ((add > 0) and (add <= max(len1,len2))):
+ common += add
+ common = min(common,len1,len2)
+ return common
+
+ # If periodic, it's a little more complicated.
+ else:
+ if len(periodic) != 2:
+ print "periodic array must be of size 2."
+ return 0
+
+ seg1.sort()
+ flen1 = seg1[1] - seg1[0]
+ len1 = flen1 - int(flen1)
+ seg2.sort()
+ flen2 = seg2[1] - seg2[0]
+ len2 = flen2 - int(flen2)
+
+ periodic.sort()
+ scale = periodic[1] - periodic[0]
+
+ if (abs(int(flen1)-int(flen2)) >= scale):
+ return min(flen1,flen2)
+
+ # Adjust for periodicity
+ seg1[0] = na.mod(seg1[0],scale) + periodic[0]
+ seg1[1] = seg1[0] + len1
+ if (seg1[1] > periodic[1]): seg1[1] -= scale
+ seg2[0] = na.mod(seg2[0],scale) + periodic[0]
+ seg2[1] = seg2[0] + len2
+ if (seg2[1] > periodic[1]): seg2[1] -= scale
+
+ # create list of non-periodic segments
+ pseg1 = []
+ if (seg1[0] >= seg1[1]):
+ pseg1.append([seg1[0],periodic[1]])
+ pseg1.append([periodic[0],seg1[1]])
+ else:
+ pseg1.append(seg1)
+ pseg2 = []
+ if (seg2[0] >= seg2[1]):
+ pseg2.append([seg2[0],periodic[1]])
+ pseg2.append([periodic[0],seg2[1]])
+ else:
+ pseg2.append(seg2)
+
+ # Add up common segments.
+ common = min(int(flen1),int(flen2))
+
+ for subseg1 in pseg1:
+ for subseg2 in pseg2:
+ common += commonSegment(subseg1,subseg2)
+
+ return common
Modified: trunk/yt/extensions/lightcone/LightCone.py
==============================================================================
--- trunk/yt/extensions/lightcone/LightCone.py (original)
+++ trunk/yt/extensions/lightcone/LightCone.py Tue Nov 11 13:58:37 2008
@@ -28,6 +28,7 @@
import numpy as na
import random as rand
import tables as h5
+from Common_nVolume import *
class LightCone(object):
def __init__(self,EnzoParameterFile,LightConeParameterFile):
@@ -232,6 +233,10 @@
self.recycleSolution = True
self.recycleRandomSeed = newSeed
+ # Keep track of fraction of volume in common between the original and recycled solution.
+ commonVolume = 0.0
+ totalVolume = 0.0
+
# Seed random number generator with new seed.
rand.seed(self.recycleRandomSeed)
@@ -245,11 +250,29 @@
newCenter = [rand.random(),rand.random(),rand.random()]
+ # Make list of rectangle corners to calculate common volume.
+ newCube = na.zeros(shape=(len(newCenter),2))
+ oldCube = na.zeros(shape=(len(newCenter),2))
+ for w in range(len(newCenter)):
+ if (w == output['ProjectionAxis']):
+ newCube[w] = [output['ProjectionCenter'][w] - 0.5 * output['DepthBoxFraction'],
+ output['ProjectionCenter'][w] + 0.5 * output['DepthBoxFraction']]
+ oldCube[w] = newCube[w]
+ else:
+ newCube[w] = [newCenter[w] - 0.5 * output['WidthBoxFraction'],
+ newCenter[w] + 0.5 * output['WidthBoxFraction']]
+ oldCube[w] = [output['ProjectionCenter'][w] - 0.5 * output['WidthBoxFraction'],
+ output['ProjectionCenter'][w] + 0.5 * output['WidthBoxFraction']]
+
+ commonVolume += commonNVolume(oldCube,newCube,periodic=na.array([[0,1],[0,1],[0,1]]))
+ totalVolume += output['DepthBoxFraction'] * output['WidthBoxFraction']**2
+
# Replaces centers for every axis except the line of sight axis.
for w in range(len(newCenter)):
if (w != self.lightConeSolution[q]['ProjectionAxis']):
self.lightConeSolution[q]['ProjectionCenter'][w] = newCenter[w]
+ mylog.info("Recycled solution has %.2e of total volume in common with the original." % (commonVolume/totalVolume))
def SaveLightConeSolution(self,file="light_cone.dat"):
"Write out a text file with information on light cone solution."
More information about the yt-svn
mailing list