def voronoi(network,
            geometry,
            offset,
              **kwargs):
    r"""
    Offset the throat vertices effectively erroding the facet by the offset distance supplied
    """    
    Nt = geometry.num_throats()
    area = sp.ndarray(Nt)
    perimeter = sp.ndarray(Nt)
    offset_verts = sp.ndarray(Nt,dtype=object)
    offset_error = sp.ndarray(Nt)
    throat_COM = sp.ndarray([Nt,3])
    for i in range(Nt):
        offset_rand = (sp.random.rand(1)-0.5)*offset            
        throat_verts=geometry["throat.vertices"][i]
        throat_normal=geometry["throat.normal"][i]
        area[i],perimeter[i],offset_verts[i],throat_COM[i],offset_error[i] = vo.get_throat_geom(throat_verts,throat_normal,offset)
    
    for i in range(Nt):
        if offset_error[i] > 0 and len(offset_verts[i]) > 0:
            offset_verts[i]=[]
    "Properties that depend on the offset vertices are the area, perimeter and the centroid or COM"
    "To speed things up we could save them all now rather than processing them individually"
    if kwargs["set_dependent"]==True:
        geometry["throat.area"]=area
        geometry["throat.perimeter"]=perimeter
        geometry["throat.centroid"]=throat_COM
    
    return offset_verts

def innerRoots(matrix):
    ris = diagRoots(matrix)
    global bb,eb,dims,wb,i,j 
    bb,eb,dims,wb,i,j = findBlocks0(matrix)
    for z in range(1,j):
        for t in range(0,j-z):
            rij = matrix[bb[t]:eb[t]+1,bb[t+z]:eb[t+z]+1]
            tempsum = summ(ris,dims[t],dims[t+z],t,z)
            if (dims[t]==1 and dims[z+t]==1):
                tnoto = rij - tempsum
                ris[bb[t]:eb[t]+1,bb[t+z]:eb[t+z]+1] = tnoto/(selectBlock(t,t,ris)+selectBlock(z+t,z+t,ris))

            elif (dims[t]==2 and dims[t+z]==1):
                tnoto = rij - tempsum
                ujj = selectBlock(z+t,z+t,ris)
                uii = selectBlock(t,t,ris)
                sysris = sp.ndarray((2,1))
                coeff = sp.array(uii)
                coeff[0,0] += ujj[0,0]
                coeff[1,1] += ujj[0,0]
                sysris = np.linalg.solve(coeff,tnoto)
                ris[bb[t]:eb[t]+1,bb[t+z]:eb[t+z]+1] = sysris

            elif (dims[t]==1 and dims[t+z]==2):
                tnoto = rij - tempsum
                ujj = selectBlock(z+t,z+t,ris)
                uii = selectBlock(t,t,ris)
                sysris = sp.ndarray((2,1))
                coeff = sp.zeros((2,2))
                coeff[0,0] = uii[0,0] + ujj[0,0]
                coeff[1,1] = uii[0,0] + ujj[1,1]
                coeff[0,1] = ujj[1,0]
                coeff[1,0] = ujj[0,1]
                sysris = np.linalg.solve(coeff,tnoto.transpose())
                print tempsum
                ris[bb[t]:eb[t]+1,bb[t+z]:eb[t+z]+1] = sysris.transpose()

            else:
                tnoto = rij - tempsum
                ujj = selectBlock(z+t,z+t,ris)
                uii = selectBlock(t,t,ris)
                coeff = sp.ndarray((4,4))
                coeff[0,0] = uii[0,0]+ujj[0,0]; coeff[0,1] = uii[0,1]; coeff[0,2] = ujj[1,0]; coeff[0,3] = 0
                coeff[1,0] = uii[1,0]; coeff[1,1] = uii[1,1]+ujj[0,0]; coeff[1,2] = 0; coeff[1,3] = ujj[1,0]
                coeff[2,0] = ujj[0,1]; coeff[2,1] = 0; coeff[2,2] = uii[0,0]+ujj[1,1]; coeff[2,3] = uii[0,1]
                coeff[3,0] = 0; coeff[3,1] = ujj[0,1]; coeff[3,2]= uii[1,0]; coeff[3,3]=uii[1,1] + ujj[1,1]
                b = rij - tempsum; cc = sp.ndarray((4,1)); cc[0] = b[0,0]; cc[1] = b[1,0]; cc[2]=b[0,1]; cc[3]=b[1,1]

                temp = sp.zeros((2,2))
                FF = np.linalg.solve(coeff,cc)
                for k in range(0,2):
                    for h in range(0,2):
                        temp[h,k] = FF[2*k+h,0];
    #                      temp[0,0] = FF[0,0]; temp[0,1] = FF[2,0]; temp[1,0] = FF[1,0];temp[1,1] = FF[3,0]

                ris[bb[t]:eb[t]+1,bb[t+z]:eb[t+z]+1] = temp
                #ris[bb[t]:eb[t],bb[t+z]:eb[t+z]+1]
    return ris

def getGenerationProgress(resumeFile):
  
  def _getProgress(splitLine, _numParams):
    err = float(splitLine[-_numParams - 1])
    try:
      uniformRange = float(splitLine[-_numParams - 2])
      logRange = float(splitLine[-_numParams - 3])
      fR = float(splitLine[-_numParams - 4])
      pR = max(uniformRange, scipy.log(logRange))
    except ValueError:
      fR = float('NaN')
      pR = float('NaN')
    return (err, fR, pR)
  
  def _getNextLine(_fIn):
    # get the next non-empty, non-comment line from file
    line = next(_fIn)[:-1].split('#')[0]
    while not line:
      line = next(_fIn)[:-1].split('#')[0]
    return line
  
  with open(resumeFile, 'r') as fIn:
    # get the number of parameters from the first line
    line = _getNextLine(fIn)
    numParams = int(line.split()[0])
    
    # skim through the parameter descriptions and other uninteresting stuff
    while 'generation history' not in line:
      line = next(fIn)
    
    # get the number of generations in the history
    numGenerations = int(line.split(None)[0])
    
    # get the generation history
    errors = scipy.ndarray(numGenerations)
    fRange = scipy.ndarray(numGenerations)
    pRange = scipy.ndarray(numGenerations)
    for genNum in range(numGenerations):
      line = _getNextLine(fIn)
      (errors[genNum], fRange[genNum], pRange[genNum]) = \
        _getProgress(line.split(None), numParams)
    
    # get the distribution of current errors
    line = _getNextLine(fIn)
    populationSize = int(line.split()[0])
    currentErrors = scipy.ndarray(populationSize)
    for n in range(populationSize):
      line = _getNextLine(fIn)
      currentErrors[n] = float(line.split()[0])
  
  return (errors, fRange, pRange, currentErrors)

def voronoi(geometry, vertices='throat.centroid', **kwargs):
    r"""
    Calculate the centroid from the mean of the throat centroids

    Parameters
    ----------
    geometry : OpenPNM Geometry object
        The Geometry object which this model is associated.  This is needed to
        access the values of the ``vertices``.

    vertices : string
        The dictionary contain of the array containing the throat vertice
        coordiantes.  The default is 'throat.centroid'

    """
    network = geometry._net
    value = _sp.ndarray([geometry.num_pores(), 3])
    value.fill(0.0)
    pore_map = geometry.map_pores(target=network,
                                  pores=geometry.pores(),
                                  return_mapping=True)
    for i, net_pore in enumerate(pore_map['target']):
        geom_pore = pore_map['source'][i]
        net_throats = network.find_neighbor_throats(net_pore)
        geom_throats = network.map_throats(target=geometry,
                                           throats=net_throats,
                                           return_mapping=True)['target']
        verts = geometry[vertices][geom_throats]
        " Ignore all zero centroids "
        verts = verts[~_sp.all(verts == 0, axis=1)]
        if len(verts) > 0:
            value[geom_pore] = _sp.mean(verts, axis=0)

    return value

 def _generate_masked_mesh(self, cell_mask=None):
     r"""
     Generates the mesh based on the cell mask provided
     """
     #
     if cell_mask is None:
         cell_mask = sp.ones(self.data_map.shape, dtype=bool)
     #
     # initializing arrays
     self._edges = sp.ones(0, dtype=str)
     self._merge_patch_pairs = sp.ones(0, dtype=str)
     self._create_blocks(cell_mask)
     #
     # building face arrays
     mapper = sp.ravel(sp.array(cell_mask, dtype=int))
     mapper[mapper == 1] = sp.arange(sp.count_nonzero(mapper))
     mapper = sp.reshape(mapper, (self.nz, self.nx))
     mapper[~cell_mask] = -sp.iinfo(int).max
     #
     boundary_dict = {
         'bottom':
             {'bottom': mapper[0, :][cell_mask[0, :]]},
         'top':
             {'top': mapper[-1, :][cell_mask[-1, :]]},
         'left':
             {'left': mapper[:, 0][cell_mask[:, 0]]},
         'right':
             {'right': mapper[:, -1][cell_mask[:, -1]]},
         'front':
             {'front': mapper[cell_mask]},
         'back':
             {'back': mapper[cell_mask]},
         'internal':
             {'bottom': [], 'top': [], 'left': [], 'right': []}
     }
     #
     # determining cells linked to a masked cell
     cell_mask = sp.where(~sp.ravel(cell_mask))[0]
     inds = sp.in1d(self._field._cell_interfaces, cell_mask)
     inds = sp.reshape(inds, (len(self._field._cell_interfaces), 2))
     inds = inds[:, 0].astype(int) + inds[:, 1].astype(int)
     inds = (inds == 1)
     links = self._field._cell_interfaces[inds]
     #
     # adjusting order so masked cells are all on links[:, 1]
     swap = sp.in1d(links[:, 0], cell_mask)
     links[swap] = links[swap, ::-1]
     #
     # setting side based on index difference
     sides = sp.ndarray(len(links), dtype='<U6')
     sides[sp.where(links[:, 1] == links[:, 0]-self.nx)[0]] = 'bottom'
     sides[sp.where(links[:, 1] == links[:, 0]+self.nx)[0]] = 'top'
     sides[sp.where(links[:, 1] == links[:, 0]-1)[0]] = 'left'
     sides[sp.where(links[:, 1] == links[:, 0]+1)[0]] = 'right'
     #
     # adding each block to the internal face dictionary
     inds = sp.ravel(mapper)[links[:, 0]]
     for side, block_id in zip(sides, inds):
         boundary_dict['internal'][side].append(block_id)
     self.set_boundary_patches(boundary_dict, reset=True)

    def find_throat_facets(self, throats=None):
        r"""
        Finds the coordinates of the Voronoi pores that define the facet or
        ridge between the pore-pairs associated with the given throat.

        Parameters
        ----------
        throats : array_like
            The throats whose facets are sought.  The given throats should be
            from the 'delaunay' network. If no throats are specified, all
            'delaunay' throats are assumed.

        Notes
        -----
        The method is not well optimized as it scans through each given throat
        inside a for-loop, so it could be slow for large networks.

        """
        if throats is None:
            throats = self.throats('delaunay')
        else:
            throats = self.filter_by_label(throats, labels='delaunay')
        if 'throat.facet_coords' not in self.keys():
            self['throat.facet_coords'] = sp.ndarray((self.Nt, ), dtype=object)
        tvals = self['throat.interconnect'].astype(int)
        am = self.create_adjacency_matrix(data=tvals, sprsfmt='lil')
        for t in throats:
            P12 = self['throat.conns'][t]
            Ps = list(set(am.rows[P12][0]).intersection(am.rows[P12][1]))
            if sp.size(Ps) > 0:
                self['throat.facet_coords'][t] = self['pore.coords'][Ps]

    def find_pore_hulls(self, pores=None):
        r"""
        Finds the coordinates of the Voronoi pores that define the convex hull
        around the given pores.

        Parameters
        ----------
        pores : array_like
            The pores whose convex hull are sought.  The given pores should be
            from the 'delaunay' network.  If no pores are given, then the hull
            is found for all 'delaunay' pores.

        Notes
        -----
        This metod is not fully optimized as it scans through each pore in a
        for-loop, so could be slow for large networks.
        """
        if pores is None:
            pores = self.pores('delaunay')
        else:
            pores = self.filter_by_label(pores, labels='delaunay')
        if 'pore.hull_coords' not in self.keys():
            self['pore.hull_coords'] = sp.ndarray((self.Np, ), dtype=object)
        tvals = self['throat.interconnect'].astype(int)
        am = self.create_adjacency_matrix(data=tvals, sprsfmt='lil')
        for p in pores:
            Ps = am.rows[p]
            if sp.size(Ps) > 0:
                self['pore.hull_coords'][p] = self['pore.coords'][Ps]

def c2c(network, geometry, pore_centroid='pore.centroid',
        throat_centroid='throat.centroid', **kwargs):
    r"""
    Calculate the centre to centre distance from centroid of pore1 to centroid
    of throat to centroid of pore2.

    Parameters
    ----------
    pore_centroid and throat_centroid : string
        The dictionary keys to the arrays containing the pore and throat
        centroid coordinates.

    Notes
    -----
    This is tricky as connections are defined at Network level but centroids
    are stored on geometry. The pore and throat map relates the geometry index
    to the network index but we must look up the index of the map to go back
    to geometry index of the connected pores. This will probably break down
    when a throat connects two different geometries.
    """
    throats = geometry.map_throats(network, geometry.throats())
    connections = network['throat.conns'][throats]
    net_pore1 = connections[:, 0]
    net_pore2 = connections[:, 1]
    pore_centroids = network[pore_centroid]
    throat_centroids = network[throat_centroid][throats]
    v1 = throat_centroids-pore_centroids[net_pore1]
    v2 = throat_centroids-pore_centroids[net_pore2]
    value = _sp.ndarray(len(connections))
    for i in range(len(connections)):
        value[i] = _sp.linalg.norm(v1[i])+_sp.linalg.norm(v2[i])
    return value

def setTransitions(dimension):
    """ setting transitions in the state space """
    space_size = np.power(2, dimension)
    transition = np.ndarray(shape=(space_size, space_size), dtype=bool)
    transition.fill(False)  # False)

    state1 = [0 for i in range(dimension)]
    state2 = state1[:]
    for i in range(dimension):
        state2[i] = 1
        transition[st2Ind(state1)][st2Ind(state2)] = True  # forward transition
        transition[st2Ind(state2)][st2Ind(state1)] = True  # backward transition
        state1 = state2[:]

    state1 = [0 for i in range(dimension)]
    state2 = state1[:]
    for i in range(dimension):
        state2[dimension-i-1] = 1
        transition[st2Ind(state1)][st2Ind(state2)] = True  # forward transition
        transition[st2Ind(state2)][st2Ind(state1)] = True  # backward transition
        state1 = state2[:]

    printTransitions(transition, dimension)

    return transition

def voronoi(network, geometry, **kwargs):
    r"""
    Update the pore vertices from the voronoi vertices
    """
    pores = geometry.map_pores(network, geometry.pores())
    value = _sp.ndarray(len(pores), dtype=object)
    for i in range(len(pores)):
        value[i] = _sp.asarray(list(network['pore.vert_index'][pores[i]].values()))
    return value

def voronoi(network, geometry, **kwargs):
    r"""
    Update the pore vertices from the voronoi vertices
    """
    throats = geometry.map_throats(network, geometry. throats())
    value = _sp.ndarray(len(throats), dtype=object)
    for i in range(len(throats)):
        value[i] = \
            _sp.asarray(list(network['throat.vert_index'][throats[i]].values()))
    return value

def voronoi(geometry,
            pore_vertices='pore.vertices',
            **kwargs):
    r"""
    Calculate the centroid of the pore from the voronoi vertices - C.O.M
    """
    verts = geometry[pore_vertices]
    value = _sp.ndarray([len(verts),3])
    for i,vert in enumerate(verts):
        value[i] = _sp.array([vert[:,0].mean(),vert[:,1].mean(),vert[:,2].mean()])
    return value

def twobytworoot(matrix):
    if (matrix.shape[0]==2):
        a = sp.sqrt(sp.linalg.eigvals(matrix)[0]).real
        ris=sp.ndarray(shape=(2,2))
        ris[0,0] = a + (1/(4*a))*(matrix[0,0] - matrix[1,1])
        ris[1,1] = a - (1/(4*a))*(matrix[0,0] - matrix[1,1])
        ris[0,1] = (1/(2*a))*matrix[0,1]
        ris[1,0] = (1/(2*a))*matrix[1,0]
    else:
        return sp.sqrt(matrix)
    return ris

def findBlocks0(matrix):
    dim = matrix.shape[0]
    i=j=0; eb = sp.ndarray(dim,sp.integer); bb = sp.ndarray(dim,sp.integer)
    dims = sp.ndarray(dim,sp.int8); wb = sp.ndarray(dim)
    while (i<dim):
        if (i<dim-1 and matrix[i+1,i]!=0):
            bb[j] = i
            eb[j] = i+1
            dims[j] = 2; wb[i] = wb[i+1]=j
            i+=1

        else:
            bb[j] = eb[j] = i
            dims[j] = 1
            wb[i] = j

        i+=1
        j+=1

    return bb,eb,dims,wb,i,j #i,wb dims maybe are useless

def centre_of_mass(geometry, vertices='throat.offset_vertices', **kwargs):
    r"""
    Calculate the centre of mass of the throat from the voronoi vertices.
    """
    Nt = geometry.num_throats()
    outer_verts = geometry['throat.vertices']
    offset_verts = geometry[vertices]
    normal = geometry['throat.normal']
    z_axis = [0, 0, 1]
    value = _sp.ndarray([Nt, 3])
    for i in range(Nt):
        if len(offset_verts[i]) > 2:
            verts = offset_verts[i]
        elif len(outer_verts[i]) > 2:
            verts = outer_verts[i]
        else:
            verts = []
        if len(verts) > 0:
            # For boundaries some facets will already be aligned with the axis -
            # if this is the case a rotation is unnecessary and could also cause
            # problems
            angle = tr.angle_between_vectors(normal[i], z_axis)
            if angle == 0.0 or angle == _sp.pi:
                "We are already aligned"
                rotate_input = False
                facet = verts
            else:
                rotate_input = True
                M = tr.rotation_matrix(tr.angle_between_vectors(normal[i], z_axis),
                                       tr.vector_product(normal[i], z_axis))
                facet = _sp.dot(verts, M[:3, :3].T)
            # Now we have a rotated facet aligned with the z axis - make 2D
            facet_2D = _sp.column_stack((facet[:, 0], facet[:, 1]))
            z = _sp.unique(_sp.around(facet[:, 2], 10))
            if len(z) == 1:
                # We need the vertices arranged in order so perform a convex hull
                hull = ConvexHull(facet_2D)
                ordered_facet_2D = facet_2D[hull.vertices]
                # Call the routine to calculate an area wighted centroid from the
                # 2D polygon
                COM_2D = vo.PolyWeightedCentroid2D(ordered_facet_2D)
                COM_3D = _sp.hstack((COM_2D, z))
                # If we performed a rotation we need to rotate back
                if (rotate_input):
                    MI = tr.inverse_matrix(M)
                    # Unrotate the offset coordinates using the inverse of the
                    # original rotation matrix
                    value[i] = _sp.dot(COM_3D, MI[:3, :3].T)
                else:
                    value[i] = COM_3D
            else:
                print('Rotation Failed: ' + str(_sp.unique(facet[:, 2])))

    return value

def findGrad(imagePixels):
	dx = ndimage.sobel(imagePixels, 0) # horizontal derivative
	dy = ndimage.sobel(imagePixels, 1) # vertical derivative
	# mag = numpy.hypot(dx, dy)  # magnitude
	# mag *= 255.0 / numpy.max(mag)  # normalize (Q&D)
	# scipy.misc.imshow(mag)
	grad = scipy.ndarray((dx.shape[0],dx.shape[1], 2), dtype=float)
	for y in range(dx.shape[0]):
		for x in range(dx.shape[1]):
			grad[y][x][0] = dx[y][x]
			grad[y][x][1] = dy[y][x]
	return grad

def blockRoot(matrix):
    if (matrix.shape[0]==2):
        a = csr3(gev(matrix)).real
        ris=sp.ndarray(shape=(2,2))
        ris[0,0] = a + (1/(4*a))*(matrix[0,0] - matrix[1,1])
        ris[1,1] = a - (1/(4*a))*(matrix[0,0] - matrix[1,1])
        ris[0,1] = (1/(2*a))*matrix[0,1]
        ris[1,0] = (1/(2*a))*matrix[1,0]
    else:
        #print(sp.sqrt(matrix))
        return sp.sqrt(matrix)
    return ris

def voronoi2(geometry,
             vertices='throat.centroid',
             **kwargs):
    r"""
    Calculate the centroid from the mean of the throat centroids
    """
    value = _sp.ndarray([geometry.num_pores(),3])
    pore_map = geometry.map_pores(geometry.pores(),geometry._net)
    for geom_pore,net_pore in pore_map:
        net_throats = geometry._net.find_neighbor_throats(net_pore)
        geom_throats = geometry._net.map_throats(net_throats,geometry)[:,1]
        verts = geometry[vertices][geom_throats]
        value[geom_pore]=_sp.array([verts[:,0].mean(),verts[:,1].mean(),verts[:,2].mean()])
    return value

    def __init__(self,KF_gen,KF_init,KF_prior,optF,upper,lower,dim):
	self.KF_gen=KF_gen
	self.KF_hyp=KF_init
	self.KF_prior=KF_prior
        self.KFs=[KF_gen(self.KF_hyp)]
        self.optF=optF
        self.upper=upper
        self.lower=lower
        self.nsam=0
        self.X=sp.matrix(sp.ndarray([0,dim]))
        self.Y=sp.matrix(sp.ndarray([0,1]))
        self.best=[None,sp.Inf]
        self.dim=dim
        self.Ky=[]
        self.KyRs=[]
        self.llks=[]
        self.renorm=1.
	self.cc=0
	self.finished=False
	self.fudgelimit=32
	self.hyp_trace=sp.matrix(KF_init)
	
        return

    def Optimize(self):
        """
        Performs the optimization.
        """
        self.result=optimize.fmin(self.wrapper, 
                                      x0=ndarray((self.num_inputs,),buffer=array(self.starting_points),offset=0,dtype=float), 
#                                      x0=ndarray( (self.num_inputs,1) ,buffer=array([0.784318808, 4.540607953, -11.919391073,-100]),dtype=float),
                                      args=[[]],  
                                      maxiter= self.max_evaluation,
                                      xtol= self.xtol,
                                      ftol= self.ftol,
                                      full_output=True
                                      )
        
        self.final_pop=[my_candidate(self.result[0],self.result[1])]