from numpy import *
import os

try:
    import ctypes
except:
    print '''
ERROR: ctypes not installed properly. (ctypes needed for opt_height_field)
        '''
    import sys
    sys.exit()

def platform_shared_library_suffix():
    import sys
    result = 'so'
    if 'win' in sys.platform.lower(): result = 'dll'
    ## No else if, because we want darwin to override win (which is a substring of darwin)
    if 'darwin' in sys.platform.lower(): result = 'dylib'
    return result

def to_ctypes_array( sequence, ctypes_type ):
    '''
    Returns a contiguous ctypes array suitable for extracting
    a pointer (with .ctypes.data_as() or ctypes_array_ref()) to pass to a C function.
    '''
    return ascontiguousarray( asarray( sequence, dtype = ctypes_type ) )
def ctypes_array_ref( ctypes_array, ctypes_type ):
    '''
    Returns a ctypes.POINTER to the internal data of 'ctypes_array' of type 'ctypes_type'.
    
    WARNING: Do not chain the calls to to_ctypes_array() and ctypes_array_ref() as in
               ctypes_array_ref( to_ctypes_array( sequence, ctypes_type ), ctypes_type ).
             ctypes_array_ref returns a pointer to the internal numpy storage of the array and
             cannot keep a reference to the array, so it will be garbage collected if these
             functions are chained together.
             Instead, keep a temporary reference to the result of to_ctypes_array():
               contig = to_ctypes_array( sequence, ctypes_type )
               ctypes_array_ref( contig, ctypes_type )
    '''
    return ctypes_array.ctypes.data_as( ctypes.POINTER( ctypes_type ) )

lib = ctypes.cdll.LoadLibrary( os.path.join( '.', 'libtoposmooth2.' + platform_shared_library_suffix() ) )

real_t = ctypes.c_double
size_t = ctypes.c_long
class c_triangle_t( ctypes.Structure ):
    _fields_ = [ ( "v", ctypes.c_long * 3 ) ]
class c_edge_t( ctypes.Structure ):
    _fields_ = [ ( "v", ctypes.c_long * 2 ) ]

__smooth_mesh2d_monotonically_decreasing = lib.smooth_mesh2d_monotonically_decreasing_d
__smooth_mesh2d_monotonically_decreasing.argtypes = \
    [
        size_t,
        ctypes.POINTER( real_t ),
        
        size_t,
        ctypes.POINTER( c_triangle_t ),
        
        ctypes.POINTER( real_t ),
        ctypes.POINTER( real_t ),
        
        size_t,
        ctypes.POINTER( c_edge_t ),
        
        real_t,
        ctypes.c_int
    ]
__smooth_mesh2d_monotonically_decreasing.restype = None
def smooth_mesh2d_monotonically_decreasing( vertices_xy, triangles, values, directed_edges, stay_close_to_original_values_weight = 0., smooth_deformation_not_smooth_function = False ):
    '''
    Smoothes values associated with vertices in a 2D triangle mesh.
    
    Vertices are ordered from 0 to N-1.
    A vertex's 2D position is given by looking up its index in the 'vertices_xy' array,
    which is an N-by-2 array.
    
    Connectivity information is given by the 'triangles' array.
    Each element in the array is a tuple containing three vertex indices.
    
    A vertex's function value is given by the corresponding value in the array
    'original_values'.
    This function returns new values guaranteed to be
    monotonically decreasing, and yet smooth and close to 'original_values'.
    
    The directed edges of the monotonicity graph are given by the 'directed_edges' array.
    Each element in the array is a tuple of two indices for the edge's start and end vertex indices.
    NOTE: Sources and sinks are considered to be the fixed nodes.
    
    The optional parameter 'stay_close_to_original_values_weight' determines how strongly
    the output should try to stay close to the original values.
    A value <= 0. means that they should not try to stay close at all, and the resulting
    values will be purely smoothed.
    The default value is 0.
    
    The optional parameter 'smooth_deformation_not_smooth_function', if true, means that
    the smoothing term attempts to find a smooth deformation of the original values to
    make them monotonic, rather than finding actually smooth values.
    The default value is False.
    '''
    
    from itertools import izip
    
    vertices_xy = to_ctypes_array( vertices_xy, real_t )
    
    c_triangles = ( c_triangle_t * len( triangles ) )()
    for c_tri, tri in izip( c_triangles, triangles ):
        c_tri.v[:] = tri[:]
    
    values = to_ctypes_array( values, real_t )
    
    c_edges = ( c_edge_t * len( directed_edges ) )()
    for c_edge, edge in izip( c_edges, directed_edges ):
        c_edge.v[:] = edge[:]
    
    new_values_out = to_ctypes_array( zeros( values.shape ), real_t )
    
    __smooth_mesh2d_monotonically_decreasing(
        
        len( vertices_xy ),
        ctypes_array_ref( vertices_xy, real_t ),
        
        len( triangles ),
        c_triangles,
        
        ctypes_array_ref( values, real_t ),
        ctypes_array_ref( new_values_out, real_t ),
        
        len( directed_edges ),
        c_edges,
        
        stay_close_to_original_values_weight,
        smooth_deformation_not_smooth_function
        )
    
    return new_values_out


__prepare_smooth_mesh2d_monotonically_decreasing_d = lib.prepare_smooth_mesh2d_monotonically_decreasing_d
__prepare_smooth_mesh2d_monotonically_decreasing_d.argtypes = \
    [
        size_t,
        ctypes.POINTER( real_t ),
        
        size_t,
        ctypes.POINTER( c_triangle_t ),
        
        ctypes.POINTER( real_t ),
        
        size_t,
        ctypes.POINTER( c_edge_t ),
        
        real_t,
        ctypes.c_int
    ]
__prepare_smooth_mesh2d_monotonically_decreasing_d.restype = ctypes.c_void_p
def prepare_smooth_mesh2d_monotonically_decreasing( vertices_xy, triangles, values, directed_edges, stay_close_to_original_values_weight = 0., smooth_deformation_not_smooth_function = False ):
    '''
    Similar to smooth_mesh2d_monotonically_decreasing(), except instead of returning the optimized
    values, returns an object suitable for passing to
    other *_smooth_mesh2d_monotonically_decreasing_d() functions.
    '''
    from itertools import izip
    
    vertices_xy = to_ctypes_array( vertices_xy, real_t )
    
    c_triangles = ( c_triangle_t * len( triangles ) )()
    for c_tri, tri in izip( c_triangles, triangles ):
        c_tri.v[:] = tri[:]
    
    values = to_ctypes_array( values, real_t )
    
    c_edges = ( c_edge_t * len( directed_edges ) )()
    for c_edge, edge in izip( c_edges, directed_edges ):
        c_edge.v[:] = edge[:]
    
    prepared = __prepare_smooth_mesh2d_monotonically_decreasing_d(
        
        len( vertices_xy ),
        ctypes_array_ref( vertices_xy, real_t ),
        
        len( triangles ),
        c_triangles,
        
        ctypes_array_ref( values, real_t ),
        
        len( directed_edges ),
        c_edges,
        
        stay_close_to_original_values_weight,
        smooth_deformation_not_smooth_function
        )
    
    return prepared


__release_smooth_mesh2d_monotonically_decreasing_d = lib.release_smooth_mesh2d_monotonically_decreasing_d
__release_smooth_mesh2d_monotonically_decreasing_d.argtypes = \
    [
        ctypes.c_void_p
    ]
__release_smooth_mesh2d_monotonically_decreasing_d.restype = None
def release_smooth_mesh2d_monotonically_decreasing( prepared ):
    '''
    Releases the object returned from prepare_smooth_mesh2d_monotonically_decreasing().
    The object may no longer be used.
    '''
    
    __release_smooth_mesh2d_monotonically_decreasing_d( prepared )


__num_dofs_smooth_mesh2d_monotonically_decreasing_d = lib.num_dofs_smooth_mesh2d_monotonically_decreasing_d
__num_dofs_smooth_mesh2d_monotonically_decreasing_d.argtypes = \
    [
        ctypes.c_void_p
    ]
__num_dofs_smooth_mesh2d_monotonically_decreasing_d.restype = ctypes.c_long
def num_dofs_smooth_mesh2d_monotonically_decreasing( prepared ):
    return int( __num_dofs_smooth_mesh2d_monotonically_decreasing_d( prepared ) )


__num_values_smooth_mesh2d_monotonically_decreasing_d = lib.num_values_smooth_mesh2d_monotonically_decreasing_d
__num_values_smooth_mesh2d_monotonically_decreasing_d.argtypes = \
    [
        ctypes.c_void_p
    ]
__num_values_smooth_mesh2d_monotonically_decreasing_d.restype = ctypes.c_long
def num_values_smooth_mesh2d_monotonically_decreasing( prepared ):
    return int( __num_values_smooth_mesh2d_monotonically_decreasing_d( prepared ) )


__initial_dofs_smooth_mesh2d_monotonically_decreasing_d = lib.initial_dofs_smooth_mesh2d_monotonically_decreasing_d
__initial_dofs_smooth_mesh2d_monotonically_decreasing_d.argtypes = \
    [
        ctypes.c_void_p
    ]
__initial_dofs_smooth_mesh2d_monotonically_decreasing_d.restype = ctypes.POINTER( real_t )

__release_initial_dofs_smooth_mesh2d_monotonically_decreasing_d = lib.release_initial_dofs_smooth_mesh2d_monotonically_decreasing_d
__release_initial_dofs_smooth_mesh2d_monotonically_decreasing_d.argtypes = \
    [
        ctypes.c_void_p,
        ctypes.POINTER( real_t )
    ]
__release_initial_dofs_smooth_mesh2d_monotonically_decreasing_d.restype = None
def initial_dofs_smooth_mesh2d_monotonically_decreasing( prepared ):
    
    initial = __initial_dofs_smooth_mesh2d_monotonically_decreasing_d( prepared )
    
    result = zeros( num_dofs_smooth_mesh2d_monotonically_decreasing( prepared ), dtype = real_t )
    result[ : len(result) ] = initial[ : len(result) ]
    
    __release_initial_dofs_smooth_mesh2d_monotonically_decreasing_d( prepared, initial )
    
    return result


__values_smooth_mesh2d_monotonically_decreasing_d = lib.values_smooth_mesh2d_monotonically_decreasing_d
__values_smooth_mesh2d_monotonically_decreasing_d.argtypes = \
    [
        ctypes.c_void_p,
        ctypes.POINTER( real_t )
    ]
__values_smooth_mesh2d_monotonically_decreasing_d.restype = ctypes.POINTER( real_t )

__release_values_smooth_mesh2d_monotonically_decreasing_d = lib.release_values_smooth_mesh2d_monotonically_decreasing_d
__release_values_smooth_mesh2d_monotonically_decreasing_d.argtypes = \
    [
        ctypes.c_void_p,
        ctypes.POINTER( real_t )
    ]
__release_values_smooth_mesh2d_monotonically_decreasing_d.restype = None
def values_smooth_mesh2d_monotonically_decreasing( prepared, x ):
    
    x = to_ctypes_array( x, real_t )
    
    values = __values_smooth_mesh2d_monotonically_decreasing_d(
        prepared,
        ctypes_array_ref( x, real_t )
        )
    
    result = zeros( num_values_smooth_mesh2d_monotonically_decreasing( prepared ), dtype = real_t )
    result[ : len(result) ] = values[ : len(result) ]
    
    __release_values_smooth_mesh2d_monotonically_decreasing_d( prepared, values )
    
    return result


__func_smooth_mesh2d_monotonically_decreasing_d = lib.func_smooth_mesh2d_monotonically_decreasing_d
__func_smooth_mesh2d_monotonically_decreasing_d.argtypes = \
    [
        ctypes.c_void_p,
        ctypes.POINTER( real_t )
    ]
__func_smooth_mesh2d_monotonically_decreasing_d.restype = real_t
def func_smooth_mesh2d_monotonically_decreasing( prepared, x ):
    
    x = to_ctypes_array( x, real_t )
    
    return float( __func_smooth_mesh2d_monotonically_decreasing_d(
        prepared,
        ctypes_array_ref( x, real_t ),
        ) )


__gradient_smooth_mesh2d_monotonically_decreasing_d = lib.gradient_smooth_mesh2d_monotonically_decreasing_d
__gradient_smooth_mesh2d_monotonically_decreasing_d.argtypes = \
    [
        ctypes.c_void_p,
        ctypes.POINTER( real_t )
    ]
__gradient_smooth_mesh2d_monotonically_decreasing_d.restype = ctypes.POINTER( real_t )

__release_gradient_smooth_mesh2d_monotonically_decreasing_d = lib.release_gradient_smooth_mesh2d_monotonically_decreasing_d
__release_gradient_smooth_mesh2d_monotonically_decreasing_d.argtypes = \
    [
        ctypes.c_void_p,
        ctypes.POINTER( real_t )
    ]
__release_gradient_smooth_mesh2d_monotonically_decreasing_d.restype = None
def gradient_smooth_mesh2d_monotonically_decreasing( prepared, x ):
    
    x = to_ctypes_array( x, real_t )
    
    g = __gradient_smooth_mesh2d_monotonically_decreasing_d(
        prepared,
        ctypes_array_ref( x, real_t ),
        )
    
    result = zeros( len(x), dtype = real_t )
    result[ : len(result) ] = g[ : len(result) ]
    
    __release_gradient_smooth_mesh2d_monotonically_decreasing_d( prepared, g )
    
    return result


class c_hess_t( ctypes.Structure ):
    _fields_ = [
        ( "num_entries", ctypes.c_long ),
        ( "row", ctypes.POINTER( ctypes.c_long ) ),
        ( "col", ctypes.POINTER( ctypes.c_long ) ),
        ( "val", ctypes.POINTER( real_t ) )
        ]
__hessian_smooth_mesh2d_monotonically_decreasing_d = lib.hessian_smooth_mesh2d_monotonically_decreasing_d
__hessian_smooth_mesh2d_monotonically_decreasing_d.argtypes = \
    [
        ctypes.c_void_p,
        ctypes.POINTER( real_t )
    ]
__hessian_smooth_mesh2d_monotonically_decreasing_d.restype = c_hess_t

__release_hessian_smooth_mesh2d_monotonically_decreasing_d = lib.release_hessian_smooth_mesh2d_monotonically_decreasing_d
__release_hessian_smooth_mesh2d_monotonically_decreasing_d.argtypes = \
    [
        ctypes.c_void_p,
        c_hess_t
    ]
__release_hessian_smooth_mesh2d_monotonically_decreasing_d.restype = None
def hessian_smooth_mesh2d_monotonically_decreasing( prepared, x ):
    
    x = to_ctypes_array( x, real_t )
    
    hess = __hessian_smooth_mesh2d_monotonically_decreasing_d(
        prepared,
        ctypes_array_ref( x, real_t ),
        )
    
    nnz = int( hess.num_entries )
    result = array( hess.row[:nnz] ), array( hess.col[:nnz] ), array( hess.val[:nnz] )
    
    __release_hessian_smooth_mesh2d_monotonically_decreasing_d( prepared, hess )
    
    return result


class smoother_mesh2d_monotonically_decreasing( object ):
    def __init__( self, *args, **kwargs ):
        self.smoother = prepare_smooth_mesh2d_monotonically_decreasing( *args, **kwargs )
    def __del__( self ):
        ## If we are called as the program is terminating, global functions seem to be set to None.
        if release_smooth_mesh2d_monotonically_decreasing is not None:
            release_smooth_mesh2d_monotonically_decreasing( self.smoother )
    
    def num_dofs( self ):
        return num_dofs_smooth_mesh2d_monotonically_decreasing( self.smoother )
    
    def num_values( self ):
        return num_values_smooth_mesh2d_monotonically_decreasing( self.smoother )
    
    def initial_dofs( self ):
        return initial_dofs_smooth_mesh2d_monotonically_decreasing( self.smoother )
    
    def values( self, dofs ):
        return values_smooth_mesh2d_monotonically_decreasing( self.smoother, dofs )
    
    def f( self, dofs ):
        return func_smooth_mesh2d_monotonically_decreasing( self.smoother, dofs )
    
    def g( self, dofs ):
        return gradient_smooth_mesh2d_monotonically_decreasing( self.smoother, dofs )
    
    def h( self, dofs ):
        return hessian_smooth_mesh2d_monotonically_decreasing( self.smoother, dofs )


############ Test code below ############

def interior_tree_on_mesh( mesh, exclude_vertices = [] ):
    '''
    Given a TriMesh 'mesh' and a collection of vertices to skip,
    returns a tree on the mesh as the tuple:
        (
            root vertex index in mesh.vs,
            list of tuples of vertices indicating directed edges,
            list of all vertex indices (into mesh.vs) contained in the tree
            ).
    The returned tree will contain only interior vertices (no boundary vertices).
    
    NOTE: If no interior tree exists, the result will be ( None, [], [] ).
    
    untested
    '''
    
    ## Build a tree.
    ## NOTE: Only use non-boundary vertices throughout.
    ## 1 Pick a root.
    ## 2 Make a queue containing the root.
    ## 3 While the queue is not empty,
    ## 3a dequeue the front vertex.
    ## 3b add all neighboring vertices to the queue, so long as they aren't already in the tree.
    ## 3c add edges from the dequeued vertex to all neighboring vertices.
    
    root = None
    tree_edges = []
    tree_vertices = []
    
    ## 1
    for v in xrange( len( mesh.vs ) ):
        if not mesh.vertex_is_boundary( v ) and v not in exclude_vertices:
            root = v
            break
    ## 2
    queue = []
    if root is not None:
        queue.append( root )
    
    ## 3
    while len( queue ) > 0:
        ## 3a
        front = queue[0]
        del queue[0]
        tree_vertices.append( front )
        
        ## 3b
        neighs = [ v for v in mesh.vertex_vertex_neighbors( front ) if not mesh.vertex_is_boundary( v ) and v not in tree and v not in exclude_vertices ]
        queue.extend( neighs )
        
        ## 3c
        tree_edges.extend( [ ( front, n ) for n in neighs ] )
    
    return root, tree_edges, tree_vertices

def simple_tree( mesh ):
    
    directed_edges = []
    parented = set()
    
    for start,end in mesh.edges:
        if mesh.vs[start][2] > mesh.vs[end][2] and end not in parented:
            directed_edges.append( ( start, end ) )
            parented.add( end )
    
    return directed_edges

def watershedding_forest( mesh ):
    '''
    Given a TriMesh 'mesh' encoding a height field on z values,
    returns a tuple containing: (
        the list of directed edges in the mesh as two-tuples containing vertex indices for start, end,
        a list of connected components as sets of vertex indices,
        a list of vertex indices that merged connected components
        )
    The directed edges make up a forest of trees containing monotonically descending edges;
    that is, the watershedding begins with the entire mesh flooded and then the water level
    decreases.
    
    untested
    '''
    
    ## 1 Start a list of connected components (lists of vertices).
    ##   (I can't use "sets of vertices" because I need to preserve the order in which vertices were added.)
    ## 2 Sort the mesh vertices by their z values in descending order.
    ## 3 Iterate over the sorted vertices:
    ## 3a If the vertex is not the neighbor of any already-seen vertex,
    ##    start a new connected component with it.
    ## 3b If a vertex is the neighbor of an already-seen vertex in a single connected component,
    ##    add the directed edge from the lowest (most-recent) already-seen vertex to the new one,
    ##    and add the new vertex to the connected component.
    ## 3c If the vertex is the neighbor of already-seen vertices in two connected components,
    ##    it's a saddle; add both directed edges and merge the connected components.
    ##    UPDATE: Don't merge the connected components, the vertex lies on a seperatrix.
    ##    UPDATE: Or bail at this point, because I don't think the optimization code can handle shared leaf nodes yet!
    ##    TODO Q: Should I add the two directed edges and then add the new vertex to a special list of
    ##            vertices not to add new edges from; i.e., don't add it to the connected components
    ##            and don't merge the two connected components.
    
    ## 1
    directed_edges = []
    connected_components = []
    seperatrices = []
    
    ## 2
    decreasing = [ (xyz[2], vi) for vi,xyz in enumerate( mesh.vs ) ]
    decreasing.sort()
    decreasing.reverse()
    decreasing = [ vi for (z,vi) in decreasing ]
    
    ## 3
    for vi in decreasing:
        neighs = frozenset( mesh.vertex_vertex_neighbors( vi ) )
        neigh_components = [ neighs.intersection( frozenset( component ) ) for component in connected_components ]
        which_comps = [ i for i, comp in enumerate( neigh_components ) if len( comp ) > 0 ]
        num_neigh_components = len( which_comps )
        
        ## 3a
        if 0 == num_neigh_components:
            connected_components.append( [ vi ] )
        
        ## 3b
        elif 1 == num_neigh_components:
            component = connected_components[ which_comps[0] ]
            for component_vi in reversed( component ):
                if component_vi in neighs:
                    directed_edges.append( ( component_vi, vi ) )
                    component.append( vi )
                    break
        
        ## 3c
        elif 1 < num_neigh_components:
            ## For now, bail
            break
    
    ## There are too many values in the connected components.
    #print 'connected_components:'
    #from pprint import pprint
    #pprint( connected_components )
    print '# vertices in connected_components:', [ len( comp ) for comp in connected_components ]
    
    return directed_edges, connected_components, seperatrices

def test_minimal():
    # result = smooth_mesh2d_monotonically_decreasing( zeros(10), ones(9), 0, 0 )
    # print result
    pass

def run_parameters_test( path, mesh_opt = None ):
    print "================ test parameters: " + path + " ================"
    
    if mesh_opt is None: mesh_opt = mesh_opt_c
    
    from trimesh import TriMesh
    mesh = TriMesh()
    
    f = file( path )
    
    num_vertices = int( f.readline() )
    for i in xrange( num_vertices ):
        sline = f.readline().strip().split()
        ## Add z values.
        sline.append( 0. )
        mesh.vs.append( tuple([ float(el) for el in sline ]) )
    
    mesh.vs = asarray( mesh.vs )
    
    num_triangles = int( f.readline() )
    for i in xrange( num_triangles ):
        sline = f.readline().strip().split()
        mesh.faces.append( tuple([ int(el) for el in sline ]) )
    
    for i in xrange( num_vertices ):
        mesh.vs[i,2] = float( f.readline() )
    
    directed_edges = []
    num_edges = int( f.readline() )
    for i in xrange( num_edges ):
        sline = f.readline().strip().split()
        directed_edges.append( tuple([ int(el) for el in sline ]) )
    
    stay_close_to_original_values_weight = float( f.readline() )
    smooth_deformation_not_smooth_function = bool(int( f.readline() ))
    
    ## Save the mesh.
    mesh.write_OBJ( path + '.obj' )
    
    #print 'mesh.faces:'
    #print asarray( mesh.faces )
    print '# mesh.faces:', len( mesh.faces )
    
    #print 'mesh.vs.xy:'
    #print mesh.vs[:,:2].round(2)
    print '# mesh.vs:', len( mesh.vs )
    
    print 'mesh.vs.z:'
    print mesh.vs[:,2].round(3)
    
    ## There are too many directed edges.
    # print 'directed edges:'
    # print directed_edges
    print '# directed edges:', len( directed_edges )
    
    undirected_edges = set( [ frozenset( edge ) for edge in directed_edges ] )
    print '# undirected edges:', len( undirected_edges )
    if len( directed_edges ) != len( undirected_edges ):
        print 'WARNING: # directed and undirected edges differ!  This means that both directions appear in the list!  You will die!'
    
    ## To inspect the data:
    #from pydb import debugger
    #debugger()
    
    new_values = mesh_opt( mesh.vs[:,:2], mesh.faces, mesh.vs[:,2], directed_edges,
        stay_close_to_original_values_weight,
        smooth_deformation_not_smooth_function
        )
    
    print 'new mesh.vs.z:'
    print asarray( new_values ).round(3)
    
    ## Update the mesh's z values.
    mesh.vs[:,2] = new_values
    
    return mesh

def run_test( mesh, test_name = '_untitled', mesh_opt = None ):
    print "================ test" + test_name + " ================"
    
    if mesh_opt is None: mesh_opt = mesh_opt_c
    
    mesh.vs = asarray( mesh.vs )
    
    #print 'mesh.faces:'
    #print asarray( mesh.faces )
    print '# mesh.faces:', len( mesh.faces )
    
    #print 'mesh.vs.xy:'
    #print mesh.vs[:,:2].round(2)
    print '# mesh.vs:', len( mesh.vs )
    
    print 'mesh.vs.z:'
    print mesh.vs[:,2].round(3)
    
    
    ## We need to build a tree.
    mname = 'regular-101-sin10x_sin10y.obj'
    pname = 'r101-edges.pickle'
    directed_edges = None
    if os.path.split( argv[1] )[1] == mname:
        try:
            import cPickle as pickle
            f = open( os.path.join( '.', pname ), 'r' )
            directed_edges = pickle.load( f )
            del f
        except IOError:
            #directed_edges = interior_tree_on_mesh( mesh )[1]
            directed_edges = watershedding_forest( mesh )[0]
            f = open( os.path.join( '.', pname ), 'w' )
            pickle.dump( directed_edges, f, pickle.HIGHEST_PROTOCOL )
            f.flush()
            f.close()
            del f
    else:
        #directed_edges = interior_tree_on_mesh( mesh )[1]
        directed_edges = watershedding_forest( mesh )[0]
    
    ## There are too many directed edges.
    # print 'directed edges:'
    # print directed_edges
    print '# directed edges:', len( directed_edges )
    
    new_values = mesh_opt( mesh.vs[:,:2], mesh.faces, mesh.vs[:,2], directed_edges,
        stay_close_to_original_values_weight = 0.
        #stay_close_to_original_values_weight = 1e4
        )
    
    print 'new mesh.vs.z:'
    print asarray( new_values ).round(3)
    
    ## Update the mesh's z values.
    mesh.vs[:,2] = new_values

def mesh_opt_c( *args, **kwargs ):
    return smooth_mesh2d_monotonically_decreasing( *args, **kwargs )

def mesh_opt_fmin_ncg( *args, **kwargs ):
    import scipy.sparse, scipy.optimize
    smoother = smoother_mesh2d_monotonically_decreasing( *args, **kwargs )
    
    ## I did a search from 1 to 10 and 5 was the fastest for regular-41-sin10x_sin10y.obj.
    hess_every_N_calls = 5
    hess_calls = [0]
    last_hess = [None]
    last_hess_call = [-hess_every_N_calls]
    def hess( xvec ):
        hess_calls[0] += 1
        if last_hess_call[0] + hess_every_N_calls > hess_calls[0]:
            return last_hess[0]
        
        i,j,v = smoother.h( xvec )
        h = scipy.sparse.coo_matrix( ( v, (i,j) ), shape = ( smoother.num_dofs(), smoother.num_dofs() ) )
        #print 'h:', h
        print 'hess evaluated with nnz:', len(v)
        result = h.todense()
        
        last_hess[0] = result
        last_hess_call[0] = hess_calls[0]
        
        return result
    
    def hess_simple( xvec ):
        i,j,v = smoother.h( xvec )
        h = scipy.sparse.coo_matrix( ( v, (i,j) ), shape = ( smoother.num_dofs(), smoother.num_dofs() ) )
        #print 'h:', h
        print 'hess evaluated with nnz:', len(v)
        result = h.todense()
        return result
    
    def grad( xvec ):
        g = smoother.g( xvec ).ravel()
        #print 'g:', g
        print '|grad|^2:', dot( g,g )
        return g
    
    def energy( xvec ):
        e = float( smoother.f( xvec ) )
        print 'e:', e
        return e
    
    def pr( foo ):
        print foo
        return foo
    
    ## TODO: Use scikits.openopt interface.
    opt_result = scipy.optimize.fmin_ncg(
        ## f
        lambda xvec: energy( xvec ),
        ## x0
        smoother.initial_dofs().ravel(),
        fprime = lambda xvec: grad( xvec ),
        fhess = lambda xvec: hess( xvec ),
        ## Don't output 'Optimization terminated successfully' messages and statistics.
        #disp = 0,
        ## The default relative tolerance is 1e-5; it that too large?
        #avextol = 1e-7,
        #callback = callback,
        ## We need full output to determine if it terminated successfully.
        full_output = True
        )
    #pprint( "opt_result" )
    #pprint( opt_result )
    ## Full output item 5 is the termination error code; if it's zero, optimization converged.
    converged = not opt_result[5]
    
    if not converged:
        pprint( "diverged!" )
        raise RuntimeError, "optimize_shell diverged!"
    
    return smoother.values( opt_result[0] )

def mesh_opt_fmin_bfgs( *args, **kwargs ):
    import scipy.sparse, scipy.optimize
    smoother = smoother_mesh2d_monotonically_decreasing( *args, **kwargs )
    
    def grad( xvec ):
        g = smoother.g( xvec ).ravel()
        print '|grad|^2:', dot( g,g )
        return g
    
    def energy( xvec ):
        e = float( smoother.f( xvec ) )
        print 'e:', e
        return e
    
    def pr( foo ):
        print foo
        return foo
    
    opt_result = scipy.optimize.fmin_bfgs(
        ## f
        lambda xvec: energy( xvec ),
        ## x0
        smoother.initial_dofs().ravel(),
        fprime = lambda xvec: grad( xvec ),
        ## Don't output 'Optimization terminated successfully' messages and statistics.
        #disp = 0,
        ## The default relative tolerance is 1e-5; it that too large?
        #avextol = 1e-7,
        #callback = callback,
        ## We need full output to determine if it terminated successfully.
        #full_output = True
        )
    #pprint( "opt_result" )
    #pprint( opt_result )
    ## Full output item 5 is the termination error code; if it's zero, optimization converged.
    #converged = not opt_result[5]
    ## It doesn't seem to provide a 'converged' output value.
    converged = True
    
    return smoother.values( opt_result[0] )

def test1():
    # run_test( zeros(10), ones(9), 0, 0, '1' )
    pass

argv = None
def main():
    global argv
    import sys
    argv = list( sys.argv )
    
    mesh_opt = None
    kSciPyArg = 'scipy'
    if kSciPyArg in argv:
        print 'Using scipy.optimize.fmin_ncg'
        mesh_opt = mesh_opt_fmin_ncg
        #mesh_opt = mesh_opt_fmin_bfgs
        del argv[ argv.index( kSciPyArg ) ]
    
    if len( argv ) >= 2:
        
        if os.path.splitext( argv[1] )[1].lower() == '.obj':
            from trimesh import TriMesh
            mesh = TriMesh.FromOBJ_FileName( argv[1] )
            
            run_test( mesh, test_name = ' ' + os.path.split( argv[1] )[1], mesh_opt = mesh_opt )
            #run_test( mesh, ' ' + os.path.split( argv[1] )[1] )
            #run_test_scipy( mesh, ' ' + os.path.split( argv[1] )[1] )
        else:
            mesh = run_parameters_test( argv[1], mesh_opt = mesh_opt )
        
        output_dir = os.path.join( '.', 'output' )
        if not os.path.exists( output_dir ): os.mkdir( output_dir )
        ## In case the mkdir() fails or there is a file named 'output_dir', bail.
        if not os.path.isdir( output_dir ): sys.exit(-1)
        
        ## Find a unique name.  Don't clobber existing output files.
        fname_no_ext = os.path.join( output_dir, os.path.splitext( os.path.split( argv[1] )[1] )[0] + '-opt_height_field' )
        suffix_count = 0
        suffix = ''
        ext = '.obj'
        while os.path.exists( fname_no_ext + suffix + ext ):
            suffix_count += 1
            suffix = ' ' + str( suffix_count )
        
        mesh.write_OBJ( fname_no_ext + suffix + ext )
        sys.exit(0)
    
    # test_minimal()
    # test1()

if __name__ == "__main__": main()
