[Yt-svn] yt-commit r941 - trunk/yt/extensions/lightcone
britton at wrangler.dreamhost.com
britton at wrangler.dreamhost.com
Sun Nov 16 23:32:13 PST 2008
Author: britton
Date: Sun Nov 16 23:32:12 2008
New Revision: 941
URL: http://yt.spacepope.org/changeset/941
Log:
Added routines to create light cone projections that minimize the amount of
sampled volume they have in common. If requested, the routine will try to
use recycled solutions before making completely new realizations.
Also, there were a number of bugs in the light cone generator that only showed
up when trying to actually make recycled projections or more than one
projection from a single light cone object. That should all be fixed now.
Added:
trunk/yt/extensions/lightcone/UniqueSolution.py
Modified:
trunk/yt/extensions/lightcone/LightCone.py
trunk/yt/extensions/lightcone/LightConeProjection.py
Modified: trunk/yt/extensions/lightcone/LightCone.py
==============================================================================
--- trunk/yt/extensions/lightcone/LightCone.py (original)
+++ trunk/yt/extensions/lightcone/LightCone.py Sun Nov 16 23:32:12 2008
@@ -25,13 +25,16 @@
from yt.extensions.lightcone import *
from yt.logger import lagosLogger as mylog
+import copy
+import os
import numpy as na
import random as rand
import tables as h5
from Common_nVolume import *
class LightCone(object):
- def __init__(self,EnzoParameterFile,LightConeParameterFile):
+ def __init__(self,EnzoParameterFile,LightConeParameterFile,verbose=True):
+ self.verbose = verbose
self.EnzoParameterFile = EnzoParameterFile
self.LightConeParameterFile = LightConeParameterFile
self.enzoParameters = {}
@@ -40,6 +43,7 @@
self.timeOutputs = []
self.allOutputs = []
self.lightConeSolution = []
+ self.masterSolution = [] # kept to compare with recycled solutions
self.projectionStack = []
self.projectionWeightFieldStack = []
@@ -54,6 +58,14 @@
self._ReadLightConeParameterFile()
self._ReadEnzoParameterFile()
+ # Create output directory.
+ if (os.path.exists(self.lightConeParameters['OutputDir'])):
+ if not(os.path.isdir(self.lightConeParameters['OutputDir'])):
+ mylog.error("Output directory exists, but is not a directory: %s." % self.lightConeParameters['OutputDir'])
+ self.lightConeParameters['OutputDir'] = './'
+ else:
+ os.mkdir(self.lightConeParameters['OutputDir'])
+
# Calculate redshifts for dt data dumps.
if (self.enzoParameters.has_key('dtDataDump')):
self._CalculateTimeDumpRedshifts()
@@ -71,7 +83,7 @@
self.recycleSolution = False
if seed is not None:
- self.lightConeParameters['RandomSeed'] = seed
+ self.lightConeParameters['RandomSeed'] = int(seed)
# Make light cone using minimum number of projections.
if (self.lightConeParameters['UseMinimumNumberOfProjections']):
@@ -94,14 +106,14 @@
while (z > output['redshift']):
output = output['previous']
if (output is None):
- mylog.error("CalculateLightConeSolution: search for data output went off the end the stack.")
- mylog.error("Could not calculate light cone solution.")
+ if self.verbose: mylog.error("CalculateLightConeSolution: search for data output went off the end the stack.")
+ if self.verbose: mylog.error("Could not calculate light cone solution.")
return
if (output['redshift'] == self.lightConeSolution[-1]['redshift']):
- mylog.error("CalculateLightConeSolution: No data dump between z = %f and %f." % \
+ if self.verbose: mylog.error("CalculateLightConeSolution: No data dump between z = %f and %f." % \
((self.lightConeSolution[-1]['redshift'] - self.lightConeSolution[-1]['deltaz']),
self.lightConeSolution[-1]['redshift']))
- mylog.error("Could not calculate light cone solution.")
+ if self.verbose: mylog.error("Could not calculate light cone solution.")
return
self.lightConeSolution.append(output)
z = self.lightConeSolution[-1]['redshift'] - self.lightConeSolution[-1]['deltaz']
@@ -123,7 +135,7 @@
self.lightConeSolution.append(nextOutput)
nextOutput = nextOutput['next']
- mylog.info("CalculateLightConeSolution: Used %d data dumps to get from z = %f to %f." % (len(self.lightConeSolution),
+ if self.verbose: mylog.info("CalculateLightConeSolution: Used %d data dumps to get from z = %f to %f." % (len(self.lightConeSolution),
self.lightConeParameters['InitialRedshift'],
self.lightConeParameters['FinalRedshift']))
@@ -145,9 +157,9 @@
self.enzoParameters['CosmologyHubbleConstantNow'] / self.enzoParameters['CosmologyComovingBoxSize']
# Simple error check to make sure more than 100% of box depth is never required.
if (self.lightConeSolution[q]['DepthBoxFraction'] > 1.0):
- mylog.error("Warning: box fraction required to go from z = %f to %f is %f" % (self.lightConeSolution[q]['redshift'],z_next,
+ if self.verbose: mylog.error("Warning: box fraction required to go from z = %f to %f is %f" % (self.lightConeSolution[q]['redshift'],z_next,
self.lightConeSolution[q]['DepthBoxFraction']))
- mylog.error("Full box delta z is %f, but it is %f to the next data dump." % (self.lightConeSolution[q]['deltaz'],
+ if self.verbose: mylog.error("Full box delta z is %f, but it is %f to the next data dump." % (self.lightConeSolution[q]['deltaz'],
self.lightConeSolution[q]['redshift']-z_next))
# Calculate fraction of box required for width corresponding to requested image size.
@@ -161,6 +173,9 @@
self.lightConeSolution[q]['ProjectionAxis'] = rand.randint(0,2)
self.lightConeSolution[q]['ProjectionCenter'] = [rand.random(),rand.random(),rand.random()]
+ # Store this as the mast solution.
+ self.masterSolution = [copy.deepcopy(q) for q in self.lightConeSolution]
+
# Clear out some stuff.
del self.allOutputs
del self.redshiftOutputs
@@ -170,6 +185,12 @@
def ProjectLightCone(self,field,weight_field=None,**kwargs):
"Create projections for light cone, then add them together."
+ # Clear projection stack.
+ self.projectionStack = []
+ self.projectionWeightFieldStack = []
+ if (self.lightConeSolution[-1].has_key('object')):
+ del self.lightConeSolution[-1]['object']
+
if not(self.lightConeParameters['OutputDir'].endswith("/")):
self.lightConeParameters['OutputDir'] += "/"
@@ -179,7 +200,7 @@
for q,output in enumerate(self.lightConeSolution):
name = "%s%s_%04d_%04d" % (self.lightConeParameters['OutputDir'],self.lightConeParameters['OutputPrefix'],
q,len(self.lightConeSolution))
- output['object'] = EnzoStaticOutput(output['filename'])
+ output['object'] = lagos.EnzoStaticOutput(output['filename'])
frb = LightConeProjection(output,field,pixels,weight_field=weight_field,
save_image=self.lightConeParameters['SaveLightConeSlices'],
name=name,**kwargs)
@@ -234,12 +255,12 @@
"""
if recycle:
- mylog.info("Recycling solution made with %s with new seed %s." % (self.lightConeParameters['RandomSeed'],
+ if self.verbose: mylog.info("Recycling solution made with %s with new seed %s." % (self.lightConeParameters['RandomSeed'],
newSeed))
- self.recycleRandomSeed = newSeed
+ self.recycleRandomSeed = int(newSeed)
else:
- mylog.info("Creating new solution with random seed %s." % newSeed)
- self.lightConeParameters['RandomSeed'] = newSeed
+ if self.verbose: mylog.info("Creating new solution with random seed %s." % newSeed)
+ self.lightConeParameters['RandomSeed'] = int(newSeed)
self.recycleRandomSeed = 0
self.recycleSolution = recycle
@@ -249,16 +270,18 @@
totalVolume = 0.0
# Seed random number generator with new seed.
- rand.seed(newSeed)
+ rand.seed(int(newSeed))
for q,output in enumerate(self.lightConeSolution):
# It is necessary to make the same number of calls to the random number generator
# so the original solution willbe produced if the same seed is given.
# Get random projection axis and center.
- # If recyclsing, axis will get thrown away since it is used in creating a unique projection object.
+ # If recycling, axis will get thrown away since it is used in creating a unique projection object.
newAxis = rand.randint(0,2)
- if not recycle:
+ if recycle:
+ output['ProjectionAxis'] = self.masterSolution[q]['ProjectionAxis']
+ else:
output['ProjectionAxis'] = newAxis
newCenter = [rand.random(),rand.random(),rand.random()]
@@ -267,20 +290,22 @@
newCube = na.zeros(shape=(len(newCenter),2))
oldCube = na.zeros(shape=(len(newCenter),2))
for w in range(len(newCenter)):
+ if (w == self.masterSolution[q]['ProjectionAxis']):
+ oldCube[w] = [self.masterSolution[q]['ProjectionCenter'][w] - 0.5 * self.masterSolution[q]['DepthBoxFraction'],
+ self.masterSolution[q]['ProjectionCenter'][w] + 0.5 * self.masterSolution[q]['DepthBoxFraction']]
+ else:
+ oldCube[w] = [self.masterSolution[q]['ProjectionCenter'][w] - 0.5 * self.masterSolution[q]['WidthBoxFraction'],
+ self.masterSolution[q]['ProjectionCenter'][w] + 0.5 * self.masterSolution[q]['WidthBoxFraction']]
+
if (w == output['ProjectionAxis']):
- oldCube[w] = [output['ProjectionCenter'][w] - 0.5 * output['DepthBoxFraction'],
- output['ProjectionCenter'][w] + 0.5 * output['DepthBoxFraction']]
- # If recycling old solution, then keep center for axis in line of sight.
if recycle:
newCube[w] = oldCube[w]
else:
- newCube[w] = [newCenter[w] - 0.5 * output['DepthBoxFraction'],
- newCenter[w] + 0.5 * output['DepthBoxFraction']]
+ newCube[w] = [newCenter[w] - 0.5 * self.masterSolution[q]['DepthBoxFraction'],
+ newCenter[w] + 0.5 * self.masterSolution[q]['DepthBoxFraction']]
else:
- oldCube[w] = [output['ProjectionCenter'][w] - 0.5 * output['WidthBoxFraction'],
- output['ProjectionCenter'][w] + 0.5 * output['WidthBoxFraction']]
- newCube[w] = [newCenter[w] - 0.5 * output['WidthBoxFraction'],
- newCenter[w] + 0.5 * output['WidthBoxFraction']]
+ newCube[w] = [newCenter[w] - 0.5 * self.masterSolution[q]['WidthBoxFraction'],
+ newCenter[w] + 0.5 * self.masterSolution[q]['WidthBoxFraction']]
commonVolume += commonNVolume(oldCube,newCube,periodic=na.array([[0,1],[0,1],[0,1]]))
totalVolume += output['DepthBoxFraction'] * output['WidthBoxFraction']**2
@@ -290,12 +315,20 @@
if not(recycle and (w == self.lightConeSolution[q]['ProjectionAxis'])):
self.lightConeSolution[q]['ProjectionCenter'][w] = newCenter[w]
- mylog.info("Fraction of total volume in common with old solution is %.2e." % (commonVolume/totalVolume))
+ if recycle:
+ if self.verbose: mylog.info("Fractional common volume between master and recycled solution is %.2e" % (commonVolume/totalVolume))
+ else:
+ if self.verbose: mylog.info("Fraction of total volume in common with old solution is %.2e." % (commonVolume/totalVolume))
+ self.masterSolution = [copy.deepcopy(q) for q in self.lightConeSolution]
+
+ def RestoreMasterSolution(self):
+ "Reset the active light cone solution to the master solution."
+ self.lightConeSolution = [copy.deepcopy(q) for q in self.masterSolution]
def SaveLightConeSolution(self,file="light_cone.dat"):
"Write out a text file with information on light cone solution."
- mylog.info("Saving light cone solution to %s." % file)
+ if self.verbose: mylog.info("Saving light cone solution to %s." % file)
f = open(file,'w')
if self.recycleSolution:
@@ -318,14 +351,15 @@
def SaveLightConeStack(self,filename=None):
"Save the light cone projection stack as a 3d array in and hdf5 file."
if (filename is None):
- filename = "%s%s_data" % (self.lightConeParameters['OutputDir'],self.lightConeParameters['OutputPrefix'])
- filename += ".h5"
+ filename = "%s/%s_data" % (self.lightConeParameters['OutputDir'],self.lightConeParameters['OutputPrefix'])
+ if not(filename.endswith('.h5')):
+ filename += ".h5"
if (len(self.projectionStack) == 0):
- mylog.error("SaveLightConeStack: no projection data loaded.")
+ if self.verbose: mylog.error("SaveLightConeStack: no projection data loaded.")
return
- mylog.info("Writing light cone data to %s." % filename)
+ if self.verbose: mylog.info("Writing light cone data to %s." % filename)
output = h5.openFile(filename, "a")
@@ -368,7 +402,7 @@
distance2 = co.ComovingRadialDistance(z2,z) * self.enzoParameters['CosmologyHubbleConstantNow']
iteration += 1
if (iteration > max_Iterations):
- mylog.error("CalculateDeltaZMax: Warning - max iterations exceeded for z = %f (delta z = %f)." % (z,na.fabs(z2-z)))
+ if self.verbose: mylog.error("CalculateDeltaZMax: Warning - max iterations exceeded for z = %f (delta z = %f)." % (z,na.fabs(z2-z)))
break
output['deltaz'] = na.fabs(z2-z)
@@ -427,7 +461,7 @@
param = param.strip()
vals = vals.strip()
except ValueError:
- mylog.error("Skipping line: %s" % line)
+ if self.verbose: mylog.error("Skipping line: %s" % line)
continue
if EnzoParameterDict.has_key(param):
t = map(EnzoParameterDict[param], vals.split())
@@ -462,7 +496,7 @@
param = param.strip()
vals = vals.strip()
except ValueError:
- mylog.error("Skipping line: %s" % line)
+ if self.verbose: mylog.error("Skipping line: %s" % line)
continue
if LightConeParameterDict.has_key(param):
t = map(LightConeParameterDict[param], vals.split())
Modified: trunk/yt/extensions/lightcone/LightConeProjection.py
==============================================================================
--- trunk/yt/extensions/lightcone/LightConeProjection.py (original)
+++ trunk/yt/extensions/lightcone/LightConeProjection.py Sun Nov 16 23:32:12 2008
@@ -24,6 +24,11 @@
"""
from yt.extensions.lightcone import *
+from yt.logger import lagosLogger as mylog
+import yt.lagos as lagos
+import yt.raven as raven
+import numpy as na
+import copy
def LightConeProjection(lightConeSlice,field,pixels,weight_field=None,save_image=False,name="",**kwargs):
"Create a single projection to be added into the light cone stack."
@@ -106,7 +111,7 @@
# 3. The Shift Problem
# Shift projection by random x and y offsets.
- offset = lightConeSlice['ProjectionCenter']
+ offset = copy.deepcopy(lightConeSlice['ProjectionCenter'])
# Delete depth coordinate.
del offset[lightConeSlice['ProjectionAxis']]
Added: trunk/yt/extensions/lightcone/UniqueSolution.py
==============================================================================
--- (empty file)
+++ trunk/yt/extensions/lightcone/UniqueSolution.py Sun Nov 16 23:32:12 2008
@@ -0,0 +1,216 @@
+"""
+Functions to generate unique light cone solutions.
+
+Author: Britton Smith <brittons at origins.colorado.edu>
+Affiliation: CASA/University of CO, Boulder
+Homepage: http://yt.enzotools.org/
+License:
+ Copyright (C) 2008 Britton Smith. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from LightCone import *
+from Common_nVolume import *
+from yt.logger import lagosLogger as mylog
+import numpy as na
+import random as rand
+import sys
+
+def ProjectUniqueLightCones(EnzoParameterFile,LightConeParameterFile,seedFile,field,save_stack=True,**kwargs):
+ "Make light cone projections using list of random seeds in a file."
+
+ seedList = _ReadSeedFile(seedFile)
+
+ lc = LightCone(EnzoParameterFile,LightConeParameterFile)
+ prefix = lc.lightConeParameters['OutputPrefix']
+ lc.CalculateLightConeSolution(seed=0)
+ lastSeed = None
+
+ for seed in seedList:
+ if (seed['master'] != lastSeed):
+ lc.RerandomizeLightConeSolution(seed['master'],recycle=False)
+ lastSeed = seed['master']
+ if (seed['recycle'] is not None):
+ lc.RerandomizeLightConeSolution(seed['recycle'],recycle=True)
+
+ lc.lightConeParameters['OutputPrefix'] = "%s_%s_%s" % (prefix,seed['master'],seed['recycle'])
+ lc.ProjectLightCone(field,**kwargs)
+ if save_stack:
+ stackFilename = "%s/%s_%s.h5" % (lc.lightConeParameters['OutputDir'],
+ lc.lightConeParameters['OutputPrefix'],field)
+ lc.SaveLightConeStack(filename=stackFilename)
+
+def FindUniqueSolutions(EnzoParameterFile,LightConeParameterFile,solutions=100,seed=None,
+ max_overlap=0.25,failures=10,recycle=True,filename='unique.dat'):
+ "Find a set of random seeds that will give light cones will minimal volume overlap."
+
+ solution1 = LightCone(EnzoParameterFile,LightConeParameterFile,verbose=False)
+ solution2 = LightCone(EnzoParameterFile,LightConeParameterFile,verbose=False)
+ solution1.CalculateLightConeSolution(seed=0)
+ solution2.CalculateLightConeSolution(seed=0)
+
+ uniqueSeeds = []
+ if recycle:
+ master = None
+ recycleFails = 0
+ newRecycleSeed = None
+ fails = 0
+
+ maxCommon = 0.0
+
+ # Need to continuall save and reset the state of the random number generator
+ # since it is being reset by the light cone generator.
+ if seed is None:
+ state = None
+ else:
+ rand.seed(seed)
+ state = rand.getstate()
+
+ while (len(uniqueSeeds) < solutions):
+ sys.stderr.write("Created %d unique solutions.\r" % len(uniqueSeeds))
+
+ # Create new random seed.
+ if (recycle and master is not None):
+ newSeed = master
+ if state is not None: rand.setstate(state)
+ newRecycleSeed = rand.randint(1,1e9)
+ state = rand.getstate()
+ else:
+ if state is not None: rand.setstate(state)
+ newSeed = rand.randint(1,1e9)
+ state = rand.getstate()
+ if recycle:
+ master = newSeed
+ recycleFails = 0
+ newRecycleSeed = None
+
+ solution1.RerandomizeLightConeSolution(newSeed,recycle=False)
+ if newRecycleSeed is not None:
+ solution1.RerandomizeLightConeSolution(newRecycleSeed,recycle=True)
+
+ # Compare with all other seeds.
+ testPass = True
+ for uniqueSeed in uniqueSeeds:
+ solution2.RerandomizeLightConeSolution(uniqueSeed['master'],recycle=False)
+ if uniqueSeed['recycle'] is not None:
+ solution2.RerandomizeLightConeSolution(uniqueSeed['recycle'],recycle=True)
+
+ common = CompareSolutions(solution1.lightConeSolution,solution2.lightConeSolution)
+
+ if (common > max_overlap):
+ testPass = False
+ break
+ else:
+ maxCommon = max(maxCommon,common)
+
+ if testPass:
+ uniqueSeeds.append({'master':newSeed,'recycle':newRecycleSeed})
+ fails = 0
+ recycleFails = 0
+ else:
+ if newRecycleSeed is None:
+ fails += 1
+ else:
+ recycleFails += 1
+ if (recycleFails >= failures):
+ fails += 1
+ print ""
+ mylog.info("Max recycled failures reached with master seed %d." % newSeed)
+ master = None
+ if (fails >= failures):
+ print ""
+ mylog.error("Max consecutive failures reached.")
+ break
+
+ mylog.info("Created %d unique solutions." % len(uniqueSeeds))
+ mylog.info("Maximum common volume is %.2e." % maxCommon)
+ _WriteSeedFile(uniqueSeeds,filename)
+ return uniqueSeeds
+
+def CompareSolutions(solution1,solution2):
+ "Calculate common volume between two light cone solutions."
+
+ if (len(solution1) != len(solution2)):
+ mylog.error("Cannot compare light cone solutions with unequal numbers of slices.")
+ return -1
+
+ commonVolume = 0.0
+ totalVolume = 0.0
+
+ # Check that solution volumes are the same.
+ if((solution1[0]['DepthBoxFraction'] * solution1[0]['WidthBoxFraction']**2) !=
+ (solution2[0]['DepthBoxFraction'] * solution2[0]['WidthBoxFraction']**2)):
+ mylog.error("Light cone solutions do not have equal volumes, will use the smaller one.")
+
+ for q in range(len(solution1)):
+ cube1 = na.zeros(shape=(len(solution1[q]['ProjectionCenter']),2))
+ volume1 = 1.0
+ for w in range(len(cube1)):
+ if (w == solution1[q]['ProjectionAxis']):
+ width = solution1[q]['DepthBoxFraction']
+ else:
+ width = solution1[q]['WidthBoxFraction']
+ volume1 *= width
+ cube1[w] = [solution1[q]['ProjectionCenter'][w] - 0.5 * width,
+ solution1[q]['ProjectionCenter'][w] + 0.5 * width]
+
+ cube2 = na.zeros(shape=(len(solution2[q]['ProjectionCenter']),2))
+ volume2 = 1.0
+ for w in range(len(cube2)):
+ if (w == solution2[q]['ProjectionAxis']):
+ width = solution2[q]['DepthBoxFraction']
+ else:
+ width = solution2[q]['WidthBoxFraction']
+ volume2 *= width
+ cube2[w] = [solution2[q]['ProjectionCenter'][w] - 0.5 * width,
+ solution2[q]['ProjectionCenter'][w] + 0.5 * width]
+
+ totalVolume += min(volume1,volume2)
+ commonVolume += commonNVolume(cube1,cube2,periodic=na.array([[0,1],[0,1],[0,1]]))
+
+ return (commonVolume/totalVolume)
+
+def _ReadSeedFile(filename):
+ "Read list of random seeds from a file."
+
+ mylog.info("Reading random seed list from %s." % filename)
+
+ seedList = []
+
+ lines = file(filename)
+ for line in lines:
+ line = line.strip()
+ onLine = line.split(',')
+ if (len(onLine) == 1):
+ seedList.append({'master':onLine[0], 'recycle':None})
+ else:
+ seedList.append({'master':onLine[0], 'recycle':onLine[1]})
+
+ return seedList
+
+def _WriteSeedFile(seedList,filename):
+ "Write list of random seeds to a file."
+
+ mylog.info("Writing random seed list to %s." % filename)
+
+ f = open(filename,'w')
+ for seed in seedList:
+ if seed['recycle'] is None:
+ f.write("%s\n" % seed['master'])
+ else:
+ f.write("%s,%s\n" % (seed['master'],seed['recycle']))
+ f.close()
More information about the yt-svn
mailing list