#!/usr/bin/env python2.7
"""
relion_it.py
============

Script for automated, on-the-fly single-particle analysis in RELION 3

Authors: Sjors H.W. Scheres, Takanori Nakane & Colin M. Palmer

Usage:
    relion_it.py  [extra_options.py [extra_options2.py ....] ] [--gui] [--continue]

To get started, go to the intended location of your RELION project directory and make sure your micrographs are
accessible from within it (e.g. in a subdirectory called `Movies/' - use a symlink if necessary). Then run this
script, providing the names of files containing options if needed. (To call the script, you'll need to enter the full
path to it, put the directory containing it on your PATH environment variable, or put a copy of the script in the
current directory.)

Run with the `--gui' option to launch a simple GUI which will set up a run from a few basic options. (The GUI can
also be used to save a complete options file that you can then edit as required.)

Once the script is running, open a normal RELION GUI to see what's happening and visualise the results.

See below for full instructions including how to handle errors. If you have any problems, please edit the script as
needed, call on your local Python expert or email the CCP-EM mailing list (https://www.jiscmail.ac.uk/ccpem).


Overview
--------

relion_it.py creates a number of RELION jobs and then runs one or more `relion_pipeliner' processes to schedule them
(exactly like using the "Schedule" button in the RELION GUI). Instructions and information are printed to the terminal
by relion_it.py as it runs.

relion_it.py uses a large number of options to control how the jobs are run. It's designed to be very flexible and so
these options can be changed in a number of ways:

- The easiest way is to use the simple GUI (enabled by passing the `--gui' argument), which allows you to set a few
  simple options. These are then used to calculate appropriate values for the complete set of options. (See "Using the
  GUI" below for more information on this.)

- For more control, options can be put into one or more Python files (with a simple "option_name = value" format or
  with more complicated calculations - see "Options files" below for more information). The names of these options
  files can passed as command line arguments to relion_it.py.

- For maximum control, you can make your own copy of this script and change the option values and the code itself
  however you want.

Before running relion_it.py, you need to make sure you're in your intended RELION project directory, and that your
movie files are accessible by relative paths within that directory (as usual for a RELION project). You could do this
by moving the files from the microscope straight into the project directory, using a symlink from your project
directory to the real location of the data, or running a script to create a new symlink to each micrograph as it is
collected.


Options files
-------------

relion_it.py uses a large number of options for controlling both the flow of the script and the parameters for
individual jobs. These options can be read from Python script files when relion_it.py is started.

The options are all listed the body of the script below, with a comment to explain each option. One way to use this
script is to copy it in its entirety into your project directory, edit the options directly in the script and then
run it (with no command line arguments). However, it's often better to keep the script in the RELION source directory
(where it can be updated easily) and use options files to configure it.

An example of a simple options file is:

angpix = 1.06

This would override the default pixel size value, but leave all other options at their defaults.

The options files are read and interpreted as Python scripts. A simple list of "option_name = value" lines is all
that is needed, though you can also use any Python commands you like to do more complex calculations. To generate
an example file containing all of the options, run "relion_it.py --gui" and then click the "Save options" button,
which will save all the current options to a file called `relion_it_options.py' in the working directory.

The options are named descriptively so you can probably understand what most of them do quite easily. For more help on
any particular option, look at the comment above its definition in this script, or search the script's code to see
how it is used.

Options files can be useful as templates. As an example, at Diamond Light Source's eBIC facility, we have a template
file called `dls_cluster_options.py' that contains the necessary settings to make relion_it.py submit most of its jobs
to run on the DLS GPU cluster. You could also set up standard templates for a particular microscope (say, voltage and
Cs settings) or for a particular project or computer configuration.

When relion_it.py starts, it reads all options files in the order they are given on the command line. Subsequent files
will override earlier ones, so the last value given for any particular option will be the value that is used.

If you start relion_it.py with the `--continue' argument, it will automatically add `relion_it_options.py' to the end
of the list of options files. This means that if you are in a project directory where the relion_it.py GUI has
previously been used, all options will be defined in the relion_it_options.py file and they will override any other
options files given on the command line. (This is very useful for restarting the script after a problem, but it would
be pointless to combine `--continue' with any options template files.)

Note that if relion_it.py finds option names that it doesn't recognise while it's reading an options file, it will
print a warning (but continue anyway). If you've been editing options files by hand, you should check the output from
relion_it.py when it starts to make sure there are no typos in the options you wanted to set. (If you're using local
variables for intermediate Python calculations in an options file, it's a good idea to use names starting with a
leading underscore so you can immediately tell them apart from warnings about genuine spelling mistakes.)


Using the GUI
-------------

The GUI provides a simple way to start new projects with relion_it.py. If you want to use it, prepare your project
directory as described above, then start the GUI with "relion_it.py --gui". (If you're using any template options
files, you can give those too, for example "relion_it.py /path/to/site/options.py --gui".)

The window that appears should be self-explanatory. Fill in the options as needed for your project, and use the check
boxes on the right to control what processing steps will be done. When you're ready, click either "Save options" or
"Save & run". The program will check the values you've entered and then use them to calculate a few extra options for
relion_it.py. The options will then be saved to a file called `relion_it_options.py', and if you clicked "Save & run"
the processing run will start immediately.

If any of the entered values are invalid (for example, if there are letters in a field which should be a number), the
GUI will display a message box with an error when you click one of the buttons. It will also display a warning if any
values appear to be incorrect (but you can choose to ignore the warning by clicking "OK").

The GUI will try to calculate some extra options from the values you enter using the following rules:

1. If a 3D reference is given, use a single pass with reference-based autopicking, minimum distance between particles
   of 0.7 times the particle size, and a batch size of 100,000 particles.

2. If no 3D reference is given, run a first pass with reference-free LoG autopicking and a batch size of 10,000, and
   then a second pass with reference-based autopicking and a batch size of 100,000.

These options should be sensible in many cases, but if you'd like to change them, save the options from the GUI using
the "Save options" button, close the GUI, and edit the `relion_it_options.py' file to change the option values as
needed. You can then start the processing run with "relion_it.py --continue".


Running the pipelines
---------------------

relion_it.py uses several different scheduling pipelines to run its jobs. While each one is running, a file is
created in the project directory called `RUNNING_PIPELINER_<name>'. A log of the jobs run by that pipeline is stored
in `pipeline_<name>.log'.

If you want to stop one of the pipelines for any reason, delete its `RUNNING_' file and within a minute or two the
pipeliner will notice that the file has been removed and stop.

relion_it.py itself uses a similar file called `RUNNING_RELION_IT', and you can delete this to stop the script (which
will not affect any pipelines that are already running). It keeps a list of all of the jobs it has submitted in a
file called `RELION_IT_SUBMITTED_JOBS'. This file can be edited manually if necessary (but not while the script is
running!) Most of the jobs are run by the `preprocessing' pipeline. This will do the following:

  1. Import movies
  2. Motion correction
  3. CTF estimation
  4. Particle auto-picking
  5. Particle extraction
  6. Batch selection

After a number of particles have been extracted (1,000 by default), a 2D classification job will be run to provide
feedback on the quality of the data collection and particle picking.

Particles are split into batches of a fixed size (default 10,000 for the first pass with no reference, or 100,000
otherwise). The first batch is special: as it grows, the 2D classification job is re-run repeatedly to provide early
feedback on the quality of the data. For subsequent batches, the script waits for each batch to be complete before
running 2D classification on it.

You can provide reference structures for auto-picking and 3D classification. (If you provide a 3D reference in the
GUI it will automatically be used for both tasks.)

If you do not provide a reference for auto-picking, reference-free LoG picking will be used. If you do not provide a
reference for classification, relion_it.py will run the preprocessing pipeline twice. In the first pass, an initial
model will be generated, and then a second pass of preprocessing will be done using the initial model as a reference
for auto-picking and classification.

relion_it.py makes an effort to try to identify a suitable reference to use from the classes produced by the
InitialModel job, but if it selects an inappropriate reference, you can change it by stopping the pipelines and
script ("rm RUNNING_*"), updating the reference filename stored in the file named `RELION_IT_2NDPASS_3DREF', deleting
the relevant jobs (`autopick2_job' and those following) from the `RELION_IT_SUBMITTED_JOBS' file, then restarting the
pipeline with "relion_it.py --continue".


Fixing problems
---------------

One-off job failure
```````````````````
Occasionally, a single job can fail with an isolated error, for example if there are temporary network problems while
working on a remote filesystem. If this happens, RELION will wait forever for the files to appear that would indicate
the job has finished. In the meantime, no new jobs will be run, which can cause a backlog of micrographs to build up.

To fix this (for a preprocessing job), you can just try to re-run the job from the RELION GUI. Select the job in the
"Running jobs" list, then click "Job actions" -> "Mark as finished". Select the job again in the "Finished jobs"
list, then click "Continue!" to re-start the job.

That approach should work for preprocessing jobs, but probably won't work for classification or inital model
generation jobs, since those cannot be continued and must instead be restarted from the beginning. The best way to do
that is to restart the job manually, outside the RELION GUI, and then when the job finishes RELION should continue as
if the job had never failed.

For example, with a failed local job:

    ps -e | grep relion                            # to check if the job is still active
    kill <process_id>                              # to stop the job
    # now re-run the commands from the job's `note.txt' file

or with a job that was submitted to an SGE cluster queue:

    qstat                                          # to check if the job is still active in the queue
    qdel <job_id>                                  # to remove the job from the queue
    qsub job_type/job_directory/run_submit.script  # to re-submit the job

The other option is to just run a new job from the RELION GUI in the normal way (select the job you want to "copy" in
the jobs list, make a "new" job by clicking on the job type in the list in the top-left of the GUI, then click
"Run!"). However, if you do this, relion_it.py will not know about the new job and will not run any further
downstream processing based on it. In this situation, you can either continue to process your data manually in RELION,
or you could edit the `RELION_IT_SUBMITTED_JOBS' file to replace the failed job with the manual one, and delete the
jobs that followed the original one. After that, if you re-run the script it should continue as normal from that
job onwards.


Repeated job failure
````````````````````
If a job fails repeatedly, it usually indicates that there is some problem with the job parameters or the files that
the job needs to access.

In favourable cases, it's possible you could fix the problem by selecting the job in the RELION GUI, changing one of
the parameters that is not greyed out, then clicking "Continue!". Often, though, the problem will be with one of the
parameters that can't be changed for a job that already exists, so the job will need to be deleted and recreated with
a different set of parameters.

To handle this situation, stop all of the pipelines and the relion_it.py script ("rm RUNNING_*"), then identify and
fix the problem. Often, the problem will be an error in one of the job parameters, which can usually be fixed by
changing one of the script options (for example by changing the settings in `relion_it_options.py', if you originally
used the GUI to start the run).

If the problem is caused by missing files from an upstream job, you might need to check the output of previous jobs
and look in the job directories to figure out what the problem is. Again, if it's an error in the parameters for a
job, you can probably fix it by editing `relion_it_options.py'.

After changing any script options, you'll need to use the RELION GUI to delete the affected job and all jobs
downstream of it, and also remove them from the list in the `RELION_IT_SUBMITTED_JOBS' file. Then you should be able
to restart the pipelines by running "relion_it.py --continue".

If you still can't get a particular job to run without errors, you can at least continue to run the upstream jobs
that are working properly. You can do this either by changing the options for relion_it.py (there are options to
switch off 2D or 3D classification, or to stop after CTF estimation), or by manually scheduling the jobs you want
using the RELION GUI. Remember that after running relion_it.py, you have a normal RELION project, so if the script
can't do what you want, you can simply stop it and then use all of RELION's normal job management and scheduling
abilities.


Advanced usage
--------------

It's possible to customise many aspects of the way relion_it.py works, but the details go beyond the scope of this
introduction. Simple customisation can be done by setting appropriate option values (see "Option files" above). For
more substantial changes, you might need to edit the script's Python code to get the behaviour you want. Most of the
important logic is in the `run_pipeline()' function so that's a good place to start. Good luck!

"""

from __future__ import print_function
from __future__ import division  # always use float division

import argparse
import glob
import inspect
import math
import os
import runpy
import time
import traceback

try:
    import Tkinter as tk
    import tkMessageBox
    import tkFileDialog
except ImportError:
    # The GUI is optional. If the user requests it, it will fail when it tries
    # to open so we can ignore the error for now.
    pass  

# Constants
PIPELINE_STAR = 'default_pipeline.star'
RUNNING_FILE = 'RUNNING_RELION_IT'
SECONDPASS_REF3D_FILE = 'RELION_IT_2NDPASS_3DREF'
SETUP_CHECK_FILE = 'RELION_IT_SUBMITTED_JOBS'
PREPROCESS_SCHEDULE_PASS1 = 'PREPROCESS'
PREPROCESS_SCHEDULE_PASS2 = 'PREPROCESS_PASS2'
OPTIONS_FILE = 'relion_it_options.py'

class RelionItOptions(object):
    """
    Options for the relion_it pipeline setup script.
    
    When initialised, this contains default values for all options. Call
    ``update_from()`` to override the defaults with a dictionary of new values.
    """
    #############################################################################
    # Change the parameters below to reflect your experiment                    #
    # Current defaults reflect cryo-ARM betagal data set of RELION-3.0 tutorial #
    #############################################################################

    ### General parameters
    # Pixel size in Angstroms in the input movies
    angpix = 0.885
    # Acceleration voltage (in kV)
    voltage = 200
    # Polara = 2.0; Talos/Krios = 2.7; some Cryo-ARM = 1.4 
    Cs = 1.4


    ### Import images (Linux wild card; movies as *.mrc, *.mrcs, *.tiff or *.tif; single-frame micrographs as *.mrc)
    import_images = 'Movies/*.tiff'
    # Are these multi-frame movies? Set to False for single-frame micrographs (and motion-correction will be skipped)
    images_are_movies = True


    ### MotionCorrection parameters
    # Dose in electrons per squared Angstrom per fraction
    motioncor_doseperframe = 1.277
    # Gain-reference image in MRC format (only necessary if input movies are not yet gain-corrected, e.g. compressed TIFFs from K2)
    motioncor_gainreference = 'Movies/gain.mrc'
    # EER upsampling (1 = 4K, 2 = 8K). If you use 8K rendering, the pixel size (angpix) MUST be the half of the physical pixel size and the motioncor_binning should be 2.
    eer_upsampling = 1
    # EER fractionation. The dose rate (motioncor_doseperframe) is e/A2/fraction after this fractionation.
    eer_grouping = 20


    ### CTF estimation parameters
    # Most cases won't need changes here...

    ### Autopick parameters
    # Use reference-free Laplacian-of-Gaussian picking (otherwise use reference-based template matching instead)
    autopick_do_LoG = True
    # Minimum and maximum diameter in Angstrom for the LoG filter
    autopick_LoG_diam_min = 150
    autopick_LoG_diam_max = 180
    # Use positive values (0-1) to pick fewer particles; use negative values (-1-0) to pick more particles
    autopick_LoG_adjust_threshold = 0.0
    autopick_LoG_upper_threshold = 999.0
    #
    # OR:
    #
    # References for reference-based picking (when autopick_do_LoG = False)
    autopick_2dreferences = ''
    # OR: provide a 3D references for reference-based picking (when autopick_do_LoG = False)
    autopick_3dreference = ''

    # Threshold for reference-based autopicking (threshold 0 will pick too many particles. Default of 0.4 is hopefully better. Ultimately, just hope classification will sort it all out...)
    autopick_refs_threshold = 0.4
    # Minimum inter-particle distance for reference-based picking (~70% of particle diameter often works well)
    autopick_refs_min_distance = 120
    #
    # For both LoG and refs:
    #
    # Use this to remove false positives from carbon edges (useful range: 1.0-1.2, -1 to switch off)
    autopick_stddev_noise = -1
    # Use this to remove false positives from carbon edges (useful range: -0.5-0.0; -999 to switch off)
    autopick_avg_noise = -999


    ### Extract parameters
    # Box size of particles in the averaged micrographs (in pixels)
    extract_boxsize = 256
    # Down-scale the particles upon extraction?
    extract_downscale = False
    # Box size of the down-scaled particles (in pixels)
    extract_small_boxsize = 64
    # In second pass, down-scale the particles upon extraction?
    extract2_downscale = False
    # In second pass, box size of the down-scaled particles (in pixels)
    extract2_small_boxsize = 128

    
    ### Now perform 2D and/or 3D classification with the extracted particles?
    do_class2d = True
    # And/or perform 3D classification?
    do_class3d = True
    # Repeat 2D and/or 3D-classification for batches of this many particles
    batch_size = 10000
    # Number of 2D classes to use
    class2d_nr_classes  = 50
    # Diameter of the mask used for 2D/3D classification (in Angstrom)
    mask_diameter = 190
    # Symmetry group (when using SGD for initial model generation, C1 may work best)
    symmetry = 'C1'
    #
    ### 3D-classification parameters
    # Number of 3D classes to use
    class3d_nr_classes = 4
    # Have initial 3D model? If not, calculate one using SGD initial model generation
    have_3d_reference = False
    # Initial reference model
    class3d_reference = ''
    # Is reference on correct greyscale?
    class3d_ref_is_correct_greyscale = False
    # Has the initial reference been CTF-corrected?
    class3d_ref_is_ctf_corrected = True
    # Initial lowpass filter on reference
    class3d_ini_lowpass = 40 


    ### Use the largest 3D class from the first batch as a 3D reference for a second pass of autopicking? (only when do_class3d is True)
    do_second_pass = True
    # Only move on to template-based autopicking if the 3D references achieves this resolution (in A)
    minimum_resolution_3dref_2ndpass = 20
    # In the second pass, perform 2D classification?
    do_class2d_pass2 = True
    # In the second pass, perform 3D classification?
    do_class3d_pass2 = False
    # Batch size in the second pass
    batch_size_pass2 = 100000
    

    ###################################################################################
    ############ Often the parameters below can be kept the same for a given set-up
    ###################################################################################

    ### Repeat settings for entire pipeline
    # Repeat the pre-processing runs this many times (or until RUNNING_PIPELINER_default_PREPROCESS file is deleted)
    preprocess_repeat_times = 999
    # Wait at least this many minutes between each repeat cycle
    preprocess_repeat_wait = 1
    ### Stop after CTF estimation? I.e., skip autopicking, extraction, 2D/3D classification, etc?
    stop_after_ctf_estimation = False
    # Check every this many minutes if enough particles have been extracted for a new batch of 2D-classification
    batch_repeat_time = 1


    ### MotionCorrection parameters
    # Use RELION's own implementation of motion-correction (CPU-only) instead of the UCSF implementation?
    motioncor_do_own = True
    # The number of threads (only for RELION's own implementation) is optimal when nr_movie_frames/nr_threads = integer
    motioncor_threads = 6
    # Exectutable of UCSF MotionCor2
    motioncor_exe = '/public/EM/MOTIONCOR2/MotionCor2'
    # On which GPU(s) to execute UCSF MotionCor2
    motioncor_gpu = '0'
    # How many MPI processes to use for running motion correction?
    motioncor_mpi = 4
    # Local motion-estimation patches for MotionCor2
    motioncor_patches_x = 4
    motioncor_patches_y = 4
    # B-factor in A^2 for downweighting of high-spatial frequencies
    motioncor_bfactor = 150
    # Use binning=2 for super-resolution movies
    motioncor_binning = 1
    # Provide a defect file for your camera if you have one
    motioncor_defectfile = ''
    # orientation of the gain-reference w.r.t your movies (if input movies are not yet gain-corrected, e.g. TIFFs)
    motioncor_gainflip = 'No flipping (0)'
    motioncor_gainrot = 'No rotation (0)'
    # Other arguments for MotionCor2
    motioncor_other_args = ''
    # Submit motion correction job to the cluster?
    motioncor_submit_to_queue = False
    

    ### CTF estimation parameters
    # Amplitude contrast (Q0)
    ampl_contrast = 0.1
    # CTFFIND-defined parameters
    ctffind_boxsize = 512
    ctffind_astigmatism = 100
    ctffind_maxres = 5
    ctffind_minres = 30
    ctffind_defocus_max = 50000
    ctffind_defocus_min = 5000
    ctffind_defocus_step = 500
    # For Gctf: ignore parameters on the 'Searches' tab?
    ctffind_do_ignore_search_params = True
    # For Gctf: perform equi-phase averaging?
    ctffind_do_EPA = True
    # Also estimate phase shifts (for VPP data)
    ctffind_do_phaseshift = False
    # Executable to Kai Zhang's Gctf
    gctf_exe = '/public/EM/Gctf/bin/Gctf'
    # On which GPU(s) to execute Gctf
    gctf_gpu = '0'
    # Use Alexis Rohou's CTFFIND4 (CPU-only) instead?
    use_ctffind_instead = True
    # Executable for Alexis Rohou's CTFFIND4
    ctffind4_exe = '/public/EM/ctffind/ctffind.exe'
    # How many MPI processes to use for running CTF estimation?
    ctffind_mpi = 8
    # Submit CTF estimation job to the cluster?
    ctffind_submit_to_queue = False


    ### Autopick parameters
    # Use GPU-acceleration for autopicking?
    autopick_do_gpu = True
    # Which GPU(s) to use for autopicking
    autopick_gpu = '0'
    # Low-pass filter for auto-picking the micrographs
    autopick_lowpass = 20
    # Shrink factor for faster picking (0 = fastest; 1 = slowest)
    autopick_shrink_factor = 0
    # How many MPI processes to use for running auto-picking?
    autopick_mpi = 1
     # Additional arguments for autopicking
    autopick_other_args = ''
    # Submit Autopick job to the cluster?
    autopick_submit_to_queue = False
    # Are the references CTF-corrected?
    autopick_refs_are_ctf_corrected = True
    # Do the references have inverted contrast wrt the micrographs?
    autopick_refs_have_inverted_contrast = True
    # Ignore CTFs until the first peak
    autopick_refs_ignore_ctf1stpeak = False
    # Diameter of mask for the references (in A; negative value for automated detection of mask diameter)
    autopick_refs_mask_diam = -1
    # In-plane angular sampling interval
    autopick_inplane_sampling = 10
    # Symmetry of the 3D reference for autopicking
    autopick_3dref_symmetry = 'C1'
    # 3D angular sampling for generating projections of the 3D reference for autopicking (30 degrees is usually enough)
    autopick_3dref_sampling = '30 degrees'
    # Pixel size in the provided 2D/3D references (negative for same as in motion-corrected movies)
    autopick_ref_angpix = -1

    ### Extract parameters
    # Diameter for background normalisation (in pixels; negative value: default is 75% box size)
    extract_bg_diameter = -1
    # How many MPI processes to use for running particle extraction?
    extract_mpi = 1
    # Submit Extract job to the cluster?
    extract_submit_to_queue = False


    ## Discard particles based on average/stddev values? (this may be important for SGD initial model generation)
    do_discard_on_image_statistics = False
    # Discard images that have average/stddev values that are more than this many sigma away from the ensemble average
    discard_sigma = 4
    # Submit discard job to the cluster?
    discard_submit_to_queue = False


    #### Common relion_refine paremeters used for 2D/3D classification and initial model generation 
    # Read all particles in one batch into memory?
    refine_preread_images = False
    # Or copy particles to scratch disk?
    refine_scratch_disk = ''
    # Number of pooled particles?
    refine_nr_pool = 10
    # Use GPU-acceleration?
    refine_do_gpu = True
    # Which GPU to use (different from GPU used for pre-processing?)
    refine_gpu = '1'
    # How many MPI processes to use
    refine_mpi = 1
    # How many threads to use
    refine_threads = 6
    # Skip padding?
    refine_skip_padding = False
    # Submit jobs to the cluster?
    refine_submit_to_queue = False
    # Use fast subsets in 2D/3D classification when batch_size is bigger than this
    refine_batchsize_for_fast_subsets = 10000

    ### 2D classification parameters
    # Wait with the first 2D classification batch until at least this many particles are extracted
    minimum_batch_size = 10000
    # Number of iterations to perform in 2D classification
    # Must be at least 20 for fast subsets
    class2d_nr_iter = 20
    # Rotational search step (in degrees)
    class2d_angle_step = 6
    # Offset search range (in pixels)
    class2d_offset_range = 5
    # Offset search step (in pixels)
    class2d_offset_step = 1
    # Option to ignore the CTFs until their first peak (try this if all particles go into very few classes) 
    class2d_ctf_ign1stpeak = False
    # Additional arguments to pass to relion-refine
    class2d_other_args = ''
    

    ### 3D classification parameters
    # Number of iterations to perform in 3D classification
    # Must be at least 20 for fast subsets
    class3d_nr_iter = 20
    # Reference mask
    class3d_reference_mask = ''
    # Option to ignore the CTFs until their first peak (try this if all particles go into very few classes) 
    class3d_ctf_ign1stpeak = False
    # Regularisation parameter (T)
    class3d_T_value = 4
    # Angular sampling step
    class3d_angle_step = '7.5 degrees'
    # Offset search range (in pixels)
    class3d_offset_range = 5
    # Offset search step (in pixels)
    class3d_offset_step = 1
    # Additional arguments to pass to relion-refine
    class3d_other_args = ''


    ## SGD initial model generation
    # Number of models to generate simulatenously (K>1 may be useful for getting rid of outliers in the particle images)
    inimodel_nr_classes = 4
    # Ignore CTFs until first peak?
    inimodel_ctf_ign1stpeak = False
    # Enforce non-negative solvent?
    inimodel_solvent_flatten = True
    # Initial angular sampling
    inimodel_angle_step = '15 degrees'
    # Initial search range (in pixels)
    inimodel_offset_range = 6 
    # Initial offset search step (in pixels)
    inimodel_offset_step = 2
    # Number of initial iterations
    inimodel_nr_iter_initial = 50
    # Number of in-between iterations
    inimodel_nr_iter_inbetween = 200
    # Number of final iterations
    inimodel_nr_iter_final = 50
    # Frequency to write out information
    inimodel_freq_writeout = 10
    # Initial resolution (in A)
    inimodel_resol_ini = 35
    # Final resolution (in A)
    inimodel_resol_final = 15
    # Initial mini-batch size
    inimodel_batchsize_ini = 100
    # Final mini-batch size
    inimodel_batchsize_final = 500
    # Increased noise variance half-life (off, i.e. -1, by default; values of ~1000 have been observed to be useful in difficult cases)
    inimodel_sigmafudge_halflife = -1
    # Additional arguments to pass to relion_refine (skip annealing to get rid of outlier particles)
    inimodel_other_args = ' --sgd_skip_anneal '


    ### Cluster submission settings
    # Name of the queue to which to submit the job
    queue_name = 'openmpi'
    # Name of the command used to submit scripts to the queue
    queue_submit_command = 'qsub'
    # The template for your standard queue job submission script
    queue_submission_template = '/public/EM/RELION/relion/bin/qsub.csh'
    # Minimum number of dedicated cores that need to be requested on each node
    queue_minimum_dedicated = 1

    ### End of options
    
    #######################################################################
    ############ typically no need to change anything below this line
    #######################################################################

    def update_from(self, other):
        """
        Update this RelionItOptions object from a dictionary.
        
        Special values (with names like '__xxx__') are removed, allowing this
        method to be given a dictionary containing the namespace from a script
        run with ``runpy``.
        """
        while len(other) > 0:
            key, value = other.popitem()
            if not (key.startswith('__') and key.endswith('__')): # exclude __name__, __builtins__ etc.
                if hasattr(self, key):
                    setattr(self, key, value)
                else:
                    print(" RELION_IT: Unrecognised option '{}'".format(key))
    
    def print_options(self, out_file=None):
        """
        Print the current options.
        
        This method prints the options in the same format as they are read,
        allowing options to be written to a file and re-used.

        Args:
            out_file: A file object (optional). If supplied, options will be
                written to this file, otherwise they will be printed to
                sys.stdout.

        Raises:
            ValueError: If there is a problem printing the options.
        """
        out_file.write("# Options file for relion_it.py\n\n")
        seen_start = False
        option_names = [key for key in dir(self) if (not (key.startswith('__') and key.endswith('__'))
                                                     and not callable(getattr(self, key)))]

        # Parse the source code for this class, and write out all comments along with option lines containing new values
        for line in inspect.getsourcelines(RelionItOptions)[0]:
            line = line.strip()
            if not seen_start:
                if line != "### General parameters":
                    # Ignore lines until this one
                    continue
                seen_start = True
            if line == "### End of options":
                # Stop here
                break
            if line.startswith('#') or len(line) == 0:
                # Print comments or blank lines as-is
                out_file.write(line + "\n");
            else:
                # Assume all other lines define an option name and value. Replace with new value.
                equals_index = line.find('=')
                if equals_index > 0:
                    option_name = line[:equals_index].strip()
                    if option_name in option_names:
                        out_file.write('{} = {}\n'.format(option_name, repr(getattr(self, option_name))))
                        option_names.remove(option_name)
                    else:
                        # This error should not occur. If it does, there is probably a programming error.
                        raise ValueError("Unrecognised option name '{}'".format(option_name))
        if len(option_names) > 0:
            # This error should not occur. If it does, there is probably a programming error.
            raise ValueError("Some options were not written to the output file: {}".format(option_names))


class RelionItGui(object):

    def __init__(self, main_window, options):
        self.main_window = main_window
        self.options = options

        # Convenience function for making file browser buttons
        def new_browse_button(master, var_to_set, filetypes=(('MRC file', '*.mrc'), ('All files', '*'))):
            def browse_command():
                chosen_file = tkFileDialog.askopenfilename(filetypes=filetypes)
                if chosen_file is not None:
                    # Make path relative if it's in the current directory
                    if chosen_file.startswith(os.getcwd()):
                        chosen_file = os.path.relpath(chosen_file)
                    var_to_set.set(chosen_file)
            return tk.Button(master, text="Browse...", command=browse_command)

        ### Create GUI

        main_frame = tk.Frame(main_window)
        main_frame.pack(fill=tk.BOTH, expand=1)

        left_frame = tk.Frame(main_frame)
        left_frame.pack(side=tk.LEFT, anchor=tk.N, fill=tk.X, expand=1)

        right_frame = tk.Frame(main_frame)
        right_frame.pack(side=tk.LEFT, anchor=tk.N, fill=tk.X, expand=1)

        ###
        
        expt_frame = tk.LabelFrame(left_frame, text="Experimental details", padx=5, pady=5)
        expt_frame.pack(padx=5, pady=5, fill=tk.X, expand=1)
        tk.Grid.columnconfigure(expt_frame, 1, weight=1)

        row = 0

        tk.Label(expt_frame, text="Voltage (kV):").grid(row=row, sticky=tk.W)
        self.voltage_entry = tk.Entry(expt_frame)
        self.voltage_entry.grid(row=row, column=1, sticky=tk.W+tk.E)
        self.voltage_entry.insert(0, str(options.voltage))

        row += 1
        
        tk.Label(expt_frame, text="Cs (mm):").grid(row=row, sticky=tk.W)
        self.cs_entry = tk.Entry(expt_frame)
        self.cs_entry.grid(row=row, column=1, sticky=tk.W+tk.E)
        self.cs_entry.insert(0, str(options.Cs))

        row += 1
        
        tk.Label(expt_frame, text="Phase plate?").grid(row=row, sticky=tk.W)
        self.phaseplate_var = tk.IntVar()
        phaseplate_button = tk.Checkbutton(expt_frame, var=self.phaseplate_var)
        phaseplate_button.grid(row=row, column=1, sticky=tk.W)
        if options.ctffind_do_phaseshift:
            phaseplate_button.select()

        row += 1

        tk.Label(expt_frame, text=u"Pixel size (\u212B):").grid(row=row, sticky=tk.W)
        self.angpix_var = tk.StringVar()  # for data binding
        self.angpix_entry = tk.Entry(expt_frame, textvariable=self.angpix_var)
        self.angpix_entry.grid(row=row, column=1, sticky=tk.W+tk.E)
        self.angpix_entry.insert(0, str(options.angpix))

        row += 1
        
        tk.Label(expt_frame, text=u"Exposure rate (e\u207B / \u212B\u00B2 / frame):").grid(row=row, sticky=tk.W)
        self.exposure_entry = tk.Entry(expt_frame)
        self.exposure_entry.grid(row=row, column=1, sticky=tk.W + tk.E)
        self.exposure_entry.insert(0, str(options.motioncor_doseperframe))

        ###

        particle_frame = tk.LabelFrame(left_frame, text="Particle details", padx=5, pady=5)
        particle_frame.pack(padx=5, pady=5, fill=tk.X, expand=1)
        tk.Grid.columnconfigure(particle_frame, 1, weight=1)

        row = 0

        tk.Label(particle_frame, text=u"Longest diameter (\u212B):").grid(row=row, sticky=tk.W)
        self.particle_max_diam_var = tk.StringVar()  # for data binding
        self.particle_max_diam_entry = tk.Entry(particle_frame, textvariable=self.particle_max_diam_var)
        self.particle_max_diam_entry.grid(row=row, column=1, sticky=tk.W+tk.E, columnspan=2)
        self.particle_max_diam_entry.insert(0, str(options.autopick_LoG_diam_max))

        row += 1

        tk.Label(particle_frame, text=u"Shortest diameter (\u212B):").grid(row=row, sticky=tk.W)
        self.particle_min_diam_entry = tk.Entry(particle_frame)
        self.particle_min_diam_entry.grid(row=row, column=1, sticky=tk.W+tk.E, columnspan=2)
        self.particle_min_diam_entry.insert(0, str(options.autopick_LoG_diam_min))

        row += 1
        
        tk.Label(particle_frame, text="3D reference (optional):").grid(row=row, sticky=tk.W)
        self.ref_3d_var = tk.StringVar()  # for data binding
        self.ref_3d_entry = tk.Entry(particle_frame, textvariable=self.ref_3d_var)
        self.ref_3d_entry.grid(row=row, column=1, sticky=tk.W+tk.E)
        self.ref_3d_entry.insert(0, str(options.autopick_3dreference))

        new_browse_button(particle_frame, self.ref_3d_var).grid(row=row, column=2)

        row += 1

        tk.Label(particle_frame, text=u"Mask diameter (\u212B):").grid(row=row, sticky=tk.W)
        self.mask_diameter_var = tk.StringVar()  # for data binding
        self.mask_diameter_entry = tk.Entry(particle_frame, textvariable=self.mask_diameter_var)
        self.mask_diameter_entry.grid(row=row, column=1, sticky=tk.W+tk.E)
        self.mask_diameter_entry.insert(0, str(options.mask_diameter))
        self.mask_diameter_px = tk.Label(particle_frame, text="= NNN px")
        self.mask_diameter_px.grid(row=row, column=2,sticky=tk.W)

        row += 1

        tk.Label(particle_frame, text="Box size (px):").grid(row=row, sticky=tk.W)
        self.box_size_var = tk.StringVar()  # for data binding
        self.box_size_entry = tk.Entry(particle_frame, textvariable=self.box_size_var)
        self.box_size_entry.grid(row=row, column=1, sticky=tk.W+tk.E)
        self.box_size_entry.insert(0, str(options.extract_boxsize))
        self.box_size_in_angstrom = tk.Label(particle_frame, text=u"= NNN \u212B")
        self.box_size_in_angstrom.grid(row=row, column=2,sticky=tk.W)

        row += 1

        tk.Label(particle_frame, text="Down-sample to (px):").grid(row=row, sticky=tk.W)
        self.extract_small_boxsize_var = tk.StringVar()  # for data binding
        self.extract_small_boxsize_entry = tk.Entry(particle_frame, textvariable=self.extract_small_boxsize_var)
        self.extract_small_boxsize_entry.grid(row=row, column=1, sticky=tk.W+tk.E)
        self.extract_small_boxsize_entry.insert(0, str(options.extract_small_boxsize))
        self.extract_angpix = tk.Label(particle_frame, text=u"= NNN \u212B/px")
        self.extract_angpix.grid(row=row, column=2,sticky=tk.W)

        row += 1

        tk.Label(particle_frame, text="Calculate for me:").grid(row=row, sticky=tk.W)
        self.auto_boxsize_var = tk.IntVar()
        auto_boxsize_button = tk.Checkbutton(particle_frame, var=self.auto_boxsize_var)
        auto_boxsize_button.grid(row=row, column=1, sticky=tk.W)
        auto_boxsize_button.select()

        ###

        project_frame = tk.LabelFrame(right_frame, text="Project details", padx=5, pady=5)
        project_frame.pack(padx=5, pady=5, fill=tk.X, expand=1)
        tk.Grid.columnconfigure(project_frame, 1, weight=1)

        row = 0

        tk.Label(project_frame, text="Project directory:").grid(row=row, sticky=tk.W)
        tk.Label(project_frame, text=os.getcwd(), anchor=tk.W).grid(row=row, column=1, sticky=tk.W, columnspan=2)

        row += 1

        tk.Label(project_frame, text="Pattern for movies:").grid(row=row, sticky=tk.W)
        self.import_images_var = tk.StringVar()  # for data binding
        self.import_images_entry = tk.Entry(project_frame, textvariable=self.import_images_var)
        self.import_images_entry.grid(row=row, column=1, sticky=tk.W+tk.E)
        self.import_images_entry.insert(0, self.options.import_images)

        import_button = new_browse_button(project_frame, self.import_images_var,
                                          filetypes=(('Image file', '{*.mrc, *.mrcs, *.tif, *.tiff}'), ('All files', '*')))
        import_button.grid(row=row, column=2)

        row += 1
        
        tk.Label(project_frame, text="Gain reference (optional):").grid(row=row, sticky=tk.W)
        self.gainref_var = tk.StringVar()  # for data binding
        self.gainref_entry = tk.Entry(project_frame, textvariable=self.gainref_var)
        self.gainref_entry.grid(row=row, column=1, sticky=tk.W+tk.E)
        self.gainref_entry.insert(0, self.options.motioncor_gainreference)

        new_browse_button(project_frame, self.gainref_var).grid(row=row, column=2)

        ###

        pipeline_frame = tk.LabelFrame(right_frame, text="Pipeline control", padx=5, pady=5)
        pipeline_frame.pack(padx=5, pady=5, fill=tk.X, expand=1)
        tk.Grid.columnconfigure(expt_frame, 1, weight=1)

        row = 0

        tk.Label(pipeline_frame, text="Stop after CTF estimation?").grid(row=row, sticky=tk.W)
        self.stop_after_ctf_var = tk.IntVar()
        stop_after_ctf_button = tk.Checkbutton(pipeline_frame, var=self.stop_after_ctf_var)
        stop_after_ctf_button.grid(row=row, column=1, sticky=tk.W)
        if options.stop_after_ctf_estimation:
            stop_after_ctf_button.select()

        row += 1

        tk.Label(pipeline_frame, text="Do 2D classification?").grid(row=row, sticky=tk.W)
        self.class2d_var = tk.IntVar()
        class2d_button = tk.Checkbutton(pipeline_frame, var=self.class2d_var)
        class2d_button.grid(row=row, column=1, sticky=tk.W)
        if options.do_class2d:
            class2d_button.select()

        row += 1

        tk.Label(pipeline_frame, text="Do 3D classification?").grid(row=row, sticky=tk.W)
        self.class3d_var = tk.IntVar()
        class3d_button = tk.Checkbutton(pipeline_frame, var=self.class3d_var)
        class3d_button.grid(row=row, column=1, sticky=tk.W)
        if options.do_class3d:
            class3d_button.select()

        row += 1

        tk.Label(pipeline_frame, text="Do second pass? (only if no 3D ref)").grid(row=row, sticky=tk.W)
        self.second_pass_var = tk.IntVar()
        second_pass_button = tk.Checkbutton(pipeline_frame, var=self.second_pass_var)
        second_pass_button.grid(row=row, column=1, sticky=tk.W)
        if options.do_second_pass:
            second_pass_button.select()

        row += 1

        tk.Label(pipeline_frame, text="Do 2D classification (2nd pass)?").grid(row=row, sticky=tk.W)
        self.class2d_pass2_var = tk.IntVar()
        class2d_pass2_button = tk.Checkbutton(pipeline_frame, var=self.class2d_pass2_var)
        class2d_pass2_button.grid(row=row, column=1, sticky=tk.W)
        class2d_pass2_button.select()
        if options.do_class2d_pass2:
            class2d_pass2_button.select()

        row += 1

        tk.Label(pipeline_frame, text="Do 3D classification (2nd pass)?").grid(row=row, sticky=tk.W)
        self.class3d_pass2_var = tk.IntVar()
        class3d_pass2_button = tk.Checkbutton(pipeline_frame, var=self.class3d_pass2_var)
        class3d_pass2_button.grid(row=row, column=1, sticky=tk.W)
        if options.do_class3d_pass2:
            class3d_pass2_button.select()

        ### Add logic to the box size boxes

        def calculate_box_size(particle_size_pixels):
            # Use box 20% larger than particle and ensure size is even
            box_size_exact = 1.2 * particle_size_pixels
            box_size_int = int(math.ceil(box_size_exact))
            return box_size_int + box_size_int % 2

        def calculate_downscaled_box_size(box_size_pix, angpix):
            for small_box_pix in (48, 64, 96, 128, 160, 192, 256, 288, 300, 320, 360,
                                  384, 400, 420, 450, 480, 512, 640, 768, 896, 1024):
                # Don't go larger than the original box
                if small_box_pix > box_size_pix:
                    return box_size_pix
                # If Nyquist freq. is better than 8.5 A, use this downscaled box, otherwise continue to next size up
                small_box_angpix = angpix * box_size_pix / small_box_pix
                if small_box_angpix < 4.25:
                    return small_box_pix
            # Fall back to a warning message
            return "Box size is too large!"

        def update_box_size_labels(*args_ignored, **kwargs_ignored):
            try:
                angpix = float(self.angpix_entry.get())
            except ValueError:
                # Can't update any of the labels without angpix
                self.mask_diameter_px.config(text="= NNN px")
                self.box_size_in_angstrom.config(text=u"= NNN \u212B")
                self.extract_angpix.config(text=u"= NNN \u212B/px")
                return
            try:
                mask_diameter = float(self.mask_diameter_entry.get())
                mask_diameter_px = mask_diameter / angpix
                self.mask_diameter_px.config(text="= {:.1f} px".format(mask_diameter_px))
            except (ValueError, ZeroDivisionError):
                self.mask_diameter_px.config(text="= NNN px")
                # Don't return - an error here doesn't stop us calculating the other labels
            try:
                box_size = float(self.box_size_entry.get())
                box_angpix = angpix * box_size
                self.box_size_in_angstrom.config(text=u"= {:.1f} \u212B".format(box_angpix))
            except ValueError:
                # Can't update these without the box size
                self.box_size_in_angstrom.config(text=u"= NNN \u212B")
                self.extract_angpix.config(text=u"= NNN \u212B/px")
                return
            try:
                extract_small_boxsize = float(self.extract_small_boxsize_entry.get())
                small_box_angpix = box_angpix / extract_small_boxsize
                self.extract_angpix.config(text=u"= {:.3f} \u212B/px".format(small_box_angpix))
            except (ValueError, ZeroDivisionError):
                # Can't update the downscaled pixel size unless the downscaled box size is valid
                self.extract_angpix.config(text=u"= NNN \u212B/px")

        def update_box_sizes(*args_ignored, **kwargs_ignored):
            # Always activate entry boxes - either we're activating them anyway, or we need to edit the text.
            # For text editing we need to activate the box first then deactivate again afterwards.
            self.mask_diameter_entry.config(state=tk.NORMAL)
            self.box_size_entry.config(state=tk.NORMAL)
            self.extract_small_boxsize_entry.config(state=tk.NORMAL)
            if self.get_var_as_bool(self.auto_boxsize_var):
                try:
                    particle_size_angstroms = float(self.particle_max_diam_entry.get())
                    mask_diameter = 1.1 * particle_size_angstroms
                    self.mask_diameter_entry.delete(0, tk.END)
                    self.mask_diameter_entry.insert(0, str(mask_diameter))
                    angpix = float(self.angpix_entry.get())
                    particle_size_pixels = particle_size_angstroms / angpix
                    box_size = calculate_box_size(particle_size_pixels)
                    self.box_size_entry.delete(0, tk.END)
                    self.box_size_entry.insert(0, str(box_size))
                    small_boxsize = calculate_downscaled_box_size(int(box_size), angpix)
                    self.extract_small_boxsize_entry.delete(0, tk.END)
                    self.extract_small_boxsize_entry.insert(0, str(small_boxsize))
                except:
                    # Ignore errors - they will be picked up if the user tries to save the options
                    pass
                self.mask_diameter_entry.config(state=tk.DISABLED)
                self.box_size_entry.config(state=tk.DISABLED)
                self.extract_small_boxsize_entry.config(state=tk.DISABLED)
            update_box_size_labels()

        self.box_size_var.trace('w', update_box_size_labels)
        self.extract_small_boxsize_var.trace('w', update_box_size_labels)

        self.angpix_var.trace('w', update_box_sizes)
        self.particle_max_diam_var.trace('w', update_box_sizes)
        auto_boxsize_button.config(command=update_box_sizes)

        ### Add logic to the check boxes

        def update_pipeline_control_state(*args_ignored, **kwargs_ignored):
            new_state = tk.DISABLED if self.stop_after_ctf_var.get() else tk.NORMAL
            class2d_button.config(state=new_state)
            class3d_button.config(state=new_state)
            self.particle_max_diam_entry.config(state=new_state)
            self.particle_min_diam_entry.config(state=new_state)
            self.ref_3d_entry.config(state=new_state)
            # Update the box size controls with care to avoid activating them when we shouldn't
            auto_boxsize_button.config(state=new_state)
            if new_state == tk.DISABLED:
                self.mask_diameter_entry.config(state=new_state)
                self.box_size_entry.config(state=new_state)
                self.extract_small_boxsize_entry.config(state=new_state)
            else:
                update_box_sizes()
            can_do_second_pass = (self.class3d_var.get()
                                  and len(self.ref_3d_var.get()) == 0
                                  and not self.stop_after_ctf_var.get())
            second_pass_button.config(state=tk.NORMAL if can_do_second_pass else tk.DISABLED)
            will_do_second_pass = can_do_second_pass and self.second_pass_var.get()
            class2d_pass2_button.config(state=tk.NORMAL if will_do_second_pass else tk.DISABLED)
            class3d_pass2_button.config(state=tk.NORMAL if will_do_second_pass else tk.DISABLED)

        stop_after_ctf_button.config(command=update_pipeline_control_state)
        class3d_button.config(command=update_pipeline_control_state)
        second_pass_button.config(command=update_pipeline_control_state)
        self.ref_3d_var.trace('w', update_pipeline_control_state)

        ###

        button_frame = tk.Frame(right_frame)
        button_frame.pack(padx=5, pady=5, fill=tk.X, expand=1)

        self.run_button = tk.Button(button_frame, text="Save & run", command=self.run_pipeline)
        self.run_button.pack(padx=5, pady=5, side=tk.RIGHT)

        self.save_button = tk.Button(button_frame, text="Save options", command=self.save_options)
        self.save_button.pack(padx=5, pady=5, side=tk.RIGHT)

        # Show initial pixel sizes
        update_box_sizes()

    def get_var_as_bool(self, var):
        """Helper function to convert a Tk IntVar (linked to a checkbox) to a boolean value"""
        return True if var.get() == 1 else False

    def fetch_options_from_gui(self):
        """
        Fetch the current values from the GUI widgets and store them in the options object.

        Returns:
            A list of warning messages about possible incorrect option values.

        Raises:
            ValueError: If an option value is invalid.
        """

        opts = self.options
        warnings = []

        opts.stop_after_ctf_estimation = self.get_var_as_bool(self.stop_after_ctf_var)
        opts.do_class2d = self.get_var_as_bool(self.class2d_var)
        opts.do_class3d = self.get_var_as_bool(self.class3d_var)
        opts.do_second_pass = self.get_var_as_bool(self.second_pass_var)
        opts.do_class2d_pass2 = self.get_var_as_bool(self.class2d_pass2_var)
        opts.do_class3d_pass2 = self.get_var_as_bool(self.class3d_pass2_var)

        try:
            opts.voltage = float(self.voltage_entry.get())
        except ValueError:
            raise ValueError("Voltage must be a number")
        if opts.voltage <= 0.0:
            warnings.append("- Voltage should be a positive number")

        try:
            opts.Cs = float(self.cs_entry.get())
        except ValueError:
            raise ValueError("Cs must be a number")

        opts.ctffind_do_phaseshift = self.get_var_as_bool(self.phaseplate_var)

        try:
            opts.angpix = float(self.angpix_entry.get())
        except ValueError:
            raise ValueError("Pixel size must be a number")
        if opts.angpix <= 0.0:
            warnings.append("- Pixel size should be a positive number")

        try:
            opts.motioncor_doseperframe = float(self.exposure_entry.get())
        except ValueError:
            raise ValueError("Exposure rate must be a number")
        if opts.motioncor_doseperframe <= 0.0:
            warnings.append("- Exposure rate should be a positive number")

        try:
            opts.autopick_LoG_diam_max = float(self.particle_max_diam_entry.get())
        except ValueError:
            if len(self.particle_max_diam_entry.get()) == 0 and opts.stop_after_ctf_estimation:
                # This was left blank and won't be used, set to zero to avoid errors in calculations later
                opts.autopick_LoG_diam_max = 0.0
            else:
                raise ValueError("Particle longest diameter must be a number")

        try:
            opts.autopick_LoG_diam_min = float(self.particle_min_diam_entry.get())
        except ValueError:
            if len(self.particle_min_diam_entry.get()) == 0 and opts.stop_after_ctf_estimation:
                # This was left blank and won't be used, set to zero to avoid errors in calculations later
                opts.autopick_LoG_diam_min = 0.0
            else:
                raise ValueError("Particle shortest diameter must be a number")

        opts.autopick_3dreference = self.ref_3d_entry.get()
        if len(opts.autopick_3dreference) > 0 and not os.path.isfile(opts.autopick_3dreference):
            warnings.append("- 3D reference file '{}' does not exist".format(opts.autopick_3dreference))

        try:
            opts.mask_diameter = float(self.mask_diameter_entry.get())
        except ValueError:
            raise ValueError("Mask diameter must be a number")
        if opts.mask_diameter <= 0:
            warnings.append("- Mask diameter should be a positive number")

        try:
            opts.extract_boxsize = int(self.box_size_entry.get())
        except ValueError:
            raise ValueError("Box size must be a number")
        if opts.extract_boxsize <= 0:
            warnings.append("- Box size should be a positive number")

        try:
            opts.extract_small_boxsize = int(self.extract_small_boxsize_entry.get())
            opts.extract2_small_boxsize = opts.extract_small_boxsize
            opts.extract_downscale = True
            opts.extract2_downscale = True
        except ValueError:
            raise ValueError("Down-sampled box size must be a number")
        if opts.extract_small_boxsize <= 0:
            warnings.append("- Down-sampled box size should be a positive number")

        opts.import_images = self.import_images_entry.get()
        if opts.import_images.startswith(('/', '..')):
            warnings.append("- Movies should be located inside the project directory")
        if '*' not in opts.import_images:
            warnings.append("- Pattern for input movies should normally contain a '*' to select more than one file")

        opts.motioncor_gainreference = self.gainref_entry.get()
        if len(opts.motioncor_gainreference) > 0 and not os.path.isfile(opts.motioncor_gainreference):
            warnings.append("- Gain reference file '{}' does not exist".format(opts.motioncor_gainreference))

        return warnings

    def calculate_full_options(self):
        """
        Update the options from the values that have been fetched from the GUI.

        This method uses the values that the user has set in the GUI to calculate a number of other options for the
        script.
        """
        opts = self.options

        # If we have a 3D reference, do a single pass with a large batch size
        if len(opts.autopick_3dreference) > 0:
            opts.autopick_do_LoG = False
            opts.autopick_refs_min_distance = opts.autopick_LoG_diam_max * 0.7
            opts.class3d_reference = opts.autopick_3dreference
            opts.do_second_pass = False
        else:
            # No 3D reference - do LoG autopicking in the first pass
            opts.autopick_do_LoG = True
            opts.class3d_reference = ''

        # Now set a sensible batch size (leaving batch_size_pass2 at its default 100,000)
        if opts.do_second_pass:
            opts.batch_size = 10000
        else:
            opts.batch_size = 100000

    def save_options(self):
        """
        Update the full set of options from the values in the GUI, and save them to a file.

        Returns:
            True if the options were valid and saved successfully, otherwise False.
        """
        try:
            warnings = self.fetch_options_from_gui()
            if len(warnings) == 0 or tkMessageBox.askokcancel("Warning", "\n".join(warnings), icon='warning',
                                                              default=tkMessageBox.CANCEL):
                self.calculate_full_options()
                print(" RELION_IT: Writing all options to {}".format(OPTIONS_FILE))
                if os.path.isfile(OPTIONS_FILE):
                    print(" RELION_IT: File {0} already exists; renaming old copy to {0}~".format(OPTIONS_FILE))
                    os.rename(OPTIONS_FILE, OPTIONS_FILE + '~')
                with open(OPTIONS_FILE, 'w') as optfile:
                    self.options.print_options(optfile)
                return True
        except Exception as ex:
            tkMessageBox.showerror("Error", ex.message)
            traceback.print_exc()
        return False

    def run_pipeline(self):
        """
        Update the full set of options from the values in the GUI, close the GUI and run the pipeline.
        """
        if self.save_options():
            self.main_window.destroy()
            run_pipeline(self.options)


def safe_load_star(filename, max_try=5, wait=10, expected=[]):
    for _ in range(max_try):
        try:
            star = load_star(filename)
            entry = star

            # make sure the expected key is present
            for key in expected:
               entry = entry[key]
            return star
        except:
            print("safe_load_star is retrying to read: ", filename, ", expected key:", expected)
            import time
            time.sleep(wait)
    assert False, "Failed to read a star file: " + filename

def load_star(filename):
    from collections import OrderedDict
    
    datasets = OrderedDict()
    current_data = None
    current_colnames = None
    
    in_loop = 0 # 0: outside 1: reading colnames 2: reading data

    for line in open(filename):
        line = line.strip()
        
        # remove comments
        comment_pos = line.find('#')
        if comment_pos > 0:
            line = line[:comment_pos]

        if line == "":
            if in_loop == 2:
                in_loop = 0
            continue

        if line.startswith("data_"):
            in_loop = 0

            data_name = line[5:]
            current_data = OrderedDict()
            datasets[data_name] = current_data

        elif line.startswith("loop_"):
            current_colnames = []
            in_loop = 1

        elif line.startswith("_"):
            if in_loop == 2:
                in_loop = 0

            elems = line[1:].split()
            if in_loop == 1:
                current_colnames.append(elems[0])
                current_data[elems[0]] = []
            else:
                current_data[elems[0]] = elems[1]

        elif in_loop > 0:
            in_loop = 2
            elems = line.split()
            assert len(elems) == len(current_colnames), ("Error in STAR file {}, number of elements in {} does not match number of column names {}"
                                                         .format(filename, elems, current_colnames))
            for idx, e in enumerate(elems):
                current_data[current_colnames[idx]].append(e)        
        
    return datasets

# Don't get stuck in infinite while True loops....
def CheckForExit():
    if not os.path.isfile(RUNNING_FILE):
        print(" RELION_IT:", RUNNING_FILE, "file no longer exists, exiting now ...")
        exit(0)

# Allow direct progressing to the second pass
def getSecondPassReference():
    if os.path.isfile(SECONDPASS_REF3D_FILE):
        with open(SECONDPASS_REF3D_FILE, 'r') as myfile:
            filename, angpix = myfile.readlines()
    else:
        filename = ''
        angpix = '0'
    return filename.replace('\n',''), angpix.replace('\n','')

def getJobName(name_in_script, done_file):
    jobname = None
    # See if we've done this job before, i.e. whether it is in the done_file
    if (os.path.isfile(done_file)):
        f = open(done_file,'r')
        for line in f:
            elems = line.split()
            if len(elems) < 3: continue 
            if elems[0] == name_in_script:
                jobname = elems[2]
                break
        f.close()

    return jobname

def addJob(jobtype, name_in_script, done_file, options, alias=None):
    jobname = getJobName(name_in_script, done_file)

    # If we hadn't done it before, add it now
    if (jobname is not None):
        already_had_it = True 
    else:
        already_had_it = False
        optionstring = ''
        for opt in options[:]:
            optionstring += opt + ';'

        command = 'relion_pipeliner --addJob ' + jobtype + ' --addJobOptions "' + optionstring + '"'
        if alias is not None:
                command += ' --setJobAlias "' + alias + '"'
        #print("DEBUG: Running " + command)
        os.system(command)

        pipeline = safe_load_star(PIPELINE_STAR, expected=['pipeline_processes', 'rlnPipeLineProcessName'])
        jobname = pipeline['pipeline_processes']['rlnPipeLineProcessName'][-1]
        
        # Now add the jobname to the done_file
        f = open(done_file,'a')
        f.write(name_in_script + ' = ' + jobname + '\n')
        f.close()

    # return the name of the job in the RELION pipeline, e.g. 'Import/job001/'
    return jobname, already_had_it

def RunJobs(jobs, repeat, wait, schedulename):

    runjobsstring = ''
    for job in jobs[:]:
        runjobsstring += job + ' '

    command = 'relion_pipeliner --schedule ' + schedulename + ' --repeat ' + str(repeat) + ' --min_wait ' + str(wait) + ' --RunJobs "' + runjobsstring + '" &' 

    os.system(command)

def WaitForJob(wait_for_this_job, seconds_wait):
    time.sleep(seconds_wait)
    print(" RELION_IT: waiting for job to finish in", wait_for_this_job)
    while True:
        pipeline = safe_load_star(PIPELINE_STAR, expected=['pipeline_processes', 'rlnPipeLineProcessName'])
        myjobnr = -1
        for jobnr in range(0,len(pipeline['pipeline_processes']['rlnPipeLineProcessName'])):
            jobname = pipeline['pipeline_processes']['rlnPipeLineProcessName'][jobnr]
            if jobname == wait_for_this_job:
                myjobnr = jobnr
        if myjobnr < 0:
            print(" ERROR: cannot find ", wait_for_this_job, " in ", PIPELINE_STAR)
            exit(1)

        status = int(pipeline['pipeline_processes']['rlnPipeLineProcessStatus'][myjobnr])
        if status == 2:
            print(" RELION_IT: job in", wait_for_this_job, "has finished now")
            return
        else:
            CheckForExit()
            time.sleep(seconds_wait)

def find_split_job_output(prefix, n, max_digits=6):
    import os.path
    for i in range(max_digits):
        filename = prefix + str(n).rjust(i, '0') + '.star'
        if os.path.isfile(filename):
            return filename
    return None

def writeManualPickingGuiFile(my_part_diam):
    if not os.path.isfile('.gui_manualpickrun.job'):
        with open('.gui_manualpickrun.job', 'w') as g:
            g.write("""job_type == 3
Pixel size (A) == -1
Black value: == 0
Blue value:  == 0
MetaDataLabel for color: == rlnParticleSelectZScore
Scale for CTF image: == 1
Particle diameter (A): == {}
Blue<>red color particles? == No
Highpass filter (A) == -1
Lowpass filter (A) == 20
Scale for micrographs: == 0.2
Red value:  == 2
Sigma contrast: == 3
White value: == 0
""".format(my_part_diam))

    return


def findBestClass(model_star_file, use_resol=True):
    model_star = safe_load_star(model_star_file)
    best_resol = 999
    best_size = 0 
    best_class = 0
    for iclass in range(0, len(model_star['model_classes']['rlnReferenceImage'])):
        mysize = float(model_star['model_classes']['rlnClassDistribution'][iclass])
        myresol = float(model_star['model_classes']['rlnEstimatedResolution'][iclass])
        if (not use_resol and (mysize > best_size or (mysize == best_size and myresol < best_resol))) \
        or (use_resol and (myresol < best_resol or (myresol == best_resol and mysize > best_size))):
            best_size = mysize
            best_class = model_star['model_classes']['rlnReferenceImage'][iclass]
            best_resol = myresol

    print(" RELION_IT: found best class:",best_class,"with class size of",best_size,"and resolution of",best_resol)
    return best_class, best_resol, model_star['model_general']['rlnPixelSize']

def findOutputModelStar(job_dir):
    found = None
    try:
        job_star = safe_load_star(job_dir + "job_pipeline.star", expected=['pipeline_output_edges', 'rlnPipeLineEdgeToNode'])
        for output_file in job_star["pipeline_output_edges"]['rlnPipeLineEdgeToNode']:
            if output_file.endswith("_model.star"):
                found = output_file
                break
    except:
        pass

    return found

def run_pipeline(opts):
    """
    Configure and run the RELION 3 pipeline with the given options.
    
    Args:
        opts: options for the pipeline, as a RelionItOptions object.
    """

    # if this really necessary? dont think so...
    if (os.path.isfile(PIPELINE_STAR) == False): 
        g = open(PIPELINE_STAR,'w')
        g.write('data_pipeline_general\n')
        g.write('_rlnPipeLineJobCounter 1\n')
        g.close()

    # Write RUNNING_RELION_IT file, when deleted, this script will stop
    with open(RUNNING_FILE, 'w'):
        pass

    # Write mainGUI project file, so GUI won't ask to set up a project
    with open('.gui_projectdir', 'w'):
        pass

    #### Set up GUI file for Manualpick job to allow easy viewing of autopick results
    if opts.autopick_do_LoG:
        my_part_diam = opts.autopick_LoG_diam_min
    else:
        my_part_diam = opts.autopick_refs_min_distance
    writeManualPickingGuiFile(my_part_diam)

    ### Prepare the list of queue arguments for later use
    queue_options = ['Submit to queue? == Yes',
                     'Queue name:  == {}'.format(opts.queue_name),
                     'Queue submit command: == {}'.format(opts.queue_submit_command),
                     'Standard submission script: == {}'.format(opts.queue_submission_template),
                     'Minimum dedicated cores per node: == {}'.format(opts.queue_minimum_dedicated)]

    # If we're only doing motioncorr and ctf estimation, then forget about the second pass and the batch processing
    if opts.stop_after_ctf_estimation:
        opts.do_class2d = False
        opts.do_class3d = False
        opts.do_second_pass = False

    if opts.do_second_pass:
        nr_passes = 2
    else:
        nr_passes = 1

    # if SECONDPASS_REF3D_FILE exists, go straight into the second pass
    first_pass = 0
    if opts.do_second_pass:
        secondpass_ref3d, secondpass_ref3d_angpix = getSecondPassReference()
        if not secondpass_ref3d == '':
            print(' RELION_IT: found', secondpass_ref3d,'with angpix=',secondpass_ref3d_angpix,'as a 3D reference for second pass in file',SECONDPASS_REF3D_FILE)
            print(' RELION_IT: if the automatic selection of the reference turned out to be unsatisfactory,')
            print(' RELION_IT: you can re-run the second pass with another reference by:')
            print(' RELION_IT:  stopping the pipeline by deleting RUNNING_*')
            print(' RELION_IT:  updating the reference filename in',SECONDPASS_REF3D_FILE)
            print(' RELION_IT:  deleting relevant jobs (autopick2_job and followings) in',SETUP_CHECK_FILE)
            print(' RELION_IT:  and restarting the pipeline.')
            first_pass = 1
            opts.autopick_3dreference = secondpass_ref3d
            opts.autopick_ref_angpix = secondpass_ref3d_angpix
            opts.autopick_2dreferences = ''
            opts.autopick_do_LoG = False
            opts.class3d_reference = secondpass_ref3d
            opts.have_3d_reference = True

    # Allow to perform two passes through the entire pipeline (PREPROCESS and CLASS2D/3D batches)
    # The second pass, a 3D reference generated in the first pass will be used for template-based autopicking
    for ipass in range(first_pass, nr_passes):


        #### Set up the Import job
        import_options = ['Raw input files: == {}'.format(opts.import_images),
                          'Import raw movies/micrographs? == Yes',
	                  'Pixel size (Angstrom): == {}'.format(opts.angpix),
	                  'Voltage (kV): == {}'.format(opts.voltage),
                          'Spherical aberration (mm): == {}'.format(opts.Cs),
                          'Amplitude contrast: == {}'.format(opts.ampl_contrast)]

        if opts.images_are_movies:
            import_options.append('Are these multi-frame movies? == Yes')
        else:
            import_options.append('Are these multi-frame movies? == No')

        import_job, already_had_it  = addJob('Import','import_job', SETUP_CHECK_FILE, import_options)

        if opts.images_are_movies:
            #### Set up the MotionCor job
            motioncorr_options = ['Input movies STAR file: == {}movies.star'.format(import_job),
                                  'MOTIONCOR2 executable: == {}'.format(opts.motioncor_exe),
                                  'Defect file: == {}'.format(opts.motioncor_defectfile),
                                  'Gain-reference image: == {}'.format(opts.motioncor_gainreference),
                                  'Gain flip: == {}'.format(opts.motioncor_gainflip),
                                  'Gain rotation: == {}'.format(opts.motioncor_gainrot),
                                  'Do dose-weighting? == Yes',
                                  'Dose per frame (e/A2): == {}'.format(opts.motioncor_doseperframe),
                                  'Number of patches X: ==  {}'.format(opts.motioncor_patches_x),
                                  'Number of patches Y: == {}'.format(opts.motioncor_patches_y),
                                  'Bfactor: ==  {}'.format(opts.motioncor_bfactor),
                                  'Binning factor: == {}'.format(opts.motioncor_binning),
                                  'Which GPUs to use: == {}'.format(opts.motioncor_gpu),
                                  'Other MOTIONCOR2 arguments == {}'.format(opts.motioncor_other_args),
                                  'Number of threads: == {}'.format(opts.motioncor_threads),
                                  'Number of MPI procs: == {}'.format(opts.motioncor_mpi),
                                  'Additional arguments: == --eer_upsampling {} --eer_grouping {}'.format(opts.eer_upsampling, opts.eer_grouping)]

            if (opts.motioncor_do_own):
                motioncorr_options.append('Use RELION\'s own implementation? == Yes')
                if opts.use_ctffind_instead:
                    motioncorr_options.append('Save sum of power spectra? == Yes')
                else:
                    motioncorr_options.append('Save sum of power spectra? == No')
            else:
                motioncorr_options.append('Use RELION\'s own implementation? == No')

            if opts.motioncor_submit_to_queue:
                motioncorr_options.extend(queue_options)

            motioncorr_job, already_had_it  = addJob('MotionCorr', 'motioncorr_job', SETUP_CHECK_FILE, motioncorr_options)


        #### Set up the CtfFind job
        ctffind_options = ['Amount of astigmatism (A): == {}'.format(opts.ctffind_astigmatism),
                           'FFT box size (pix): == {}'.format(opts.ctffind_boxsize),
                           'Maximum defocus value (A): == {}'.format(opts.ctffind_defocus_max),
                           'Minimum defocus value (A): == {}'.format(opts.ctffind_defocus_min),
                           'Defocus step size (A): == {}'.format(opts.ctffind_defocus_step),
                           'Maximum resolution (A): == {}'.format(opts.ctffind_maxres),
                           'Minimum resolution (A): == {}'.format(opts.ctffind_minres),
                           'Gctf executable: == {}'.format(opts.gctf_exe),
                           'Which GPUs to use: == {}'.format(opts.gctf_gpu),
                           'CTFFIND-4.1 executable: == {}'.format(opts.ctffind4_exe),
                           'Number of MPI procs: == {}'.format(opts.ctffind_mpi)]

        if opts.images_are_movies:
            ctffind_options.append('Input micrographs STAR file: == {}{}'.format(motioncorr_job, 'corrected_micrographs.star'))
        else:
            ctffind_options.append('Input micrographs STAR file: == {}{}'.format(import_job, 'micrographs.star'))

        if opts.use_ctffind_instead:
            ctffind_options.append('Use CTFFIND-4.1? == Yes')
            ctffind_options.append('Use Gctf instead? == No')
            ctffind_options.append('Use power spectra from MotionCorr job? == Yes')
        else:
            ctffind_options.append('Use CTFFIND-4.1? == No')
            ctffind_options.append('Use Gctf instead? == Yes')
            ctffind_options.append('Use power spectra from MotionCorr job? == No')
            if (opts.ctffind_do_ignore_search_params):
                ctffind_options.append('Ignore \'Searches\' parameters? == Yes')
            else:
                ctffind_options.append('Ignore \'Searches\' parameters? == No')
            if (opts.ctffind_do_EPA):
                 ctffind_options.append('Perform equi-phase averaging? == Yes')
            else:
                 ctffind_options.append('Perform equi-phase averaging? == No')

        if opts.ctffind_do_phaseshift:
            ctffind_options.append('Estimate phase shifts? == Yes')
        else:
            ctffind_options.append('Estimate phase shifts? == No')

        if opts.ctffind_submit_to_queue:
            ctffind_options.extend(queue_options)

        ctffind_job, already_had_it  = addJob('CtfFind', 'ctffind_job', SETUP_CHECK_FILE, ctffind_options)

        runjobs = [import_job]
        if opts.images_are_movies:
            runjobs.append(motioncorr_job)
        runjobs.append(ctffind_job)

        # There is an option to stop on-the-fly processing after CTF estimation
        if not opts.stop_after_ctf_estimation:
            autopick_options = ['Input micrographs for autopick: == {}micrographs_ctf.star'.format(ctffind_job),
                                'Min. diameter for LoG filter (A) == {}'.format(opts.autopick_LoG_diam_min),
                                'Max. diameter for LoG filter (A) == {}'.format(opts.autopick_LoG_diam_max),
                                'Maximum resolution to consider (A) == {}'.format(opts.autopick_lowpass),
                                'Adjust default threshold (stddev): == {}'.format(opts.autopick_LoG_adjust_threshold),
                                'Upper threshold (stddev): == {}'.format(opts.autopick_LoG_upper_threshold),
                                '2D references: == {}'.format(opts.autopick_2dreferences),
                                '3D reference: == {}'.format(opts.autopick_3dreference),
                                'Symmetry: == {}'.format(opts.autopick_3dref_symmetry),
                                'Pixel size in references (A) == {}'.format(opts.autopick_ref_angpix),
                                '3D angular sampling: == {}'.format(opts.autopick_3dref_sampling),
                                'In-plane angular sampling (deg) == {}'.format(opts.autopick_inplane_sampling),
                                'Picking threshold: == {}'.format(opts.autopick_refs_threshold),
                                'Minimum inter-particle distance (A): == {}'.format(opts.autopick_refs_min_distance),
                                'Mask diameter (A) == {}'.format(opts.autopick_refs_mask_diam),
                                'Maximum stddev noise: == {}'.format(opts.autopick_stddev_noise),
                                'Minimum avg noise: == {}'.format(opts.autopick_avg_noise),
                                'Shrink factor: == {}'.format(opts.autopick_shrink_factor),
                                'Which GPUs to use: == {}'.format(opts.autopick_gpu),
                                'Additional arguments: == {}'.format(opts.autopick_other_args),
                                'Number of MPI procs: == {}'.format(opts.autopick_mpi)]

            if not opts.autopick_3dreference == '':
                autopick_options.append('OR: provide a 3D reference? == Yes')
            else:
                autopick_options.append('OR: provide a 3D reference? == No')

            if opts.autopick_do_LoG:
                autopick_options.append('OR: use Laplacian-of-Gaussian? == Yes')
            else:
                autopick_options.append('OR: use Laplacian-of-Gaussian? == No')

            if opts.autopick_refs_are_ctf_corrected:
                autopick_options.append('Are References CTF corrected? == Yes')
            else:
                autopick_options.append('Are References CTF corrected? == No')

            if opts.autopick_refs_have_inverted_contrast:
                autopick_options.append('References have inverted contrast? == Yes')
            else:
                autopick_options.append('References have inverted contrast? == No')

            if opts.autopick_refs_ignore_ctf1stpeak:
                autopick_options.append('Ignore CTFs until first peak? == Yes')
            else:
                autopick_options.append('Ignore CTFs until first peak? == No')

            if opts.autopick_do_gpu and (not opts.autopick_do_LoG):
                autopick_options.append('Use GPU acceleration? == Yes')
            else:
                autopick_options.append('Use GPU acceleration? == No')

            if opts.autopick_submit_to_queue:
                autopick_options.extend(queue_options)
            
            if ipass == 0:
                autopick_job_name = 'autopick_job'
                autopick_alias = 'pass 1'
            else:
                autopick_job_name = 'autopick2_job'
                autopick_alias = 'pass 2'

            autopick_job, already_had_it  = addJob('AutoPick', autopick_job_name, SETUP_CHECK_FILE, autopick_options, alias=autopick_alias)
            runjobs.append(autopick_job)

            #### Set up the Extract job
            extract_options = ['Input coordinates:  == {}coords_suffix_autopick.star'.format(autopick_job),
                               'micrograph STAR file:  == {}micrographs_ctf.star'.format(ctffind_job),
                               'Diameter background circle (pix):  == {}'.format(opts.extract_bg_diameter),
                               'Particle box size (pix): == {}'.format(opts.extract_boxsize),
                               'Number of MPI procs: == {}'.format(opts.extract_mpi)]

            if ipass == 0:
                if opts.extract_downscale:
                    extract_options.append('Rescale particles? == Yes')
                    extract_options.append('Re-scaled size (pixels):  == {}'.format(opts.extract_small_boxsize))
            else:
                if opts.extract2_downscale:
                    extract_options.append('Rescale particles? == Yes')
                    extract_options.append('Re-scaled size (pixels):  == {}'.format(opts.extract2_small_boxsize))

            if opts.extract_submit_to_queue:
                extract_options.extend(queue_options)

            if ipass == 0:
                extract_job_name = 'extract_job'
                extract_alias = 'pass 1'
            else:
                extract_job_name = 'extract2_job'
                extract_alias = 'pass 2'

            extract_job, already_had_it  = addJob('Extract', extract_job_name, SETUP_CHECK_FILE, extract_options, alias=extract_alias)
            runjobs.append(extract_job)

            if (ipass == 0 and (opts.do_class2d or opts.do_class3d)) or (ipass == 1 and (opts.do_class2d_pass2 or opts.do_class3d_pass2)):
                #### Set up the Select job to split the particle STAR file into batches
                split_options = ['OR select from particles.star: == {}particles.star'.format(extract_job),
                                 'OR: split into subsets? == Yes',
                                 'OR: number of subsets:  == -1']

                if ipass == 0:
                    split_job_name = 'split_job'
                    split_options.append('Subset size:  == {}'.format(opts.batch_size))
                    split_alias = 'into {}'.format(opts.batch_size) 
                else:
                    split_job_name = 'split2_job'
                    split_options.append('Subset size:  == {}'.format(opts.batch_size_pass2))
                    split_alias = 'into {}'.format(opts.batch_size_pass2)

                split_job, already_had_it = addJob('Select', split_job_name, SETUP_CHECK_FILE, split_options, alias=split_alias)

                # Now start running stuff
                runjobs.append(split_job)

        # Now execute the entire preprocessing pipeliner
        if ipass == 0:
            preprocess_schedule_name = PREPROCESS_SCHEDULE_PASS1
        else:
            preprocess_schedule_name = PREPROCESS_SCHEDULE_PASS2
        RunJobs(runjobs, opts.preprocess_repeat_times, opts.preprocess_repeat_wait, preprocess_schedule_name)
        print(' RELION_IT: submitted',preprocess_schedule_name,'pipeliner with', opts.preprocess_repeat_times,'repeats of the preprocessing jobs')
        print(' RELION_IT: this pipeliner will run in the background of your shell. You can stop it by deleting the file RUNNING_PIPELINER_'+preprocess_schedule_name)


        ########## From now on, process extracted particles in batches for 2D or 3D classification, only perform SGD inimodel for first batch and if no 3D reference is available

        # There is again an option to stop here...
        if (ipass == 0 and (opts.do_class2d or opts.do_class3d)) or (ipass == 1 and (opts.do_class2d_pass2 or opts.do_class3d_pass2)):

            ### If necessary, rescale the 3D reference in the second pass!
            # TODO: rescale initial reference if different from movies?
            if ipass == 1 and (opts.extract_downscale or opts.extract2_downscale):
                particles_angpix = opts.angpix
                if opts.images_are_movies:
                    particles_angpix = particles_angpix * opts.motioncor_binning
                if opts.extract2_downscale:
                    particles_angpix = particles_angpix * opts.extract_boxsize / opts.extract2_small_boxsize
                    particles_boxsize = opts.extract2_small_boxsize
                else:
                    particles_boxsize = opts.extract_boxsize
                if abs(float(particles_angpix) - float(opts.autopick_ref_angpix)) > 0.01:
                    # Now rescale the reference for 3D classification
                    opts.class3d_reference = opts.autopick_3dreference.replace('.mrc','_rescaled.mrc')
                    print(' RELION_IT: rescaling the 3D reference from pixel size',opts.autopick_ref_angpix,'to',particles_angpix,'and saving the new reference as',opts.class3d_reference)
                    command = 'relion_image_handler --i ' + opts.autopick_3dreference + ' --o ' + opts.class3d_reference + ' --angpix ' + str(opts.autopick_ref_angpix) + ' --rescale_angpix ' + str(particles_angpix) + ' --new_box ' + str(particles_boxsize)
                    os.system(command)


            print(' RELION_IT: now entering an infinite loop for batch-processing of particles. You can stop this loop by deleting the file', RUNNING_FILE)
            
            # It could be that this is a restart, so check previous_batch1_size in the output directory.
            # Also check the presence of class2d_job_batch_001 in case the first job was not submitted yet.
            first_split_file = find_split_job_output(split_job + 'particles_split', 1)
            if getJobName("class2d_job_batch_001", SETUP_CHECK_FILE) is not None and \
               first_split_file is not None:
                batch1 = safe_load_star(first_split_file, expected=['particles', 'rlnMicrographName'])
                previous_batch1_size = len(batch1['particles']['rlnMicrographName'])
            else:
                previous_batch1_size = 0

            continue_this_pass = True
            while continue_this_pass:

                have_new_batch = False
                nr_batches = len(glob.glob(split_job + "particles_split*.star"))
                for ibatch  in range(0, nr_batches):
                    iibatch = ibatch + 1
                    batch_name = find_split_job_output(split_job + "particles_split", iibatch)

                    batch = safe_load_star(batch_name, expected=['particles', 'rlnMicrographName'])
                    batch_size = len(batch['particles']['rlnMicrographName'])
                    rerun_batch1 = False
                    if ( iibatch == 1 and batch_size > previous_batch1_size and batch_size > opts.minimum_batch_size ):
                        previous_batch1_size = batch_size
                        rerun_batch1 = True

                    particles_star_file = batch_name

                    # The first batch is special: perform 2D classification with smaller batch size (but at least minimum_batch_size) and keep overwriting in the same output directory
                    if ( rerun_batch1 or batch_size == opts.batch_size):


                        # Discard particles with odd average/stddev values
                        if opts.do_discard_on_image_statistics:

                            #### Run a Select job to get rid of particles with outlier average/stddev values...
                            discard_options = ['OR select from particles.star: == {}'.format(batch_name),
                                               'OR: select on image statistics? == Yes',
                                               'Sigma-value for discarding images: == {}'.format(opts.discard_sigma),
                                               'Metadata label for images: == rlnImageName']

                            if ipass == 0:
                                discard_job_name = 'discard_job'
                            else:
                                discard_job_name = 'discard2_job'

                            if opts.discard_submit_to_queue:
                                discard_options.extend(queue_options)

                            discard_job, already_had_it = addJob('Select', discard_job_name, SETUP_CHECK_FILE, discard_options)

                            if ((not already_had_it) or rerun_batch1):
                                have_new_batch = True
                                RunJobs([discard_job], 1, 1, 'DISCARD')
                                print(" RELION_IT: submitted job to discard based on image statistics for", batch_size ,"particles in", batch_name)

                                # Wait here until this Discard job is finished. Check every thirty seconds
                                WaitForJob(discard_job, 30)

                            particles_star_file = discard_job + 'particles.star'


                        # 2D classification
                        if (ipass == 0 and opts.do_class2d) or (ipass == 1 and opts.do_class2d_pass2):

                            class2d_options = ['Input images STAR file: == {}'.format(particles_star_file),
                                               'Number of classes: == {}'.format(opts.class2d_nr_classes),
                                               'Mask diameter (A): == {}'.format(opts.mask_diameter),
                                               'Number of iterations: == {}'.format(opts.class2d_nr_iter),
                                               'Angular search range - psi (deg): == {}'.format(opts.class2d_angle_step), 
                                               'Offset search range (pix): == {}'.format(opts.class2d_offset_range),
                                               'Offset search step (pix): == {}'.format(opts.class2d_offset_step),
                                               'Number of pooled particles: == {}'.format(opts.refine_nr_pool),
                                               'Which GPUs to use: == {}'.format(opts.refine_gpu),
                                               'Number of MPI procs: == {}'.format(opts.refine_mpi),
                                               'Number of threads: == {}'.format(opts.refine_threads),
                                               'Copy particles to scratch directory: == {}'.format(opts.refine_scratch_disk),
                                               'Additional arguments: == {}'.format(opts.class2d_other_args)]

                            if batch_size > opts.refine_batchsize_for_fast_subsets:
                                class2d_options.append('Use fast subsets (for large data sets)? == Yes')
                            else:
                                class2d_options.append('Use fast subsets (for large data sets)? == No')

                            if opts.refine_do_gpu:
                                class2d_options.append('Use GPU acceleration? == Yes')
                            else:
                                class2d_options.append('Use GPU acceleration? == No')

                            if opts.class2d_ctf_ign1stpeak:
                                class2d_options.append('Ignore CTFs until first peak? == Yes')
                            else:
                                class2d_options.append('Ignore CTFs until first peak? == No')

                            if opts.refine_preread_images:
                                class2d_options.append('Pre-read all particles into RAM? == Yes')
                            else:
                                class2d_options.append('Pre-read all particles into RAM? == No')

                            if opts.refine_submit_to_queue:
                                class2d_options.extend(queue_options)

                            if ipass == 0:
                                jobname = 'class2d_job_batch_{:03d}'.format(iibatch)
                                alias = 'pass1_batch_{:03d}'.format(iibatch)
                            else:
                                jobname = 'class2d_pass2_job_batch_{:03d}'.format(iibatch)
                                alias = 'pass2_batch_{:03d}'.format(iibatch)

                            class2d_job, already_had_it = addJob('Class2D', jobname, SETUP_CHECK_FILE, class2d_options, alias=alias)
 
                            if ((not already_had_it) or rerun_batch1):
                                have_new_batch = True
                                RunJobs([class2d_job], 1, 1, 'CLASS2D')
                                print(" RELION_IT: submitted 2D classification with", batch_size ,"particles in", class2d_job)

                                # Wait here until this Class2D job is finished. Check every thirty seconds
                                WaitForJob(class2d_job, 30)

                    # Perform 3D classification
                    if (ipass == 0 and opts.do_class3d) or (ipass == 1 and opts.do_class3d_pass2):

                        # Do SGD initial model generation only in the first pass, when no reference is provided AND only for the first (complete) batch, for subsequent batches use that model
                        if (not opts.have_3d_reference) and ipass == 0 and iibatch == 1 and batch_size == opts.batch_size:

                            inimodel_options = ['Input images STAR file: == {}'.format(particles_star_file),
                                                'Symmetry: == {}'.format(opts.symmetry),
                                                'Mask diameter (A): == {}'.format(opts.mask_diameter),
                                                'Number of classes: == {}'.format(opts.inimodel_nr_classes),
                                                'Initial angular sampling: == {}'.format(opts.inimodel_angle_step), 
                                                'Offset search range (pix): == {}'.format(opts.inimodel_offset_range),
                                                'Offset search step (pix): == {}'.format(opts.inimodel_offset_step),
                                                'Number of initial iterations: == {}'.format(opts.inimodel_nr_iter_initial), 
                                                'Number of in-between iterations: == {}'.format(opts.inimodel_nr_iter_inbetween), 
                                                'Number of final iterations: == {}'.format(opts.inimodel_nr_iter_final),
                                                'Write-out frequency (iter): == {}'.format(opts.inimodel_freq_writeout), 
                                                'Initial resolution (A): == {}'.format(opts.inimodel_resol_ini), 
                                                'Final resolution (A): == {}'.format(opts.inimodel_resol_final), 
                                                'Initial mini-batch size: == {}'.format(opts.inimodel_batchsize_ini), 
                                                'Final mini-batch size: == {}'.format(opts.inimodel_batchsize_final), 
                                                'Increased noise variance half-life: == {}'.format(opts.inimodel_sigmafudge_halflife), 
                                                'Number of pooled particles: == 1',
                                                'Which GPUs to use: == {}'.format(opts.refine_gpu),
                                                'Number of MPI procs: == {}'.format(opts.refine_mpi),
                                                'Number of threads: == {}'.format(opts.refine_threads),
                                                'Copy particles to scratch directory: == {}'.format(opts.refine_scratch_disk),
                                                'Additional arguments: == {}'.format(opts.inimodel_other_args)]

                            if opts.inimodel_solvent_flatten:
                                inimodel_options.append('Flatten and enforce non-negative solvent? == Yes')
                            else:    
                                inimodel_options.append('Flatten and enforce non-negative solvent? == No')

                            if opts.refine_skip_padding:
                                inimodel_options.append('Skip padding? == Yes')
                            else:    
                                inimodel_options.append('Skip padding? == No')

                            if opts.refine_do_gpu:
                                inimodel_options.append('Use GPU acceleration? == Yes')
                            else:
                                inimodel_options.append('Use GPU acceleration? == No')

                            if opts.inimodel_ctf_ign1stpeak:
                                inimodel_options.append('Ignore CTFs until first peak? == Yes')
                            else:
                                inimodel_options.append('Ignore CTFs until first peak? == No')

                            if opts.refine_preread_images:
                                inimodel_options.append('Pre-read all particles into RAM? == Yes')
                            else:
                                inimodel_options.append('Pre-read all particles into RAM? == No')

                            if opts.refine_submit_to_queue:
                                inimodel_options.extend(queue_options)

                            inimodel_job, already_had_it = addJob('InitialModel', 'inimodel', SETUP_CHECK_FILE, inimodel_options)              

                            if (not already_had_it):
                                have_new_batch = True
                                RunJobs([inimodel_job], 1, 1, 'INIMODEL')
                                print(" RELION_IT: submitted initial model generation with", batch_size ,"particles in", inimodel_job)

                                # Wait here until this inimodel job is finished. Check every thirty seconds
                                WaitForJob(inimodel_job, 30)

                            sgd_model_star = findOutputModelStar(inimodel_job)
                            if sgd_model_star is None:
                                print(" RELION_IT: Initial model generation " + inimodel_job + " does not contain expected output maps.")
                                print(" RELION_IT: This job should have finished, but you may continue it from the GUI.")
                                raise Exception("ERROR!! quitting the pipeline.") # TODO: MAKE MORE ROBUST

                            # Use the model of the largest class for the 3D classification below
                            total_iter = opts.inimodel_nr_iter_initial + opts.inimodel_nr_iter_inbetween + opts.inimodel_nr_iter_final
                            best_inimodel_class, best_inimodel_resol, best_inimodel_angpix = findBestClass(sgd_model_star, use_resol=True)
                            opts.class3d_reference = best_inimodel_class
                            opts.class3d_ref_is_correct_greyscale = True
                            opts.class3d_ref_is_ctf_corrected = True
                            opts.have_3d_reference = True


                        if opts.have_3d_reference:
                            # Now perform the actual 3D classification
                            class3d_options = ['Input images STAR file: == {}'.format(particles_star_file),
                                               'Reference map: == {}'.format(opts.class3d_reference),
                                               'Initial low-pass filter (A): == {}'.format(opts.class3d_ini_lowpass),
                                               'Symmetry: == {}'.format(opts.symmetry),
                                               'Regularisation parameter T: == {}'.format(opts.class3d_T_value),
                                               'Reference mask (optional): == {}'.format(opts.class3d_reference_mask),
                                               'Number of classes: == {}'.format(opts.class3d_nr_classes),
                                               'Mask diameter (A): == {}'.format(opts.mask_diameter),
                                               'Number of iterations: == {}'.format(opts.class3d_nr_iter),
                                               'Angular sampling interval: == {}'.format(opts.class3d_angle_step), 
                                               'Offset search range (pix): == {}'.format(opts.class3d_offset_range),
                                               'Offset search step (pix): == {}'.format(opts.class3d_offset_step),
                                               'Number of pooled particles: == {}'.format(opts.refine_nr_pool),
                                               'Which GPUs to use: == {}'.format(opts.refine_gpu),
                                               'Number of MPI procs: == {}'.format(opts.refine_mpi),
                                               'Number of threads: == {}'.format(opts.refine_threads),
                                               'Copy particles to scratch directory: == {}'.format(opts.refine_scratch_disk),
                                               'Additional arguments: == {}'.format(opts.class3d_other_args)]

                            if batch_size > opts.refine_batchsize_for_fast_subsets:
                                class3d_options.append('Use fast subsets (for large data sets)? == Yes')
                            else:
                                class3d_options.append('Use fast subsets (for large data sets)? == No')
                            
                            if opts.class3d_ref_is_correct_greyscale:
                                class3d_options.append('Ref. map is on absolute greyscale? == Yes')
                            else:
                                class3d_options.append('Ref. map is on absolute greyscale? == No')

                            if opts.class3d_ref_is_ctf_corrected:
                                class3d_options.append('Has reference been CTF-corrected? == Yes')
                            else:
                                class3d_options.append('Has reference been CTF-corrected? == No')

                            if opts.refine_skip_padding:
                                class3d_options.append('Skip padding? == Yes')
                            else:    
                                class3d_options.append('Skip padding? == No')

                            if opts.refine_do_gpu:
                                class3d_options.append('Use GPU acceleration? == Yes')
                            else:
                                class3d_options.append('Use GPU acceleration? == No')

                            if opts.class3d_ctf_ign1stpeak:
                                class3d_options.append('Ignore CTFs until first peak? == Yes')
                            else:
                                class3d_options.append('Ignore CTFs until first peak? == No')

                            if opts.refine_preread_images:
                                class3d_options.append('Pre-read all particles into RAM? == Yes')
                            else:
                                class3d_options.append('Pre-read all particles into RAM? == No')

                            if opts.refine_submit_to_queue:
                                class3d_options.extend(queue_options)

                            if ipass == 0:
                                jobname = 'class3d_job_batch_{:03d}'.format(iibatch)
                                alias = 'pass1_batch_{:03d}'.format(iibatch)
                            else:
                                jobname = 'class3d2_job_batch_{:03d}'.format(iibatch)
                                alias = 'pass2_batch_{:03d}'.format(iibatch)

                            class3d_job, already_had_it = addJob('Class3D', jobname, SETUP_CHECK_FILE, class3d_options, alias=alias) 

                            if ((not already_had_it) or rerun_batch1):
                                have_new_batch = True
                                RunJobs([class3d_job], 1, 1, 'CLASS3D')
                                print(' RELION_IT: submitted 3D classification with', batch_size ,'particles in', class3d_job)

                                # Wait here until this Class2D job is finished. Check every thirty seconds
                                WaitForJob(class3d_job, 30)

                            class3d_model_star = findOutputModelStar(class3d_job)
                            if class3d_model_star is None:
                                print(" RELION_IT: 3D Classification " + class3d_job + " does not contain expected output maps.")
                                print(" RELION_IT: This job should have finished, but you may continue it from the GUI.")
                                raise Exception("ERROR!! quitting the pipeline.") # TODO: MAKE MORE ROBUST

                            best_class3d_class, best_class3d_resol, best_class3d_angpix = findBestClass(class3d_model_star, use_resol=True)

                            # Once the first batch in the first pass is completed: move on to the second pass
                            if (ipass == 0 and opts.do_second_pass and iibatch == 1 and best_class3d_resol < opts.minimum_resolution_3dref_2ndpass):
                                opts.autopick_3dreference = best_class3d_class
                                opts.autopick_ref_angpix = best_class3d_angpix
                                opts.autopick_2dreferences = ''
                                opts.autopick_do_LoG = False
                                opts.class3d_reference = best_class3d_class
                                opts.have_3d_reference = True
                                opts.autopick_3dref_symmetry = opts.symmetry

                                # Stop the PREPROCESS pipeliner of the first pass by removing its RUNNING file
                                filename_to_remove = 'RUNNING_PIPELINER_'+preprocess_schedule_name
                                if os.path.isfile(filename_to_remove):
                                    print(' RELION_IT: removing file',filename_to_remove,'to stop the pipeliner from the first pass')
                                    os.remove(filename_to_remove)

                                # Generate a file to indicate we're in the second pass, so that restarts of the python script will be smooth
                                g = open(SECONDPASS_REF3D_FILE,'w')
                                g.write(str(best_class3d_class)+'\n'+str(best_class3d_angpix)+'\n')
                                g.close()

                                # Move out of this ipass of the passes loop....
                                ibatch = nr_batches+1
                                continue_this_pass = False
                                print(' RELION_IT: moving on to the second pass using',opts.autopick_3dreference,'for template-based autopicking')
                                # break out of the for-loop over the batches
                                break


                if not have_new_batch:
                    CheckForExit()
                    # The following prevents checking the particles.star file too often
                    time.sleep(60*opts.batch_repeat_time)


def main():
    """
    Run the RELION 3 pipeline.
    
    Options files given as command line arguments will be opened in order and
    used to update the default options.
    """
    # Start by parsing arguments
    # (If --help is given, the program will print a usage message and exit)
    parser = argparse.ArgumentParser()
    parser.add_argument("extra_options", nargs="*", metavar="extra_options.py",
                        help="Python files containing options for relion_it.py")
    parser.add_argument("--gui", action="store_true", help="launch a simple GUI to set options")
    parser.add_argument("--continue", action="store_true", dest="continue_",
                        help="continue a previous run by loading options from ./relion_it_options.py")
    args = parser.parse_args()

    print(' RELION_IT: -------------------------------------------------------------------------------------------------------------------')
    print(' RELION_IT: script for automated, on-the-fly single-particle analysis in RELION (>= 3.1)')
    print(' RELION_IT: authors: Sjors H.W. Scheres, Takanori Nakane & Colin M. Palmer')
    print(' RELION_IT: ')
    print(' RELION_IT: usage: ./relion_it.py [extra_options.py [extra_options2.py ....] ] [--gui] [--continue]')
    print(' RELION_IT: ')
    print(' RELION_IT: this script will check whether processes are still running using files with names starting with RUNNING')
    print(' RELION_IT:   you can restart this script after stopping previous processes by deleting all RUNNING files')
    print(' RELION_IT: this script keeps track of already submitted jobs in a filed called',SETUP_CHECK_FILE)
    print(' RELION_IT:   upon a restart, jobs present in this file will be continued (for preprocessing), or ignored when already finished')
    print(' RELION_IT: if you would like to re-do a specific job from scratch (e.g. because you changed its parameters)')
    print(' RELION_IT:   remove that job, and those that depend on it, from the',SETUP_CHECK_FILE)
    print(' RELION_IT: -------------------------------------------------------------------------------------------------------------------')
    print(' RELION_IT: ')
    
    # Make sure no other version of this script are running...
    if os.path.isfile(RUNNING_FILE):
        print(" RELION_IT: ERROR:", RUNNING_FILE, "is already present: delete this file and make sure no other copy of this script is running. Exiting now ...")
        exit(0)

    # Also make sure the preprocessing pipeliners are stopped before re-starting this script
    for checkfile in ('RUNNING_PIPELINER_'+PREPROCESS_SCHEDULE_PASS1, 'RUNNING_PIPELINER_'+PREPROCESS_SCHEDULE_PASS2):
        if os.path.isfile(checkfile):
            print(" RELION_IT: ERROR:", checkfile, "is already present: delete this file and make sure no relion_pipeliner job is still running. Exiting now ...")
            exit(0)

    if args.continue_:
        print(' RELION_IT: continuing a previous run. Options will be loaded from ./relion_it_options.py')
        args.extra_options.append(OPTIONS_FILE)

    opts = RelionItOptions()
    for user_opt_file in args.extra_options:
        print(' RELION_IT: reading options from {}'.format(user_opt_file))
        user_opts = runpy.run_path(user_opt_file)
        opts.update_from(user_opts)

    if args.gui:
        print(' RELION_IT: launching GUI...')
        tk_root = tk.Tk()
        tk_root.title("relion_it.py setup")
        RelionItGui(tk_root, opts)
        tk_root.mainloop()
    else:
        run_pipeline(opts)


if __name__ == "__main__":
    main()