223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

changed 'range' keyword to 'limit', now using proper ASCII table (with column 'intensity') for linearODF 

hybridIO_linODFsampling simplified
This commit is contained in:
Martin Diehl 2015-04-15 14:37:46 +00:00
parent e29628b459
commit f5762209dc
3 changed files with 165 additions and 189 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
code/IO.f90
View File

@ -611,7 +611,7 @@ function IO_hybridIA(Nast,ODFfileName)
read(FILEUNIT,'(a1024)') line read(FILEUNIT,'(a1024)') line
positions = IO_stringPos(line,7_pInt) positions = IO_stringPos(line,7_pInt)
select case ( IO_lc(IO_StringValue(line,positions,1_pInt,.true.)) ) select case ( IO_lc(IO_StringValue(line,positions,1_pInt,.true.)) )
case ('range') case ('limit')
gotRange = .true. gotRange = .true.
do j = 2_pInt,6_pInt,2_pInt do j = 2_pInt,6_pInt,2_pInt
select case (IO_lc(IO_stringValue(line,positions,j))) select case (IO_lc(IO_stringValue(line,positions,j)))

5
code/config/example.linearODF
View File

@ -1,7 +1,8 @@
3 header 4 header
range phi1 90 PHI 90 phi2 90 limit phi1 90 PHI 90 phi2 90
delta phi1 5 PHI 5 phi2 5 delta phi1 5 PHI 5 phi2 5
origin voxelBoundary origin voxelBoundary
intensity
0.905 0.905
0.4025 0.4025
0.1925 0.1925

347
processing/pre/hybridIA_linODFsampling.py
View File

@ -3,11 +3,12 @@
from optparse import OptionParser from optparse import OptionParser
import damask import damask
import os,sys,math,re,random,string import os,sys,math,re,random,string
import numpy as np
scriptID = string.replace('$Id$','\n','\\n') scriptID = string.replace('$Id$','\n','\\n')
scriptName = scriptID.split()[1] scriptName = scriptID.split()[1]
random.seed() random.seed(1)
# --- helper functions --- # --- helper functions ---
@ -25,10 +26,9 @@ def binsAsBin(bins,intervals):
def EulersAsBins(Eulers,intervals,deltas,center): def EulersAsBins(Eulers,intervals,deltas,center):
""" return list of Eulers translated into 3D bins list """ """ return list of Eulers translated into 3D bins list """
return [\ return [int((euler+(0.5-center)*delta)//delta)%interval \
int((euler+(0.5-center)*delta)//delta)%interval \ for euler,delta,interval in zip(Eulers,deltas,intervals) \
for euler,delta,interval in zip(Eulers,deltas,intervals) \ ]
]
def binAsEulers(bin,intervals,deltas,center): def binAsEulers(bin,intervals,deltas,center):
""" compound bin number translated into list of Eulers """ """ compound bin number translated into list of Eulers """
@ -46,24 +46,17 @@ def directInvRepetitions(probability,scale):
return nDirectInv return nDirectInv
# --- sampling methods --- # ---------------------- sampling methods -----------------------------------------------------------------------
# ----- efficient algorithm --------- # ----- efficient algorithm ---------
def directInversion (ODF,nSamples): def directInversion (ODF,nSamples):
""" ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians) """ """ ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians) """
nBins = ODF['intervals'][0]*ODF['intervals'][1]*ODF['intervals'][2] nOptSamples = max(ODF['nNonZero'],nSamples) # random subsampling if too little samples requested
deltas = [limit/intervals for limit,intervals in zip(ODF['limits'],ODF['intervals'])]
# calculate repetitions of each orientation
if re.search(r'hybrid',sys.argv[0],re.IGNORECASE): # my script's name contains "hybrid"
nOptSamples = max(ODF['nNonZero'],nSamples) # random subsampling if too little samples requested
else: # blunt integer approximation
nOptSamples = nSamples
nInvSamples = 0 nInvSamples = 0
repetition = [None]*nBins repetition = [None]*ODF['nBins']
probabilityScale = nOptSamples # guess probabilityScale = nOptSamples # guess
scaleLower = 0.0 scaleLower = 0.0
@ -89,27 +82,27 @@ def directInversion (ODF,nSamples):
nInvSamples = nInvSamplesUpper nInvSamples = nInvSamplesUpper
scale = scaleUpper scale = scaleUpper
print 'created set of',nInvSamples,'samples (',float(nInvSamples)/nOptSamples-1.0,') with scaling',scale,'delivering',nSamples print 'created set of',nInvSamples,'samples (',float(nInvSamples)/nOptSamples-1.0,') with scaling',scale,'delivering',nSamples
repetition = [None]*nBins # preallocate and clear repetition = [None]*ODF['nBins'] # preallocate and clear
for bin in range(nBins): # loop over bins for bin in range(ODF['nBins']): # loop over bins
repetition[bin] = int(round(ODF['dV_V'][bin]*scale)) # calc repetition repetition[bin] = int(round(ODF['dV_V'][bin]*scale)) # calc repetition
# build set # build set
set = [None]*nInvSamples set = [None]*nInvSamples
i = 0 i = 0
for bin in range(nBins): for bin in range(ODF['nBins']):
set[i:i+repetition[bin]] = [bin]*repetition[bin] # fill set with bin, i.e. orientation set[i:i+repetition[bin]] = [bin]*repetition[bin] # fill set with bin, i.e. orientation
i += repetition[bin] # advance set counter i += repetition[bin] # advance set counter
orientations = [None]*nSamples orientations = [None]*nSamples
reconstructedODF = [0.0]*nBins reconstructedODF = [0.0]*ODF['nBins']
unitInc = 1.0/nSamples unitInc = 1.0/nSamples
for j in range(nSamples): for j in range(nSamples):
if (j == nInvSamples-1): ex = j if (j == nInvSamples-1): ex = j
else: ex = int(round(random.uniform(j+0.5,nInvSamples-0.5))) else: ex = int(round(random.uniform(j+0.5,nInvSamples-0.5)))
bin = set[ex] bin = set[ex]
bins = binAsBins(bin,ODF['intervals']) bins = binAsBins(bin,ODF['interval'])
Eulers = binAsEulers(bin,ODF['intervals'],deltas,ODF['center']) Eulers = binAsEulers(bin,ODF['interval'],ODF['delta'],ODF['center'])
orientations[j] = '%g\t%g\t%g' %( math.degrees(Eulers[0]),math.degrees(Eulers[1]),math.degrees(Eulers[2]) ) orientations[j] = '%g\t%g\t%g' %( math.degrees(Eulers[0]),math.degrees(Eulers[1]),math.degrees(Eulers[2]) )
reconstructedODF[bin] += unitInc reconstructedODF[bin] += unitInc
set[ex] = set[j] # exchange orientations set[ex] = set[j] # exchange orientations
@ -124,10 +117,8 @@ def MonteCarloEulers (ODF,nSamples):
countMC = 0 countMC = 0
maxdV_V = max(ODF['dV_V']) maxdV_V = max(ODF['dV_V'])
nBins = ODF['intervals'][0]*ODF['intervals'][1]*ODF['intervals'][2]
deltas = [limit/intervals for limit,intervals in zip(ODF['limits'],ODF['intervals'])]
orientations = [None]*nSamples orientations = [None]*nSamples
reconstructedODF = [0.0]*nBins reconstructedODF = [0.0]*ODF['nBins']
unitInc = 1.0/nSamples unitInc = 1.0/nSamples
for j in range(nSamples): for j in range(nSamples):
@ -136,9 +127,9 @@ def MonteCarloEulers (ODF,nSamples):
while MC > ODF['dV_V'][bin]: while MC > ODF['dV_V'][bin]:
countMC += 1 countMC += 1
MC = maxdV_V*random.random() MC = maxdV_V*random.random()
Eulers = [limit*random.random() for limit in ODF['limits']] Eulers = [limit*random.random() for limit in ODF['limit']]
bins = EulersAsBins(Eulers,ODF['intervals'],deltas,ODF['center']) bins = EulersAsBins(Eulers,ODF['interval'],ODF['delta'],ODF['center'])
bin = binsAsBin(bins,ODF['intervals']) bin = binsAsBin(bins,ODF['interval'])
orientations[j] = '%g\t%g\t%g' %( math.degrees(Eulers[0]),math.degrees(Eulers[1]),math.degrees(Eulers[2]) ) orientations[j] = '%g\t%g\t%g' %( math.degrees(Eulers[0]),math.degrees(Eulers[1]),math.degrees(Eulers[2]) )
reconstructedODF[bin] += unitInc reconstructedODF[bin] += unitInc
@ -150,10 +141,8 @@ def MonteCarloBins (ODF,nSamples):
countMC = 0 countMC = 0
maxdV_V = max(ODF['dV_V']) maxdV_V = max(ODF['dV_V'])
nBins = ODF['intervals'][0]*ODF['intervals'][1]*ODF['intervals'][2]
deltas = [limit/intervals for limit,intervals in zip(ODF['limits'],ODF['intervals'])]
orientations = [None]*nSamples orientations = [None]*nSamples
reconstructedODF = [0.0]*nBins reconstructedODF = [0.0]*ODF['nBins']
unitInc = 1.0/nSamples unitInc = 1.0/nSamples
for j in range(nSamples): for j in range(nSamples):
@ -162,8 +151,8 @@ def MonteCarloBins (ODF,nSamples):
while MC > ODF['dV_V'][bin]: while MC > ODF['dV_V'][bin]:
countMC += 1 countMC += 1
MC = maxdV_V*random.random() MC = maxdV_V*random.random()
bin = int(nBins * random.random()) bin = int(ODF['nBins'] * random.random())
Eulers = binAsEulers(bin,ODF['intervals'],deltas,ODF['center']) Eulers = binAsEulers(bin,ODF['interval'],ODF['delta'],ODF['center'])
orientations[j] = '%g\t%g\t%g' %( math.degrees(Eulers[0]),math.degrees(Eulers[1]),math.degrees(Eulers[2]) ) orientations[j] = '%g\t%g\t%g' %( math.degrees(Eulers[0]),math.degrees(Eulers[1]),math.degrees(Eulers[2]) )
reconstructedODF[bin] += unitInc reconstructedODF[bin] += unitInc
@ -173,10 +162,8 @@ def MonteCarloBins (ODF,nSamples):
def TothVanHoutteSTAT (ODF,nSamples): def TothVanHoutteSTAT (ODF,nSamples):
""" ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians) """ """ ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians) """
nBins = ODF['intervals'][0]*ODF['intervals'][1]*ODF['intervals'][2]
deltas = [limit/intervals for limit,intervals in zip(ODF['limits'],ODF['intervals'])]
orientations = [None]*nSamples orientations = [None]*nSamples
reconstructedODF = [0.0]*nBins reconstructedODF = [0.0]*ODF['nBins']
unitInc = 1.0/nSamples unitInc = 1.0/nSamples
selectors = [random.random() for i in range(nSamples)] selectors = [random.random() for i in range(nSamples)]
@ -186,10 +173,10 @@ def TothVanHoutteSTAT (ODF,nSamples):
cumdV_V = 0.0 cumdV_V = 0.0
countSamples = 0 countSamples = 0
for bin in range(nBins) : for bin in range(ODF['nBins']) :
cumdV_V += ODF['dV_V'][bin] cumdV_V += ODF['dV_V'][bin]
while indexSelector < nSamples and selectors[indexSelector] < cumdV_V: while indexSelector < nSamples and selectors[indexSelector] < cumdV_V:
Eulers = binAsEulers(bin,ODF['intervals'],deltas,ODF['center']) Eulers = binAsEulers(bin,ODF['interval'],ODF['delta'],ODF['center'])
orientations[countSamples] = '%g\t%g\t%g' %( math.degrees(Eulers[0]),math.degrees(Eulers[1]),math.degrees(Eulers[2]) ) orientations[countSamples] = '%g\t%g\t%g' %( math.degrees(Eulers[0]),math.degrees(Eulers[1]),math.degrees(Eulers[2]) )
reconstructedODF[bin] += unitInc reconstructedODF[bin] += unitInc
countSamples += 1 countSamples += 1
@ -203,188 +190,176 @@ def TothVanHoutteSTAT (ODF,nSamples):
# -------------------------------------------------------------------- # --------------------------------------------------------------------
# MAIN # MAIN
# -------------------------------------------------------------------- # --------------------------------------------------------------------
identifiers = {
'limit': ['phi1','phi','phi2'],
'delta': ['phi1','phi','phi2'],
}
mappings = {
'limit': lambda x: math.radians(float(x)),
'delta': lambda x: math.radians(float(x)),
'origin': lambda x: str(x),
}
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
Transform linear binned data into Euler angles. Transform linear binned data into Euler angles.
""", version=string.replace(scriptID,'\n','\\n')
)
parser.add_option('-f', '--file', dest='file', type='string', metavar = 'string', \ """, version = scriptID)
help='file name, the input file is generated by the script "OIMlinear2linearODF.py"')
parser.add_option('-o', '--output', dest='output', type='string', metavar = 'string', \ parser.add_option('-n', '--number', dest='number', type='int', metavar = 'int',
help='the prefix of output files name.')
parser.add_option('-n', '--number', dest='number', type='int', metavar = 'int', \
help='the number of orientations needed to be generated, the default is [%default]') help='the number of orientations needed to be generated, the default is [%default]')
parser.add_option('-a','--algorithm', dest='algorithm', type='string', metavar = 'string', \ parser.add_option('-a','--algorithm', dest='algorithm', type='string', metavar = 'string',
help='''The algorithm adopted, three algorithms are provided, help='''The algorithm adopted, three algorithms are provided,
that is: that is:
[IA]: direct inversion, [IA]: direct inversion,
[STAT]: Van Houtte, [STAT]: Van Houtte,
[MC]: Monte Carlo. [MC]: Monte Carlo.
the default is [%default].''') the default is [%default].''') #make (multiple) choice
parser.add_option('-p','--phase', dest='phase', type='int', metavar = 'int', parser.add_option('-p','--phase', dest='phase', type='int', metavar = 'int',
help='phase index to be used [%default]') help='phase index to be used [%default]')
parser.add_option('--crystallite', dest='crystallite', type='int', metavar = 'int', parser.add_option('--crystallite', dest='crystallite', type='int', metavar = 'int',
help='crystallite index to be used [%default]') help='crystallite index to be used [%default]')
parser.add_option('-c', '--configuration', dest='config', action='store_true',
help='output material configuration [%default]')
parser.set_defaults(number = 500) parser.set_defaults(number = 500)
parser.set_defaults(output = 'texture')
parser.set_defaults(algorithm = 'IA') parser.set_defaults(algorithm = 'IA')
parser.set_defaults(phase = 1) parser.set_defaults(phase = 1)
parser.set_defaults(crystallite = 1) parser.set_defaults(crystallite = 1)
(options,filenames) = parser.parse_args()
options = parser.parse_args()[0]
# check usage
if not os.path.exists(options.file):
parser.error('binnedODF file does not exist'); sys.exit()
nameBinnedODF = options.file
nameSampledODF = options.output
nSamples = options.number nSamples = options.number
methods = [options.algorithm] methods = [options.algorithm]
#open binned ODF #--- setup file handles ---------------------------------------------------------------------------
try: files = []
fileBinnedODF = open(nameBinnedODF,'r') if filenames == []:
except: files.append({'name':'STDIN',
print 'unable to open binnedODF:', nameBinnedODF; 'input':sys.stdin,
sys.exit(1); 'output':sys.stdout,
'croak':sys.stderr,
# process header info })
ODF = {}
ODF['limits'] = [math.radians(float(limit)) for limit in fileBinnedODF.readline().split()[0:3]]
ODF['deltas'] = [math.radians(float(delta)) for delta in fileBinnedODF.readline().split()[0:3]]
#ODF['intervals'] = [int(interval) for interval in fileBinnedODF.readline().split()[0:3]]
ODF['intervals'] = [int(round(limit/delta)) for limit,delta in zip(ODF['limits'],ODF['deltas'])]
nBins = ODF['intervals'][0]*ODF['intervals'][1]*ODF['intervals'][2]
print 'Limit: ', [math.degrees(limit) for limit in ODF['limits']]
print 'Delta: ', [math.degrees(delta) for delta in ODF['deltas']]
print 'Interval: ', [interval + 1 for interval in ODF['intervals']]
centering = fileBinnedODF.readline()
if re.search('cell',centering.lower()):
ODF['center'] = 0.5
print 'cell-centered data (offset %g)'%ODF['center']
else: else:
ODF['center'] = 0.0 for name in filenames:
print 'vertex-centered data (offset %g)'%ODF['center'] if os.path.exists(name):
files.append({'name':name,
'input':open(name),
'output':open(name+'_tmp','w'),
'croak':sys.stdout,
})
fileBinnedODF.readline() # skip blank delimiter #--- loop over input files ------------------------------------------------------------------------
for file in files:
file['croak'].write('\033[1m' + scriptName + '\033[0m: ' + (file['name'] if file['name'] != 'STDIN' else '') + '\n')
# read linear binned data ODF = {
linesBinnedODF = fileBinnedODF.readlines() 'limit': np.empty(3,'d'),
fileBinnedODF.close() 'delta': np.empty(3,'d'),
'interval':np.empty(3,'i'),
'origin': ''
}
if len(linesBinnedODF) != nBins: table = damask.ASCIItable(file['input'],file['output'],buffered = False)
print 'expecting', nBins, 'values but got', len(linesBinnedODF) table.head_read()
sys.exit(1)
# build binnedODF array for header in table.info:
print 'reading',nBins,'values' headitems = map(str.lower,header.split())
sumdV_V = 0.0 if len(headitems) == 0: continue
ODF['dV_V'] = [None]*nBins if headitems[0] in mappings.keys():
ODF['nNonZero'] = 0 if headitems[0] in identifiers.keys():
dg = ODF['deltas'][0]*2*math.sin(ODF['deltas'][1]/2.0)*ODF['deltas'][2] for i in xrange(len(identifiers[headitems[0]])):
for bin in range(nBins): ODF[headitems[0]][i] = \
ODF['dV_V'][bin] = \ mappings[headitems[0]](headitems[headitems.index(identifiers[headitems[0]][i])+1])
max(0.0,float(linesBinnedODF[bin])) * dg * \ else:
math.sin(((bin//ODF['intervals'][2])%ODF['intervals'][1]+ODF['center'])*ODF['deltas'][1]) ODF[headitems[0]] = mappings[headitems[0]](headitems[1])
if ODF['dV_V'][bin] > 0.0:
sumdV_V += ODF['dV_V'][bin]
ODF['nNonZero'] += 1
for bin in range(nBins): ODF['dV_V'][bin] /= sumdV_V # normalize dV/V ODF['interval'] = np.array([int(round(limit/delta)) for limit,delta in zip(ODF['limit'],ODF['delta'])])
ODF['nBins'] = ODF['interval'].prod()
if re.search('boundary',ODF['origin'].lower()):
ODF['center'] = 0.5
else:
ODF['center'] = 0.0
binnedODF=table.data_readArray([table.labels.index('intensity')])
print 'non-zero fraction:', float(ODF['nNonZero'])/nBins,'(%i/%i)'%(ODF['nNonZero'],nBins) if binnedODF[0] != ODF['nBins']:
print 'Volume integral of ODF:', sumdV_V print 'expecting', ODF['nBins'], 'values but got', len(linesBinnedODF)
print 'Reference Integral:', ODF['limits'][0]*ODF['limits'][2]*(1-math.cos(ODF['limits'][1])) sys.exit(1)
# call methods # build binnedODF array
Functions = {'IA': 'directInversion', 'STAT': 'TothVanHoutteSTAT', 'MC': 'MonteCarloBins'} sumdV_V = 0.0
Orientations = {} ODF['dV_V'] = [None]*ODF['nBins']
ReconstructedODF = {} ODF['nNonZero'] = 0
for method in methods: dg = ODF['delta'][0]*2*math.sin(ODF['delta'][1]/2.0)*ODF['delta'][2]
Orientations[method], ReconstructedODF[method] = (globals()[Functions[method]])(ODF,nSamples) for bin in range(ODF['nBins']):
ODF['dV_V'][bin] = \
# calculate accuracy of sample max(0.0,table.data[bin,0]) * dg * \
squaredDiff = {} math.sin(((bin//ODF['interval'][2])%ODF['interval'][1]+ODF['center'])*ODF['delta'][1])
squaredRelDiff = {}
mutualProd = {}
indivSum = {}
indivSquaredSum = {}
for method in ['orig']+methods:
squaredDiff[method] = 0.0
squaredRelDiff[method] = 0.0
mutualProd[method] = 0.0
indivSum[method] = 0.0
indivSquaredSum[method] = 0.0
for bin in range(nBins):
for method in methods:
squaredDiff[method] += (ODF['dV_V'][bin] - ReconstructedODF[method][bin])**2
if ODF['dV_V'][bin] > 0.0: if ODF['dV_V'][bin] > 0.0:
squaredRelDiff[method] += (ODF['dV_V'][bin] - ReconstructedODF[method][bin])**2/ODF['dV_V'][bin]**2 sumdV_V += ODF['dV_V'][bin]
mutualProd[method] += ODF['dV_V'][bin]*ReconstructedODF[method][bin] ODF['nNonZero'] += 1
indivSum[method] += ReconstructedODF[method][bin]
indivSquaredSum[method] += ReconstructedODF[method][bin]**2 for bin in range(ODF['nBins']): ODF['dV_V'][bin] /= sumdV_V # normalize dV/V
print 'non-zero fraction:', float(ODF['nNonZero'])/ODF['nBins'],'(%i/%i)'%(ODF['nNonZero'],ODF['nBins'])
print 'Volume integral of ODF:', sumdV_V
print 'Reference Integral:', ODF['limit'][0]*ODF['limit'][2]*(1-math.cos(ODF['limit'][1]))
# call methods
Functions = {'IA': 'directInversion', 'STAT': 'TothVanHoutteSTAT', 'MC': 'MonteCarloBins'}
method = Functions[options.algorithm]
Orientations, ReconstructedODF = (globals()[method])(ODF,nSamples)
# calculate accuracy of sample
squaredDiff = {'orig':0.0,method:0.0}
squaredRelDiff = {'orig':0.0,method:0.0}
mutualProd = {'orig':0.0,method:0.0}
indivSum = {'orig':0.0,method:0.0}
indivSquaredSum = {'orig':0.0,method:0.0}
for bin in range(ODF['nBins']):
squaredDiff[method] += (ODF['dV_V'][bin] - ReconstructedODF[bin])**2
if ODF['dV_V'][bin] > 0.0:
squaredRelDiff[method] += (ODF['dV_V'][bin] - ReconstructedODF[bin])**2/ODF['dV_V'][bin]**2
mutualProd[method] += ODF['dV_V'][bin]*ReconstructedODF[bin]
indivSum[method] += ReconstructedODF[bin]
indivSquaredSum[method] += ReconstructedODF[bin]**2
indivSum['orig'] += ODF['dV_V'][bin] indivSum['orig'] += ODF['dV_V'][bin]
indivSquaredSum['orig'] += ODF['dV_V'][bin]**2 indivSquaredSum['orig'] += ODF['dV_V'][bin]**2
print 'sqrt(N*)RMSD of ODFs:\t', [math.sqrt(nSamples*squaredDiff[method]) for method in methods] print 'sqrt(N*)RMSD of ODFs:\t', math.sqrt(nSamples*squaredDiff[method])
print 'RMSrD of ODFs:\t', [math.sqrt(squaredRelDiff[method]) for method in methods] print 'RMSrD of ODFs:\t', math.sqrt(squaredRelDiff[method])
print 'rMSD of ODFs:\t', [squaredDiff[method]/indivSquaredSum['orig'] for method in methods] print 'rMSD of ODFs:\t', squaredDiff[method]/indivSquaredSum['orig']
#print 'correlation slope:\t', [(nBins*mutualProd[method]-indivSum['orig']*indivSum[method])/(nBins*indivSquaredSum['orig']-indivSum['orig']**2) for method in ('IA','STAT','MC')] print 'nNonZero correlation slope:\t', (ODF['nNonZero']*mutualProd[method]-indivSum['orig']*indivSum[method])/\
#print 'correlation confidence:\t', [(mutualProd[method]-indivSum['orig']*indivSum[method]/nBins)/\ (ODF['nNonZero']*indivSquaredSum['orig']-indivSum['orig']**2)
# (nBins*math.sqrt((indivSquaredSum['orig']/nBins-(indivSum['orig']/nBins)**2)*(indivSquaredSum[method]/nBins-(indivSum[method]/nBins)**2))) for method in ('IA','STAT','MC')] print 'nNonZero correlation confidence:\t',\
print 'nNonZero correlation slope:\t', [(ODF['nNonZero']*mutualProd[method]-indivSum['orig']*indivSum[method])/(ODF['nNonZero']*indivSquaredSum['orig']-indivSum['orig']**2) for method in methods] (mutualProd[method]-indivSum['orig']*indivSum[method]/ODF['nNonZero'])/\
print 'nNonZero correlation confidence:\t', [(mutualProd[method]-indivSum['orig']*indivSum[method]/ODF['nNonZero'])/\ (ODF['nNonZero']*math.sqrt((indivSquaredSum['orig']/ODF['nNonZero']-(indivSum['orig']/ODF['nNonZero'])**2)*\
(ODF['nNonZero']*math.sqrt((indivSquaredSum['orig']/ODF['nNonZero']-(indivSum['orig']/ODF['nNonZero'])**2)*(indivSquaredSum[method]/ODF['nNonZero']-(indivSum[method]/ODF['nNonZero'])**2))) for method in methods] (indivSquaredSum[method]/ODF['nNonZero']-(indivSum[method]/ODF['nNonZero'])**2)))
for method in methods:
if method == 'IA' and nSamples < ODF['nNonZero']: if method == 'IA' and nSamples < ODF['nNonZero']:
strOpt = '(%i)'%ODF['nNonZero'] strOpt = '(%i)'%ODF['nNonZero']
else:
strOpt = ''
try:
fileSampledODF = open(nameSampledODF+'.'+method+'sampled_'+str(nSamples)+strOpt, 'w')
fileSampledODF.write('%i\n'%nSamples)
fileSampledODF.write('\n'.join(Orientations[method])+'\n')
fileSampledODF.close()
except:
print 'unable to write sampledODF:', nameSampledODF+'.'+method+'sampled_'+str(nSamples)+strOpt
try:
fileRegressionODF = open(nameSampledODF+'.'+method+'regression_'+str(nSamples)+strOpt, 'w')
fileRegressionODF.write('\n'.join([a+'\t'+b for (a,b) in zip(map(str,ReconstructedODF[method]),map(str,ODF['dV_V']))])+'\n')
fileRegressionODF.close()
except:
print 'unable to write RegressionODF:', nameSampledODF+'.'+method+'regression_'+str(nSamples)+strOpt
if options.config: formatwidth = 1
try: file['output'].write('#-------------------#')
fileConfig = open(nameSampledODF+'.'+method+str(nSamples)+'.config', 'w') # write config file file['output'].write('\n<microstructure>\n')
formatwidth = 1 file['output'].write('#-------------------#\n')
fileConfig.write('#-------------------#')
fileConfig.write('\n<microstructure>\n') for i,ID in enumerate(xrange(nSamples)):
fileConfig.write('#-------------------#\n') file['output'].write('[Grain%s]\n'%(str(ID+1).zfill(formatwidth)) + \
'crystallite %i\n'%options.crystallite + \
for i,ID in enumerate(xrange(nSamples)): '(constituent) phase %i texture %s fraction 1.0\n'%(options.phase,str(ID+1).rjust(formatwidth)))
fileConfig.write('[Grain%s]\n'%(str(ID+1).zfill(formatwidth)) + \
'crystallite %i\n'%options.crystallite + \ file['output'].write('\n#-------------------#')
'(constituent) phase %i texture %s fraction 1.0\n'%(options.phase,str(ID+1).rjust(formatwidth))) file['output'].write('\n<texture>\n')
file['output'].write('#-------------------#\n')
fileConfig.write('\n#-------------------#') for ID in xrange(nSamples):
fileConfig.write('\n<texture>\n') eulers = re.split(r'[\t]', Orientations[ID].strip())
fileConfig.write('#-------------------#\n')
for ID in xrange(nSamples): file['output'].write('[Grain%s]\n'%(str(ID+1).zfill(formatwidth)) + \
eulers = re.split(r'[\t]', Orientations[method][ID].strip()) '(gauss) phi1 %10.5f Phi %10.5f phi2 %10.6f scatter 0.0 fraction 1.0\n'\
%(float(eulers[0]),float(eulers[1]),float(eulers[2])))
fileConfig.write('[Grain%s]\n'%(str(ID+1).zfill(formatwidth)) + \ #--- output finalization --------------------------------------------------------------------------
'(gauss) phi1 %10.5f Phi %10.5f phi2 %10.6f scatter 0.0 fraction 1.0\n'\ if file['name'] != 'STDIN':
%(float(eulers[0]),float(eulers[1]),float(eulers[2]))) file['output'].close()
fileConfig.close() os.rename(file['name']+'_tmp',
except: os.path.splitext(file['name'])[0] +'_'+method+'%s'%('_material.config'))
print 'unable to write material.config file:', nameSampledODF+'.'+method+str(nSamples)+'.config'