some bug fixes and updated functionality, imporved plots, now also rsplit and mutl

reduction_tools/partialator.py

@@ -15,9 +15,11 @@ python partialator.py -s <path-to-stream-file>
                        -c <path-to-cell-file>
                        -b number of resolution bins - must be > 20
                        -r high-res limt. Needs a default. Default set to 1.3
+                       -a max-adu. Default = 12000
 
 # output
 - scaled/merged files
+- an mtz file
 - useful plots
 - useful summerized .dat files
 """
@@ -60,7 +62,7 @@ def wait_for_jobs( job_ids, total_jobs ):
             job_ids.difference_update(completed_jobs)
             time.sleep(2)
 
-def run_partialator( proc_dir, name, stream, pointgroup, model, iterations, cell, shells, part_h_res ):
+def run_partialator( proc_dir, name, stream, pointgroup, model, iterations, cell, shells, part_h_res, adu ):
 
     # partialator file name
     part_run_file = "{0}/partialator_{1}.sh".format( proc_dir, name )
@@ -74,13 +76,14 @@ def run_partialator( proc_dir, name, stream, pointgroup, model, iterations, cell
     part_sh.write( "  -o merged_{0}.hkl \\\n".format( name ) )
     part_sh.write( "  -y {0} \\\n".format( pointgroup ) )
     part_sh.write( "  --model={0} \\\n".format( model ) )
+    part_sh.write( "  --max-adu={0} \\\n".format( adu ) )
     part_sh.write( "  --iterations={0}\n\n".format( iterations ) )
-    part_sh.write( "check_hkl --shell-file=mult.dat *.hkl -p {0} --nshells={1} --highres={2} >  check_hkl.log\n".format( cell, shells, part_h_res ) )
+    part_sh.write( "check_hkl --shell-file=mult.dat *.hkl -p {0} --nshells={1} --highres={2} &> check_hkl.log\n".format( cell, shells, part_h_res ) )
-    part_sh.write( "check_hkl --ltest --ignore-negs --shell-file=ltest.dat *.hkl -p {0} --nshells={1} --highres={2} >  ltest.log\n".format( cell, shells, part_h_res ) )
+    part_sh.write( "check_hkl --ltest --ignore-negs --shell-file=ltest.dat *.hkl -p {0} --nshells={1} --highres={2} &>  ltest.log\n".format( cell, shells, part_h_res ) )
-    part_sh.write( "check_hkl --wilson --shell-file=wilson.dat *.hkl -p {0} --nshells={1} --highres={2} >  wilson.log\n".format( cell, shells, part_h_res ) )
+    part_sh.write( "check_hkl --wilson --shell-file=wilson.dat *.hkl -p {0} --nshells={1} --highres={2} &> wilson.log\n".format( cell, shells, part_h_res ) )
-    part_sh.write( "compare_hkl --fom=Rsplit --shell-file=Rsplit.dat *.hkl1 *hkl2 -p {0} --nshells={1} --highres={2} >  Rsplit.log\n".format( cell, shells, part_h_res ) )
+    part_sh.write( "compare_hkl --fom=Rsplit --shell-file=Rsplit.dat *.hkl1 *hkl2 -p {0} --nshells={1} --highres={2} &>  Rsplit.log\n".format( cell, shells, part_h_res ) )
-    part_sh.write( "compare_hkl --fom=cc --shell-file=cc.dat *.hkl1 *hkl2  -p {0} --nshells={1} --highres={2} >  cc.log\n".format( cell, shells, part_h_res ) )
+    part_sh.write( "compare_hkl --fom=cc --shell-file=cc.dat *.hkl1 *hkl2  -p {0} --nshells={1} --highres={2} &> cc.log\n".format( cell, shells, part_h_res ) )
-    part_sh.write( "compare_hkl --fom=ccstar --shell-file=ccstar.dat *.hkl1 *hkl2  -p {0} --nshells={1} --highres={2} >  ccstar.log\n".format( cell, shells, part_h_res ) )
+    part_sh.write( "compare_hkl --fom=ccstar --shell-file=ccstar.dat *.hkl1 *hkl2  -p {0} --nshells={1} --highres={2} &>  ccstar.log\n".format( cell, shells, part_h_res ) )
     part_sh.close()
 
     # make file executable
@@ -97,7 +100,7 @@ def make_process_dir( dir ):
         if e.errno != errno.EEXIST:
             raise
 
-def summary_stats( cc_dat, ccstar_dat, mult_dat ):
+def summary_stats( cc_dat, ccstar_dat, mult_dat, rsplit_dat ):
 
     # read all files into pd
     # function to sort out different column names
@@ -111,22 +114,28 @@ def summary_stats( cc_dat, ccstar_dat, mult_dat ):
         elif var == "mult":
             cols = [ "d(nm)", "nref", "poss", "comp", "obs",
                     "mult", "snr", "I", "d", "min", "max" ]
+        elif var == "rsplit":
+            cols = [ "d(nm)", "rsplit", "nref", "d", "min", "max"  ]
         
         df = pd.read_csv( dat, names=cols, skiprows=1, sep="\s+" )
 
+        print(df)
         return df
         
+
     # make df
     cc_df = read_dat( cc_dat, "cc" )
     ccstar_df = read_dat( ccstar_dat, "ccstar" )
     mult_df = read_dat( mult_dat, "mult" )
+    rsplit_df = read_dat( rsplit_dat, "rsplit" )
 
     # remove unwanted cols
     cc_df = cc_df[ [ "cc" ] ]
     ccstar_df = ccstar_df[ [ "ccstar" ] ]
+    rsplit_df = rsplit_df[ [ "rsplit" ] ]
     
     # merge dfs
-    stats_df = pd.concat( [ mult_df, cc_df, ccstar_df], axis=1, join="inner" )
+    stats_df = pd.concat( [ mult_df, cc_df, ccstar_df, rsplit_df ], axis=1, join="inner" )
 
     # make 1/d, 1/d^2 column
     stats_df[ "1_d" ] = 1 / stats_df.d
@@ -136,7 +145,7 @@ def summary_stats( cc_dat, ccstar_dat, mult_dat ):
     stats_df = stats_df[ [  "1_d", "1_d2", "d", "min",
                             "max", "nref", "poss",
                             "comp", "obs", "mult",
-                            "snr", "I", "cc", "ccstar"] ]
+                            "snr", "I", "cc", "ccstar", "rsplit" ] ]
 
     # change nan to 0 
     stats_df = stats_df.fillna(0)
@@ -170,31 +179,82 @@ def get_metric( d2_series, cc_series, cut_off ):
 def summary_fig( stats_df ):
 
     # plot results
-    cc_fig, (ax1, ax2) = plt.subplots(1, 2)
+    cc_fig, axs = plt.subplots(2, 2)
     cc_fig.suptitle( "cc and cc* vs resolution" )
 
-    # indexed images plot
+    # cc plot
     color = "tab:red"
-    ax1.set_xlabel( "1/d2 (1/A)" )
+    axs[0,0].set_xlabel( "1/d (1/A)" )
-    ax1.set_ylabel("CC" )
+    axs[0,0].set_ylabel("CC" )
-    ax1.axhline(y = 0.3, color="black", linestyle = "dashed")
+    axs[0,0].set_ylim( 0, 1 )
-    ax1.plot(stats_df[ "1_d" ], stats_df.cc, color=color)
+    axs[0,0].axhline(y = 0.3, color="black", linestyle = "dashed")
-    ax1.tick_params(axis="y", labelcolor=color)
+    axs[0,0].plot(stats_df[ "1_d" ], stats_df.cc, color=color)
+    axs[0,0].tick_params(axis="y", labelcolor=color)
 
-    # label color
+    # cc* plot
     color = "tab:blue"
-    ax2.set_xlabel( "1/d (1/A)" )
+    axs[0,1].set_xlabel( "1/d (1/A)" )
-    ax2.set_ylabel("CC*", color=color)
+    axs[0,1].set_ylabel("CC*", color=color)
-    ax2.axhline(y = 0.7, color="black", linestyle = "dashed")
+    axs[0,1].set_ylim( 0, 1 )
-    ax2.plot(stats_df[ "1_d" ], stats_df.ccstar, color=color)
+    axs[0,1].axhline(y = 0.7, color="black", linestyle = "dashed")
-    ax2.tick_params(axis='y', labelcolor=color)
+    axs[0,1].plot(stats_df[ "1_d" ], stats_df.ccstar, color=color)
-    plt.show()
+    axs[0,1].tick_params(axis='y', labelcolor=color)
+
+    # rsplit plot
+    color = "tab:green"
+    axs[1,0].set_xlabel( "1/d (1/A)" )
+    axs[1,0].set_ylabel("Rsplit", color=color)
+    axs[1,0].plot(stats_df[ "1_d" ], stats_df.rsplit, color=color)
+    axs[1,0].tick_params(axis='y', labelcolor=color)
+
+
+    # rsplit plot
+    color = "tab:purple"
+    axs[1,1].set_xlabel( "1/d (1/A)" )
+    axs[1,1].set_ylabel("Multiplicity", color=color)
+    axs[1,1].plot(stats_df[ "1_d" ], stats_df.mult, color=color)
+    axs[1,1].tick_params(axis='y', labelcolor=color)
 
     # save figure
     plt.tight_layout()
-    plt.savefig("cc_curves.png")
+    plt.savefig("plots.png")
 
-def main( cwd, name, stream, pointgroup, model, iterations, cell, shells, part_h_res ):
+def get_mean_cell( stream ):
+
+    # get uc values from stream file
+    # example - Cell parameters 7.71784 7.78870 3.75250 nm, 90.19135 90.77553 90.19243 deg
+    # scrub clen and return - else nan
+    try:
+        pattern = r"Cell\sparameters\s(\d+\.\d+)\s(\d+\.\d+)\s(\d+\.\d+)\snm,\s(\d+\.\d+)\s(\d+\.\d+)\s(\d+\.\d+)\sdeg"
+        cell_lst = re.findall( pattern, stream )
+        xtals = len( cell_lst )
+    except AttributeError:
+        return np.nan
+    
+    cols = [ "a", "b", "c", "alpha", "beta", "gamma" ]
+    cell_df = pd.DataFrame( cell_lst, columns=cols )
+
+    mean_a = round( cell_df.a.mean()*10, 3 )
+    mean_b = round( cell_df.b.mean()*10, 3 )
+    mean_c = round( cell_df.c.mean()*10, 3 )
+    mean_alpha = round( cell_df.alpha.mean(), 3 )
+    mean_beta = round( cell_df.beta.mean(), 3 )
+    mean_gamma = round( cell_df.gamma.mean(), 3 )
+
+    return mean_a, mean_b, mean_c, mean_alpha, mean_beta, mean_gamma
+
+def make_mtz( hkl, mtz ):
+
+    # write file
+    mtz_sh = open( mtz_run_file, "w" )
+    mtz_sh.write( "#!/bin/sh\n\n" )
+    mtz_sh.write( "module purge\n" )
+    mtz_sh.write( "module load ccp4/8.0\n" )
+    mtz_sh.write( "f2mtz HKLIN {0} HKLOUT {1}\n",format( hkl, mtz ) )
+    mtz_sh.write( "TITLE Reflections from " )
+    mtz_sh.close()
+
+def main( cwd, name, stream, pointgroup, model, iterations, cell, shells, part_h_res, adu ):
 
     print( "begin job" )
     # submitted job set
@@ -209,7 +269,7 @@ def main( cwd, name, stream, pointgroup, model, iterations, cell, shells, part_h
 
     print( "making partialator files" )
     # make partialator run file
-    part_run_file = run_partialator( part_dir, name, stream, pointgroup, model, iterations, cell, shells, part_h_res )
+    part_run_file = run_partialator( part_dir, name, stream, pointgroup, model, iterations, cell, shells, part_h_res, adu )
 
     # submit job
     job_id = submit_job( part_run_file )
@@ -225,9 +285,10 @@ def main( cwd, name, stream, pointgroup, model, iterations, cell, shells, part_h
     cc_dat = "cc.dat"
     ccstar_dat = "ccstar.dat"
     mult_dat = "mult.dat"
+    rsplit_dat = "rsplit.dat"
 
     # make summary data table
-    stats_df = summary_stats( cc_dat, ccstar_dat, mult_dat )
+    stats_df = summary_stats( cc_dat, ccstar_dat, mult_dat, rsplit_dat )
     print( stats_df.to_string() )
     print_df = stats_df[ [  "1_d", "d", "min",
                             "max", "nref", "poss",
@@ -252,21 +313,23 @@ if __name__ == "__main__":
     parser.add_argument(
         "-n",
         "--name",
-        help="name of partialator run, also name of folder where data will be processed. Default = partialator",
+        help="name of partialator run, also name of folder where data will be processed.",
         type=str,
-        default="partialator"
+        required=True
     )
     parser.add_argument(
         "-s",
         "--stream_file",
         help="path to stream file",
-        type=os.path.abspath
+        type=os.path.abspath,
+        required=True
     )
     parser.add_argument(
         "-p",
         "--pointgroup",
         help="pointgroup used by CrystFEL for partialator run",
-        type=str
+        type=str,
+        required=True
     )
     parser.add_argument(
         "-m",
@@ -286,7 +349,8 @@ if __name__ == "__main__":
         "-c",
         "--cell_file",
         help="path to CrystFEL cell file for partialator",
-        type=os.path.abspath
+        type=os.path.abspath,
+        required=True
     )
     parser.add_argument(
         "-b",
@@ -302,11 +366,19 @@ if __name__ == "__main__":
         type=float,
         default=1.3
     )
+    parser.add_argument(
+        "-a",
+        "--max_adu",
+        help="maximum detector counts to allow. Default is 12000",
+        type=int,
+        default=12000
+    )
     args = parser.parse_args()
     # run main
     cwd = os.getcwd()
     print( "top working directory = {0}".format( cwd ) )
-    main( cwd, args.name, args.stream_file, args.pointgroup, args.model, args.iterations, args.cell_file, args.bins, args.resolution )
+    main( cwd, args.name, args.stream_file, args.pointgroup, args.model, args.iterations, args.cell_file, args.bins, args.resolution, args.max_adu )
+