copied some stcopied some scripts from the Code folder that are of interest for DAMASK
This commit is contained in:
Martin Diehl
parent
70eaa9f8b8
commit
f1df6cf40f
|
|
@ -0,0 +1,67 @@
|
|||
#!/usr/bin/env python
|
||||
|
||||
import os,sys,re
|
||||
|
||||
try:
|
||||
inName = sys.argv[1]
|
||||
outName = os.path.splitext(inName)[0]+'.linearODF'
|
||||
nPhi1,nPHI,nPhi2 = map(int,sys.argv[2:5])
|
||||
except:
|
||||
print "\nusage:",sys.argv[0],"file nPhi1 nPHI nPhi2\n"
|
||||
sys.exit(1)
|
||||
|
||||
N = (nPhi1-1)*(nPHI-1)*(nPhi2-1)
|
||||
|
||||
try:
|
||||
inFile = open(inName,'r')
|
||||
content = inFile.readlines()
|
||||
except:
|
||||
print 'unable to read:',inName
|
||||
sys.exit(1)
|
||||
try:
|
||||
outFile = open(outName,'w')
|
||||
except:
|
||||
print 'unable to write:',outName
|
||||
sys.exit(1)
|
||||
|
||||
ODF = [[[[None] for k in range(nPhi2)] for j in range(nPHI)] for i in range(nPhi1)]
|
||||
linear = [None]*N
|
||||
line = 0
|
||||
|
||||
while (content[line].startswith('#')): # skip comments at start of file
|
||||
line += 1
|
||||
|
||||
for iPhi1 in range(nPhi1):
|
||||
for iPHI in range(nPHI):
|
||||
for iPhi2 in range(nPhi2):
|
||||
words = content[line].split()
|
||||
ODF[iPhi1][iPHI][iPhi2] = float(words[3]) # extract intensity (in column 4)
|
||||
line += 1
|
||||
|
||||
for iPhi1 in range(nPhi1-1):
|
||||
for iPHI in range(nPHI-1):
|
||||
for iPhi2 in range(nPhi2-1):
|
||||
linear[iPhi1*(nPHI-1)*(nPhi2-1)+iPHI*(nPhi2-1)+iPhi2] = (\
|
||||
ODF[iPhi1 ][iPHI ][iPhi2 ] + \
|
||||
ODF[iPhi1 ][iPHI ][iPhi2+1] + \
|
||||
ODF[iPhi1 ][iPHI+1][iPhi2 ] + \
|
||||
ODF[iPhi1 ][iPHI+1][iPhi2+1] + \
|
||||
ODF[iPhi1+1][iPHI ][iPhi2 ] + \
|
||||
ODF[iPhi1+1][iPHI ][iPhi2+1] + \
|
||||
ODF[iPhi1+1][iPHI+1][iPhi2 ] + \
|
||||
ODF[iPhi1+1][iPHI+1][iPhi2+1] \
|
||||
) / 8.0
|
||||
|
||||
|
||||
inFile.close()
|
||||
|
||||
outFile.write('rangePhi1?\trangePHI?\trangePhi2?\n')
|
||||
outFile.write('%i\t%i\t%i needs to be converted to angular steps\n'%(nPhi1-1,nPHI-1,nPhi2-1))
|
||||
outFile.write('cell-centered data\n')
|
||||
outFile.write('\n')
|
||||
|
||||
for i in range(N):
|
||||
outFile.write('%g\n'%(linear[i]))
|
||||
|
||||
outFile.close()
|
||||
|
||||
|
|
@ -0,0 +1,105 @@
|
|||
#!/usr/bin/env python
|
||||
|
||||
import os,sys,math,re
|
||||
|
||||
# --- helper functions ---
|
||||
|
||||
def binAsBins(bin,intervals):
|
||||
""" explode compound bin into 3D bins list """
|
||||
bins = [0]*3
|
||||
bins[0] = (bin//(intervals[1] * intervals[2])) % intervals[0]
|
||||
bins[1] = (bin//intervals[2]) % intervals[1]
|
||||
bins[2] = bin % intervals[2]
|
||||
return bins
|
||||
|
||||
def binsAsBin(bins,intervals):
|
||||
""" implode 3D bins into compound bin """
|
||||
return (bins[0]*intervals[1] + bins[1])*intervals[2] + bins[2]
|
||||
|
||||
def EulersAsBins(Eulers,intervals,deltas,center):
|
||||
""" return list of Eulers translated into 3D bins list """
|
||||
return [\
|
||||
int((euler+(0.5-center)*delta)//delta)%interval \
|
||||
for euler,delta,interval in zip(Eulers,deltas,intervals) \
|
||||
]
|
||||
|
||||
def binAsEulers(bin,intervals,deltas,center):
|
||||
""" compound bin number translated into list of Eulers """
|
||||
Eulers = [0.0]*3
|
||||
Eulers[2] = (bin%intervals[2] + center)*deltas[2]
|
||||
Eulers[1] = (bin//intervals[2]%intervals[1] + center)*deltas[1]
|
||||
Eulers[0] = (bin//(intervals[2]*intervals[1]) + center)*deltas[0]
|
||||
return Eulers
|
||||
|
||||
|
||||
|
||||
# check usage
|
||||
try:
|
||||
inName = sys.argv[1]
|
||||
outLimits = sys.argv[2:5]
|
||||
except:
|
||||
print "usage:",sys.argv[0],"nameLinearODF limitPhi1 limitPHI limitPhi2"
|
||||
sys.exit(1);
|
||||
|
||||
|
||||
#open binned ODF
|
||||
try:
|
||||
inFile = open(inName,'r')
|
||||
except:
|
||||
print 'unable to open binnedODF:', inName;
|
||||
sys.exit(1);
|
||||
|
||||
# process header info
|
||||
ODF = {}
|
||||
inLimits = [math.radians(int(float(limit))) for limit in inFile.readline().split()]
|
||||
outLimits = [math.radians(int(float(limit))) for limit in outLimits]
|
||||
deltas = [math.radians(float(delta)) for delta in inFile.readline().split()]
|
||||
inIntervals = [int(limit/delta) for limit,delta in zip(inLimits,deltas)]
|
||||
outIntervals = [int(limit/delta) for limit,delta in zip(outLimits,deltas)]
|
||||
inBins = inIntervals[0]*inIntervals[1]*inIntervals[2]
|
||||
|
||||
print 'Limit:', [math.degrees(limit) for limit in inLimits]
|
||||
print 'Crop:', [math.degrees(limit) for limit in outLimits]
|
||||
print 'Delta:', [math.degrees(delta) for delta in deltas]
|
||||
print 'Interval:', inIntervals
|
||||
print 'Interval:', outIntervals
|
||||
|
||||
centering = inFile.readline()
|
||||
if re.search('cell',centering.lower()):
|
||||
ODF['center'] = 0.5
|
||||
print 'cell-centered data (offset %g)'%ODF['center']
|
||||
else:
|
||||
ODF['center'] = 0.0
|
||||
print 'vertex-centered data (offset %g)'%ODF['center']
|
||||
|
||||
inFile.readline() # skip blank delimiter
|
||||
|
||||
# read linear binned data
|
||||
inODF = map(float,inFile.readlines())
|
||||
inFile.close()
|
||||
|
||||
if len(inODF) != inBins:
|
||||
print 'expecting', inBins, 'values but got', len(inODF)
|
||||
sys.exit(1)
|
||||
|
||||
try:
|
||||
outName = os.path.splitext(inName)[0]+'_%ix%ix%i'%(outIntervals[0],outIntervals[1],outIntervals[2])+'.linearODF'
|
||||
outFile = open(outName,'w')
|
||||
except:
|
||||
print 'unable to write:',outName
|
||||
sys.exit(1)
|
||||
|
||||
outFile.write('%g\t%g\t%g\n'%(\
|
||||
math.degrees(outIntervals[0]*deltas[0]),\
|
||||
math.degrees(outIntervals[1]*deltas[1]),\
|
||||
math.degrees(outIntervals[2]*deltas[2]) ))
|
||||
outFile.write('%i\t%i\t%i\n'%(math.degrees(deltas[0]),math.degrees(deltas[1]),math.degrees(deltas[2])))
|
||||
outFile.write('%s-centered data\n'%{True:'vertex',False:'cell'}[ODF['center']==0.0])
|
||||
outFile.write('\n')
|
||||
|
||||
for phi1 in range(outIntervals[0]):
|
||||
for Phi in range(outIntervals[1]):
|
||||
for phi2 in range(outIntervals[2]):
|
||||
outFile.write('%g\n'%(inODF[((phi1*inIntervals[1])+Phi)*inIntervals[2]+phi2]))
|
||||
|
||||
outFile.close()
|
||||
|
|
@ -0,0 +1,331 @@
|
|||
#!/usr/bin/env python
|
||||
|
||||
import os,sys,math,re,random
|
||||
random.seed()
|
||||
|
||||
|
||||
# --- helper functions ---
|
||||
|
||||
def binAsBins(bin,intervals):
|
||||
""" explode compound bin into 3D bins list """
|
||||
bins = [0]*3
|
||||
bins[0] = (bin//(intervals[1] * intervals[2])) % intervals[0]
|
||||
bins[1] = (bin//intervals[2]) % intervals[1]
|
||||
bins[2] = bin % intervals[2]
|
||||
return bins
|
||||
|
||||
def binsAsBin(bins,intervals):
|
||||
""" implode 3D bins into compound bin """
|
||||
return (bins[0]*intervals[1] + bins[1])*intervals[2] + bins[2]
|
||||
|
||||
def EulersAsBins(Eulers,intervals,deltas,center):
|
||||
""" return list of Eulers translated into 3D bins list """
|
||||
return [\
|
||||
int((euler+(0.5-center)*delta)//delta)%interval \
|
||||
for euler,delta,interval in zip(Eulers,deltas,intervals) \
|
||||
]
|
||||
|
||||
def binAsEulers(bin,intervals,deltas,center):
|
||||
""" compound bin number translated into list of Eulers """
|
||||
Eulers = [0.0]*3
|
||||
Eulers[2] = (bin%intervals[2] + center)*deltas[2]
|
||||
Eulers[1] = (bin//intervals[2]%intervals[1] + center)*deltas[1]
|
||||
Eulers[0] = (bin//(intervals[2]*intervals[1]) + center)*deltas[0]
|
||||
return Eulers
|
||||
|
||||
def directInvRepetitions(probability,scale):
|
||||
""" calculate number of samples drawn by direct inversion """
|
||||
nDirectInv = 0
|
||||
for bin in range(len(probability)): # loop over bins
|
||||
nDirectInv += int(round(probability[bin]*scale)) # calc repetition
|
||||
return nDirectInv
|
||||
|
||||
|
||||
# --- sampling methods ---
|
||||
|
||||
# ----- efficient algorithm ---------
|
||||
|
||||
def directInversion (ODF,nSamples):
|
||||
""" 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'])]
|
||||
|
||||
# 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
|
||||
repetition = [None]*nBins
|
||||
probabilityScale = nOptSamples # guess
|
||||
|
||||
scaleLower = 0.0
|
||||
nInvSamplesLower = 0
|
||||
scaleUpper = float(nOptSamples)
|
||||
incFactor = 1.0
|
||||
nIter = 0
|
||||
nInvSamplesUpper = directInvRepetitions(ODF['dV_V'],scaleUpper)
|
||||
while (\
|
||||
(scaleUpper-scaleLower > scaleUpper*1e-15 or nInvSamplesUpper < nOptSamples) and \
|
||||
nInvSamplesUpper != nOptSamples \
|
||||
): # closer match required?
|
||||
if nInvSamplesUpper < nOptSamples:
|
||||
scaleLower,scaleUpper = scaleUpper,scaleUpper+incFactor*(scaleUpper-scaleLower)/2.0
|
||||
incFactor *= 2.0
|
||||
nInvSamplesLower,nInvSamplesUpper = nInvSamplesUpper,directInvRepetitions(ODF['dV_V'],scaleUpper)
|
||||
else:
|
||||
scaleUpper = (scaleLower+scaleUpper)/2.0
|
||||
incFactor = 1.0
|
||||
nInvSamplesUpper = directInvRepetitions(ODF['dV_V'],scaleUpper)
|
||||
nIter += 1
|
||||
print '%i:(%12.11f,%12.11f) %i <= %i <= %i'%(nIter,scaleLower,scaleUpper,nInvSamplesLower,nOptSamples,nInvSamplesUpper)
|
||||
nInvSamples = nInvSamplesUpper
|
||||
scale = scaleUpper
|
||||
print 'created set of',nInvSamples,'samples (',float(nInvSamples)/nOptSamples-1.0,') with scaling',scale,'delivering',nSamples
|
||||
repetition = [None]*nBins # preallocate and clear
|
||||
|
||||
for bin in range(nBins): # loop over bins
|
||||
repetition[bin] = int(round(ODF['dV_V'][bin]*scale)) # calc repetition
|
||||
|
||||
# build set
|
||||
set = [None]*nInvSamples
|
||||
i = 0
|
||||
for bin in range(nBins):
|
||||
set[i:i+repetition[bin]] = [bin]*repetition[bin] # fill set with bin, i.e. orientation
|
||||
i += repetition[bin] # advance set counter
|
||||
|
||||
orientations = [None]*nSamples
|
||||
reconstructedODF = [0.0]*nBins
|
||||
unitInc = 1.0/nSamples
|
||||
for j in range(nSamples):
|
||||
if (j == nInvSamples-1): ex = j
|
||||
else: ex = int(round(random.uniform(j+0.5,nInvSamples-0.5)))
|
||||
bin = set[ex]
|
||||
bins = binAsBins(bin,ODF['intervals'])
|
||||
Eulers = binAsEulers(bin,ODF['intervals'],deltas,ODF['center'])
|
||||
orientations[j] = '%g\t%g\t%g' %( math.degrees(Eulers[0]),math.degrees(Eulers[1]),math.degrees(Eulers[2]) )
|
||||
reconstructedODF[bin] += unitInc
|
||||
set[ex] = set[j] # exchange orientations
|
||||
|
||||
return orientations, reconstructedODF
|
||||
|
||||
|
||||
# ----- trial and error algorithms ---------
|
||||
|
||||
def MonteCarloEulers (ODF,nSamples):
|
||||
""" ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians) """
|
||||
|
||||
countMC = 0
|
||||
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
|
||||
reconstructedODF = [0.0]*nBins
|
||||
unitInc = 1.0/nSamples
|
||||
|
||||
for j in range(nSamples):
|
||||
MC = maxdV_V*2.0
|
||||
bin = 0
|
||||
while MC > ODF['dV_V'][bin]:
|
||||
countMC += 1
|
||||
MC = maxdV_V*random.random()
|
||||
Eulers = [limit*random.random() for limit in ODF['limits']]
|
||||
bins = EulersAsBins(Eulers,ODF['intervals'],deltas,ODF['center'])
|
||||
bin = binsAsBin(bins,ODF['intervals'])
|
||||
orientations[j] = '%g\t%g\t%g' %( math.degrees(Eulers[0]),math.degrees(Eulers[1]),math.degrees(Eulers[2]) )
|
||||
reconstructedODF[bin] += unitInc
|
||||
|
||||
return orientations, reconstructedODF, countMC
|
||||
|
||||
|
||||
def MonteCarloBins (ODF,nSamples):
|
||||
""" ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians) """
|
||||
|
||||
countMC = 0
|
||||
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
|
||||
reconstructedODF = [0.0]*nBins
|
||||
unitInc = 1.0/nSamples
|
||||
|
||||
for j in range(nSamples):
|
||||
MC = maxdV_V*2.0
|
||||
bin = 0
|
||||
while MC > ODF['dV_V'][bin]:
|
||||
countMC += 1
|
||||
MC = maxdV_V*random.random()
|
||||
bin = int(nBins * random.random())
|
||||
Eulers = binAsEulers(bin,ODF['intervals'],deltas,ODF['center'])
|
||||
orientations[j] = '%g\t%g\t%g' %( math.degrees(Eulers[0]),math.degrees(Eulers[1]),math.degrees(Eulers[2]) )
|
||||
reconstructedODF[bin] += unitInc
|
||||
|
||||
return orientations, reconstructedODF
|
||||
|
||||
|
||||
def TothVanHoutteSTAT (ODF,nSamples):
|
||||
""" 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
|
||||
reconstructedODF = [0.0]*nBins
|
||||
unitInc = 1.0/nSamples
|
||||
|
||||
selectors = [random.random() for i in range(nSamples)]
|
||||
selectors.sort()
|
||||
indexSelector = 0
|
||||
|
||||
cumdV_V = 0.0
|
||||
countSamples = 0
|
||||
|
||||
for bin in range(nBins) :
|
||||
cumdV_V += ODF['dV_V'][bin]
|
||||
while indexSelector < nSamples and selectors[indexSelector] < cumdV_V:
|
||||
Eulers = binAsEulers(bin,ODF['intervals'],deltas,ODF['center'])
|
||||
orientations[countSamples] = '%g\t%g\t%g' %( math.degrees(Eulers[0]),math.degrees(Eulers[1]),math.degrees(Eulers[2]) )
|
||||
reconstructedODF[bin] += unitInc
|
||||
countSamples += 1
|
||||
indexSelector += 1
|
||||
|
||||
print 'created set of',countSamples,'when asked to deliver',nSamples
|
||||
|
||||
return orientations, reconstructedODF
|
||||
|
||||
|
||||
|
||||
|
||||
# check usage
|
||||
try:
|
||||
nameBinnedODF = sys.argv[1]
|
||||
nSamples = int(float(sys.argv[2]))
|
||||
except:
|
||||
print "\nusage:",os.path.basename(sys.argv[0]),"nameLinearODF nSamples [nameSampledODF]\n"
|
||||
sys.exit(1);
|
||||
|
||||
methods = ['IA','STAT']
|
||||
|
||||
#open binned ODF
|
||||
try:
|
||||
fileBinnedODF = open(nameBinnedODF,'r')
|
||||
except:
|
||||
print 'unable to open binnedODF:', nameBinnedODF;
|
||||
sys.exit(1);
|
||||
|
||||
# process header info
|
||||
ODF = {}
|
||||
ODF['limits'] = [math.radians(float(limit)) for limit in fileBinnedODF.readline().split()]
|
||||
ODF['deltas'] = [math.radians(float(delta)) for delta in fileBinnedODF.readline().split()]
|
||||
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:', ODF['intervals']
|
||||
|
||||
centering = fileBinnedODF.readline()
|
||||
if re.search('cell',centering.lower()):
|
||||
ODF['center'] = 0.5
|
||||
print 'cell-centered data (offset %g)'%ODF['center']
|
||||
else:
|
||||
ODF['center'] = 0.0
|
||||
print 'vertex-centered data (offset %g)'%ODF['center']
|
||||
|
||||
fileBinnedODF.readline() # skip blank delimiter
|
||||
|
||||
# read linear binned data
|
||||
linesBinnedODF = fileBinnedODF.readlines()
|
||||
fileBinnedODF.close()
|
||||
|
||||
if len(linesBinnedODF) != nBins:
|
||||
print 'expecting', nBins, 'values but got', len(linesBinnedODF)
|
||||
sys.exit(1)
|
||||
|
||||
# build binnedODF array
|
||||
print 'reading',nBins,'values'
|
||||
sumdV_V = 0.0
|
||||
ODF['dV_V'] = [None]*nBins
|
||||
ODF['nNonZero'] = 0
|
||||
dg = ODF['deltas'][0]*2*math.sin(ODF['deltas'][1]/2.0)*ODF['deltas'][2]
|
||||
for bin in range(nBins):
|
||||
ODF['dV_V'][bin] = \
|
||||
max(0.0,float(linesBinnedODF[bin])) * dg * \
|
||||
math.sin(((bin//ODF['intervals'][2])%ODF['intervals'][1]+ODF['center'])*ODF['deltas'][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
|
||||
|
||||
print 'non-zero fraction:', float(ODF['nNonZero'])/nBins,'(%i/%i)'%(ODF['nNonZero'],nBins)
|
||||
print 'Volume integral of ODF:', sumdV_V
|
||||
print 'Reference Integral:', ODF['limits'][0]*ODF['limits'][2]*(1-math.cos(ODF['limits'][1]))
|
||||
|
||||
# call methods
|
||||
Functions = {'IA': 'directInversion', 'STAT': 'TothVanHoutteSTAT', 'MC': 'MonteCarloBins'}
|
||||
Orientations = {}
|
||||
ReconstructedODF = {}
|
||||
for method in methods:
|
||||
Orientations[method], ReconstructedODF[method] = (globals()[Functions[method]])(ODF,nSamples)
|
||||
|
||||
# calculate accuracy of sample
|
||||
squaredDiff = {}
|
||||
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:
|
||||
squaredRelDiff[method] += (ODF['dV_V'][bin] - ReconstructedODF[method][bin])**2/ODF['dV_V'][bin]**2
|
||||
mutualProd[method] += ODF['dV_V'][bin]*ReconstructedODF[method][bin]
|
||||
indivSum[method] += ReconstructedODF[method][bin]
|
||||
indivSquaredSum[method] += ReconstructedODF[method][bin]**2
|
||||
indivSum['orig'] += ODF['dV_V'][bin]
|
||||
indivSquaredSum['orig'] += ODF['dV_V'][bin]**2
|
||||
|
||||
print 'sqrt(N*)RMSD of ODFs:\t', [math.sqrt(nSamples*squaredDiff[method]) for method in methods]
|
||||
print 'RMSrD of ODFs:\t', [math.sqrt(squaredRelDiff[method]) for method in methods]
|
||||
print 'rMSD of ODFs:\t', [squaredDiff[method]/indivSquaredSum['orig'] for method in methods]
|
||||
#print 'correlation slope:\t', [(nBins*mutualProd[method]-indivSum['orig']*indivSum[method])/(nBins*indivSquaredSum['orig']-indivSum['orig']**2) for method in ('IA','STAT','MC')]
|
||||
#print 'correlation confidence:\t', [(mutualProd[method]-indivSum['orig']*indivSum[method]/nBins)/\
|
||||
# (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 slope:\t', [(ODF['nNonZero']*mutualProd[method]-indivSum['orig']*indivSum[method])/(ODF['nNonZero']*indivSquaredSum['orig']-indivSum['orig']**2) for method in methods]
|
||||
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)*(indivSquaredSum[method]/ODF['nNonZero']-(indivSum[method]/ODF['nNonZero'])**2))) for method in methods]
|
||||
|
||||
# write result
|
||||
try:
|
||||
nameSampledODF = sys.argv[3]
|
||||
except:
|
||||
sys.exit(0) # that's it folks
|
||||
|
||||
for method in methods:
|
||||
if method == 'IA' and nSamples < 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
|
||||
|
|
@ -0,0 +1,121 @@
|
|||
#!/usr/bin/env python
|
||||
|
||||
import os,sys,math,re
|
||||
from optparse import OptionParser
|
||||
|
||||
def integerFactorization(i):
|
||||
|
||||
j = int(math.floor(math.sqrt(float(i))))
|
||||
while (j>1 and int(i)%j != 0):
|
||||
j -= 1
|
||||
return j
|
||||
|
||||
def positiveRadians(angle):
|
||||
|
||||
angle = math.radians(float(angle))
|
||||
while angle < 0.0:
|
||||
angle += 2.0*math.pi
|
||||
|
||||
return angle
|
||||
|
||||
|
||||
def getHeader(sizeX,sizeY,step):
|
||||
|
||||
return [ \
|
||||
'# TEM_PIXperUM 1.000000', \
|
||||
'# x-star 0.509548', \
|
||||
'# y-star 0.795272', \
|
||||
'# z-star 0.611799', \
|
||||
'# WorkingDistance 18.000000', \
|
||||
'#', \
|
||||
'# Phase 1', \
|
||||
'# MaterialName Al', \
|
||||
'# Formula Fe', \
|
||||
'# Info', \
|
||||
'# Symmetry 43', \
|
||||
'# LatticeConstants 2.870 2.870 2.870 90.000 90.000 90.000', \
|
||||
'# NumberFamilies 4', \
|
||||
'# hklFamilies 1 1 0 1 0.000000 1', \
|
||||
'# hklFamilies 2 0 0 1 0.000000 1', \
|
||||
'# hklFamilies 2 1 1 1 0.000000 1', \
|
||||
'# hklFamilies 3 1 0 1 0.000000 1', \
|
||||
'# Categories 0 0 0 0 0 ', \
|
||||
'#', \
|
||||
'# GRID: SquareGrid', \
|
||||
'# XSTEP: ' + str(step), \
|
||||
'# YSTEP: ' + str(step), \
|
||||
'# NCOLS_ODD: ' + str(sizeX), \
|
||||
'# NCOLS_EVEN: ' + str(sizeX), \
|
||||
'# NROWS: ' + str(sizeY), \
|
||||
'#', \
|
||||
'# OPERATOR: ODFsammpling', \
|
||||
'#', \
|
||||
'# SAMPLEID: ', \
|
||||
'#', \
|
||||
'# SCANID: ', \
|
||||
'#', \
|
||||
]
|
||||
|
||||
|
||||
parser = OptionParser(usage='%prog [options] datafile(s)')
|
||||
parser.add_option("-c", "--column", type="int",\
|
||||
dest="column",\
|
||||
help="starting column of Euler triplet")
|
||||
parser.add_option("-s", "--skip", type="int",\
|
||||
dest="skip",\
|
||||
help="skip this many lines of heading info [%default]")
|
||||
|
||||
parser.set_defaults (column = 1)
|
||||
parser.set_defaults (skip = 0)
|
||||
|
||||
(options, files) = parser.parse_args()
|
||||
options.column -= 1
|
||||
|
||||
if files == []:
|
||||
parser.error('no input file specified...')
|
||||
sys.exit(1)
|
||||
if options.column < 0:
|
||||
parser.error('column needs to be 1,2,...')
|
||||
sys.exit(1)
|
||||
|
||||
while len(files) > 0:
|
||||
textureFilename = files.pop()
|
||||
baseName = os.path.splitext(textureFilename)[0]
|
||||
|
||||
# open texture file and read content
|
||||
textureFile = open(textureFilename)
|
||||
content = textureFile.readlines()
|
||||
textureFile.close()
|
||||
|
||||
m = re.match('(\d+)\s+head',content[0],re.I)
|
||||
if m != None and options.skip == 0:
|
||||
options.skip = int(m.group(1))+1
|
||||
|
||||
# extract orientation angles
|
||||
angles = [map(positiveRadians,line.split()[options.column:options.column+3]) for line in content[options.skip:]]
|
||||
|
||||
nPoints = len(angles)
|
||||
sizeY = integerFactorization(nPoints)
|
||||
sizeX = nPoints / sizeY
|
||||
print '%s: %i * %i = %i (== %i)'%(baseName,sizeX,sizeY,sizeX*sizeY,nPoints)
|
||||
# write ang file
|
||||
|
||||
try:
|
||||
|
||||
# write header
|
||||
angFile = open(baseName + '.ang','w')
|
||||
for line in getHeader(sizeX,sizeY,1.0):
|
||||
angFile.write(line + '\n')
|
||||
|
||||
# write data
|
||||
counter = 0
|
||||
for point in angles:
|
||||
angFile.write(''.join(['%10.5f'%angle for angle in point])+
|
||||
''.join(['%10.5f'%coord for coord in [counter%sizeX,counter//sizeX]])+
|
||||
' 100.0 1.0 0 1 1.0\n')
|
||||
counter += 1
|
||||
|
||||
angFile.close()
|
||||
|
||||
except:
|
||||
print 'unable to write',baseName
|
||||
Loading…
Reference in New Issue