diff --git a/burnman/eos/slb.py b/burnman/eos/slb.py index 1eb4ad8c8..525e15243 100644 --- a/burnman/eos/slb.py +++ b/burnman/eos/slb.py @@ -126,7 +126,7 @@ def volume(self, pressure, temperature, params): args = (pressure,temperature,V_0,T_0,Debye_0,n,a1_ii,a2_iikk,b_iikk,b_iikkmm) try: sol = bracket(_delta_pressure, params['V_0'], 1.e-2*params['V_0'], args) - except: + except ValueError: raise Exception('Cannot find a volume, perhaps you are outside of the range of validity for the equation of state?') return opt.brentq(_delta_pressure, sol[0], sol[1] ,args=args) diff --git a/misc/pyrolite_uncertainty.py b/misc/pyrolite_uncertainty.py index 985a39d78..a358cd8ea 100644 --- a/misc/pyrolite_uncertainty.py +++ b/misc/pyrolite_uncertainty.py @@ -60,15 +60,19 @@ def realize_pyrolite(): pv_fe_num = 0.07 fp_fe_num = 0.2 - mg_perovskite = minerals.SLB_2011_ZSB_2013.mg_perovskite(); realize_mineral(mg_perovskite) - fe_perovskite = minerals.SLB_2011_ZSB_2013.fe_perovskite(); realize_mineral(fe_perovskite) - wuestite = minerals.SLB_2011_ZSB_2013.wuestite(); realize_mineral(wuestite) - periclase = minerals.SLB_2011_ZSB_2013.periclase(); realize_mineral(periclase) + mg_perovskite = minerals.SLB_2011_ZSB_2013.mg_perovskite() + realize_mineral(mg_perovskite) + fe_perovskite = minerals.SLB_2011_ZSB_2013.fe_perovskite() + realize_mineral(fe_perovskite) + wuestite = minerals.SLB_2011_ZSB_2013.wuestite() + realize_mineral(wuestite) + periclase = minerals.SLB_2011_ZSB_2013.periclase() + realize_mineral(periclase) perovskite = HelperSolidSolution( [ mg_perovskite, fe_perovskite], [1.0-pv_fe_num, pv_fe_num]) ferropericlase = HelperSolidSolution( [ periclase, wuestite], [1.0-fp_fe_num, fp_fe_num]) - pyrolite = burnman.Composite( [ (perovskite, x_pv), (ferropericlase, x_fp) ] ) + pyrolite = burnman.Composite([perovskite,ferropericlase], [x_pv, x_fp]) pyrolite.set_method('slb3') anchor_temperature = normal(loc = 1935.0, scale = 200.0) @@ -77,54 +81,53 @@ def realize_pyrolite(): def output_rock( rock, file_handle ): - for ph in rock.staticphases: - if( isinstance(ph.mineral, burnman.minerals_base.helper_solid_solution) ): - for min in ph.mineral.endmembers: + for ph in rock.phases: + if( isinstance(ph, HelperSolidSolution) ): + for min in ph.endmembers: file_handle.write( '\t' + min.to_string() + '\n') for key in min.params: file_handle.write('\t\t' + key + ': ' + str(min.params[key]) + '\n') else: - file_handle.write( '\t' + ph.mineral.to_string() + '\n' ) - for key in ph.mineral.params: - file_handle.write('\t\t' + key + ': ' + str(ph.mineral.params[key]) + '\n') + file_handle.write( '\t' + ph.to_string() + '\n' ) + for key in ph.params: + file_handle.write('\t\t' + key + ': ' + str(ph.params[key]) + '\n') def realization_to_array(rock, anchor_t): arr = [anchor_t] names = ['anchor_T'] - for ph in rock.staticphases: - if( isinstance(ph.mineral, burnman.minerals_base.helper_solid_solution) ): - for min in ph.mineral.endmembers: + for ph in rock.phases: + if isinstance(ph, HelperSolidSolution): + for min in ph.endmembers: for key in min.params: if key != 'equation_of_state': arr.append(min.params[key]) names.append(min.to_string()+'.'+key) else: - for key in ph.mineral.params: + for key in ph.params: if key != 'equation_of_state': arr.append(ph.mineral.params[key]) names.append(ph.mineral.to_string()+'.'+key) return arr, names + def array_to_rock(arr, names): rock, _ = realize_pyrolite() anchor_t = arr[0] idx = 1 - for ph in rock.staticphases: - if( isinstance(ph.mineral, burnman.minerals_base.helper_solid_solution) ): - for min in ph.mineral.endmembers: - for key in min.params: - if key != 'equation_of_state': - assert(names[idx]==min.to_string()+'.'+key) - min.params[key] = arr[idx] + for phase in rock.phases: + if isinstance(phase, HelperSolidSolution): + for mineral in phase.endmembers: + while mineral.to_string() in names[idx]: + key = names[idx].split('.')[-1] + if key != 'equation_of_state' and key != 'F_0' and key != 'T_0' and key != 'P_0': + assert(mineral.to_string() in names[idx]) + mineral.params[key] = arr[idx] idx += 1 else: - for key in ph.mineral.params: - if key != 'equation_of_state': - assert(names[idx]==ph.mineral.to_string()+'.'+key) - ph.mineral.params[key] = arr[idx] - idx += 1 + raise Exception("unknown type") return rock, anchor_t + #set up the seismic model seismic_model = burnman.seismic.PREM() npts = 10 @@ -139,21 +142,17 @@ def array_to_rock(arr, names): whattodo = "" +dbname = "" -goodfits = []; -names = []; +goodfits = [] +names = [] -if len(sys.argv)<3: - print("options:") - print("run ") - print("plot") - print("plotgood ...") - print("plotone 1") -else: +if len(sys.argv)>=2: whattodo = sys.argv[1] +if len(sys.argv)>=3: dbname = sys.argv[2] -if whattodo=="plotgood": +if whattodo=="plotgood" and len(sys.argv)>2: files=sys.argv[2:] print("files:",files) names = pickle.load(open(files[0]+".names","rb")); @@ -166,7 +165,7 @@ def array_to_rock(arr, names): a = pickle.load(open(f,"rb")) allfits.extend(a) b = a - b = [i for i in a if i[erridx]<3e-5] + # b = [i for i in a if i[erridx]<3e-5] -- filter, need to adjust error value print("adding %d out of %d"%(len(b),len(a))) goodfits.extend(b) @@ -192,7 +191,7 @@ def array_to_rock(arr, names): figure=plt.figure(dpi=150,figsize=figsize) plt.subplots_adjust(hspace=0.3) for name in names: - if name.endswith(".n") or name.endswith(".V_0") or name.endswith(".molar_mass"): + if name.endswith(".n") or name.endswith(".V_0") or name.endswith(".molar_mass") or name.endswith(".F_0") or name.endswith(".P_0"): i+=1 continue plt.subplot(5,8,idx) @@ -218,6 +217,7 @@ def array_to_rock(arr, names): for entry in goodfits: trace.append(entry[i]) hist,bins = np.histogram(trace) + n, bins, patches = plt.hist(np.array(trace), 20, facecolor='green', alpha=0.75, normed=True) (mu, sigma) = norm.fit(np.array(trace)) y = mlab.normpdf( bins, mu, sigma) @@ -228,14 +228,15 @@ def array_to_rock(arr, names): i+=1 plt.savefig('good.png') + print("Writing good.png") #plt.show() figsize=(8,6) figure=plt.figure(dpi=150,figsize=figsize) for fit in goodfits: - print(fit) - print(names) + #print(fit) + #print(names) rock, anchor_t = array_to_rock(fit, names) temperature = burnman.geotherm.adiabatic(pressure, anchor_t, rock) @@ -263,9 +264,10 @@ def array_to_rock(arr, names): plt.savefig('goodones.png') + print("Writing goodones.png") plt.show() -if whattodo=="plotone": +elif whattodo=="plotone": figsize=(6,5) figure=plt.figure(dpi=100,figsize=figsize) @@ -447,11 +449,8 @@ def array_to_rock(arr, names): rock.set_averaging_scheme(burnman.averaging_schemes.HashinShtrikmanAverage()) rho, vs, vphi = rock.evaluate(['rho','v_s','v_phi'], pressure, temperature) - err_vs, err_vphi, err_rho = burnman.compare_l2(depths/np.mean(depths), vs/np.mean(seis_vs), vphi/np.mean(seis_vphi), \ - rho/np.mean(seis_rho), seis_vs/np.mean(seis_vs), seis_vphi/np.mean(seis_vphi), seis_rho/np.mean(seis_rho)) - error = np.sum([err_rho, err_vphi, err_vs]) - - print(error, err_rho, err_vphi, err_vs) + err_vs, err_vphi, err_rho = burnman.compare_l2(depths, [vs, vphi, rho], [seis_vs, seis_vphi, seis_rho]) + error = np.sum([err_rho/np.mean(seis_rho), err_vphi/np.mean(seis_vphi), err_vs/np.mean(seis_vs)]) figsize=(6,5) prop={'size':12} @@ -496,7 +495,7 @@ def array_to_rock(arr, names): #plt.show() -if whattodo=="run": +elif whattodo=="run" and len(sys.argv)>2: outfile = open(fname, 'w') outfile.write("#pressure\t Vs \t Vp \t rho \n") best_fit_file = open('output_pyrolite_closest_fit.txt', 'w') @@ -512,21 +511,22 @@ def array_to_rock(arr, names): print("realization", i+1) try: #create the ith model - pyrolite, anchor_temperature= realize_pyrolite() + pyrolite, anchor_temperature = realize_pyrolite() temperature = burnman.geotherm.adiabatic(pressure, anchor_temperature, pyrolite) #calculate the seismic observables - rock.set_averaging_scheme(burnman.averaging_schemes.HashinShtrikmanAverage()) + pyrolite.set_averaging_scheme(burnman.averaging_schemes.HashinShtrikmanAverage()) rho, vs, vphi = pyrolite.evaluate(['rho','v_s','v_phi'], pressure, temperature) #estimate the misfit with the seismic model - err_rho, err_vphi, err_vs = burnman.compare_l2(depths/np.mean(depths), vs/np.mean(seis_vs), vphi/np.mean(seis_vphi), \ - rho/np.mean(seis_rho), seis_vs/np.mean(seis_vs), seis_vphi/np.mean(seis_vphi), seis_rho/np.mean(seis_rho)) - error = np.sum([err_rho, err_vphi, err_vs]) + err_vs, err_vphi, err_rho = burnman.compare_l2(depths, [vs, vphi, rho], [seis_vs, seis_vphi, seis_rho]) + + error = np.sum([err_rho/np.mean(seis_rho), err_vphi/np.mean(seis_vphi), err_vs/np.mean(seis_vs)]) + if error < min_error: min_error = error - print(error) + print("new best realization with error", error) best_fit_file.write('Current best fit : '+str(error) + '\n' ) output_rock(pyrolite, best_fit_file) @@ -555,6 +555,7 @@ def array_to_rock(arr, names): outfile.close() best_fit_file.close() + elif whattodo=="error": values = [1957.020221991886, 1.6590112209181886, 249335164670.39246, 170883524675.03842, 0.8922515920546608, 4.083536182853109, 1.4680357687136616, 2.445e-05, 907.6618871363347, 0.1, 5, 1.4575168081960164, 1.3379195339709193, 260344929478.3809, 138077598973.27307, 0.17942226498091196, 1.3948903373340595, 1.436924855529012, 2.549e-05, 881.2532665499875, 0.1319, 5, 3.1204661890247394, 2.1411938868468483, 164407523972.7836, 131594720803.07439, 1.855224221011796, 3.867545309505681, 1.2953203656315155, 1.124e-05, 769.8199298156555, 0.0403, 2, 2.8860489779521985, 1.4263617489128713, 177341125271.45096, 59131041052.46985, 2.352310980469468, 5.1279202520952545, 1.6021924873676925, 1.226e-05, 440.13042122457716, 0.0718, 2, -1.6065263588976038, 7.5954915681374134e-05, 9.6441602176002807e-07, 4.4326026287552629e-05, 3.0664473372061482e-05] names = ['anchor_T', "'burnman.minerals.SLB_2011_ZSB_2013.mg_perovskite'.Gprime_0", "'burnman.minerals.SLB_2011_ZSB_2013.mg_perovskite'.K_0", "'burnman.minerals.SLB_2011_ZSB_2013.mg_perovskite'.G_0", "'burnman.minerals.SLB_2011_ZSB_2013.mg_perovskite'.q_0", "'burnman.minerals.SLB_2011_ZSB_2013.mg_perovskite'.Kprime_0", "'burnman.minerals.SLB_2011_ZSB_2013.mg_perovskite'.grueneisen_0", "'burnman.minerals.SLB_2011_ZSB_2013.mg_perovskite'.V_0", "'burnman.minerals.SLB_2011_ZSB_2013.mg_perovskite'.Debye_0", "'burnman.minerals.SLB_2011_ZSB_2013.mg_perovskite'.molar_mass", "'burnman.minerals.SLB_2011_ZSB_2013.mg_perovskite'.n", "'burnman.minerals.SLB_2011_ZSB_2013.mg_perovskite'.eta_s_0", "'burnman.minerals.SLB_2011_ZSB_2013.fe_perovskite'.Gprime_0", "'burnman.minerals.SLB_2011_ZSB_2013.fe_perovskite'.K_0", "'burnman.minerals.SLB_2011_ZSB_2013.fe_perovskite'.G_0", "'burnman.minerals.SLB_2011_ZSB_2013.fe_perovskite'.q_0", "'burnman.minerals.SLB_2011_ZSB_2013.fe_perovskite'.Kprime_0", "'burnman.minerals.SLB_2011_ZSB_2013.fe_perovskite'.grueneisen_0", "'burnman.minerals.SLB_2011_ZSB_2013.fe_perovskite'.V_0", "'burnman.minerals.SLB_2011_ZSB_2013.fe_perovskite'.Debye_0", "'burnman.minerals.SLB_2011_ZSB_2013.fe_perovskite'.molar_mass", "'burnman.minerals.SLB_2011_ZSB_2013.fe_perovskite'.n", "'burnman.minerals.SLB_2011_ZSB_2013.fe_perovskite'.eta_s_0", "'burnman.minerals.SLB_2011_ZSB_2013.periclase'.Gprime_0", "'burnman.minerals.SLB_2011_ZSB_2013.periclase'.K_0", "'burnman.minerals.SLB_2011_ZSB_2013.periclase'.G_0", "'burnman.minerals.SLB_2011_ZSB_2013.periclase'.q_0", "'burnman.minerals.SLB_2011_ZSB_2013.periclase'.Kprime_0", "'burnman.minerals.SLB_2011_ZSB_2013.periclase'.grueneisen_0", "'burnman.minerals.SLB_2011_ZSB_2013.periclase'.V_0", "'burnman.minerals.SLB_2011_ZSB_2013.periclase'.Debye_0", "'burnman.minerals.SLB_2011_ZSB_2013.periclase'.molar_mass", "'burnman.minerals.SLB_2011_ZSB_2013.periclase'.n", "'burnman.minerals.SLB_2011_ZSB_2013.periclase'.eta_s_0", "'burnman.minerals.SLB_2011_ZSB_2013.wuestite'.Gprime_0", "'burnman.minerals.SLB_2011_ZSB_2013.wuestite'.K_0", "'burnman.minerals.SLB_2011_ZSB_2013.wuestite'.G_0", "'burnman.minerals.SLB_2011_ZSB_2013.wuestite'.q_0", "'burnman.minerals.SLB_2011_ZSB_2013.wuestite'.Kprime_0", "'burnman.minerals.SLB_2011_ZSB_2013.wuestite'.grueneisen_0", "'burnman.minerals.SLB_2011_ZSB_2013.wuestite'.V_0", "'burnman.minerals.SLB_2011_ZSB_2013.wuestite'.Debye_0", "'burnman.minerals.SLB_2011_ZSB_2013.wuestite'.molar_mass", "'burnman.minerals.SLB_2011_ZSB_2013.wuestite'.n", "'burnman.minerals.SLB_2011_ZSB_2013.wuestite'.eta_s_0", 'err', 'err_rho', 'err_vphi', 'err_vs'] @@ -566,10 +567,10 @@ def array_to_rock(arr, names): rho, vs, vphi = rock.evaluate(['rho','v_s','v_phi'], pressure, temperature) - err_rho, err_vphi, err_vs = burnman.compare_l2(depths/np.mean(depths), vs/np.mean(seis_vs), vphi/np.mean(seis_vphi), \ - rho/np.mean(seis_rho), seis_vs/np.mean(seis_vs), seis_vphi/np.mean(seis_vphi), seis_rho/np.mean(seis_rho)) + err_vs, err_vphi, err_rho = burnman.compare_l2(depths, [vs, vphi, rho], [seis_vs, seis_vphi, seis_rho]) error = np.sum([err_rho, err_vphi, err_vs]) + print(error, err_rho, err_vphi, err_vs) @@ -645,4 +646,15 @@ def array_to_rock(arr, names): fig.set_size_inches(6.0, 6.0) if "RUNNING_TESTS" not in globals(): fig.savefig("pyrolite_uncertainty.pdf",bbox_inches='tight', dpi=100) + print("Writing pyrolite_uncertainty.pdf") plt.show() + +else: + print("Options:") + print(" run -- run realizations and write into given database name") + print(" plot -- plot given database") + print(" plotgood ... -- aggregate databases and plot") + print(" plotone -- plot a single hardcoded nice result") + print(" error -- testing, compute errors of a single hardcoded realization") + + diff --git a/misc/ref/colors.py.out b/misc/ref/colors.py.out new file mode 100644 index 000000000..e69de29bb diff --git a/misc/ref/performance.py.out b/misc/ref/performance.py.out new file mode 100644 index 000000000..5b2c1582a --- /dev/null +++ b/misc/ref/performance.py.out @@ -0,0 +1 @@ +8251.24506997 diff --git a/misc/ref/pyrolite_uncertainty.py.out b/misc/ref/pyrolite_uncertainty.py.out new file mode 100644 index 000000000..30584daf1 --- /dev/null +++ b/misc/ref/pyrolite_uncertainty.py.out @@ -0,0 +1,6 @@ +Options: + run -- run realizations and write into given database name + plot -- plot given database + plotgood ... -- aggregate databases and plot + plotone -- plot a single hardcoded nice result + error -- testing, compute errors of a single hardcoded realization diff --git a/test.sh b/test.sh index e9d254a28..154ffaade 100755 --- a/test.sh +++ b/test.sh @@ -107,8 +107,12 @@ cd .. echo "checking misc/ ..." cd misc -for test in `ls paper*.py` +for test in `ls *.py` do + [ $test == "gen_doc.py" ] && echo " *** skipping $test !" && continue + [ $test == "helper_solid_solution.py" ] && echo " *** skipping $test !" && continue + [ $test == "__init__.py" ] && echo " *** skipping $test !" && continue + [ $test == "table.py" ] && echo " *** skipping $test !" && continue [ $test == "paper_opt_pv_old.py" ] && echo " *** skipping $test !" && continue testit $test $fulldir