<jupyter_start><jupyter_text>Munch MuseumSearch using experimental tabulated phases in DatabaseExp Created using MunchMuzeum Notebook.I found the hematite manualy and fitted its gamma.  Imports<jupyter_code>from XRDXRFutils import Phase,DatabaseXRD, DataXRD, SpectraXRD, GaussNewton, PhaseList, PhaseMap, PhaseSearch, PhaseMapSave
from XRDXRFutils import GammaMap,ChiMap

import os
import pickle


from joblib import Parallel, delayed
import h5py
from sklearn.linear_model import LinearRegression
from scipy.optimize import curve_fit, least_squares

from numpy import linspace,concatenate,sqrt,log,histogram,array
from matplotlib.pyplot import sca,vlines,show,fill_between,sca,legend,imshow,subplots,plot,xlim,ylim,xlabel,ylabel,cm,title,scatter,colorbar,figure,vlines
from sklearn.cluster import KMeans,MiniBatchKMeans

from multiprocessing import Pool

from PIL import Image

def f_linear(x,a,b):
    return a*x + b
def f_loss(x,t,y):
    return (x[0]*t + x[1]) - y<jupyter_output><empty_output><jupyter_text>Define Paths and Spectra Parameters<jupyter_code>path_xrd = '/home/shared/dataXRDXRF/MunchMuseum/M491/ProfiloXRD/'
path_database = '/home/shared/DatabaseXRD'

path_data = 'data/' # data of intermediate results, for fast loading
path_figures = 'figures/' # figures generated by the script
path_results = 'results/' # results generated by the script: raw data, tif maps

min_theta = 20
max_theta = 53
min_intensity = 0.1 # among the tabulated peaks, selects only the ones above this threshold of intensity (scale between 0 and 1)
first_n_peaks = None # selects the first n most intense peaks (if None, leaves all the peaks)

sigma = 0.2<jupyter_output><empty_output><jupyter_text>Read XRD Datafrom xrd.h5<jupyter_code>try:
    data = DataXRD().load_h5(path_xrd + 'xrd.h5')
except:
    print('Reading from raw data.')
    data = DataXRD().read_params(path_xrd + 'Scanning_Parameters.txt').read(path_xrd).calibrate_from_file(path_xrd + 'calibration.ini').remove_background(std = 5).save_h5(path_xrd + 'xrd.h5')
    
print("a: %.1f s: %.1f beta: %.3f"%(data.opt[0],data.opt[1],data.opt[2]))

figure(figsize=(6,4))
im = imshow(data.data.sum(axis=2))
show()<jupyter_output>Loading: /home/shared/dataXRDXRF/MunchMuseum/M491/ProfiloXRD/xrd.h5
a: -1214.7 s: 1952.7 beta: 51.663<jupyter_text>Read database Define PhasesThis is for simplification. Phases can be selected iteratively from database using 'Tab'<jupyter_code>database = DatabaseXRD().read_cifs(path_database)
databaseExp = DatabaseXRD().read_cifs('DatabaseExp/')

print('Phases in database:',len(database))
print('Phases in databaseEXP:',len(databaseExp))

lazurite = database['Lazurite'][0]
hydrocerussite = database['Hydrocerussite'][0]
cinnabar = database['Cinnabar'][1]
barite = database['Barite'][0]
spinel = database['Spinel'][0]
calcite = database['Calcite'][0]
hematite = database['Hematite'][4]

huntite = database['Huntite'][0]
as4 = database['As4 O6'][0]

orpiment = database['Orpiment'][0]
cuprorivaite = database['Cuprorivaite'][0]

hydrocerussite = databaseExp['Hydrocerussite'][0]
cinnabar = databaseExp['Cinnabar'][0]
lazurite = databaseExp['Lazurite'][0]
hematite = databaseExp['Hematite'][1]

phases_a_s = PhaseList([hydrocerussite,cinnabar,lazurite,hematite])

phases_a_s.get_theta(min_intensity=min_intensity,
                 min_theta = min_theta,
                 max_theta = max_theta,
                first_n_peaks = first_n_peaks)

if 'pmax_a' in locals():
    data.opt[0] = pmax_a
    data.opt[1] = pmax_s
pme = ChiMap().from_data(data,phases_a_s,sigma = sigma)
%%time
pme = pme.search()
L1loss, MSEloss, overlap_area = pme.metrics()
chi = pme.chi()
fig,ax = subplots(len(pme.phases),1,figsize=(12,10))

for i,phase in enumerate(pme.phases):
    ax[i].set_title(phase.label)
    p = ax[i].imshow(chi[...,i],vmin=0,vmax=1.2)
    colorbar(p,ax = ax[i])
    
show()
fig,ax = subplots(len(pme.phases),1,figsize=(12,10))

rescaling_chi = pme.chi() * data.rescaling**0.5

for i,phase in enumerate(pme.phases):
    ax[i].set_title(phase.label)
    p = ax[i].imshow(rescaling_chi[...,i],vmin=0,vmax=20)
    colorbar(p,ax = ax[i])
show()<jupyter_output><empty_output><jupyter_text>Histogram of $a$If $a$ is spead over too large area it might be that the phases are to right or a phase is missing<jupyter_code>%%time
opt = pme.opt()
a = opt[...,0]
s = opt[...,1]

vmin = -1220
vmax = -1190

h,b = histogram(a,bins=512)

figure(figsize=(8,4))

plot(b[:-1],h)
xlim(b[0],b[-1])
ylim(0,h.max())

vlines(vmin,0,h.max(),'k',ls='--',lw=1)
vlines(vmax,0,h.max(),'k',ls='--',lw=1)

xlabel('$a$')
ylabel(r'count($a$)')
title(r'Histogram of $a$')

figure(figsize=(8,8))
title('Distribution map of $a$')
im = imshow(a,cmap='Spectral',vmin=vmin,vmax=vmax)
colorbar(im,fraction=0.011)<jupyter_output>CPU times: user 99.7 ms, sys: 1.88 ms, total: 102 ms
Wall time: 99.2 ms<jupyter_text>Plotting the $a,s$ dependenceThere is a slight notion of a second $as$ dependence but it is weak.<jupyter_code>%matplotlib inline

opt,var = curve_fit(f_linear,a.flatten(),s.flatten())
res = least_squares(f_loss,x0=opt,args=(a.flatten(),s.flatten()),loss='cauchy')

linear_y = f_linear(a.flatten(),*opt)
cauchy_y = f_linear(a.flatten(),*res['x'])

print('Linear:',opt)
print('Cauchy:',res['x'])

plot(a.flatten(),s.flatten(),'.',alpha=0.01)

x = linspace(a.min(),a.max(),10)
plot(x,f_linear(x,*opt),'-.',lw=2,label='fit linear')
plot(x,f_linear(x,*res['x']),'--',lw=2,label='fit cauchy')
plot(data.opt[0],data.opt[1],'k+',ms=12,label='inital fit')

print(a.mean(),s.mean())

legend(frameon=False)

xlabel(r'$a$')
ylabel(r'$s$')

pmax_a = b[h.argmax()]
pmax_s = f_linear(pmax_a, *res['x'])

print(pmax_a,pmax_s)

plot(pmax_a,pmax_s,'r+',ms=12,label='most likely')

show()<jupyter_output>Linear: [ -1.32903602 341.16690989]
Cauchy: [-8.55064888e-01  9.13784409e+02]
-1206.0355826708758 1944.0316415604664
-1209.8693463895806 1948.3012066019305
{'fraction_non_alphanumeric': 0.16228144627397725, 'fraction_numerical': 0.04311661857070107, 'mean_word_length': 6.310173697270471, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 22, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '91216', 'license': 'MIT', 'n_tokens_mistral': 2515, 'n_tokens_neox': 2367, 'n_words': 432, 'path': 'notebooks/MunchMuzeumExp.ipynb', 'repo_name': 'zpreisler/XRDXRFutils', 'size': '5891'}