def test_wwinp_to_h5m():
    thisdir = os.path.dirname(__file__)
    wwinp = os.path.join(thisdir, "files_test_wwinp_to_h5m/wwinp_test.e")
    output = os.path.join(os.getcwd(), "wwinp_mesh.h5m")

    if output in os.listdir("."):
        os.remove(output)

    wwinp_to_h5m.cartesian(wwinp, output)

    with open(output) as f:
        written = f.read()

    expected_sm = ScdMesh.fromFile(os.path.join(thisdir, "files_test_wwinp_to_h5m/expected_wwinp_mesh.h5m"))
    written_sm = ScdMesh.fromFile(output)

    for x in range(0, 5):
        for y in range(0, 6):
            for z in range(0, 10):
                expected_voxel = expected_sm.getHex(x, y, z)
                expected = expected_sm.imesh.getTagHandle("ww_n_group_001")[expected_voxel]
                written_voxel = written_sm.getHex(x, y, z)
                written = written_sm.imesh.getTagHandle("ww_n_group_001")[written_voxel]
                assert_equal(written, expected)

    os.remove(output)

def test_wwinp_to_h5m():
    wwinp = 'files_test_wwinp_to_h5m/wwinp_test.txt'
    output = 'wwinp_mesh.h5m'
    
    if output in os.listdir('.'):
        os.remove(output)
    
    wwinp_to_h5m.cartesian(wwinp, output)
    
    with open(output) as f:
        written = f.read()
    
    expected_sm = ScdMesh.fromFile('files_test_wwinp_to_h5m/expected_wwinp_mesh.h5m')
    written_sm = ScdMesh.fromFile(output)
    
    for x in range(0,5):
        for y in range(0,6):
            for z in range(0,10):
                expected_voxel = expected_sm.getHex(x,y,z)
                expected = expected_sm.imesh.getTagHandle('ww_n_group_001')[expected_voxel]
                written_voxel = written_sm.getHex(x,y,z)
                written = written_sm.imesh.getTagHandle('ww_n_group_001')[written_voxel]
                assert_equal(written, expected)
    
    os.remove(output)

    def test_get_divs(self):
        x = [1, 2.5, 4, 6.9]
        y = [-12, -10, -.5]
        z = [100, 200]

        sm = ScdMesh( x, y, z )

        self.assertEqual( sm.getDivisions('x'), x )
        self.assertEqual( sm.getDivisions('y'), y )
        self.assertEqual( sm.getDivisions('z'), z )

    def test_large_iterator(self):

        print "building large mesh"
        big = ScdMesh(range(1,100), range(101,200), range(201,300))
        print "iterating (1)"
        for i in big.iterateHex():
            pass
        print "iterating (2)"
        for i in big.iterateHex( 'yzx' ):
            pass

    def test_get_hex(self):
        # mesh with valid i values 0-4, j values 0-3, k values 0-2
        sm = ScdMesh( range(11,16), range(21,25), range(31,34) )
        def check( e ):
            self.assertTrue( isinstance(e, iBase.Entity) )
        check(sm.getHex(0, 0, 0))
        check(sm.getHex(1, 1, 1))
        check(sm.getHex(3, 0, 0))
        check(sm.getHex(3, 2, 1))

        self.assertRaises(ScdMeshError, sm.getHex, -1,-1,-1)
        self.assertRaises(ScdMeshError, sm.getHex, 4, 0, 0)
        self.assertRaises(ScdMeshError, sm.getHex, 0, 3, 0)
        self.assertRaises(ScdMeshError, sm.getHex, 0, 0, 2)

    def test_get_vtx(self):
        # mesh with valid i values 0-4, j values 0-3, k values 0-2
        x_range = numpy.array(range(10,15),dtype=numpy.float64)
        y_range = numpy.array(range(21,24),dtype=numpy.float64)
        z_range = numpy.array(range(31,33),dtype=numpy.float64)

        sm = ScdMesh( x_range, y_range, z_range, self.mesh ) 

        for i,x in enumerate(x_range):
            for j,y in enumerate(y_range):
                for k,z in enumerate(z_range):
                    vtx = sm.getVtx(i,j,k)
                    vcoord = self.mesh.getVtxCoords( vtx )
                    self.assertTrue( all( vcoord == [x,y,z]) )

def main():
    """Method defines an option parser and handles command-line
    usage of this module.

    Notes
    -----
    Requires command line arguments to be passed - otherwise prints help
    information.
    """

    usage = "usage: %prog input-h5m-file [options] arg"
    parser = OptionParser(usage)
    
    # Input and mesh file names
    parser.add_option("-o","--output",action="store",dest="output", \
            default="gammas",help="Option specifies the name of the 'gammas'" \
            "file. Default: %default")
    #
   
    (options, args) = parser.parse_args()

    # Create ScdMesh object, which also loads 'meshfile' into mesh.
    sm = ScdMesh.fromFile(args[0])

    gen_gammas_file_from_h5m(sm, options.output)

    return 1

def main():
    """ACTION: Method defines an option parser and handles command-line
    usage of this module.
    REQUIRES: command line arguments to be passed - otherwise prints help
    information.
    """

    usage = "usage: %prog ENERGY_FILE MESH_FILE [options] \n\n" \
            "ENERGY_FILE is a list of the energy bins for each photon energy " \
            "group, with a single energy per line. MESH_FILE is the MOAB " \
            "mesh that will store the contents of ENERGY_FILE in the tag " \
            "'PHTN_ERGS'."
    parser = OptionParser(usage)

    (options, args) = parser.parse_args()

    fr = open(args[0])
    sm = ScdMesh.fromFile(args[1])

    # Call the method to read fr and tag mesh
    read_and_tag_phtn_ergs(fr, sm)

    mesh.save(args[1])
    fr.close()

    return 1

def main( arguments = None ):

    #Instatiate options parser
    parser = OptionParser\
             (usage='%prog <ww mesh> [options]')

    parser.add_option('-o', dest='output_name', default='wwinp.out',\
        help='Name of WWINP output file, default=%default')

    parser.add_option('-t', action='store_true', dest='totals_bool',\
        default=False, \
        help='If multiple energy groups exist, only use Total \
         default=%default')

    (opts, args) = parser.parse_args( arguments )

    if len(args) != 1:
        parser.error\
        ( '\nNeed exactly 1 argument: <ww mesh>' )

    # load mesh
    ww_mesh = ScdMesh.fromFile(args[0])


    write_wwinp(ww_mesh, opts.totals_bool, opts.output_name)

    print "\tWrote WWINP file '{0}'".format(opts.output_name)

    print "Complete"

def test_h5m_to_wwinp_3D_n():
    thisdir = os.path.dirname(__file__)
    ww_sm_filename = os.path.join(thisdir, 'files_test_wwinp_to_h5m/expected_ww_mesh_3D_n.h5m')
    ww_sm = ScdMesh.fromFile(ww_sm_filename)
    output = os.path.join(os.getcwd(), 'test.e')
    expected_output = os.path.join(thisdir, 'files_test_wwinp_to_h5m/3D_n.e')
    
    if output in os.listdir('.'):
        os.remove(output)
    
    totals_bool = False 
    h5m_to_wwinp.write_wwinp(ww_sm, totals_bool, output)
    
    with open(output) as f:
        written = f.readlines()

    with open(expected_output) as f:
        expected = f.readlines()

    # check to make sure the first line is the same except for the date
    assert_equal(written[0].split()[:-2], expected[0].split()[:-2])

    # check to make sure files are the same length
    assert_equal(len(written), len(expected))

    # check to make sure the rest of the lines have the same values
    # since the number formats are different, float comparisons are used
    for i in range(1, len(expected)):
        for j in range(0, len(expected[i].split())):
            assert_equal(float(written[i].split()[j]), float(expected[i].split()[j]))

    os.remove(output)

def test_wwinp_to_h5m_3D_n():
    thisdir = os.path.dirname(__file__)
    wwinp = os.path.join(thisdir, 'files_test_wwinp_to_h5m/3D_n.e')
    expected_file = os.path.join(thisdir, 'files_test_wwinp_to_h5m/expected_ww_mesh_3D_n.h5m')

    written_sm = wwinp_to_h5m.cartesian(wwinp)
    expected_sm = ScdMesh.fromFile(expected_file)
    
    #verify weight window lower bounds are the same
    for x in range(0,14):
        for y in range(0,8):
            for z in range(0,6):
                for e_group in range(1, 8):
                    expected_voxel = expected_sm.getHex(x,y,z)
                    expected = expected_sm.imesh.getTagHandle('ww_n_group_00{0}'.format(e_group))[expected_voxel]
                    written_voxel = written_sm.getHex(x,y,z)
                    written = written_sm.imesh.getTagHandle('ww_n_group_00{0}'.format(e_group))[written_voxel]
                    assert_equal(written, expected)

    #verify correct particle identifier
    assert_equal(written_sm.imesh.getTagHandle('particle')[written_sm.imesh.rootSet], 1)

    #verify correct energy upper bounds
    expected_E = [1E-9, 1E-8, 1E-7, 1E-6, 1E-5, 1E-4, 1E-3]
    written_E = written_sm.imesh.getTagHandle('E_upper_bounds')[written_sm.imesh.rootSet]
    
    for i in range(0, len(expected_E)):
        assert_equal(written_E[i], expected_E[i])

def test_wwinp_to_h5m_1D_p():
    thisdir = os.path.dirname(__file__)
    wwinp = os.path.join(thisdir, 'files_test_wwinp_to_h5m/1D_p.e')
    expected_file = os.path.join(thisdir, 'files_test_wwinp_to_h5m/expected_ww_mesh_1D_p.h5m')

    written_sm = wwinp_to_h5m.cartesian(wwinp)
    expected_sm = ScdMesh.fromFile(expected_file)
 
    #verify weight window lower bounds are the same
    for x in range(0,1):
        for y in range(0,1):
            for z in range(0,9):
                expected_voxel = expected_sm.getHex(x,y,z)
                expected = expected_sm.imesh.getTagHandle('ww_n_group_001')[expected_voxel]
                written_voxel = written_sm.getHex(x,y,z)
                written = written_sm.imesh.getTagHandle('ww_n_group_001')[written_voxel]
                assert_equal(written, expected)

    #verify correct particle identifier
    assert_equal(written_sm.imesh.getTagHandle('particle')[written_sm.imesh.rootSet], 1)

    #verify correct energy upper bounds
    expected_E = 100
    written_E = written_sm.imesh.getTagHandle('E_upper_bounds')[written_sm.imesh.rootSet]

    assert_equal(written_E, expected_E)

def main( arguments = None ) :

    # Instantiate option parser
    parser = OptionParser\
             (usage='%prog <meshtal_file> <normalization_factor> [options]')

    parser.add_option('-o', dest='mesh_output', default='flux_mesh.h5m',\
                      help = 'Name of mesh output file, default=%default.\
                             For meshtal files with multiple tallies,\
                             if the -o flag is used all tallies must be named,\
                             with file names seperated by commas and no spaces\
                             (e.g. "tally14.h5m,tally24.h5m,tally34.h5m")')
    parser.add_option('-n', dest='norm', default=None,
                      help='Normalization factor, default=%default,\
                            For meshtal files with multiple tallies, if the -n\
                            flag is used, a normalization factor must be\
                            specified for all tallies, seperated by commas but \
                            not spaces (eg. -n 1.1,2.2,3.3) ')
    parser.add_option('-m', dest='smesh_filename', default=None,
                      help='Preexisting mesh on which to tag fluxes')
                         

    (opts, args) = parser.parse_args(arguments)
    
    #if len(args) != 2 :
     #   parser.error('\nNeed 1 argument: meshtal file')
    print "\n\nRunning read_meshtal.py"
    tally_numbers, tally_lines = find_tallies(args[1])
    print "Number of tallies found: {0}\nTally number(s): {1}" \
                                     .format(len(tally_numbers), tally_numbers)

    # Parse input from options parser, generate default values
    if opts.norm :
        norm = opts.norm.split(',')
    else :
        norm = [1]*len(tally_numbers)

    if opts.mesh_output !='flux_mesh.h5m' :
        mesh_output = opts.mesh_output.split(',')
    else:
        mesh_output = []
        for n in range(0, len(tally_numbers)) :
            if len(tally_numbers) == 1 :
                mesh_output.append('flux_mesh.h5m')
            else :
                mesh_output.append('flux_mesh_tally{0}.h5m'.format(tally_numbers[n]))

    # Convert each tally to h5m and name accordingly
    for n in range(0,len(tally_numbers)) :
        print "\nNow parsing tally number {0}".format(tally_numbers[n])
        if opts.smesh_filename:
            alt_sm = ScdMesh.fromFile(opts.smesh_filename)
            sm = read_meshtal(args[1], tally_lines[n], float(norm[n]), smesh=alt_sm)
        else:
            sm = read_meshtal(args[1], tally_lines[n],float(norm[n]))
        sm.scdset.save(mesh_output[n])

        print "\tSaved tally {0} as {1}".format(tally_numbers[n], mesh_output[n])
    print "\nStructured mesh tagging complete\n\n"

def create_ww_mesh(flux_mesh, e_group_names):
    """
    This function reads a flux mesh and returns a weight window mesh tagged with 
    all zeros in every energy group. This is used for the intial MAGIC weight
    window generation, for which no preexisting weight window mesh is given.

    Parameters
    ----------
    flux_mesh : ScdMesh
        A ScdMesh tagged with fluxes in the form X_group_YYY and or
        X_group_total. Addition required tags are "particle" (1 for n, 2 for p) 
        and E_group_bounds (vector of energy upper bounds).
    e_group_names : vector of energy names
        In the form X_group_YYY or X_group_total.
        
    """

    ww_mesh = ScdMesh(flux_mesh.getDivisions('x'),\
                       flux_mesh.getDivisions('y'),\
                       flux_mesh.getDivisions('z'))

    # create ww tags
    for e_group_name in e_group_names:
            ww_tag = ww_mesh.imesh.createTag('ww_{0}'.format(e_group_name), 1, float)
            for voxel in ww_mesh.iterateHex('xyz'):
                ww_tag[voxel] = 0

    # create e_upper_bound tags
    e_upper_bounds = \
        flux_mesh.imesh.getTagHandle("E_upper_bounds")[flux_mesh.imesh.rootSet]
    if isinstance(e_upper_bounds, float):
        e_upper_bounds = [e_upper_bounds]

    e_tag = ww_mesh.imesh.createTag("E_upper_bounds", len(e_upper_bounds), float)
    e_tag[ww_mesh.imesh.rootSet] = \
        flux_mesh.imesh.getTagHandle("E_upper_bounds")[flux_mesh.imesh.rootSet]

    # create  particle tag
    particle = flux_mesh.imesh.getTagHandle("particle")[flux_mesh.imesh.rootSet]
    particle_tag = ww_mesh.imesh.createTag("particle", 1, int)
    particle_tag[ww_mesh.imesh.rootSet] = particle

    return ww_mesh

    def setUp(self):
        self.mesh = iMesh.Mesh()
        self.sm = ScdMesh( range(10,15), # i = 0,1,2,3
                           range(21,25), # j = 0,1,2
                           range(31,34), # k = 0,1
                           self.mesh )

        self.I = range(0,4)
        self.J = range(0,3)
        self.K = range(0,2)

def write_magic(flux_mesh_filename, ww_inp_mesh_filename, totals_bool, \
                null_value, output_mesh, tolerance):
    """
    This function takes the filename of the flux mesh and optional wieght window
    mesh, as well necessary parameter, sends them to the magic funtion and
    writes the resulting wieght window mesh.

    Parameters
    ----------
    flux_mesh_filename: ScdMesh file name.
        A ScdMesh tagged with fluxes in the form X_group_YYY and or
        X_group_total. Addition required tags are "particle" (1 for n, 2 for p) 
        and E_group_bounds (vector of energy upper bounds).
    totals_bool : True or False
        Determines whether magic will be applied to individual energy group
        (False), or the total group (True)
    null_value : float
        The weight window lower bound value that is assigned to voxels where the
        relative error on flux exceeds the tolerance. This is only done for 
        initial weight window lower bound generation, not subsequent iterations.
    output_mesh : string
        Filename of output mesh
    tolerance: float
        The maximum relative error allowable for the MAGIC algorithm to create
        a weight window lower bound for for a given voxel for the intial weight
        window lower bound generation, or overwrite preexisting weight window
        lower bounds for subsequent iterations.  
    ww_inp_mesh_filename : ScdMesh file name
        A preexisting weight window mesh to apply MAGIC to.      
    """

    flux_mesh = ScdMesh.fromFile(flux_mesh_filename)

    if ww_inp_mesh_filename != None:
        ww_inp_mesh = ScdMesh.fromFile(ww_inp_mesh_filename)

    else:
        ww_inp_mesh = None
   

    ww_mesh = magic(flux_mesh, totals_bool, null_value, tolerance, ww_inp_mesh)
    ww_mesh.scdset.save(output_mesh)
    print "\tWrote WW mesh file '{0}'".format(output_mesh)

def magic(flux_h5m, ww_mesh, total_bool, null_value, output, output_mesh, tolerance):
    """Runs magic.py from as a module
    """
    flux_mesh = ScdMesh.fromFile(flux_h5m)

    ww_mesh, e_groups = magic_wwinp(flux_mesh, ww_mesh, total_bool, null_value, tolerance)

    if output_mesh != 'None':
        ww_mesh.scdset.save(output_mesh)

    write_wwinp(ww_mesh, e_groups, output)

def main():
    """Method defines an option parser and handles command-line
    usage of this module.
    
    Notes
    -----
    Requires command line arguments to be passed - otherwise prints help
    information.
    """

    usage = "usage: %prog [options] arg"
    parser = OptionParser(usage)
    
    # Input and mesh file names
    parser.add_option("-p","--phtn",action="store",dest="phtnsrcfile", \
            default=False,help="The photon source strengths are read from" \
            "FILENAME.")
    parser.add_option("-m","--mesh",action="store",dest="meshfile", \
            default="",help="file to write source information to, or" \
            " file name for saving a modified mesh.")
    # Other options
    parser.add_option("-i","--isotope",action="store",dest="isotope", \
            default="TOTAL",help="The isotope string identifier or 'TOTAL'. "\
            "Default: %default")
    parser.add_option("-c","--coolingstep",action="store",dest="coolingstep", \
            default=0,help="The cooling step number or string identifier. " \
            "(0 is first cooling step)  Default: %default")
    parser.add_option("-r","--retag",action="store_true",dest="retag", \
            default=False,help="Option enables retagging of .h5m meshes. " \
            "Default: %default")
    parser.add_option("-t","--totals",action="store_true",dest="totals", \
            default=False,help="Option enables adding the total photon " \
            "source strength for all energy groups as a tag for each voxel. " \
            "Default: %default")

    (options, args) = parser.parse_args()

    # Open an ScdMesh and then call read_to_h5m
    try:
        mesh = ScdMesh.fromFile(options.meshfile)
    except ScdMeshError:
        mesh = iMesh.Mesh()
        mesh.load(options.meshfile)

    read_to_h5m( \
                options.phtnsrcfile, mesh, options.isotope, \
                options.coolingstep, options.retag, options.totals)

    if isinstance(mesh, ScdMesh):
        mesh.imesh.save(options.meshfile)
    else:
        mesh.save(options.meshfile)

    return 1

def main( arguments = None ) :

   #Instantiate option parser
    parser = OptionParser\
             (usage='%prog <meshtal_file> <normalization_factor> [options]')

    parser.add_option('-o', dest='mesh_output', default=None,\
                      help = 'Name of mesh output file, default=%default')
    parser.add_option('-n', dest='norm', default=None,
                      help = 'Normalization factor, default=%default')
    parser.add_option('-m', dest='smesh_filename', default=None,
                      help='Preexisting mesh on which to tag fluxes')
                         

    (opts, args) = parser.parse_args(arguments)
    
    #if len(args) != 2 :
     #   parser.error('\nNeed 1 argument: meshtal file')
    print "\n\nRunning read_meshtal.py"
    tally_numbers, tally_lines = find_tallies(args[1])
    print 'Number of tallies found: {0}\nTally number(s): {1}'\
                                     .format(len(tally_numbers), tally_numbers)

    # Parse input from options parser, generate default values
    if opts.norm :
        norm = opts.norm.split(',')
    else :
        norm=[1]*len(tally_numbers)

    if opts.mesh_output :
        mesh_output = opts.mesh_output.split(',')
    else:
        mesh_output=[]
        for n in range(0, len(tally_numbers)) :
            if len(tally_numbers) == 1 :
                mesh_output.append('flux_mesh.h5m')
            else :
                mesh_output.append('flux_mesh_tally{0}.h5m'.format(tally_numbers[n]))
    

    # Convert each tally to h5m and name accordingly
    for n in range(0,len(tally_numbers)) :
        print "\nNow parsing tally number {0}".format(tally_numbers[n])
        if opts.smesh_filename:
            alt_sm = ScdMesh.fromFile(opts.smesh_filename)
            sm = read_meshtal(args[1], tally_lines[n], float(norm[n]), smesh=alt_sm)
        else:
            sm = read_meshtal(args[1], tally_lines[n],float(norm[n]))

        sm.scdset.save(mesh_output[n])

        print "\tSaved tally {0} as {1}".format(tally_numbers[n], mesh_output[n])
    print '\nStructured mesh tagging complete\n\n'

def test_magic_it_0_1_group():
    thisdir = os.path.dirname(__file__)
    flux_sm_filename = os.path.join(thisdir, "files_test_magic/iteration_0_flux_1_group.h5m")
    flux_sm = ScdMesh.fromFile(flux_sm_filename)
    expected_sm_filename = os.path.join(thisdir, "files_test_magic/iteration_0_magic_1_group.h5m")
    expected_sm = ScdMesh.fromFile(expected_sm_filename)

    totals_bool = False
    null_value = 1e-3
    tolerance = 0.2

    written_sm = magic.magic(flux_sm, totals_bool, null_value, tolerance)

    # verify weight window lower bounds are the same
    for x in range(0, 3):
        for y in range(0, 3):
            for z in range(0, 3):
                expected_voxel = expected_sm.getHex(x, y, z)
                expected = expected_sm.imesh.getTagHandle("ww_n_group_001")[expected_voxel]
                written_voxel = written_sm.getHex(x, y, z)
                written = written_sm.imesh.getTagHandle("ww_n_group_001")[written_voxel]
                assert_equal(written, expected)