diff --git a/reduction_tools/partialator.py b/reduction_tools/partialator.py index 4d04d0d..bd59125 100644 --- a/reduction_tools/partialator.py +++ b/reduction_tools/partialator.py @@ -15,9 +15,11 @@ python partialator.py -s -c -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 --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( "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=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( "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 --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=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.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)" ) - ax1.set_ylabel("CC" ) - ax1.axhline(y = 0.3, color="black", linestyle = "dashed") - ax1.plot(stats_df[ "1_d" ], stats_df.cc, color=color) - ax1.tick_params(axis="y", labelcolor=color) + axs[0,0].set_xlabel( "1/d (1/A)" ) + axs[0,0].set_ylabel("CC" ) + axs[0,0].set_ylim( 0, 1 ) + axs[0,0].axhline(y = 0.3, color="black", linestyle = "dashed") + 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)" ) - ax2.set_ylabel("CC*", color=color) - ax2.axhline(y = 0.7, color="black", linestyle = "dashed") - ax2.plot(stats_df[ "1_d" ], stats_df.ccstar, color=color) - ax2.tick_params(axis='y', labelcolor=color) - plt.show() + axs[0,1].set_xlabel( "1/d (1/A)" ) + axs[0,1].set_ylabel("CC*", color=color) + axs[0,1].set_ylim( 0, 1 ) + axs[0,1].axhline(y = 0.7, color="black", linestyle = "dashed") + axs[0,1].plot(stats_df[ "1_d" ], stats_df.ccstar, color=color) + 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 ) +