## This program provides a prototype library for simplicial complexes.
## Includes input, output, skeleton and boundary evaluation, linear
## extrusion, boundary and coboundary operators, and linear combination
## of chains.
## Author: Alberto Paoluzzi (paoluzzi@dia.uniroma3.it)
## Copyright (C) 2009,2010 Dipartimento Informatica e Automazione,
## Università Roma Tre, Rome, Italy.

## This library is free software; you can redistribute it and/or
## modify it under the terms of the GNU Lesser General Public
## License as published by the Free Software Foundation; either
## version 2.1 of the License, or (at your option) any later version.

## This library is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## Lesser General Public License for more details.

## You should have received a copy of the GNU Lesser General Public
## License along with this library; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

from numpy import array
from scipy import sparse,mat,eye,linalg
from pyplasm import *
import operator,copy,datetime,integr,scipy,csv
from integr import T3


def myround(x):
    return eval('%.7f' % round(x,10))

def round_or_zero (x):
    xx = myround(x)
    if abs(xx) < 1E-05: return 0.0
    else: return xx

__debug = True

## --------------------------------------------------
## --Utility functions-------------------------------
## --------------------------------------------------    

def __evalprint__(string):
    print "\n" + string + " => ", eval(string)


def printer (args):
    """ To make an 'internal printing' of a (possibly nested) list in plasm format.

    Return a string, with squaare brackets substituted by angle brackets
    """
    def one (list_of_strings):
        """To catenate a flat list of string.

        Return a single string in plasm format.
        """
        return '<' + ', '.join(list_of_strings) + '>'
    
    if isinstance(args, list):
        return one(map(printer, args))
    else:
        return str(args)


def progressive_sum(List):
    """ To compute the coordinates associated to a list of (Num) differences.
    Starting from zero.

    Return a list of Num
    """
    return [sum(List[0:i]) for i in range(len(List)+1)]

##if __name__ == "__main__" and __debug:
##    print "\n## --progressive sum of a list------------------------"
##    __evalprint__(""" progressive_sum([10,20,10,20]) """)

def _init_CellComplex (verts, d_simplices):
    
        vertices = PointSet(verts)

        def identify(d_simplices):
            new_simplices = []
            for simplex in d_simplices:
                new_simplex = []
                for v in simplex:
##                    print type(v),v
                    if type(v) == list: v = v[0]
##                    print type(v),v
                    point = map(lambda x: round_or_zero(x),
                                vertices.points[v])
                    index = vertices.dict[str(point)]
                    new_simplex = new_simplex + [index]
                new_simplices += [new_simplex]
            return new_simplices
        
        d_simplices = identify(d_simplices)
        return (map(eval, vertices.ind.values()), d_simplices)

## --------------------------------------------------
## --Skeletons of a complex--------------------------
## --------------------------------------------------

def dimension(cells):
    """ To compute the (intrinsic) dimension of a cell complex.

    Return the last non zero cardinality of k-skeletons, with
    k <= dim.
    """
    k = map(len, cells)
    dim = len(cells) - 1
    while k[dim] == 0 and dim > 0:
        dim = dim - 1
    return dim

def ispair (arg):
    return len(arg) == 2

def issingleton (arg):
    return len(arg) == 1

def traversal(tree, path, store):
    
    path = path + [tree[0]]
    
    if ispair(tree):
        for item in tree[1]:
            traversal(item, path, store)
            
    elif issingleton(tree):
        store.append(path)
        path = []


def remove_duplicates (hasDupes):
    """ To remove the duplicates from a list of cells.
    Consider as duplicate the cell with opposite orientation (using the
    canonical representation)
    The mapping on a dictionary is used, with key the repr of the elements.

    Returns a new list without duplicates.
    """
    def revert(cell):
        """ Change the sign of a permutation by exchanging the first
        and the last elements
        """
        if len(cell) > 1:
            cell = [cell[-1]] + cell[1:-1] + [cell[0]]
        return cell

    noDupes = {}
    [operator.setitem(noDupes,repr(cell), cell) for cell in hasDupes \
         if not (noDupes.has_key(repr(cell))  or \
                 noDupes.has_key(repr(revert(cell))))]
    return noDupes.values()



def skeleton (h,cells):
    """ To compute the h-skeleton of a complex (given as a list of lists of 'cells').

    Return the list of h-cells, without duplicates.
    """
    if h < len(cells):
        A = []
        for cell in cells[h]:
            A.extend(facets(cell))
    return remove_duplicates(A)



def cell_complex (cells):
    """ To compute and store all the skeletons of a cell complex.
    The input data structure is only required to contain the 0- and the
    highest dimensional skeletons (0- and n-, with dim == n)

    Return the (possibly upgraded) list of lists of k-dimensional cells,
    with 0 <= k <= dim.
    """    
    dim = dimension(cells)
    k = map(len, cells)
    for h in range(dim,1,-1):
        cells[h-1] = skeleton(h,cells)
        k[h-1] = len(cells[h-1])
    return cells


def mktables (Complex):
    """ To build the cell database of a cell complex.

    Return a list of d+1 dictionaries, where d is the dimension of the complex.
    Key is the 'repr' of the cells; values are integers in range(0,k[d]) where
    k[d] is the cardinality of the d-skeleton.
    """
    # resetting dictionaries
    # (so that each call to mktables results in a new set of dictionaries)
    dictos = []
    
    def rotate(cell):
        """ Rotate a permutation into another in the same permutation class
        """
        m = min(cell)
        local = list(cell)
        while local[0] != m:
            local = local[1:] + [local[0]]
        cell = local
        return cell
    
    def revert(cell):
        """ Change the sign of a permutation by exchanging the first and the last elements
        """
        if len(cell) > 1:
            cell = [cell[-1]] + cell[1:-1] + [cell[0]]
        return cell

    d = 0
    for skeleton in Complex[:]:
        dictos.append({})
        [operator.setitem(dictos[d], repr(cell), skeleton.index(cell)) \
            for cell in skeleton \
                if not (dictos[d].has_key(repr(cell)) or \
                         dictos[d].has_key(repr(revert(cell))))]
        d += 1
    return dictos

def facets (cell):
    """ To compute the (non oriented) boundary of a d-simplex.
    Uses the standard method of algebraic topology:
    elimination of a vertex index to generate every face.
    Both the simplex and its facets are represented as an ordered list of vertices.

    Return a list of facets, i.e. of (d-1)-dimensional simplexes.
    """
    def rotate(cell):
        """ Rotate a permutation into another in the same permutation class
        """
        m = min(cell)
        local = list(cell)
        while local[0] != m:
            local = local[1:] + [local[0]]
        cell = local
        return cell

    def revert(cell):
        """ Change the sign of a permutation by exchanging the first and the last elements
        """
        if len(cell) > 1:
            cell = [cell[-1]] + cell[1:-1] + [cell[0]]
        return cell

    faces = []
    # extracts facets of odd index ------------
    for i in range(1, len(cell), 2):
        v = cell[i]
        face = list(cell)
        face.remove(v)
        faces.append(rotate(face))


    # extracts facets of even index ------------
    for i in range(0, len(cell), 2):
        v = cell[i]
        face = list(cell)
        face.remove(v)
        faces.append(revert(rotate(face)))
    return faces

if __name__ == "__main__" and __debug:
    print "\n --boundary operator on a single cell-----------"
    __evalprint__(""" facets([1,2]) """)
    __evalprint__(""" facets([1,2,3]) """)
    __evalprint__(""" facets([1,2,3,4]) """)
    __evalprint__(""" facets([1,2,3,4,5]) """)


def homology_maps (dictos):
    """ Compute the homology map of a (simplicial) cell complex.
    The map is given as d+1 lists (4 lists in this prototype implementation)
    [b_0, b_1, b_2, b_3] corresponding to the boundary operators

    b_d: K_d -> K_(d-1).  Notice that b_0 == [] by definition.
    Each map is given as an (ordered) list of pairs.
    
    
    Return a list of lists of pairs (tuples) of positive integer indices.
    The first element is the index of a d-cell; the second element is the
    index of an incident (d-1)-facet.
    """
    def revert(cell):
        """ Change the sign of a permutation by exchanging the first and the last elements
        """
        if len(cell) > 1:
            cell = [cell[-1]] + cell[1:-1] + [cell[0]]
        return cell

    global homology
    homology = [[],[]]
    skeleton = dictos[1]
    for cell in skeleton:
        homology[1].extend([(skeleton[cell], facet[0], ) \
           for facet in facets(eval(cell))])
    d = 1
    for skeleton in dictos[2:]:
        homology.append([])
        d += 1
        for cell in skeleton:
            for facet in facets(eval(cell)):
                try:
                    key = dictos[d-1][repr(facet)]
                    homology[d].extend([(skeleton[cell], key)])
                except KeyError:
                    key = dictos[d-1][repr(revert(facet))]
                    homology[d].extend([(skeleton[cell], key)])
    return homology


#########################################################
class PointSet(object):
#########################################################
    """ The data type to represent a set of points,
        supporting a cell complex. """
    
    ## -- __init__ Method -------------------------------
    def __init__(self, points):
                            
        points = map(AA(float), points)
        self.points = points
        self.m = len(points)
        self.dim = len(points[0])
        self.dict = {}
        k = -1
        for point in points:
            point = map(lambda x: round_or_zero(x), point)
            if str(point) not in self.dict:
                k += 1
                operator.setitem(self.dict,str(point), k)
                
        self.ind = dict([[v,k] for k,v in self.dict.items()])
            

    ## -- __repr__ Method -------------------------------
    def __repr__(self):
        '''Return a string representation of this CellComplex in the form of a list of cells.'''
        return 'Pointset: %s' % self.ind

    ## -- __str__ Method -------------------------------
    def __str__(self):
        '''Return a string representation of this CellComplex in the form of a list of cells.'''
        return 'Pointset: {0} x {1}: {2}'.format(self.m, self.dim, self.points)

    ## -- copy Method -------------------------------
    def copy(self):
        """ To make a copy of a cellComplex object. """
        return copy.deepcopy(self)

    ## -- embed Method -------------------------------
    def embed (self, n):
        added_coords = n * [0.0]
        self.points = map(lambda p: p + added_coords, self.points)
        self.dim += n

    def translate (self, vect):
        return PointSet((mat(self.points) + mat(vect)).tolist())
    
    def scale (self, vect):
        dim = len(vect)
        diag = eye(dim,dim)
        for i in range(dim):
           diag[i,i] = vect[i]
        return PointSet((mat(self.points) * diag).tolist())

    def rotate (self, axis1, axis2, angle):
        axis1 += -1
        axis2 += -1
        dim = self.dim
        rotation = eye(dim,dim)
        c = math.cos(angle)
        s = math.sin(angle)
        rotation[axis1,axis1] = c
        rotation[axis2,axis2] = c
        rotation[axis1,axis2] = s
        rotation[axis2,axis1] = -s
        return PointSet((mat(self.points) * rotation).tolist())

    def max (self):
        return map(max, array(self.points).T)

    def min (self):
        return map(min, array(self.points).T)

    def med (self):
        Max = array(map(max, array(self.points).T))
        Min = array(map(min, array(self.points).T))
        return (Max + Min)/2

    def normalize (self):
        Min = array(map(min, array(self.points).T))
        Max = array(map(max, array(self.points).T))
        translation = -Min
        size = Max-Min
        scaling = array([1/s if (s != 0) else 1.0 for s in size])
        self = self.translate(translation.tolist())
        self = self.scale(scaling.tolist())
        return self


if __name__ == "__main__" and __debug:
    print "\n --class PointSet------------------------"
    __evalprint__(""" PointSet([[0],[1],[0.5]]) """)
    __evalprint__(""" PointSet([[0,0],[1,0],[0,1],[1,1]]) """)
    __evalprint__(""" PointSet([[0,0],[1,0],[0,1],[1,1]]).translate([1,1]) """)
    __evalprint__(""" PointSet([[0,0],[1,0],[0,1],[1,1]]).scale([1,0.5]) """)
    __evalprint__(""" PointSet([[0,0],[1,0],[0,1],[1,1]]).rotate(0,1,math.pi/3).ind""")
    __evalprint__(""" PointSet([[0],[1],[0.5]]).max() """)
    __evalprint__(""" PointSet([[0,0,2],[7,0,3],[0,10,0],[1,1,1]]).max() """)
    __evalprint__(""" PointSet([[0,0,2],[7,0,3],[0,10,0],[1,1,1]]).min() """)
    __evalprint__(""" PointSet([[0,0,2],[7,0,3],[0,10,0],[1,1,1]]).med() """)
    __evalprint__(""" PointSet([[0,0,-2],[-7,0,3],[0,10,0],[1,1,1]]).normalize() """)



def embed (pointset, n):
    added_coords = n * [0.0]
    points = map(lambda p: p + added_coords, pointset.points)
    return PointSet(points)

if __name__ == "__main__" and __debug:
    print "\n -- embedding a PointSet in a higher-dimensional space -----"
    __evalprint__(""" embed (PointSet([[1.0, 1.0], [0.0, 0.0], [0.0, 1.0]]), 2) """)



def subcomplex (d,List):
    """ To compute the list of adjacent d-tuples obtained by shifting the
    list List by one element.
    
    Return the list of n d-lists, with n = len(List)-d+1
    """
    return [List[i:i+d] for i in range(len(List)-d+1)]


if __name__ == "__main__" and __debug:
    print "\n --subcomplex function------------------------"
    __evalprint__(""" subcomplex(3,range(10)) """)



def shift (n, listoflists):
    """ To increment by (int) 'n' all the (int) elements of a list of lists.

    Returns a list of lists of int.
    """
    return [map(lambda x: x+n, list) for list in listoflists]



if __name__ == "__main__" and __debug:
    print "\n --shifting a list of lists of int-------------"
    __evalprint__(""" subcomplex(3,range(10)) """)
    __evalprint__(""" shift(5,subcomplex(3,range(10))) """)


def extrude (cells, hlist):
    """ To complute the multiple linear extrusion of a d-complex.
    Map R^d -> R^(d+1), according to: Ferrucci & Paoluzzi, CAD 1991.
    'cells' is a simplicial complex, given as a list of lists;
    'hlist' is a list of heights in the added dimension.

    Return the (d+1)-dimensional simplicial complex.
    Only the 0- and (d+1)- skeletons are computed.
    In order to get the remaining skeletons, the 'cell-complex' function
    must be invoked.
    """
    dim = dimension(cells)
    if dim == 0:
        cells = [[],[]]
        lastcoords = progressive_sum(hlist)
        cells[0] = AA(LIST)(lastcoords)
        cells[1] = [[i,i+1] for i in range(0,len(hlist))]
    else:
        verts = cells[0]
        simplexes = cells[dim]
        nverts = len(verts)
        lastcoords = progressive_sum(hlist)
        nsteps = len(lastcoords)
        sliced_vertices = nsteps*[verts]
        coords_distribute = lambda x: \
            COMP([ CAT, AA(COMP([ AA(AR), DISTR ])) ])(x)
        cells[0] = coords_distribute(TRANS([nsteps*[verts],lastcoords]))
        extruded_simplexes = []
        for cell in simplexes:
            I = cell + map(lambda x: x+nverts, cell)
            LL = subcomplex(dim+2,I)
            extruded_simplexes += LL
        final_simplices = []
        for i in range(nsteps-1):
            simplex_layer = shift(nverts*i,extruded_simplexes)
            final_simplices += simplex_layer
        cells.append(final_simplices)
    return cell_complex(cells)


## --------------------------------------------------
## --Format conversion-------------------------------
## --------------------------------------------------

def complex2surface (obj):
    """ Transforms the 2-skeleton of a simplicial complex in
    a triangulated surface.

    Return a list of triples of surface vertices (nD points).
    """
    cells = obj.cells
    verts = obj.vertices.points
    surface = [[verts[i],verts[j],verts[k]] for [i,j,k] in cells[2]]
    return surface


#########################################################
class CellComplex:
#########################################################
    """ The data type to represent a simplicial complex. """


    ## -- __init__ Method -------------------------------
    def __init__(self, vertices, d_simplices):
        """ A new simplicial complex with given 0-cells and d-cells. """

        vertices, d_simplices = _init_CellComplex (vertices, d_simplices)
        
        self.vertices = PointSet(vertices)
        
        self.n = len(vertices[0])
        
        self.dim = len(d_simplices[0]) - 1
        self.cells = (self.dim + 1)*[[]]
        self.cells[0] = map(lambda x: [x], self.vertices.ind.keys())
        self.cells[self.dim] = d_simplices
        
        self.dictos = mktables(cell_complex(self.cells))
        self.inv_dict = [
            dict([[v,k] for k,v in self.dictos[i].items()])
                for i in range(self.dim + 1) ]
        self.homology = homology_maps(self.dictos)
        self.cells[-1] = map(eval, self.inv_dict[-1].values())

    ## -- __repr__ Method -------------------------------
    def __repr__(self):
        '''Return a string representation of this CellComplex in the form of a list of cells.'''
        return 'CellComplex: %s' % self.inv_dict

    ## -- __str__ Method -------------------------------
    def __str__(self):
        '''Return a string representation of this CellComplex in the form of a list of cells.'''
        return ('CellComplex: \n\t.n: {0}'+ 
                        '\n\t.dim: {1}'+ 
                        '\n\t.vertices: {2}'+ 
                        '\n\t.cells: {3}'+ 
                        '\n\t.dictos: {4}'+ 
        '').format(self.n, self.dim, self.vertices, self.cells, self.dictos)

    ## -- copy Method -------------------------------
    def copy(self):
        """ To make a copy of a cellComplex object. """
        return copy.deepcopy(self)

    ## -- view Method -------------------------------
    def view(self):
        """ To give Plasm a CellComplex to visualize.

        Return 'None'.
        """
        def plasm(d, verts, cells):
            """ To print the PLASM definition of a 'd'-skeleton.

            Returns a string, to be written on the output .psm file. The skeleton name is generated
            starting from the input file 'datafile'(.dat)
            """
            Plasm.View(Plasm.mkpol(d, verts, cells))
            
        cells = self.cells
        vertices = CAT(self.vertices.points)    
        dim = self.dim
        m = self.n
        for d in range(0,dim+1):
            plasm(m, vertices, cells[d])
            

    ## -- HPC constructor -------------------------------
    def hpolc(self):
        """ To transform a CellComplex into a Hierarchical Polyhedral Complex.

        Return an object of Hpc class.
        """            
        cells = self.cells
        vertices = CAT(self.vertices.points)    
        dim = self.dim
        m = self.n
        return Plasm.mkpol(m, vertices, cells[dim])
            

    ## -- write Method -------------------------------
    def write(self, file_name):
        """ To write a '.psm' output file, with plasm syntax.

        Return 'None'.
        """
        def plasm_def(name, d, verts, cells):
            """ To print the PLASM definition of a 'd'-skeleton.

            Returns a string, to be written on the output .psm file. The skeleton name is generated
            starting from the input file 'datafile'(.dat)
            """
            if cells != []:
                cells = [ map(int,c) for c in cells ]
                definition = "DEF " + file_name + "_K" + str(d) + " = MKPOL:" \
                         + printer([ verts, cells, "< 1 .. " + str(len(cells)) + ">" ]) \
                         + ";\n\n"
            else: definition = ""
            return definition

        cells = self.cells
        vertices = self.vertices.points

        output_file = open(file_name + ".psm", "w")
        output_file.write("% "+ str(datetime.date.today()) + \
                          " PLASM converted file " + " %\n\n" )
        output_file.write("DEF points = " + printer(vertices) + ";\n\n" )
        output_file.write("DEF " + file_name + "_K0 = " +
            "MKPOL:<points, (aa:list~intsto~len):points, (aa:list~intsto~len):points>;\n\n")

        dim = self.dim
        if dim > 0:
            for d in range(1,dim+1):
                output_file.write(plasm_def( file_name, d, "points", cells[d] ))       
            output_file.write("\nVIEW: " + file_name + "_K1;\n")
        output_file.close()


    ## -- boundary Method -------------------------------
    def boundary (self):
        """ To compute the boundary of a simplicial d-complex.

        Return the closed (d-1)-complex triangulating the boundary of the input complex.
        """
        obj = self.copy()
        cells = obj.cells
        vertices = obj.vertices.points
        d = obj.dim
        dictos = obj.dictos
        h = obj.homology
        a = array(h[d])

        V = array( len(a)*[1] )
        J = a[:,0]
        I = a[:,1]

        A = sparse.coo_matrix((V,(I,J)), shape=(max(I)+1, max(J)+1)).tocsr()

        boundary_indices = [i for i in range(A.shape[0]) if A[i].sum() == 1]
        inv_dict = dict([[v,k] for k,v in dictos[d-1].items()])
        faces = [eval(inv_dict[k]) for k in boundary_indices]
        
        return CellComplex(vertices, faces)

    ## -- extrude Method -------------------------------
    def extrude (self, hlist):
        """ To complute the multiple linear extrusion of a d-complex.
        Map R^d -> R^(d+1), according to: Ferrucci & Paoluzzi, CAD 1991.
        'cells' is a simplicial complex, given as a list of lists;
        'hlist' is a list of heights in the added dimension.

        Return the (d+1)-dimensional simplicial complex.
        Only the 0- and (d+1)- skeletons are computed.
        In order to get the remaining skeletons, the 'cell-complex' function
        must be invoked.
        """

        cells = self.cells
        dim = self.dim
        vertices = self.vertices.points
        
        if dim == 0:
            cells = [[],[]]
            lastcoords = progressive_sum(hlist)
            vertices = AA(LIST)(lastcoords)
            cells[1] = [[i,i+1] for i in range(1,len(hlist)+1)]
        else:
            verts = vertices
            simplexes = cells[dim]
            nverts = len(verts)
            lastcoords = progressive_sum(hlist)
            nsteps = len(lastcoords)
            sliced_vertices = nsteps*[verts]
            coords_distribute = lambda x: \
                COMP([ CAT, AA(COMP([ AA(AR), DISTR ])) ])(x)
            vertices = coords_distribute(TRANS([nsteps*[verts],lastcoords]))
            extruded_simplexes = []
            for cell in simplexes:
                I = cell + map(lambda x: x+nverts, cell)
                LL = subcomplex(dim+2,I)
                extruded_simplexes += LL
            final_simplices = []
            for i in range(nsteps-1):
                simplex_layer = shift(nverts*i,extruded_simplexes)
                final_simplices += simplex_layer
            cells.append(final_simplices)
            
        return CellComplex(vertices, cells[-1])


    ## -- II Method -------------------------------
    def II(self,alpha,beta,gamma):
        """ To compute surface integrals of a three-variate monomial on the complex boundary.

        Return a number.
        """
        if self.n != 3:
            print "Integration error: 'n' != 3"
        elif self.dim == 3:
            self = self.boundary()
        elif self.dim == 2:
            pass
        else:
            print "warning: obj.dim != 3"

        surface = complex2surface(self)
        w = 0.0
        for triangle in surface:
            w += T3(triangle,alpha,beta,gamma)
        return w


    ## -- III Method -------------------------------
    def III(self,alpha,beta,gamma):
        """ To compute volume integrals of a three-variate monomial on the complex boundary.

        Return a number.
        """
        if self.n != 3:
            print "Integration error: 'n' != 3"
        elif self.dim == 3:
            self = self.boundary()
        elif self.dim == 2:
            pass
        else:
            print "warning: obj.dim != 3"

        surface = complex2surface(self)
        
        def magnitude(vect):
            return math.sqrt(sum(x*x for x in vect))

        w = 0.0
        for triangle in surface:
            t = mat(triangle)
            a = (t[1] - t[0]).tolist()[0]
            b = (t[2] - t[0]).tolist()[0]
            c = [a[1]*b[2] - a[2]*b[1], a[2]*b[0] - a[0]*b[2], a[0]*b[1] - a[1]*b[0]]
            w += (c[0] / magnitude(c)) * T3(triangle,alpha+1,beta,gamma)
        return w/(alpha + 1)


    ## -- multi-resolution Method -------------------------------
    def multires (self, C2):
        """ To compute a multiresolution cellcomplex generated by embedding
        the C2 cellcomplex within self.
        C2 must be a SUBCOMPLEX of self.

        Return a new CellComplex.
        """
        
        def pointset (simplices):
            """ To compute the list of (indices of) vertices of a simplex.

            Return a list of lists of integers.
            """
            points = []
            for simplex in simplices:
                points += simplex
            return set(points)

        def swap(args):
            """ To swap the two (first) elements of a list.

            Return a list of two elements.
            """
            return [args[1], args[0]]

        C1 = self.copy()
        
        # splitta il complesso interno
        S2 = C2.split()
        C2 = CellComplex(S2.vertices.points, C2.cells[2])

        # calcola il bordo interno iniziale e l'insieme dei suoi punti
        C23 = C2.boundary()
        C23_points = pointset(C23.cells[-1])
        
        # calcola il bordo interno splittato e l'insieme dei suoi punti
        S23 = C23.split()
        C2_points = pointset(C2.cells[-1])
      
        # calcola le celle (da splittare) adiacenti al bordo
        C3_cells = [s for s in C1.cells[-1]
                    if len(set(s).difference(C2_points)) == 1 ]

        if C3_cells != []:

            C3 = CellComplex(S23.vertices.points, C3_cells)
            C3_sets = map(set, C3_cells)

            # splitta il complesso adiacente al bordo interno
            k = -1
            S3 = []
            for cell in C23.cells[-1]:
                print "cell = ", cell
                v = [c.difference(cell) for c in C3_sets
                     if len(c.difference(cell)) == 1 ]
                k += 1
                if v != []:
                    v = v[0]
                    print "v = ", v
                    S3 += [ S23.cells[-1][2*k] + list(v), S23.cells[-1][2*k+1] + list(v) ]

            C2 = C2.add(C3)
            C1 = C1.subtract(C2)

            return CellComplex (S23.vertices.points, C1.cells[-1] + S2.cells[-1] + S3)        
        else:
            return CellComplex (S23.vertices.points, S2.cells[-1])





    ## -- trivial split Method -------------------------------
    def trivial_split (self):

        """ To compute the baricentric decomposition of a d-dimensional
        obj CellComplex.

        Return a simplicial CellComplex object with k_d * (d+1)! d-cells, where
        k_d is the cardinality of the d-skeleton of obj. 
        """
        # initialize
        obj = self.copy()
        verts = obj.vertices.points
        dim = obj.dim
        n = obj.n
        cells = obj.cells
        dictos = obj.dictos
        inv_dict = obj.inv_dict

        def centroid (face):
            """ To compute the centroid od a d-face.

            Return a n-point, convex combination of d+1 n-points.
            """
            simplex = [ verts[v]  for v in face ]
            d = len(simplex)
            A = mat(simplex, dtype=scipy.float32)
            C = mat(d*[1.0/d], dtype=scipy.float32)
            point = (C * A).tolist()[0]
            return point

        def facenode (face):
            """ To compute the index of the node associated to a face.

            Return an integer number.
            """
            point = centroid(face)
            nodekey = inv_dict[0][ str(point) ]
            return nodekey

        def split1_(face):
            """ To barycentrically split a face and its boundary.

            Return a list of nested pairs.
            """
            node = facenode(face)
            return [ [node] + f for f in facets(face)]

        def sign (face):
            """ To compute the sign of a d-simplex in a d-complex.

            Return either 1 or -1 according with the permutation class of the face.
            """
            vertices = [verts[v] + [1.0] for v in face]
            return cmp(linalg.det(vertices),0)
            

        # initialize
        obj = self.copy()
        verts = obj.vertices.points
        dim = obj.dim
        n = obj.n
        cells = obj.cells
        dictos = obj.dictos
        inv_dict = obj.inv_dict
        
        # add new vertices
        base=len(verts)-1
        k = 0
        for cell in cells[-1]:
            p = centroid(cell)
            verts.append(p)
            k += 1
            operator.setitem(inv_dict[0], str(p), base+k)
            operator.setitem(dictos[0], str([base+k]), base+k )

        # build the new d-cells
        ncells = dim+1 #number of splitted cells per d-cell
        oldcells = list(cells[-1])
        cells[-1] = []
        for cell in oldcells:
            cells[-1] += split1_(cell)
                        
         # coherent orientation of the d-cells
        for i in range(len(cells[-1])):
            cell = cells[-1][i]
            if dim == n == len(cell) - 1:
                if sign(cell) < 0:
                    cells[-1][i] = [cell[-1]] + cell[1:-1] + [cell[0]]            

        return CellComplex(verts, cells[-1])


   ## -- split Method -------------------------------
    def split (self):
        """ To compute the baricentric decomposition of a d-dimensional
        obj CellComplex.

        Return a simplicial CellComplex object with k_d * (d+1)! d-cells, where
        k_d is the cardinality of the d-skeleton of obj. 
        """
        def centroid (face):
            """ To compute the centroid od a d-face.

            Return a n-point, convex combination of d+1 n-points.
            """
            simplex = [ verts[v]  for v in face ]
            d = len(simplex)
            A = mat(simplex, dtype=scipy.float32)
            C = mat(d*[1.0/d], dtype=scipy.float32)
            point = (C * A).tolist()[0]
            return point

        def facenode (face):
            """ To compute the index of the node associated to a face.

            Return an integer number.
            """
            point = centroid(face)
            nodekey = inv_dict[0][ str(point) ]
            return nodekey

        def split_(face):   
            if len(face) > 1:
                return ([ facenode(face), AA(split_)(facets(face)) ])
            return face

        def sign (face):
            """ To compute the sign of a d-simplex in a d-complex.

            Return either 1 or -1 according with the permutation class of the face.
            """
            vertices = [verts[v] + [1.0] for v in face]
            return cmp(linalg.det(vertices),0)

        # initialize
        obj = self.copy()
        verts = obj.vertices.points
        dim = obj.dim
        n = obj.n
        cells = obj.cells
        dictos = obj.dictos
        inv_dict = obj.inv_dict
        
        # add new vertices
        base=len(verts)-1
        k = 0
        for skeleton in cells[1:]:
            for face in skeleton:
                p = centroid(face)
                verts.append(p)
                k += 1
                operator.setitem(inv_dict[0], str(p), base+k)
                operator.setitem(dictos[0], str([base+k]), base+k )

        # build the new d-cells
        ncells = PROD(INTSTO(dim+1)) #number of splitted cells per d-cell
        oldcells = list(cells[-1])
        cells[-1] = []
        for cell in oldcells:
            store = []
            traversal(split_(cell), [], store)
            cells[-1] += store

        # coherent orientation of the d-cells
        for i in range(len(cells[-1])):
            cell = cells[-1][i]
            if dim == n == len(face) - 1:
                if sign(cell) < 0:
                    cells[-1][i] = [cell[-1]] + cell[1:-1] + [cell[0]]            
            

        return CellComplex(verts, cells[-1])



    ## -- offspring Method -------------------------------
    def offspring(self, germs):
        """ To compute a subcomplex of self, with d-cells incident on some 0-cell
        in the list of face lists 'germs'.
        
        Return a new cellComplex with the 0-cells of self and a subset of d-cells.
        """
        subcomplex = self.copy()
        points = subcomplex.vertices.points
        cells = subcomplex.cells
        
        germ_verts = []
        for germ in germs:
            germ_verts += germ
            
        germ_verts = set(germ_verts)
        offsimplices = [s for s in cells[-1]
                        if germ_verts.intersection(s) != set([]) ]    
        return CellComplex(points, offsimplices)


    ## -- subtract Method -------------------------------
    def subtract(self, scomplex):
        """ To compute the subtraction of d-cells of 'scomplex' from self. 'scomplex'
        must be a subcomplex of self.

        Return a new subcomplex of self.
        """
        newsub = self.copy()
        points = newsub.vertices.points
        cells = newsub.cells

        arg1 = cells[-1]
        arg2 = scomplex.cells[-1]
        new_d_cells = [d_cell for d_cell in arg1 if d_cell not in arg2]
        return CellComplex(points, new_d_cells)



    ## -- add Method -------------------------------
    def add(self, scomplex):
        """ To compute the addition of d-cells of 'scomplex' to self. 'scomplex'
        must be a subcomplex of self.

        Return a new subcomplex of self.
        """
        newsub = self.copy()
        points = newsub.vertices.points
        cells = newsub.cells

        arg1 = cells[-1]
        arg2 = scomplex.cells[-1]
        arg2 = [d_cell for d_cell in arg2 if d_cell not in arg1]
        new_d_cells = arg1 + arg2
        return CellComplex(points, new_d_cells)


    ## -- normalize Method -------------------------------
    def normalize(self):
        """ To map the CellComplex so that its containment box is mapped to the
        standard unit cube.

        Return a new subcomplex of self.
        """
        newsub = self.copy()
        newsub.vertices = newsub.vertices.normalize()
        return newsub


    ## -- translate Method -------------------------------
    def translate(self, vect):
        """ To translate the CellComplex by the vect as side effect.

        Return 'None'.
        """
        self.vertices = self.vertices.translate(vect)


    ## -- scale Method -------------------------------
    def scale(self, vect):
        """ To scale the CellComplex by the vect as side effect.

        Return 'None'.
        """
        self.vertices = self.vertices.scale(vect)


    ## -- rotate Method -------------------------------
    def rotate (self, axis1, axis2, angle):
        """ To rotate the CellComplex by the 'angle' as side effect.
        The changed coords are 'axis1', 'axis2'

        Return 'None'.
        """
        self.vertices = self.vertices.rotate(axis1, axis2, angle)



def grid (listofquotes):
    """ To generate a grid of d-simplices.

    Return the complex as a 'cells' list of lists  of simplices.
    """
    zero_simplex = [[]]
    cells = zero_simplex
    for quote in listofquotes:
        cells = extrude(cells,quote)
    return CellComplex(cells[0],cells[-1])


if __name__ == "__main__" and __debug:
    print "\n## -- Format conversion (complex -> surface) --------------"
    __evalprint__(""" complex2surface(CellComplex([[0,0],[0,1],[1,1]],[[0,1,2]])) """)
    __evalprint__(""" complex2surface(grid([2*[1.0],2*[1.0]])) """)

    
if __name__ == "__main__" and __debug:
    print "\n## --class CellComplex------------------------"
    __evalprint__(""" CellComplex([[0.0],[1.0],[0.5]],[[0,2],[2,1]]) """)
    __evalprint__(""" repr(CellComplex([[0.0],[1.0],[0.5]],[[0,2],[2,1]])) """)
    __evalprint__(""" CellComplex([[0.0],[1.0],[0.5]],[[0,2],[2,1]]).copy() """)
    __evalprint__(""" CellComplex([[0,0],[0,1],[1,1]],[[0,1,2]]).view() """)
    __evalprint__(""" CellComplex([[0,0],[0,1],[1,1]],[[0,1,2]]).write("triangle") """)
    __evalprint__(""" CellComplex([[0,0],[0,1],[1,1]],[[0,1,2]]).boundary() """)
    __evalprint__(""" CellComplex([[0,0],[0,1],[1,1]],[[0,1,2]]).extrude([1.0]) """)
    __evalprint__(""" CellComplex(
        embed(PointSet([[0,0],[0,1],[1,1]]),1).points,
        [[0,1,2]]).II(0,0,0) """)


if __name__ == "__main__" and __debug:
    print "\n## --dimension independent grid operator-------------"
    __evalprint__(""" grid([ 5*[1.0] ]) """)
    __evalprint__(""" grid([ 5*[1.0], 2*[1.0] ]) """)
    __evalprint__(""" grid([ [1.0], [1.0], [1.0] ]) """)
    __evalprint__(""" grid([[1.0],[1.0],[1.0]]).III(0,0,0) """)
    __evalprint__(""" grid([[1.0],[1.0],[1.0]]).II(0,0,0) """)
    __evalprint__(""" grid([[1.0],[1.0],[1.0]]).III(1,0,0) """)
    __evalprint__(""" grid([ [2.0], [3.0], [4.0] ]).normalize() """)


if __name__ == "__main__" :#and __debug:
    cuboid = grid([[1],[1]])
    print "\n## --linear transformation methods-------------"
    __evalprint__(""" cuboid """)
    __evalprint__(""" cuboid.translate([10,0]) """)
    __evalprint__(""" cuboid.vertices """)
    __evalprint__(""" cuboid.scale([1/10,1]) """)
    __evalprint__(""" cuboid.vertices """)
    __evalprint__(""" cuboid.rotate(1,2,math.pi/3) """)
    __evalprint__(""" cuboid.vertices """)


if __name__ == "__main__":# and __debug:
    cuboid = grid([[1],[1],[1]])
    print "\n## --chain complex maps----------------------"
    __evalprint__(""" cuboid.boundary().II(0,0,0) """)
    __evalprint__(""" cuboid.boundary().III(0,0,0) """)
    __evalprint__(""" cuboid.boundary().III(1,0,0) """)
    __evalprint__(""" cuboid.boundary().III(0,1,0) """)
    __evalprint__(""" cuboid.boundary().III(0,0,1) """)
    __evalprint__(""" cuboid.boundary().III(2,0,0) """)
    __evalprint__(""" cuboid.boundary().III(0,2,0) """)
    __evalprint__(""" cuboid.boundary().III(0,0,2) """)
    __evalprint__(""" cuboid.boundary().III(0,1,1) """)
    __evalprint__(""" cuboid.boundary().III(1,0,1) """)
    __evalprint__(""" cuboid.boundary().III(1,1,0) """)
    __evalprint__(""" cuboid.II(0,0,0) """)
    __evalprint__(""" cuboid.III(0,0,0) """)
    __evalprint__(""" cuboid.III(1,0,0) """)
    __evalprint__(""" cuboid.III(0,1,0) """)
    __evalprint__(""" cuboid.III(0,0,1) """)
    __evalprint__(""" cuboid.III(2,0,0) """)
    __evalprint__(""" cuboid.III(0,2,0) """)
    __evalprint__(""" cuboid.III(0,0,2) """)
    __evalprint__(""" cuboid.III(0,1,1) """)
    __evalprint__(""" cuboid.III(1,0,1) """)
    __evalprint__(""" cuboid.III(1,1,0) """)


if __name__ == "__main__" and __debug:
    mygrid = grid([10*[1.0],10*[1.0]])
    portion = mygrid.offspring(mygrid.cells[2][13:15])
    mydiff = mygrid.subtract(portion)
    print "\n## --subcomplex extraction & splitting-------------"
    __evalprint__(""" grid([2*[1.0],2*[1.0]]).trivial_split().view() """)
    __evalprint__(""" grid([2*[1.0],2*[1.0]]).split().view() """)
    __evalprint__(""" portion.view() """)
    __evalprint__(""" mygrid.multires(portion).view() """)
    __evalprint__(""" portion.trivial_split().view() """)
    __evalprint__(""" mydiff.view() """)
    __evalprint__(""" mydiff.add(portion).view() """)



def chain_complex (obj):
    """ To compute d adjacency matrices that implement the boundary maps.

    Return a list of sparse matrices in CSR (Compressed Sparse Row) format.
    """
    def incidence_coeff (i,j):
        face = obj.inv_dict[d-1][i]
        faces = map( repr, facets(eval(obj.inv_dict[d][j])) )
        if d == 1:
            if face == faces[0]: value = 1
            elif face == faces[1]: value = -1
        elif face in faces:
            value = 1
        else:
            value = -1
        return value

    cells = obj.cells
    homology = obj.homology
    matrices = []
 
    d = 0
    for pair_set in homology[1:]:
        d += 1
        a = array(pair_set)
        # implicit transposition of pairs
        J = a[:,0]
        I = a[:,1]
##        V = array( len(a)*[1] )
        V = array([ incidence_coeff(I[h], J[h]) for h in range(len(a)) ])

        A = sparse.coo_matrix((V,(I,J)), shape=(max(I)+1, max(J)+1))
        matrices.append(A.tocsr())
    return matrices

## --------------------------------------------------
## --file reading------------------------------------
## --------------------------------------------------


def complex_read(file_name):
    """ To read a '.dat' input file, with csv syntax.
    Each cell is read from a single line.

    Return the complex as a list of lists 'cells' of simplices.
    """
    file_name = str(file_name)
    input_file = open(file_name + ".dat", "r")
    reader = csv.reader(input_file)
    cells = [[]]
    ## ---sequential reading---------
    for row in reader:
        if row[0][0] != "#":
            d = len(cells)
            cell_dim = int(row[0][0])
            n = cell_dim - d + 1
            if row[0] == "0D" and row[-1] == "0D" :
                pass
            elif row[0][0] == "0":
                for i in range(n): cells.append([])
                cells[cell_dim].append(map(float,row[1:]))
            elif row[0][1] == "D":
                for i in range(n): cells.append([])
                cells[cell_dim].append(map(int,row[1:]))
    ## ---skeleton extraction--------
    input_file.close()
    return CellComplex(cells[0],cells[-1])




if __name__ == "__main__" and __debug:
    cuboid = grid([[1],[1],[1]])
    print "\n## --chain complex maps----------------------"
    __evalprint__(""" [mat.todense() for mat in chain_complex(cuboid) ]""")
    __evalprint__(""" [mat.todense() for mat in chain_complex(cuboid.boundary()) ]""")



def Map (mapping, pol):
    """ To Apply the "mappping" to the vertices of the CellComplex pol.
    Be careful to write "Map", with the first letter in Upper case.

    Return a new CellComplex.
    """
    newpol = copy.copy(pol)
    points = newpol.vertices.points
    newpoints = map(mapping, points)
    return CellComplex(newpoints, newpol.cells[-1])
    

if __name__ == "__main__":# and __debug:
    cuboid = grid([[1],[1],[1]])
    print "\n## --chain complex maps----------------------"
    __evalprint__(""" Map(ID, grid([[1.0],[1.0],[1.0]])) """)


def merge_complex(c1,c2):
    c = CellComplex(c1.vertices.points + c2.vertices.points, c1.cells[-1])
    ind = [c.vertices.dict[repr(v)] for v in c2.vertices.points]
    def code (k): return ind[k]
    verts = c.vertices.points
    newcells = remove_duplicates(
        c1.cells[-1] + [map(code, cell) for cell in c2.cells[-1]] )
    return CellComplex(verts,newcells)