Merge pull request geodynamics#224 from tjhei/test_and_fix_misc

Test and fix misc. Cheers!

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)
 

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 <dbname>")
-    print("plot")
-    print("plotgood <dbname1> <dbname2> ...")
-    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 <dbname>    -- run realizations and write into given database name")
+    print("  plot <dbname>   -- plot given database")
+    print("  plotgood <dbname1> <dbname2> ...   -- aggregate databases and plot")
+    print("  plotone  -- plot a single hardcoded nice result")
+    print("  error    -- testing, compute errors of a single hardcoded realization")
+
+

misc/ref/performance.py.out

@@ -0,0 +1 @@
+8251.24506997

misc/ref/pyrolite_uncertainty.py.out

@@ -0,0 +1,6 @@
+Options:
+  run <dbname>    -- run realizations and write into given database name
+  plot <dbname>   -- plot given database
+  plotgood <dbname1> <dbname2> ...   -- aggregate databases and plot
+  plotone  -- plot a single hardcoded nice result
+  error    -- testing, compute errors of a single hardcoded realization

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