text:
r"""
Groebner Fans

Sage provides much of the functionality of gfan, which is a
software package whose main function is to enumerate all reduced
Groebner bases of a polynomial ideal. The reduced Groebner bases
yield the maximal cones in the Groebner fan of the ideal. Several
subcomputations can be issued and additional tools are included.
Among these the highlights are:

-  Commands for computing tropical varieties.

-  Interactive walks in the Groebner fan of an ideal.

-  Commands for graphical renderings of Groebner fans and monomial
   ideals.

AUTHORS:

- Anders Nedergaard Jensen: Wrote the gfan C++ program, which
  implements algorithms many of which were invented by Jensen, Komei
  Fukuda, and Rekha Thomas. All the underlying hard work of the
  Groebner fans functionality of Sage depends on this C++ program.

- William Stein (2006-04-20): Wrote first version of the Sage code
  for working with Groebner fans.

- Tristram Bogart: the design of the Sage interface
  to gfan is joint work with Tristram Bogart, who also supplied
  numerous examples.

- Marshall Hampton (2008-03-25): Rewrote various functions to use
  gfan-0.3. This is still a work in progress, comments are
  appreciated on sage-devel@googlegroups.com (or personally at
  hamptonio@gmail.com).

EXAMPLES::

    sage: x,y = QQ['x,y'].gens()
    sage: i = ideal(x^2 - y^2 + 1)
    sage: g = i.groebner_fan()
    sage: g.reduced_groebner_bases()
    [[x^2 - y^2 + 1], [-x^2 + y^2 - 1]]

TESTS::

    sage: x,y = QQ['x,y'].gens()
    sage: i = ideal(x^2 - y^2 + 1)
    sage: g = i.groebner_fan()
    sage: g == loads(dumps(g))
    True

REFERENCES:

- Anders N. Jensen; Gfan, a software system for Groebner fans;
  available at
  http://www.math.tu-berlin.de/~jensen/software/gfan/gfan.html
"""
from __future__ import print_function
from __future__ import absolute_import

import string
import pexpect
from subprocess import PIPE, Popen

from sage.misc.sage_eval import sage_eval

from sage.structure.sage_object import SageObject
from sage.interfaces.gfan import gfan
from .multi_polynomial_ideal import is_MPolynomialIdeal
from .polynomial_ring_constructor import PolynomialRing
from sage.rings.rational_field import QQ
from sage.rings.integer import Integer
from sage.rings.integer_ring import ZZ
from sage.modules.free_module_element import vector
from sage.plot.all import line, Graphics, polygon
from sage.plot.plot3d.shapes2 import line3d
from sage.geometry.polyhedron.constructor import Polyhedron
from sage.geometry.fan import Fan
from sage.matrix.constructor import matrix
from sage.misc.misc_c import prod

def prefix_check(str_list):
    """
    Checks if any strings in a list are prefixes of another string in
    the list.

    EXAMPLES::

        sage: from sage.rings.polynomial.groebner_fan import prefix_check
        sage: prefix_check(['z1','z1z1'])
        False
        sage: prefix_check(['z1','zz1'])
        True
    """
    for index1 in range(len(str_list)):
        for index2 in range(len(str_list)):
            string1 = str_list[index1]
            string2 = str_list[index2]
            if index1 != index2 and str_list[index1][0:len(string2)].find(string2) != -1:
                return False
    return True

def max_degree(list_of_polys):
    """
    Computes the maximum degree of a list of polynomials

    EXAMPLES::

        sage: from sage.rings.polynomial.groebner_fan import max_degree
        sage: R.<x,y> = PolynomialRing(QQ,2)
        sage: p_list = [x^2-y,x*y^10-x]
        sage: max_degree(p_list)
        11.0
    """
    return max([float(qf.degree()) for qf in list_of_polys])

def _cone_parse(fan_dict_cone):
    """
    Utility function that parses cone information into a dict indexed by
    dimension.

    INPUT:

    -  ``fan_dict_cone`` - the value of a fan_dict with
       key 'CONES'

    EXAMPLES::

        sage: R.<x,y,z,w> = PolynomialRing(QQ,4)
        sage: I = R.ideal([x^2+y^2+z^2-1,x^4-y^2-z-1,x+y+z,w+x+y])
        sage: GF = I.groebner_fan()
        sage: TI = GF.tropical_intersection()
        sage: from sage.rings.polynomial.groebner_fan import _cone_parse
        sage: _cone_parse(TI.fan_dict['CONES'])
        {1: [[0], [1], [2], [3], [4]], 2: [[2, 3]]}
    """
    cone_dict = {}
    cur_dim = 0
    for item in fan_dict_cone:
        temp = item.split('{')[-1]
        if temp.split('}') == '':
            temp = ['']
        else:
            temp = temp.split('}')[0]
            temp = temp.split(' ')
        if item.find('Dimension') != -1:
            cur_dim = Integer(item.split(' ')[-1])
            if cur_dim > 0:
                cone_dict[cur_dim] = [[Integer(q) for q in temp if q != '']]
        else:
            if cur_dim > 0: cone_dict[cur_dim] += [[Integer(q) for q in temp if q != '']]
    return cone_dict

class PolyhedralCone(SageObject):

    def __init__(self, gfan_polyhedral_cone, ring = QQ):
        """
        Converts polymake/gfan data on a polyhedral cone into a sage class.
        Currently (18-03-2008) needs a lot of work.

        EXAMPLES::

            sage: R3.<x,y,z> = PolynomialRing(QQ,3)
            sage: gf = R3.ideal([x^8-y^4,y^4-z^2,z^2-2]).groebner_fan()
            sage: a = gf[0].groebner_cone()
            sage: a.facets()
            [[0, 0, 1], [0, 1, 0], [1, 0, 0]]
        """
        cone_keys = ['AMBIENT_DIM','DIM','IMPLIED_EQUATIONS', 'LINEALITY_DIM', 'LINEALITY_SPACE','FACETS', 'RELATIVE_INTERIOR_POINT']
        poly_lines = gfan_polyhedral_cone.split('\n')
        self.cone_dict = {}
        key_ind = 0
        cur_key = None
        for ting in poly_lines:
            if cone_keys.count(ting) > 0:
                cur_key = ting
                self.cone_dict[cur_key] = []
            elif cur_key and ting != '':
                self.cone_dict[cur_key].append(ting)
        self._facets = []
        for facet in self.cone_dict['FACETS']:
            temp_facet = facet.split('\t')[0]
            temp_facet = temp_facet.split(' ')
            temp_facet = [int(x) for x in temp_facet]
            self._facets.append(temp_facet)
        self._ambient_dim = int(self.cone_dict['AMBIENT_DIM'][0])
        self._dim = int(self.cone_dict['DIM'][0])
        self._lineality_dim = int(self.cone_dict['LINEALITY_DIM'][0])
        rel_int_pt_str = self.cone_dict['RELATIVE_INTERIOR_POINT'][0]
        self._relative_interior_point = [int(q) for q in rel_int_pt_str.split(' ')]

    def _repr_(self):
        """
        Returns a basic description of the polyhedral cone.

        EXAMPLES::

            sage: R3.<x,y,z> = PolynomialRing(QQ,3)
            sage: gf = R3.ideal([x^8-y^4,y^4-z^2,z^2-2]).groebner_fan()
            sage: a = gf[0].groebner_cone()
            sage: a # indirect doctests
            Polyhedral cone in 3 dimensions of dimension 3
        """
        return "Polyhedral cone in %s dimensions of dimension %s"%(str(self.ambient_dim()), str(self.dim()))

    def facets(self):
        """
        Returns the inward facet normals of the Groebner cone.

        EXAMPLES::

            sage: R3.<x,y,z> = PolynomialRing(QQ,3)
            sage: gf = R3.ideal([x^8-y^4,y^4-z^2,z^2-2]).groebner_fan()
            sage: a = gf[0].groebner_cone()
            sage: a.facets()
            [[0, 0, 1], [0, 1, 0], [1, 0, 0]]
        """
        return self._facets

    def ambient_dim(self):
        """
        Returns the ambient dimension of the Groebner cone.

        EXAMPLES::

            sage: R3.<x,y,z> = PolynomialRing(QQ,3)
            sage: gf = R3.ideal([x^8-y^4,y^4-z^2,z^2-2]).groebner_fan()
            sage: a = gf[0].groebner_cone()
            sage: a.ambient_dim()
            3
        """
        return self._ambient_dim

    def dim(self):
        """
        Returns the dimension of the Groebner cone.

        EXAMPLES::

            sage: R3.<x,y,z> = PolynomialRing(QQ,3)
            sage: gf = R3.ideal([x^8-y^4,y^4-z^2,z^2-2]).groebner_fan()
            sage: a = gf[0].groebner_cone()
            sage: a.dim()
            3
        """
        return self._dim

    def lineality_dim(self):
        """
        Returns the lineality dimension of the Groebner cone. This is
        just the difference between the ambient dimension and the dimension
        of the cone.

        EXAMPLES::

            sage: R3.<x,y,z> = PolynomialRing(QQ,3)
            sage: gf = R3.ideal([x^8-y^4,y^4-z^2,z^2-2]).groebner_fan()
            sage: a = gf[0].groebner_cone()
            sage: a.lineality_dim()
            0
        """
        return self._lineality_dim

    def relative_interior_point(self):
        """
        Returns a point in the relative interior of the Groebner cone.

        EXAMPLES::

            sage: R3.<x,y,z> = PolynomialRing(QQ,3)
            sage: gf = R3.ideal([x^8-y^4,y^4-z^2,z^2-2]).groebner_fan()
            sage: a = gf[0].groebner_cone()
            sage: a.relative_interior_point()
            [1, 1, 1]
        """
        return self._relative_interior_point

class PolyhedralFan(SageObject):
    def __init__(self, gfan_polyhedral_fan, parameter_indices = []):
        """
        Converts polymake/gfan data on a polyhedral fan into a sage class.

        INPUT:

        -  ``gfan_polyhedral_fan`` - output from gfan of a
           polyhedral fan.

        EXAMPLES::

            sage: R.<x,y,z> = PolynomialRing(QQ,3)
            sage: i2 = ideal(x*z + 6*y*z - z^2, x*y + 6*x*z + y*z - z^2, y^2 + x*z + y*z)
            sage: gf2 = i2.groebner_fan(verbose = False)
            sage: pf = gf2.polyhedralfan()
            sage: pf.rays()
            [[-1, 0, 1], [-1, 1, 0], [1, -2, 1], [1, 1, -2], [2, -1, -1]]
        """
        fan_keys = ['AMBIENT_DIM','DIM','LINEALITY_DIM','RAYS','N_RAYS',
                    'LINEALITY_SPACE','ORTH_LINEALITY_SPACE','F_VECTOR',
                    'CONES','MAXIMAL_CONES','PURE','SIMPLICIAL','MULTIPLICITIES']
        poly_lines = gfan_polyhedral_fan.split('\n')
        self.fan_dict = {}
        key_ind = 0
        cur_key = None
        for ting in poly_lines:
            if fan_keys.count(ting) > 0:
                cur_key = ting
                self.fan_dict[cur_key] = []
            elif cur_key and ting != '':
                self.fan_dict[cur_key].append(ting)
        self._ambient_dim = int(self.fan_dict['AMBIENT_DIM'][0])
        self._dim = int(self.fan_dict['DIM'][0])
        self._lineality_dim = int(self.fan_dict['LINEALITY_DIM'][0])
        self._rays = []
        for ray in self.fan_dict['RAYS']:
            temp_ray = ray.split('\t')[0]
            temp_ray = temp_ray.split(' ')
            temp_ray = [int(x) for x in temp_ray]
            if parameter_indices != []:
                for q in parameter_indices:
                    temp_ray = temp_ray[0:q] + [0] + temp_ray[q:]
            self._rays.append(temp_ray)
        self._cone_dict = _cone_parse(self.fan_dict['CONES'])
        self._maximal_cone_dict = _cone_parse(self.fan_dict['MAXIMAL_CONES'])
        self._str = gfan_polyhedral_fan

    def _repr_(self):
        """
        Returns a basic description of the polyhedral fan.

        EXAMPLES::

            sage: R3.<x,y,z> = PolynomialRing(QQ,3)
            sage: gf = R3.ideal([x^8-y^4,y^4-z^2,z^2-2]).groebner_fan()
            sage: pf = gf.polyhedralfan()
            sage: pf # indirect doctest
            Polyhedral fan in 3 dimensions of dimension 3
        """
        return "Polyhedral fan in %s dimensions of dimension %s"%(str(self.ambient_dim()), str(self.dim()))

    def _str_(self):
        r"""
        Returns the raw output of gfan as a string. This should only be
        needed internally as all relevant output is converted to sage
        objects.

        EXAMPLES::

            sage: R3.<x,y,z> = PolynomialRing(QQ,3)
            sage: gf = R3.ideal([x^8-y^4,y^4-z^2,z^2-2]).groebner_fan()
            sage: pf = gf.polyhedralfan()
            sage: pf._str_()
            '_application fan\n_version 2.2\n_type SymmetricFan\n\nAMBIENT_DIM\n3\n\nDIM\n3\n\nLINEALITY_DIM\n0\n\nRAYS\n0 0 1\t# 0\n0 1 0\t# 1\n1 0 0\t# 2\n\nN_RAYS\n3\n\nLINEALITY_SPACE\n\nORTH_LINEALITY_SPACE\n1 0 0\n0 1 0\n0 0 1\n\nF_VECTOR\n1 3 3 1\n\nSIMPLICIAL\n1\n\nPURE\n1\n\nCONES\n{}\t# Dimension 0\n{0}\t# Dimension 1\n{1}\n{2}\n{0 1}\t# Dimension 2\n{0 2}\n{1 2}\n{0 1 2}\t# Dimension 3\n\nMAXIMAL_CONES\n{0 1 2}\t# Dimension 3\n'
        """
        return self._str

    def ambient_dim(self):
        """
        Returns the ambient dimension of the Groebner fan.

        EXAMPLES::

            sage: R3.<x,y,z> = PolynomialRing(QQ,3)
            sage: gf = R3.ideal([x^8-y^4,y^4-z^2,z^2-2]).groebner_fan()
            sage: a = gf.polyhedralfan()
            sage: a.ambient_dim()
            3
        """
        return self._ambient_dim

    def dim(self):
        """
        Returns the dimension of the Groebner fan.

        EXAMPLES::

            sage: R3.<x,y,z> = PolynomialRing(QQ,3)
            sage: gf = R3.ideal([x^8-y^4,y^4-z^2,z^2-2]).groebner_fan()
            sage: a = gf.polyhedralfan()
            sage: a.dim()
            3
        """
        return self._dim

    def lineality_dim(self):
        """
        Returns the lineality dimension of the fan. This is the
        dimension of the largest subspace contained in the fan.

        EXAMPLES::

            sage: R3.<x,y,z> = PolynomialRing(QQ,3)
            sage: gf = R3.ideal([x^8-y^4,y^4-z^2,z^2-2]).groebner_fan()
            sage: a = gf.polyhedralfan()
            sage: a.lineality_dim()
            0
        """
        return self._lineality_dim

    def rays(self):
        """
        A list of rays of the polyhedral fan.

        EXAMPLES::

            sage: R.<x,y,z> = PolynomialRing(QQ,3)
            sage: i2 = ideal(x*z + 6*y*z - z^2, x*y + 6*x*z + y*z - z^2, y^2 + x*z + y*z)
            sage: gf2 = i2.groebner_fan(verbose = False)
            sage: pf = gf2.polyhedralfan()
            sage: pf.rays()
            [[-1, 0, 1], [-1, 1, 0], [1, -2, 1], [1, 1, -2], [2, -1, -1]]
        """
        return sorted(self._rays)

    def cones(self):
        """
        A dictionary of cones in which the keys are the cone dimensions.
        For each dimension, the value is a list of the cones,
        where each element consists of a list of ray indices.

        EXAMPLES::

            sage: R.<x,y,z> = QQ[]
            sage: f = 1+x+y+x*y
            sage: I = R.ideal([f+z*f, 2*f+z*f, 3*f+z^2*f])
            sage: GF = I.groebner_fan()
            sage: PF = GF.tropical_intersection()
            sage: PF.cones()
            {1: [[0], [1], [2], [3], [4], [5]], 2: [[0, 1], [0, 2], [0, 3], [0, 4], [1, 2], [1, 3], [2, 4], [3, 4], [1, 5], [2, 5], [3, 5], [4, 5]]}
        """
        return self._cone_dict

    def maximal_cones(self):
        """
        A dictionary of the maximal cones in which the keys are the
        cone dimensions.  For each dimension, the value is a list of
        the maximal cones, where each element consists of a list of ray indices.

        EXAMPLES::

            sage: R.<x,y,z> = QQ[]
            sage: f = 1+x+y+x*y
            sage: I = R.ideal([f+z*f, 2*f+z*f, 3*f+z^2*f])
            sage: GF = I.groebner_fan()
            sage: PF = GF.tropical_intersection()
            sage: PF.maximal_cones()
            {2: [[0, 1], [0, 2], [0, 3], [0, 4], [1, 2], [1, 3], [2, 4], [3, 4], [1, 5], [2, 5], [3, 5], [4, 5]]}
        """
        return self._maximal_cone_dict

    def f_vector(self):
        """
        The f-vector of the fan.

        EXAMPLES::

            sage: R.<x,y,z> = QQ[]
            sage: f = 1+x+y+x*y
            sage: I = R.ideal([f+z*f, 2*f+z*f, 3*f+z^2*f])
            sage: GF = I.groebner_fan()
            sage: PF = GF.tropical_intersection()
            sage: PF.f_vector()
            [1, 6, 12]
        """
        str_data = self.fan_dict['F_VECTOR'][0]
        fv = [Integer(x) for x in str_data.split(' ')]
        return fv

    def is_simplicial(self):
        """
        Whether the fan is simplicial or not.

        EXAMPLES::

            sage: R.<x,y,z> = QQ[]
            sage: f = 1+x+y+x*y
            sage: I = R.ideal([f+z*f, 2*f+z*f, 3*f+z^2*f])
            sage: GF = I.groebner_fan()
            sage: PF = GF.tropical_intersection()
            sage: PF.is_simplicial()
            True
        """
        return bool(int(self.fan_dict['SIMPLICIAL'][0]))

    def to_RationalPolyhedralFan(self):
        """
        Converts to the RationalPolyhedralFan class, which is more actively
        maintained.  While the information in each class is essentially the
        same, the methods and implementation are different.

        EXAMPLES::

            sage: R.<x,y,z> = QQ[]
            sage: f = 1+x+y+x*y
            sage: I = R.ideal([f+z*f, 2*f+z*f, 3*f+z^2*f])
            sage: GF = I.groebner_fan()
            sage: PF = GF.tropical_intersection()
            sage: fan = PF.to_RationalPolyhedralFan()
            sage: [tuple(q.facet_normals()) for q in fan]
            [(M(0, -1, 0), M(-1, 0, 0)), (M(0, 0, -1), M(-1, 0, 0)), (M(0, 0, 1), M(-1, 0, 0)), (M(0, 1, 0), M(-1, 0, 0)), (M(0, 0, -1), M(0, -1, 0)), (M(0, 0, 1), M(0, -1, 0)), (M(0, 1, 0), M(0, 0, -1)), (M(0, 1, 0), M(0, 0, 1)), (M(1, 0, 0), M(0, -1, 0)), (M(1, 0, 0), M(0, 0, -1)), (M(1, 0, 0), M(0, 0, 1)), (M(1, 0, 0), M(0, 1, 0))]

        Here we use the RationalPolyhedralFan's Gale_transform method on a tropical
        prevariety.

        .. link

        ::

            sage: fan.Gale_transform()
            [ 1  0  0  0  0  1 -2]
            [ 0  1  0  0  1  0 -2]
            [ 0  0  1  1  0  0 -2]

        """
        try:
            return self._fan
        except AttributeError:
            cdnt = []
            cones = self.cones()
            for x in cones:
                if x > 1:
                    cdnt = cdnt + cones[x]
            fan = Fan(cones=cdnt, rays = self.rays(), discard_faces=True)
            self._fan = fan
            return self._fan

class InitialForm(SageObject):
    def __init__(self, cone, rays, initial_forms):
        """
        A system of initial forms from a polynomial system.
        To each form is associated a cone and a list of
        polynomials (the initial form system itself).

        This class is intended for internal use inside of the
        TropicalPrevariety class.

        EXAMPLES::

            sage: from sage.rings.polynomial.groebner_fan import InitialForm
            sage: R.<x,y> = QQ[]
            sage: inform = InitialForm([0], [[-1, 0]], [y^2 - 1, y^2 - 2, y^2 - 3])
            sage: inform._cone
            [0]
        """
        self._cone = cone
        self._rays = rays
        self._initial_forms = initial_forms

    def cone(self):
        """
        The cone associated with the initial form system.

        EXAMPLES::

            sage: R.<x,y> = QQ[]
            sage: I = R.ideal([(x+y)^2-1,(x+y)^2-2,(x+y)^2-3])
            sage: GF = I.groebner_fan()
            sage: PF = GF.tropical_intersection()
            sage: pfi0 = PF.initial_form_systems()[0]
            sage: pfi0.cone()
            [0]
        """
        return self._cone

    def rays(self):
        """
        The rays of the cone associated with the initial form system.

        EXAMPLES::

            sage: R.<x,y> = QQ[]
            sage: I = R.ideal([(x+y)^2-1,(x+y)^2-2,(x+y)^2-3])
            sage: GF = I.groebner_fan()
            sage: PF = GF.tropical_intersection()
            sage: pfi0 = PF.initial_form_systems()[0]
            sage: pfi0.rays()
            [[-1, 0]]
        """
        return self._rays

    def internal_ray(self):
        """
        A ray internal to the cone associated with the initial form
        system.

        EXAMPLES::

            sage: R.<x,y> = QQ[]
            sage: I = R.ideal([(x+y)^2-1,(x+y)^2-2,(x+y)^2-3])
            sage: GF = I.groebner_fan()
            sage: PF = GF.tropical_intersection()
            sage: pfi0 = PF.initial_form_systems()[0]
            sage: pfi0.internal_ray()
            (-1, 0)
        """
        return sum([vector(q) for q in self.rays()])

    def initial_forms(self):
        """
        The initial forms (polynomials).

        EXAMPLES::

            sage: R.<x,y> = QQ[]
            sage: I = R.ideal([(x+y)^2-1,(x+y)^2-2,(x+y)^2-3])
            sage: GF = I.groebner_fan()
            sage: PF = GF.tropical_intersection()
            sage: pfi0 = PF.initial_form_systems()[0]
            sage: pfi0.initial_forms()
            [y^2 - 1, y^2 - 2, y^2 - 3]
        """
        return self._initial_forms

def verts_for_normal(normal, poly):
    """
    Returns the exponents of the vertices of a newton polytope
    that make up the supporting hyperplane for the given outward
    normal.

    EXAMPLES::

        sage: from sage.rings.polynomial.groebner_fan import verts_for_normal
        sage: R.<x,y,z> = PolynomialRing(QQ,3)
        sage: f1 = x*y*z - 1
        sage: f2 = f1*(x^2 + y^2 + 1)
        sage: verts_for_normal([1,1,1],f2)
        [(3, 1, 1), (1, 3, 1)]
    """
    exps = [tuple(x) for x in poly.exponents()]
    expmat = matrix(exps)
    vals = expmat*vector(QQ,normal)
    maxval = max(vals)
    outverts = []
    for i in range(len(exps)):
        if vals[i] == maxval:
            outverts.append(exps[i])
    return outverts

class TropicalPrevariety(PolyhedralFan):

    def __init__(self, gfan_polyhedral_fan, polynomial_system, poly_ring, parameters = []):
        """
        This class is a subclass of the PolyhedralFan class,
        with some additional methods for tropical prevarieties.

        INPUT:

        -  ``gfan_polyhedral_fan`` - output from gfan of a
           polyhedral fan.
        -  ``polynomial_system`` - a list of polynomials
        -  ``poly_ring`` - the polynomial ring of the list of polynomials
        -  ``parameters`` (optional) - a list of variables to be considered as parameters

        EXAMPLES::

            sage: R.<x,y,z> = QQ[]
            sage: I = R.ideal([(x+y+z)^2-1,(x+y+z)-x,(x+y+z)-3])
            sage: GF = I.groebner_fan()
            sage: TI = GF.tropical_intersection()
            sage: TI._polynomial_system
            [x^2 + 2*x*y + y^2 + 2*x*z + 2*y*z + z^2 - 1, y + z, x + y + z - 3]
        """
        parameter_indices = []
        if parameters != []:
            allvars = poly_ring.gens()
            parameter_indices = [allvars.index(q) for q in parameters]
        PolyhedralFan.__init__(self, gfan_polyhedral_fan, parameter_indices = parameter_indices)
        self._polynomial_system = polynomial_system
        self._parameters = parameters

    def initial_form_systems(self):
        """
        Returns a list of systems of initial forms for each cone
        in the tropical prevariety.

        EXAMPLES::

            sage: R.<x,y> = QQ[]
            sage: I = R.ideal([(x+y)^2-1,(x+y)^2-2,(x+y)^2-3])
            sage: GF = I.groebner_fan()
            sage: PF = GF.tropical_intersection()
            sage: pfi = PF.initial_form_systems()
            sage: for q in pfi:
            ....:     print(q.initial_forms())
            [y^2 - 1, y^2 - 2, y^2 - 3]
            [x^2 - 1, x^2 - 2, x^2 - 3]
            [x^2 + 2*x*y + y^2, x^2 + 2*x*y + y^2, x^2 + 2*x*y + y^2]
        """
        try:
            return self._initial_form_systems
        except AttributeError:
            initial_form_systems = []
            pvars = self._polynomial_system[0].parent().gens()
            nvars = len(pvars)
            for dcone in self.cones():
                for acone in self.cones()[dcone]:
                    rays = [self.rays()[i] for i in acone]
                    repray = sum([vector(q) for q in rays])
                    iforms = []
                    for poly in self._polynomial_system:
                        verts = verts_for_normal(repray, poly)
                        nform = 0
                        for x in verts:
                            factorlist = [pvars[i]**x[i] for i in range(nvars)]
                            temp_monomial = prod(factorlist)
                            nform += poly.monomial_coefficient(temp_monomial)*temp_monomial
                        iforms.append(nform)
                    initial_form_systems.append(InitialForm(acone, rays, iforms))
            self._initial_form_systems = initial_form_systems
            return self._initial_form_systems

def ring_to_gfan_format(input_ring):
    """
    Converts a ring to gfan's format.

    EXAMPLES::

        sage: R.<w,x,y,z> = QQ[]
        sage: from sage.rings.polynomial.groebner_fan import ring_to_gfan_format
        sage: ring_to_gfan_format(R)
        'Q[w, x, y, z]'
        sage: R2.<x,y> = GF(2)[]
        sage: ring_to_gfan_format(R2)
        'Z/2Z[x, y]'
    """
    if input_ring.base_ring() == QQ:
        ring_str = 'Q' + str(input_ring.gens()).replace('(','[').replace(')',']')
    elif input_ring.base_ring() == ZZ:
        ring_str = 'Z' + str(input_ring.gens()).replace('(','[').replace(')',']')
    else:
        ring_str = 'Z/' + str(input_ring.characteristic()) + 'Z'+ str(input_ring.gens()).replace('(','[').replace(')',']')
    return ring_str

def ideal_to_gfan_format(input_ring, polys):
    """
    Return the ideal in gfan's notation.

        EXAMPLES::

            sage: R.<x,y,z> = PolynomialRing(QQ,3)
            sage: polys = [x^2*y - z, y^2*z - x, z^2*x - y]
            sage: from sage.rings.polynomial.groebner_fan import ideal_to_gfan_format
            sage: ideal_to_gfan_format(R, polys)
            'Q[x, y, z]{x^2*y-z,y^2*z-x,x*z^2-y}'

        TESTS:

        Test that #20146 is fixed::

            sage: P = PolynomialRing(QQ,"x11,x12,x13,x14,x15,x21,x22,x23,x24,x25,x31,x32,x33,x34,x35"); x = P.gens(); M = Matrix(3,x)
            sage: I = P.ideal(M.minors(2))
            sage: ideal_to_gfan_format(P,I.gens())
            'Q[x11, x12, x13, x14, x15, x21, x22, x23, x24, x25, x31, x32, x33, x34, x35]{-x12*x21+x11*x22,-x13*x21+x11*x23,-x14*x21+x11*x24,-x15*x21+x11*x25,-x13*x22+x12*x23,-x14*x22+x12*x24,-x15*x22+x12*x25,-x14*x23+x13*x24,-x15*x23+x13*x25,-x15*x24+x14*x25,-x12*x31+x11*x32,-x13*x31+x11*x33,-x14*x31+x11*x34,-x15*x31+x11*x35,-x13*x32+x12*x33,-x14*x32+x12*x34,-x15*x32+x12*x35,-x14*x33+x13*x34,-x15*x33+x13*x35,-x15*x34+x14*x35,-x22*x31+x21*x32,-x23*x31+x21*x33,-x24*x31+x21*x34,-x25*x31+x21*x35,-x23*x32+x22*x33,-x24*x32+x22*x34,-x25*x32+x22*x35,-x24*x33+x23*x34,-x25*x33+x23*x35,-x25*x34+x24*x35}'
    """
    ideal_gen_str = "{" + ",".join(str(poly).replace(" ","").replace("'","") for poly in polys) + "}"
    ring_str = ring_to_gfan_format(input_ring)
    output = ring_str + ideal_gen_str
    return output

class GroebnerFan(SageObject):

    def __init__(self, I, is_groebner_basis=False, symmetry=None, verbose=False):
        """
        This class is used to access capabilities of the program Gfan.  In
        addition to computing Groebner fans, Gfan can compute other things in
        tropical geometry such as tropical prevarieties.

        INPUT:

        -  ``I`` - ideal in a multivariate polynomial ring

        -  ``is_groebner_basis`` - bool (default False). if
           True, then I.gens() must be a Groebner basis with respect to the
           standard degree lexicographic term order.

        -  ``symmetry`` - default: None; if not None, describes
           symmetries of the ideal

        -  ``verbose`` - default: False; if True, printout
           useful info during computations

        EXAMPLES::

            sage: R.<x,y,z> = QQ[]
            sage: I = R.ideal([x^2*y - z, y^2*z - x, z^2*x - y])
            sage: G = I.groebner_fan(); G
            Groebner fan of the ideal:
            Ideal (x^2*y - z, y^2*z - x, x*z^2 - y) of Multivariate Polynomial Ring in x, y, z over Rational Field

        Here is an example of the use of the tropical_intersection command, and then using the RationalPolyhedralFan
        class to compute the Stanley-Reisner ideal of the tropical prevariety::

            sage: R.<x,y,z> = QQ[]
            sage: I = R.ideal([(x+y+z)^3-1,(x+y+z)^3-x,(x+y+z)-3])
            sage: GF = I.groebner_fan()
            sage: PF = GF.tropical_intersection()
            sage: PF.rays()
            [[-1, 0, 0], [0, -1, 0], [0, 0, -1], [1, 1, 1]]
            sage: RPF = PF.to_RationalPolyhedralFan()
            sage: RPF.Stanley_Reisner_ideal(PolynomialRing(QQ,4,'A, B, C, D'))
            Ideal (A*B, A*C, B*C*D) of Multivariate Polynomial Ring in A, B, C, D over Rational Field

        """
        self.__is_groebner_basis = is_groebner_basis
        self.__symmetry = symmetry
        if symmetry:
            print("WARNING! Symmetry option not yet implemented!!")
        self.__verbose = verbose
        if not is_MPolynomialIdeal(I):
            raise TypeError("I must be a multivariate polynomial ideal")
        if prefix_check([str(R_gen) for R_gen in I.ring().gens()]) != True:
            raise RuntimeError("Ring variables cannot contain each other as prefixes")
        S = I.ring()
        R = S.base_ring()
        # todo: add support for ZZ, which only works for bases computation, not tropical intersections
        if not R.is_field():
            raise NotImplementedError("Groebner fan computation only implemented over fields")
        if not (R is QQ or (R.is_finite() and R.is_prime_field() and R.order() <= 32749)):
            # sage: previous_prime (2^15)
            # 32749
            raise NotImplementedError("Groebner fan computation only implemented over Q or GF(p) for p <= 32749.")
        if S.ngens() > 52:
            raise NotImplementedError("Groebner fan computation only implemented for rings in at most 52 variables.")

        self.__ideal = I
        self.__ring = S

    def _repr_(self):
        """
        Describes the Groebner fan and its corresponding ideal.

        EXAMPLES::

            sage: R.<q,u> = PolynomialRing(QQ,2)
            sage: gf = R.ideal([q-u,u^2-1]).groebner_fan()
            sage: gf # indirect doctest
            Groebner fan of the ideal:
            Ideal (q - u, u^2 - 1) of Multivariate Polynomial Ring in q, u over Rational Field
        """
        return "Groebner fan of the ideal:\n%s"%self.__ideal

    def __eq__(self,right):
        """
        Tests equality of Groebner fan objects.

        EXAMPLES::

            sage: R.<q,u> = PolynomialRing(QQ,2)
            sage: gf = R.ideal([q^2-u,u^2-q]).groebner_fan()
            sage: gf2 = R.ideal([u^2-q,q^2-u]).groebner_fan()
            sage: gf.__eq__(gf2)
            True
        """
        return type(self) is type(right) and self.ideal() == right.ideal()

    def ideal(self):
        """
        Return the ideal the was used to define this Groebner fan.

        EXAMPLES::

            sage: R.<x1,x2> = PolynomialRing(QQ,2)
            sage: gf = R.ideal([x1^3-x2,x2^3-2*x1-2]).groebner_fan()
            sage: gf.ideal()
            Ideal (x1^3 - x2, x2^3 - 2*x1 - 2) of Multivariate Polynomial Ring in x1, x2 over Rational Field
        """
        return self.__ideal

    def _gfan_maps(self):
        """
        INPUT: none

        OUTPUT:

        - map from Sage ring to gfan ring

        - map from gfan ring to Sage ring

        EXAMPLES::

            sage: R.<x,y,z> = PolynomialRing(QQ,3)
            sage: G = R.ideal([x^2*y - z, y^2*z - x, z^2*x - y]).groebner_fan()
            sage: G._gfan_maps()
            (Ring morphism:
              From: Multivariate Polynomial Ring in x, y, z over Rational Field
              To:   Multivariate Polynomial Ring in a, b, c over Rational Field
              Defn: x |--> a
                    y |--> b
                    z |--> c,
             Ring morphism:
              From: Multivariate Polynomial Ring in a, b, c over Rational Field
              To:   Multivariate Polynomial Ring in x, y, z over Rational Field
              Defn: a |--> x
                    b |--> y
                    c |--> z)
        """
        try:
            return self.__gfan_maps
        except AttributeError:
            S = self.__ring
            n = S.ngens()

            # Define a polynomial ring in n variables
            # that are named a,b,c,d, ..., z, A, B, C, ...
            R = S.base_ring()
            T = PolynomialRing(R, n, string.ascii_letters[:n])

            # Define the homomorphism that sends the
            # generators of S to the generators of T.
            phi = S.hom(T.gens())

            # Define the homomorphism that sends the
            # generators of T to the generators of S.
            psi = T.hom(S.gens())
            self.__gfan_maps = (phi, psi)
            return self.__gfan_maps

    def _gfan_ring(self):
        """
        Return the ring in gfan's notation

        EXAMPLES::

            sage: R.<x,y,z> = PolynomialRing(QQ,3)
            sage: G = R.ideal([x^2*y - z, y^2*z - x, z^2*x - y]).groebner_fan()
            sage: G._gfan_ring()
            'Q[x, y, z]'
        """
        return ring_to_gfan_format(self.ring())

    def _gfan_ideal(self):
        """
        Return the ideal in gfan's notation.

        EXAMPLES::

            sage: R.<x,y,z> = PolynomialRing(QQ,3)
            sage: G = R.ideal([x^2*y - z, y^2*z - x, z^2*x - y]).groebner_fan()
            sage: G._gfan_ideal()
            'Q[x, y, z]{x^2*y-z,y^2*z-x,x*z^2-y}'
        """
        try:
            return self.__gfan_ideal
        except AttributeError:
            self.__gfan_ideal = ideal_to_gfan_format(self.ring(),self.__ideal.gens())
            return self.__gfan_ideal

    def weight_vectors(self):
        """
        Returns the weight vectors corresponding to the reduced Groebner
        bases.

        EXAMPLES::

            sage: r3.<x,y,z> = PolynomialRing(QQ,3)
            sage: g = r3.ideal([x^3+y,y^3-z,z^2-x]).groebner_fan()
            sage: g.weight_vectors()
            [(3, 7, 1), (5, 1, 2), (7, 1, 4), (5, 1, 4), (1, 1, 1), (1, 4, 8), (1, 4, 10)]
            sage: r4.<x,y,z,w> = PolynomialRing(QQ,4)
            sage: g4 = r4.ideal([x^3+y,y^3-z,z^2-x,z^3 - w]).groebner_fan()
            sage: len(g4.weight_vectors())
            23
        """
        gfan_processes = Popen(['gfan','_weightvector','-m'], stdin = PIPE, stdout=PIPE, stderr=PIPE)
        ans, err = gfan_processes.communicate(input = self.gfan())
        ans = sage_eval(ans.replace('{','').replace('}','').replace('\n',''))
        ans = [vector(QQ,x) for x in ans]
        return ans

    def ring(self):
        """
        Return the multivariate polynomial ring.

        EXAMPLES::

            sage: R.<x1,x2> = PolynomialRing(QQ,2)
            sage: gf = R.ideal([x1^3-x2,x2^3-x1-2]).groebner_fan()
            sage: gf.ring()
            Multivariate Polynomial Ring in x1, x2 over Rational Field
        """
        return self.__ring

    def _gfan_reduced_groebner_bases(self):
        """
        A string of the reduced Groebner bases of the ideal as output by
        gfan.

        EXAMPLES::

            sage: R.<a,b> = PolynomialRing(QQ,2)
            sage: gf = R.ideal([a^3-b^2,b^2-a-1]).groebner_fan()
            sage: gf._gfan_reduced_groebner_bases()
            'Q[a,b]{{b^6-1+2*b^2-3*b^4,a+1-b^2},{b^2-1-a,a^3-1-a}}'
        """
        try:
            return self.__gfan_reduced_groebner_bases
        except AttributeError:
            B = self.gfan(cmd='bases')
            B = B.replace('\n','')
            self.__gfan_reduced_groebner_bases = B
            return B

    def characteristic(self):
        """
        Return the characteristic of the base ring.

        EXAMPLES::

            sage: R.<x,y,z> = PolynomialRing(QQ,3)
            sage: i1 = ideal(x*z + 6*y*z - z^2, x*y + 6*x*z + y*z - z^2, y^2 + x*z + y*z)
            sage: gf = i1.groebner_fan()
            sage: gf.characteristic()
            0
        """
        return self.__ring.characteristic()

    def reduced_groebner_bases(self):
        """
        EXAMPLES::

            sage: R.<x,y,z> = PolynomialRing(QQ, 3, order='lex')
            sage: G = R.ideal([x^2*y - z, y^2*z - x, z^2*x - y]).groebner_fan()
            sage: X = G.reduced_groebner_bases()
            sage: len(X)
            33
            sage: X[0]
            [z^15 - z, y - z^11, x - z^9]
            sage: X[0].ideal()
            Ideal (z^15 - z, y - z^11, x - z^9) of Multivariate Polynomial Ring in x, y, z over Rational Field
            sage: X[:5]
            [[z^15 - z, y - z^11, x - z^9],
            [-y + z^11, y*z^4 - z, y^2 - z^8, x - z^9],
            [-y^2 + z^8, y*z^4 - z, y^2*z^3 - y, y^3 - z^5, x - y^2*z],
            [-y^3 + z^5, y*z^4 - z, y^2*z^3 - y, y^4 - z^2, x - y^2*z],
            [-y^4 + z^2, y^6*z - y, y^9 - z, x - y^2*z]]
            sage: R3.<x,y,z> = PolynomialRing(GF(2477),3)
            sage: gf = R3.ideal([300*x^3-y,y^2-z,z^2-12]).groebner_fan()
            sage: gf.reduced_groebner_bases()
            [[z^2 - 12, y^2 - z, x^3 + 933*y],
            [-y^2 + z, y^4 - 12, x^3 + 933*y],
            [z^2 - 12, -300*x^3 + y, x^6 - 1062*z],
            [-828*x^6 + z, -300*x^3 + y, x^12 + 200]]
        """
        try:
            return self.__reduced_groebner_bases
        except AttributeError:
            G = self._gfan_reduced_groebner_bases()
            if G.find(']') != -1:
                G = G.split(']')[1]
            G = G.replace('{{','').replace('}}','').split('},{')
            S = self.__ring
            X = [ReducedGroebnerBasis(self, [S(f) for f in G[i].split(',')], G[i]) for i in range(len(G))]
            self.__reduced_groebner_bases = X
            return X

    def _gfan_mod(self):
        """
        Return the extra options to the gfan command that are used by this
        object to account for working modulo a prime or in the presence of
        extra symmetries.

        EXAMPLES::

            sage: R.<x,y> = PolynomialRing(QQ,2)
            sage: gf = R.ideal([x^3-y,y^3-x-1]).groebner_fan()
            sage: gf._gfan_mod()
            ''
        """
        try:
            return self.__gfan_mod
        except AttributeError:
            mod = ''
            #p = self.characteristic()
            #if p != 0:
            #    mod += ' --mod %s'%p
            #else:
            #    mod += ''

            if self.__is_groebner_basis:
                mod += ' -g'

            if self.__symmetry:
                mod += ' --symmetry'

            self.__gfan_mod = mod
            return self.__gfan_mod

    def gfan(self, cmd='bases', I=None, format=True):
        r"""
        Returns the gfan output as a string given an input cmd; the default
        is to produce the list of reduced Groebner bases in gfan format.

        EXAMPLES::

            sage: R.<x,y> = PolynomialRing(QQ,2)
            sage: gf = R.ideal([x^3-y,y^3-x-1]).groebner_fan()
            sage: gf.gfan()
            'Q[x,y]\n{{\ny^9-1-y+3*y^3-3*y^6,\nx+1-y^3}\n,\n{\ny^3-1-x,\nx^3-y}\n,\n{\ny-x^3,\nx^9-1-x}\n}\n'
        """
        if I is None:
            I = self._gfan_ideal()
        # todo -- put something in here (?) when self.__symmetry isn't None...
        cmd += self._gfan_mod()
        s = gfan(I, cmd, verbose=self.__verbose, format=format)
        if s.strip() == '{':
            raise RuntimeError("Error running gfan command %s on %s"%(cmd, self))
        return s

    def __iter__(self):
        """
        Returns an iterator for the reduced Groebner bases.

        EXAMPLES::

            sage: R.<x,y> = PolynomialRing(QQ,2)
            sage: gf = R.ideal([x^3-y,y^3-x-1]).groebner_fan()
            sage: a = gf.__iter__()
            sage: next(a)
            [y^9 - 3*y^6 + 3*y^3 - y - 1, -y^3 + x + 1]
        """
        for x in self.reduced_groebner_bases():
            yield x

    def __getitem__(self, i):
        """
        Gets a reduced groebner basis

        EXAMPLES;

        ::

            sage: R4.<w1,w2,w3,w4> = PolynomialRing(QQ,4)
            sage: gf = R4.ideal([w1^2-w2,w2^3-1,2*w3-w4^2,w4^2-w1]).groebner_fan()
            sage: gf[0]
            [w4^12 - 1, -1/2*w4^2 + w3, -w4^4 + w2, -w4^2 + w1]
        """
        return self.reduced_groebner_bases()[i]

    def buchberger(self):
        """
        Computes and returns a lexicographic reduced Groebner basis for the
        ideal.

        EXAMPLES::

            sage: R.<x,y,z> = PolynomialRing(QQ,3)
            sage: G = R.ideal([x - z^3, y^2 - x + x^2 - z^3*x]).groebner_fan()
            sage: G.buchberger()
            [-z^3 + y^2, -z^3 + x]
        """
        try:
            return self.__buchberger
        except AttributeError:
            B = self.gfan(cmd='buchberger')
            if B.find(']') != -1:
                B = B.split(']')[1]
            B = B.replace('}','').replace('{','')
            S = self.__ring
            B = [S(f) for f in B.split(',')]
            self.__buchberger = B
            return B

    def polyhedralfan(self):
        """
        Returns a polyhedral fan object corresponding to the reduced
        Groebner bases.

        EXAMPLES::

            sage: R3.<x,y,z> = PolynomialRing(QQ,3)
            sage: gf = R3.ideal([x^8-y^4,y^4-z^2,z^2-1]).groebner_fan()
            sage: pf = gf.polyhedralfan()
            sage: pf.rays()
            [[0, 0, 1], [0, 1, 0], [1, 0, 0]]
        """
        try:
            return self.__polyhedralfan
        except AttributeError:
            f = self.gfan(cmd='topolyhedralfan', I=self._gfan_reduced_groebner_bases())
            self.__polyhedralfan = PolyhedralFan(f)
            return PolyhedralFan(f)

    def homogeneity_space(self):
        """
        Return the homogeneity space of a the list of polynomials that
        define this Groebner fan.

        EXAMPLES::

            sage: R.<x,y> = PolynomialRing(QQ,2)
            sage: G = R.ideal([y^3 - x^2, y^2 - 13*x]).groebner_fan()
            sage: H = G.homogeneity_space()
        """
        try:
            return self.__homogeneity_space
        except AttributeError:
            h = self.gfan(cmd='homogeneityspace', format=False)
            self.__homogeneity_space = h
            return h

    def render(self, file = None, larger=False, shift=0, rgbcolor = (0,0,0), polyfill = max_degree, scale_colors = True):
        """
        Render a Groebner fan as sage graphics or save as an xfig file.

        More precisely, the output is a drawing of the Groebner fan
        intersected with a triangle. The corners of the triangle are
        (1,0,0) to the right, (0,1,0) to the left and (0,0,1) at the top.
        If there are more than three variables in the ring we extend these
        coordinates with zeros.

        INPUT:

        -  ``file`` - a filename if you prefer the output
           saved to a file. This will be in xfig format.

        -  ``shift`` - shift the positions of the variables in
           the drawing. For example, with shift=1, the corners will be b
           (right), c (left), and d (top). The shifting is done modulo the
           number of variables in the polynomial ring. The default is 0.

        -  ``larger`` - bool (default: False); if True, make
           the triangle larger so that the shape of of the Groebner region
           appears. Affects the xfig file but probably not the sage graphics
           (?)

        -  ``rgbcolor`` - This will not affect the saved xfig
           file, only the sage graphics produced.

        -  ``polyfill`` - Whether or not to fill the cones with
           a color determined by the highest degree in each reduced Groebner
           basis for that cone.

        -  ``scale_colors`` - if True, this will normalize
           color values to try to maximize the range

        EXAMPLES::

            sage: R.<x,y,z> = PolynomialRing(QQ,3)
            sage: G = R.ideal([y^3 - x^2, y^2 - 13*x,z]).groebner_fan()
            sage: test_render = G.render()

        ::

            sage: R.<x,y,z> = PolynomialRing(QQ,3)
            sage: G = R.ideal([x^2*y - z, y^2*z - x, z^2*x - y]).groebner_fan()
            sage: test_render = G.render(larger=True)

        TESTS:

        Testing the case where the number of generators is < 3. Currently,
        this should raise a ``NotImplementedError`` error.

        ::

            sage: R.<x,y> = PolynomialRing(QQ, 2)
            sage: R.ideal([y^3 - x^2, y^2 - 13*x]).groebner_fan().render()
            Traceback (most recent call last):
            ...
            NotImplementedError
        """
        S = self.__ring
        if S.ngens() < 3:
            print("For 2-D fan rendering the polynomial ring must have 3 variables (or more, which are ignored).")
            raise NotImplementedError
        cmd = 'render'
        if shift:
            cmd += ' --shiftVariables %s'%shift
        if larger:
            cmd += ' -L'
        s = self.gfan(cmd, I=self._gfan_reduced_groebner_bases().replace(' ',','), format=False)
        if file is not None:
            open(file,'w').write(s)
        sp = s.split('\n')
        sp2 = []
        for x in sp[9:]:
            xs = x.split(' ')
            y = []
            if x[0:3] != '2 3' and len(xs) > 1:
                for q in xs:
                    if q != '':
                        y.append(q)
                sp2.append(y)
        sp3 = []
        for j in range(len(sp2)):
            temp = []
            for i in range(0,len(sp2[j])-1,2):
                temp.append([float(sp2[j][i])/1200.0, float(sp2[j][i+1])/1200.0])
            sp3.append(temp)
        r_lines = Graphics()
        for x in sp3:
            r_lines = r_lines + line(x, rgbcolor = rgbcolor)
        if polyfill:
            vals = [polyfill(q) for q in self.reduced_groebner_bases()]
            if isinstance(vals[0], list) or isinstance(vals[0], tuple):
                if scale_colors:
                    vmins = [min([q[i] for q in vals]) for i in (0,1,2)]
                    vmaxs = [max([q[i] for q in vals]) for i in (0,1,2)]
                    for i in (0,1,2):
                        if vmaxs[i] == vmins[i]:
                            vmaxs[i] = vmins[i] + .01
                    for index in range(len(sp3)):
                        col = [1-(vals[index][i]-vmins[i])/(vmaxs[i]-vmins[i]) for i in (0,1,2)]
                        r_lines = r_lines + polygon(sp3[index], rgbcolor = col)
                else:
                    for index in range(len(sp3)):
                        r_lines = r_lines + polygon(sp3[index], rgbcolor = vals[i])
            else:
                if scale_colors:
                    vmin = min(vals)
                    vmax = max(vals)
                    if vmin == vmax:
                        vmax = vmin + .01
                    for index in range(len(sp3)):
                        r_lines = r_lines + polygon(sp3[index], hue = .1 + .6*(vals[index]-vmin)/(vmax-vmin))
                else:
                    for index in range(len(sp3)):
                        r_lines = r_lines + polygon(sp3[index], hue = vals[i])
        return r_lines

    def _cone_to_ieq(self, facet_list):
        """
        A simple utility function for converting a facet normal to an
        inequality form.

        EXAMPLES::

            sage: R.<x,y> = PolynomialRing(QQ,2) # dummy stuff to get a gfan object
            sage: gf = R.ideal([x^2+y,y^2]).groebner_fan()
            sage: gf._cone_to_ieq([[1,2,3,4]])
            [[0, 1, 2, 3, 4]]
        """
        ieq_list = []
        for q in facet_list:
            ieq_list.append([0] + q)
        return ieq_list

    def _embed_tetra(self, fpoint):
        """
        Takes a 4-d vector and projects it onto the plane perpendicular to
        (1,1,1,1). Stretches by a factor of 2 as well, since this is only
        for graphical display purposes.

        INPUT:

        -  ``fpoint`` - a list of four numbers

        EXAMPLES::

            sage: R.<x> = PolynomialRing(QQ,1) # dummy stuff to get a gfan object
            sage: gf = R.ideal([x^2]).groebner_fan()
            sage: gf._embed_tetra([1/2,1/2,1/2,1/2])
            [0, 0, 0]
        """
        v1 = [1,1,-1,-1]
        v2 = [1,-1,1,-1]
        v3 = [-1,1,1,-1]
        x1 = sum([fpoint[ind]*v1[ind] for ind in range(4)])
        x2 = sum([fpoint[ind]*v2[ind] for ind in range(4)])
        x3 = sum([fpoint[ind]*v3[ind] for ind in range(4)])
        return [x1,x2,x3]

    def _4d_to_3d(self, polyhedral_data):
        """
        A utility function that takes a list of 4d polytopes, projects them
        to 3d, and returns a list of edges.

        INPUT:

        - ``polyhedral_data`` -- an object with 4d vertex and adjacency
          information

        OUTPUT:

        - ``edges`` -- a list of edges in 3d - each list item is a pair of
          points

        EXAMPLES::

            sage: R4.<w,x,y,z> = PolynomialRing(QQ,4)
            sage: gf = R4.ideal([w^2-x,x^2-y,y^2-z,z^2-1]).groebner_fan()
            sage: g_cone = gf[0].groebner_cone()
            sage: g_cone_facets = g_cone.facets()
            sage: g_cone_ieqs = gf._cone_to_ieq(g_cone_facets)
            sage: cone_data = Polyhedron(ieqs = g_cone_ieqs, eqns = [[1,-1,-1,-1,-1]])
            sage: cone_lines = gf._4d_to_3d(cone_data)
            sage: cone_lines
            [[[1, 1, -1], [1, -1/3, 1/3]],
            [[1, 1, -1], [-1/7, 3/7, 5/7]],
            [[1, 1, -1], [-3/5, -1/3, -1/5]],
            [[1, -1/3, 1/3], [-1/7, 3/7, 5/7]],
            [[1, -1/3, 1/3], [-3/5, -1/3, -1/5]],
            [[-1/7, 3/7, 5/7], [-3/5, -1/3, -1/5]]]
        """
        fpoints = polyhedral_data.vertices()
        tpoints = [self._embed_tetra(q) for q in fpoints]
        edges = []
        for vertex in polyhedral_data.vertices():
            i = vertex.index()
            for adjacent_vertex in vertex.adjacent():
                j = adjacent_vertex.index()
                if j>i:
                    try:
                        edges.append([tpoints[i],tpoints[j]])
                    except Exception:
                        print(adj)
                        print('tpoints: ' + str(tpoints))
                        print('fpoints: ' + str(fpoints))
                        print(adjacent_vertex)
                        print(polyhedral_data.ieqs())
                        raise RuntimeError(adj)
        return edges

    def render3d(self, verbose = False):
        """
        For a Groebner fan of an ideal in a ring with four variables, this
        function intersects the fan with the standard simplex perpendicular
        to (1,1,1,1), creating a 3d polytope, which is then projected into
        3 dimensions. The edges of this projected polytope are returned as
        lines.

        EXAMPLES::

            sage: R4.<w,x,y,z> = PolynomialRing(QQ,4)
            sage: gf = R4.ideal([w^2-x,x^2-y,y^2-z,z^2-x]).groebner_fan()
            sage: three_d = gf.render3d()

        TESTS:

        Now test the case where the number of generators is not 4. Currently,
        this should raise a ``NotImplementedError`` error.

        ::

            sage: P.<a,b,c> = PolynomialRing(QQ, 3, order="lex")
            sage: sage.rings.ideal.Katsura(P, 3).groebner_fan().render3d()
            Traceback (most recent call last):
            ...
            NotImplementedError
        """
        S = self.__ring
        if S.ngens() != 4:
            print("For 3-D fan rendering the polynomial ring must have 4 variables")
            raise NotImplementedError
        g_cones = [q.groebner_cone() for q in self.reduced_groebner_bases()]
        g_cones_facets = [q.facets() for q in g_cones]
        g_cones_ieqs = [self._cone_to_ieq(q) for q in g_cones_facets]
        # Now the cones are intersected with a plane:
        cone_info = [Polyhedron(ieqs = q, eqns = [[1,-1,-1,-1,-1]]) for q in g_cones_ieqs]
        #This is really just for debugging
        if verbose:
            for x in cone_info:
                print(x.inequalities() + ([1,1,0,0,0],[1,0,1,0,0],[1,0,0,1,0],[1,0,0,0,1]))
                print(x.equations())
                print()
        cone_info = [Polyhedron(ieqs = x.inequalities() +
                                ([1,1,0,0,0],[1,0,1,0,0],[1,0,0,1,0],[1,0,0,0,1]),
                                eqns = x.equations()) for x in cone_info]
        all_lines = []
        for cone_data in cone_info:
            try:
                cone_lines = self._4d_to_3d(cone_data)
            except Exception:
                print(cone_data._rays)
                raise RuntimeError
            for a_line in cone_lines:
                all_lines.append(a_line)
        return sum([line3d(a_line) for a_line in all_lines])

    def _gfan_stats(self):
        """
        Return various statistics about this Groebner fan.

        EXAMPLES::

            sage: R.<x,y> = PolynomialRing(QQ)
            sage: G = R.ideal([y^3 - x^2, y^2 - 13*x]).groebner_fan()
            sage: G._gfan_stats()
            {'Dimension of homogeneity space': 0,
             'Maximal number of polynomials in Groebner basis': 3,
             'Maximal number of terms in Groebner basis': 6,
             'Maximal total degree of a Groebner basis': 4,
             'Minimal total degree of a Groebner basis': 2,
             'Number of reduced Groebner bases': 3,
             'Number of variables': 2}
        """
        try:
            return self.__stats
        except AttributeError:
            s = self.gfan(cmd='stats', I=self._gfan_reduced_groebner_bases().replace(' ',','), format=False)
            d = {}
            for v in s.split('\n'):
                if len(v) > 0:
                    a,b = v.split(':')
                    d[a] = ZZ(b)
            self.__stats = d
            return d

    def dimension_of_homogeneity_space(self):
        """
        Return the dimension of the homogeneity space.

        EXAMPLES::

            sage: R.<x,y> = PolynomialRing(QQ,2)
            sage: G = R.ideal([y^3 - x^2, y^2 - 13*x]).groebner_fan()
            sage: G.dimension_of_homogeneity_space()
            0
        """
        return self._gfan_stats()['Dimension of homogeneity space']

    def maximal_total_degree_of_a_groebner_basis(self):
        """
        Return the maximal total degree of any Groebner basis.

        EXAMPLES::

            sage: R.<x,y> = PolynomialRing(QQ,2)
            sage: G = R.ideal([y^3 - x^2, y^2 - 13*x]).groebner_fan()
            sage: G.maximal_total_degree_of_a_groebner_basis()
            4
        """
        return self._gfan_stats()['Maximal total degree of a Groebner basis']

    def minimal_total_degree_of_a_groebner_basis(self):
        """
        Return the minimal total degree of any Groebner basis.

        EXAMPLES::

            sage: R.<x,y> = PolynomialRing(QQ,2)
            sage: G = R.ideal([y^3 - x^2, y^2 - 13*x]).groebner_fan()
            sage: G.minimal_total_degree_of_a_groebner_basis()
            2
        """
        return self._gfan_stats()['Minimal total degree of a Groebner basis']

    def number_of_reduced_groebner_bases(self):
        """
        Return the number of reduced Groebner bases.

        EXAMPLES::

            sage: R.<x,y> = PolynomialRing(QQ,2)
            sage: G = R.ideal([y^3 - x^2, y^2 - 13*x]).groebner_fan()
            sage: G.number_of_reduced_groebner_bases()
            3
        """
        return self._gfan_stats()['Number of reduced Groebner bases']

    def number_of_variables(self):
        """
        Return the number of variables.

        EXAMPLES::

            sage: R.<x,y> = PolynomialRing(QQ,2)
            sage: G = R.ideal([y^3 - x^2, y^2 - 13*x]).groebner_fan()
            sage: G.number_of_variables()
            2

        ::

            sage: R = PolynomialRing(QQ,'x',10)
            sage: R.inject_variables(globals())
            Defining x0, x1, x2, x3, x4, x5, x6, x7, x8, x9
            sage: G = ideal([x0 - x9, sum(R.gens())]).groebner_fan()
            sage: G.number_of_variables()
            10
        """
        return self.__ring.ngens()

    def tropical_basis(self, check=True, verbose = False):
        """
        Return a tropical basis for the tropical curve associated to this
        ideal.

        INPUT:

        -  ``check`` - bool (default: True); if True raises a
           ValueError exception if this ideal does not define a tropical curve
           (i.e., the condition that R/I has dimension equal to 1 + the
           dimension of the homogeneity space is not satisfied).

        EXAMPLES::

            sage: R.<x,y,z> = PolynomialRing(QQ,3, order='lex')
            sage: G = R.ideal([y^3-3*x^2, z^3-x-y-2*y^3+2*x^2]).groebner_fan()
            sage: G
            Groebner fan of the ideal:
            Ideal (-3*x^2 + y^3, 2*x^2 - x - 2*y^3 - y + z^3) of Multivariate Polynomial Ring in x, y, z over Rational Field
            sage: G.tropical_basis()
            [-3*x^2 + y^3, 2*x^2 - x - 2*y^3 - y + z^3, 3/4*x + y^3 + 3/4*y - 3/4*z^3]
        """
        try:
            return self.__tropical_basis
        except AttributeError:
            pass
        cmd = 'tropicalbasis'

        I = self.ideal()
        if not I.is_homogeneous():
            cmd += ' -h'
        if check:
            if I.dimension() != 1 + self.dimension_of_homogeneity_space():
                raise ValueError("The ideal does not define a tropical curve.")

        B = self.gfan(cmd)
        if B.find(']') != -1:
                B = B.split(']')[1]
        S = self.__ring
        B = B.replace('\n','')
        B = B.replace('{','').replace('}','').split(',')
        if verbose: print(S, B)
        X = [S(f) for f in B]
        self.__tropical_basis = X
        return X

    def interactive(self, *args, **kwds):
        """
        See the documentation for self[0].interactive(). This does not work
        with the notebook.

        EXAMPLES::

            sage: print("This is not easily doc-testable; please write a good one!")
            This is not easily doc-testable; please write a good one!
        """
        self[0].interactive(*args, **kwds)

    def tropical_intersection(self, parameters = [], symmetry_generators = [], *args, **kwds):
        """
        Returns information about the tropical intersection of the
        polynomials defining the ideal.  This is the common refinement
        of the outward-pointing normal fans of the Newton polytopes of
        the generators of the ideal. Note that some people use the
        inward-pointing normal fans.

        INPUT:

        - ``parameters`` (optional) - a list of variables to be considered as parameters
        - ``symmetry_generators`` (optional) - generators of the symmetry group

        OUTPUT: a TropicalPrevariety object

        EXAMPLES::

            sage: R.<x,y,z> = PolynomialRing(QQ,3)
            sage: I = R.ideal(x*z + 6*y*z - z^2, x*y + 6*x*z + y*z - z^2, y^2 + x*z + y*z)
            sage: gf = I.groebner_fan()
            sage: pf = gf.tropical_intersection()
            sage: pf.rays()
            [[-2, 1, 1]]

            sage: R.<x,y,z> = PolynomialRing(QQ,3)
            sage: f1 = x*y*z - 1
            sage: f2 = f1*(x^2 + y^2 + z^2)
            sage: f3 = f2*(x + y + z - 1)
            sage: I = R.ideal([f1,f2,f3])
            sage: gf = I.groebner_fan()
            sage: pf = gf.tropical_intersection(symmetry_generators = '(1,2,0),(1,0,2)')
            sage: pf.rays()
            [[-2, 1, 1], [1, -2, 1], [1, 1, -2]]

            sage: R.<x,y,z> = QQ[]
            sage: I = R.ideal([(x+y+z)^2-1,(x+y+z)-x,(x+y+z)-3])
            sage: GF = I.groebner_fan()
            sage: TI = GF.tropical_intersection()
            sage: TI.rays()
            [[-1, 0, 0], [0, -1, -1], [1, 1, 1]]
            sage: GF = I.groebner_fan()
            sage: TI = GF.tropical_intersection(parameters=[y])
            sage: TI.rays()
            [[-1, 0, 0]]
        """
        try:
            return self.__tropical_intersection
        except AttributeError:
            cmd = 'tropicalintersection'
            id_str = self._gfan_ideal()
            if parameters != []:
                allvars = self.ring().gens()
                truevars = [q for q in allvars if not q in parameters]
                base_ring = self.ring().base_ring()
                new_ring = PolynomialRing(base_ring, len(truevars), ",".join(str(q) for q in truevars))
                old_polys = self.ideal().gens()
                new_polys = []
                sub = dict([[v,1] for v in parameters])
                for apoly in old_polys:
                    mons = apoly.monomials()
                    mons = [m.subs(sub) for m in mons]
                    new_polys.append(sum(mons))
                id_str = ideal_to_gfan_format(new_ring, new_polys)
            if symmetry_generators != []:
                cmd = cmd + ' --symmetryExploit'
                id_str = id_str + '{' + symmetry_generators + '}'
            f = self.gfan(cmd=cmd, I = id_str)
            pf = TropicalPrevariety(f,self.ideal().gens(), self.ring(), parameters = parameters)
            pf._gfan_output = f
            self.__tropical_intersection = pf
            return pf

    def mixed_volume(self):
        """
        Returns the mixed volume of the generators of this ideal (i.e. this is
        not really an ideal property, it can depend on the generators used).
        The generators must give a square system (as many polynomials as variables).

        EXAMPLES::

            sage: R.<x,y,z> = QQ[]
            sage: example_ideal = R.ideal([x^2-y-1,y^2-z-1,z^2-x-1])
            sage: gf = example_ideal.groebner_fan()
            sage: mv = gf.mixed_volume()
            sage: mv
            8

            sage: R2.<x,y> = QQ[]
            sage: g1 = 1 - x + x^7*y^3 + 2*x^8*y^4
            sage: g2 = 2 + y + 3*x^7*y^3 + x^8*y^4
            sage: example2 = R2.ideal([g1,g2])
            sage: example2.groebner_fan().mixed_volume()
            15
        """
        if len(self.ring().variable_names()) != self.ideal().ngens():
            raise UserWarning('not a square system')
        return Integer(self.gfan(cmd='mixedvolume'))

class ReducedGroebnerBasis(SageObject, list):
    def __init__(self, groebner_fan, gens, gfan_gens):
        """
        A class for representing reduced Groebner bases as produced by
        gfan.

        INPUT:

        -  ``groebner_fan`` - a GroebnerFan object from an
           ideal

        -  ``gens`` - the generators of the ideal

        -  ``gfan_gens`` - the generators as a gfan string

        EXAMPLES::

            sage: R.<a,b> = PolynomialRing(QQ,2)
            sage: gf = R.ideal([a^2-b^2,b-a-1]).groebner_fan()
            sage: from sage.rings.polynomial.groebner_fan import ReducedGroebnerBasis
            sage: ReducedGroebnerBasis(gf,gf[0],gf[0]._gfan_gens())
            [b - 1/2, a + 1/2]
        """
        self.__groebner_fan = groebner_fan
        list.__init__(self, gens)
        self.__gfan_gens = '{' + gfan_gens.replace(' ',',') + '}'
        self.__ring = groebner_fan._gfan_ring()

    def _repr_(self):
        """
        Returns the reduced Groebner basis as a string.

        EXAMPLES::

            sage: R.<z1,zz1> = PolynomialRing(QQ,2)
            sage: gf = R.ideal([z1^2*zz1-1,zz1-2]).groebner_fan()
            sage: rgb1 = gf.reduced_groebner_bases()[0]
            sage: rgb1 # indirect doctest
            [zz1 - 2, z1^2 - 1/2]
        """
        return list.__repr__(self)

    def _gfan_gens(self):
        """
        Returns the reduced Groebner basis as a string in gfan format.

        EXAMPLES::

            sage: R.<z1,zz1> = PolynomialRing(QQ,2)
            sage: gf = R.ideal([z1^2*zz1-1,zz1-2]).groebner_fan()
            sage: rgb1 = gf.reduced_groebner_bases()[0]
            sage: rgb1._gfan_gens()
            '{zz1-2,z1^2-1/2}'
        """
        return self.__gfan_gens

    def _gfan(self):
        """
        Returns a description of the Groebner fan this basis was derived
        from.

        EXAMPLES::

            sage: R.<z1,zz1> = PolynomialRing(QQ,2)
            sage: gf = R.ideal([z1^2*zz1-1,zz1-2]).groebner_fan()
            sage: rgb1 = gf.reduced_groebner_bases()[0]
            sage: rgb1._gfan()
            Groebner fan of the ideal:
            Ideal (z1^2*zz1 - 1, zz1 - 2) of Multivariate Polynomial Ring in z1, zz1 over Rational Field
        """
        return self.__groebner_fan

    def interactive(self, latex=False, flippable=False, wall=False,
                    inequalities=False, weight=False):
        """
        Do an interactive walk of the Groebner fan starting at this reduced
        Groebner basis.

        EXAMPLES::

            sage: R.<x,y> = PolynomialRing(QQ,2)
            sage: G = R.ideal([y^3 - x^2, y^2 - 13*x]).groebner_fan()
            sage: G[0].interactive()      # not tested
            Initializing gfan interactive mode
            *********************************************
            *     Press control-C to return to Sage     *
            *********************************************
            ....
        """
        cmd = 'gfan_interactive'
        if latex:
            cmd += ' -L'
        if flippable:
            cmd += ' -f'
        if wall:
            cmd += ' -w'
        if inequalities:
            cmd += ' -i'
        if weight:
            cmd += ' -W'
        cmd += self.__groebner_fan._gfan_mod()
        E = pexpect.spawn(cmd)
        print("Initializing gfan interactive mode")
        #E.sendline(self._gfan_ideal())
        E.sendline(self.__gfan_gens)
        print("*" * 45)
        print("*     Press control-C to return to Sage     *")
        print("*" * 45)
        try:
            E.interact()
        except OSError:
            print("Returning to Sage.")

    def groebner_cone(self, restrict=False):
        """
        Return defining inequalities for the full-dimensional Groebner cone
        associated to this marked minimal reduced Groebner basis.

        INPUT:

        -  ``restrict`` - bool (default: False); if True, add
           an inequality for each coordinate, so that the cone is restricted
           to the positive orthant.

        OUTPUT: tuple of integer vectors

        EXAMPLES::

            sage: R.<x,y> = PolynomialRing(QQ,2)
            sage: G = R.ideal([y^3 - x^2, y^2 - 13*x]).groebner_fan()
            sage: poly_cone = G[1].groebner_cone()
            sage: poly_cone.facets()
            [[-1, 2], [1, -1]]
            sage: [g.groebner_cone().facets() for g in G]
            [[[0, 1], [1, -2]], [[-1, 2], [1, -1]], [[-1, 1], [1, 0]]]
            sage: G[1].groebner_cone(restrict=True).facets()
            [[-1, 2], [1, -1]]
        """
        try:
            return self.__groebner_cone[restrict]
        except AttributeError:
            self.__groebner_cone = {}
        except KeyError:
            pass
        cmd = 'groebnercone'
        if restrict:
            cmd += ' --restrict'
        gf = self.__groebner_fan
        c = gf.gfan(cmd=cmd, I=self.__ring + self.__gfan_gens)
        return PolyhedralCone(c)

    def ideal(self):
        """
        Return the ideal generated by this basis.

        EXAMPLES::

            sage: R.<x,y,z> = PolynomialRing(QQ,3)
            sage: G = R.ideal([x - z^3, y^2 - 13*x]).groebner_fan()
            sage: G[0].ideal()
            Ideal (-13*z^3 + y^2, -z^3 + x) of Multivariate Polynomial Ring in x, y, z over Rational Field
        """
        return self.__groebner_fan.ring().ideal(self)

max_stars_repo_name: defeo/sage

max_issues_repo_name: defeo/sage

max_forks_repo_name: defeo/sage

reason: import sage,from sage


text:
#
# Stephen Vondenstein, Matthew Buckley
# 10/09/2018
#
import tensorflow as tf
import numpy as np

def iou(prediction, mask, name):
    # Compute the argmax for the output to match mask shape
    prediction = tf.expand_dims(tf.argmax(prediction, axis=-1), axis=-1)
    # Thresholds as specified by competition
    thresholds = np.arange(0.5, 0.95, 0.05, dtype='float32')
    # Compute precion at thresholds
    precision = tf.metrics.precision_at_thresholds(mask, prediction, thresholds)
    # Compute the mean for this step
    mean_iou = tf.reduce_mean(precision, name=name)

    return mean_iou

def cross_entropy(prediction, mask, num_classes, name):
    # Encode mask as one_hot to match prediction
    mask = tf.one_hot(tf.squeeze(mask), num_classes, dtype=tf.float32)
    # Compute binary cross_entropy
    loss = tf.keras.backend.binary_crossentropy(mask, prediction)
    # Compute the mean for this step
    loss = tf.reduce_mean(loss, name=name)

    return loss

ext: py

max_stars_repo_name: neurologists/ionic

max_issues_repo_name: neurologists/ionic

max_forks_repo_name: KaggleBuoys/ionic

reason: import numpy


text:
from numpy import *
from kaplot import *
import mab.parallelize
import mab.utils.numpy

class Test(object):
	def __init__(self, modelpath, dfgrid, light_model, profile_model, solution, storage_3d, orbitweights):
		self.modelpath = modelpath
		self.dfgrid = dfgrid
		self.profile_model = profile_model
		self.light_model = light_model
		self.solution = solution
		self.storage_3d = storage_3d
		self.orbitweights = orbitweights

	def run(self, args):
		self.orbitweights.load()
		self.storage_3d.load()

		ini, scope = mab.gd.schw.configure.loadini(self.solution.modelpath, objectname="galaxy")
		galaxy = scope["galaxy"]
		jeans = galaxy.jeans()
		logr = arange(-3, 3, 0.1)
		r = 10**logr
		#r = arange(0, 1.25, 0.01)
		sigmar = jeans.sigmar(r)
		document(size="30cm,25cm")
		mozaic(3,3,box)
		select(0,0)
		graph(logr, sigmar)
		graph(logr, sigmar*sqrt(1+0.5))
		select(1,0)
		orbitweights = self.orbitweights.orbitweights.reshape(self.dfgrid.n_I2, self.dfgrid.n_I1)
		df = self.orbitweights.df.reshape(self.dfgrid.n_I2, self.dfgrid.n_I1)
		indexedimage(df)

		select(2,1)
		for i in range(self.dfgrid.n_I1):#[::5]:
			red = float(i)/self.dfgrid.n_I1
			for j in range(self.dfgrid.n_I2):#[::4]:
				index = self.dfgrid.dof_index(i, j, 0)
				#print self.storage_3d.moments3d.shape
				#if (i == 15) and (j == 1):
				if 1:
					y = self.storage_3d.moments3d[index,4]
					n = self.storage_3d.moments3d[index,0]
					#y[n>0] /= n[n>0]
					#graph(self.storage_3d.x, y)#/(10**self.storage_3d.x)**3)
				#self.storage_3d.moments[
		select(1,1)
		Epot = self.profile_model.potentialr(r)
		varr_tot = sigmar * 0
		vart_tot = sigmar * 0
		dens_tot = sigmar * 0
		varr_tot_check = self.storage_3d.moments3d[index,0] * 0
		counts = varr_tot_check * 0
		Lmaxs = [self.profile_model.Lmax_at_E(E) for E in self.dfgrid.subgrid.E_borders]
		fs = []
		#self.modelpath = galaxy.modelpath
		filename = os.path.join(self.modelpath, "df", "fE.npy")
		#logger.info("loading f1(E) as: %s" % filename)
		dff = load(filename)

		filename = os.path.join(self.modelpath, "df", "fE-E.npy")
		#logger.info("loading corresponding energy as: %s" % filename)
		dfEs = load(filename)

		for i in range(self.dfgrid.n_I1): #[::5]:
			red = float(i)/self.dfgrid.n_I1
			print red
			for j in range(self.dfgrid.n_I2): #[::4]:
				green = float(j)/self.dfgrid.n_I1
				dither = self.dfgrid.dither
				I1 = i
				I2 = j
				E1 = self.dfgrid.E_borders[I1]
				E2 = self.dfgrid.E_borders[I1+1]
				L11 = self.dfgrid.l_borders[I2] * self.profile_model.Lmax_at_E(E1)
				L12 = self.dfgrid.l_borders[I2+1]*  self.profile_model.Lmax_at_E(E1)
				L21 = self.dfgrid.l_borders[I2] * self.profile_model.Lmax_at_E(E2)
				L22 = self.dfgrid.l_borders[I2+1]*  self.profile_model.Lmax_at_E(E2)
				for d1 in range(dither):
					for d2 in range(dither):
						Lmax1 = Lmaxs[i*dither+d1]
						Lmax2 = Lmaxs[(i)*dither+d1+1]
						E1 = self.dfgrid.subgrid.E_borders[i*dither+d1]
						E2 = self.dfgrid.subgrid.E_borders[(i)*dither+d1+1]
						L11 = self.dfgrid.subgrid.l_borders[j*dither+d2] * Lmax1
						L12 = self.dfgrid.subgrid.l_borders[(j)*dither+d2+1]*  Lmax1
						L21 = self.dfgrid.subgrid.l_borders[j*dither+d2] * Lmax2
						L22 = self.dfgrid.subgrid.l_borders[(j)*dither+d2+1]*  Lmax2

						def int_dens(r, L, E):
							return -2./3. * pi * (2.*E - 2. * Epot-L**2/r**2)**(3/2.) #* r**3 #* r**2#/r**3
							#return 2 * pi * L/r**2 / sqrt(2.*E - 2. * Epot-L**2/r**2)
						def int_varr(r, L, E):
							return -2./15. * pi * (2.*E - 2. * Epot-L**2/r**2)**(5/2.) #* r**3 #* r**2#/r**3
							#return 2 * pi * L/r**2 * sqrt(2.*E - 2. * Epot-L**2/r**2)
						def int_vart(r, L, E):
							return 2./15. * pi * (-4.*E + 4. * Epot - 3*L**2/r**2) * (2.*E - 2. * Epot-L**2/r**2)**(3/2.)
							#return 2 * pi * L/r**2 * sqrt(2.*E - 2. * Epot-L**2/r**2)

						Lmax = (Lmax1 + Lmax2)/2
						L1 = (L11 + L12)/2
						L2 = (L21 + L22)/2

						E = (E1+E2)/2
						L = (L1+L2)/2
						dE = E2 - E1
						dL = L2 - L1

						#mask =
						#continue
						varr = ((int_varr(r, L2, E2) - int_varr(r, L1, E2)) - (int_varr(r, L2, E1) - int_varr(r, L1, E1)))
						#varr = int_varr(r, L, E) * dE * dL
						varr[isnan(varr)] = 0
						vart = ((int_vart(r, L2, E2) - int_vart(r, L1, E2)) - (int_vart(r, L2, E1) - int_vart(r, L1, E1)))
						#vart = int_vart(r, L, E) * dE * dL
						vart[isnan(vart)] = 0
						dens = (int_dens(r, L2, E2) - int_dens(r, L1, E2)) - (int_dens(r, L2, E1) - int_dens(r, L1, E1))
						#dens = int_dens(r, L, E) * dE * dL
						dens[isnan(dens)] = 0
						#if all(varr == 0):
						#	import pdb
						#	pdb.set_trace()
						#print ">", varr
						#if (i == 15) and (j == 1):
						#g = self.dfgrid.profile_model.gEL(E, L)
						if 1:
							y = (varr) * r ** 3  / (dE * dL)  # *r**2#* self.light_model.densityr(r)*r**3
							#graph(logr, y/y.max(), color=Color(green, 0, 0))
							#graph(logr, y, color=Color(green, 0, 0))
							index = self.dfgrid.dof_index(i, j, 0)
							#y = self.storage_3d.moments3d[index,4]/self.storage_3d.moments3d[index,0]
							y = self.storage_3d.moments3d[index,4]
							y = y#/(10**self.storage_3d.x)**3
							#n = self.storage_3d.moments3d[index,0]
							#y[n>0] /= n[n>0]
							#graph(self.storage_3d.x, y/y.max(), color="red")

							#print varr
						Ei = argmin(abs(dfEs-E))
						#print Ei,
						ddE = (dfEs[Ei+1] - dfEs[Ei])
						w = dff[Ei] * galaxy.fL(L) / (ddE) #/ Lmax # /ddE * dE #* L ** (-2*self.anisotropy_beta)
						#w = df[j,i] /(dE * dL)
						varr_tot += varr * w #df[j,i] /(dE * dL)
						vart_tot += vart * w
						dens_tot += dens * w #df[j,i] /(dE * dL)
						#varr_tot += varr * orbitweights[j,i] / (dE * dL * g) #df[j,i] #/ (dE * dL * g)
						#dens_tot += dens * orbitweights[j,i] / (dE * dL * g)#df[j,i]# / (dE * dL * g)
						fs.append(df[j,i])#/(dE * dL))
						index = self.dfgrid.index_to_orbitnr(i, j, 0)
						index = self.dfgrid.n_I2 * i + j
						index = self.dfgrid.n_I1 * j + i
						#print index
						#varr_tot_check += self.storage_3d.moments3d[index,4] * self.orbitweights.orbitweights[index]# * orbitweights[j,i]
						#counts += self.storage_3d.moments3d[index,0] * self.orbitweights.orbitweights[index]#orbitweights[j,i]
		select(0,1)
		#dens_tot /= sum(fs)
		#varr_tot /= sum(fs)
		print "tot f", sum(fs), sqrt(sum(fs)), sum(df)
		nu = self.light_model.densityr(r)/self.light_model.light_profile.M
		y = varr_tot * 1.0
		y[dens_tot>0] /= (dens_tot * r ** 0)[dens_tot>0]
		#y[dens_tot>0] /= (dens_tot * r ** 0)[dens_tot>0]
		yr = y
		graph(logr, (y)**0.5)

		y = vart_tot * 1.0
		y[dens_tot>0] /= (dens_tot * r ** 0)[dens_tot>0]
		yt = y
		yt /= 2
		print y, vart_tot
		graph(logr, (y)**0.5, color="blue")
		#graph(logr, (varr_tot/(nu))**0.5)

		#y = varr_tot_check
		#y[counts>0] /= counts[counts>0]
		#print y**0.5, self.storage_3d.x

		select(0,0)
		graph(logr, (yt)**0.5, color="green")
		graph(logr, (yr)**0.5, color="green")

		#print self.orbitweights.orbitweights.shape, self.storage_3d.moments3d.shape
		m = tensordot(self.orbitweights.orbitweights, self.storage_3d.moments3d, axes=[(0,), (0,)])
		m0 = m[0,:]
		m2r = m[4,:]
		m2r[m0>0] /= m0[m0>0]
		m2t1 = m[5,:]
		m2t1[m0>0] /= m0[m0>0]
		m2t2 = m[6,:]
		m2t2[m0>0] /= m0[m0>0]
		m2t = (m2t1 + m2t2)/2.
		#print m2
		graph(self.storage_3d.x, m2r**0.5, color="blue")
		graph(self.storage_3d.x, m2t**0.5, color="blue")
		#graph(self.storage_3d.x, m2t2**0.5, color="blue", linestyle="dash")
		select(2,0)
		graph(self.storage_3d.x, m0/m0.max(), color="blue")
		graph(self.storage_3d.x, counts/counts.max(), color="red")
		y = dens_tot*r**3
		graph(logr, y/y.max(), color="green")
		beta = 1-yt/yr
		mask = ~isnan(beta)
		graph(logr[mask], beta[mask])
		ylim(-2,1)
		hline(-0.5)
		print dens_tot

		r2 = 10**self.storage_3d.x
		y = self.light_model.densityr(r2)*r2**3
		graph(self.storage_3d.x, y/y.max())

		select(0,2)
		graph(logr, dens_tot)
		select(1,2)
		graph(logr, nu, color="red")
		select(2,2)
		graph(logr, dens_tot/dens_tot.max())
		graph(logr, nu/nu.max(), color="red")

		#solution.solve()

		draw()

	def run_(self, args):
		def pot(r):
			return self.profile_model.potentialr(r)
		def f(r, E, L):
			a = sqrt(2*E - L**2/r**2 - 2*pot(r))
			return 1/r**2 * (E * r * arctan2(L * r * a , (-L**2 + 2*E*r**2 - 2*r**2*pot(r)) ) \
			+ r * arctan2(r * a,  L) * pot(r) \
			 + 0.5 * L * a)
		def f1(r, E, L):
			a = sqrt(2*E - L**2/r**2 - 2*pot(r))
			return 1/r**2 * 1 / a
		def dens1(r, E, L):
			return f1(r, E, L)
		def dens(r, E1, E2, L11, L12, L21, L22):
			L1 = 0.5 * (L11 + L21)
			L2 = 0.5 * (L12 + L22)
			L1 = (L11)
			L2 = (L12)
			return (f(r, E2, L2) - f(r, E2, L1)) - (f(r, E1, L2) - f(r, E1, L1))
			#return (f(r, E2, L12) - f(r, E2, L11)) - (f(r, E1, L12) - f(r, E1, L11))
			#return (f(r, E2, L22) - f(r, E2, L21)) - (f(r, E1, L12) - f(r, E1, L11))
			#return (f(r, E2, L22) - f(r, E1, L21)) - (f(r, E2, L12) - f(r, E1, L11))
		def dens(r, E1, E2, L11, L12, L21, L22):
			dL1 = L12 - L11
			dL2 = L22 - L21
			L1 = L11
			def t(e):
				return ((2*e-L1**2/r**2-2*pot(r))**(3./2) - (-(dL2**2+2*dL2*L1+L1**2-2*e*r**2+2*r**2*pot(r))/r**2)**(3./2)) / (3 * dL2)
			return (t(E2)-t(E1)) * 2 * pi * (dL1 + (dL2-dL1))
		# / r**2
			#dL1 = 0.5 * (L11 + L21)
			#L2 = 0.5 * (L12 + L22)
			#L1 = (L11)
			#L2 = (L12)
			#return (f(r, E2, L2) - f(r, E2, L1)) - (f(r, E1, L2) - f(r, E1, L1))
		print "args", args
		box()
		#vsplit(box)
		self.storage_3d.load()
		if 0:
			I1 = int(args[0])
			I2 = int(args[1])

			E1 = self.dfgrid.E_borders[I1]
			E2 = self.dfgrid.E_borders[I1+1]
			L11 = self.dfgrid.l_borders[I2] * self.profile_model.Lmax_at_E(E1)
			L12 = self.dfgrid.l_borders[I2+1]*  self.profile_model.Lmax_at_E(E1)
			L21 = self.dfgrid.l_borders[I2] * self.profile_model.Lmax_at_E(E2)
			L22 = self.dfgrid.l_borders[I2+1]*  self.profile_model.Lmax_at_E(E2)
			#print E1, E2, L1, L2
			#logrs = arange(-3, 3, 0.1)
			#logrs = arange(-3, 3, 0.1)
			#rs = 10**logrs
			rs = self.storage_3d.x
			rho_total = rs * 0
			dither = self.dfgrid.dither
			for d1 in range(dither):
				for d2 in range(dither):
					E = self.dfgrid.subgrid.Es[I1*dither+d1]
					L = self.dfgrid.subgrid.ls[I2*dither+d2] * self.profile_model.Lmax_at_E(E)
					#rho = dens(rs, E1, E2, L11, L12, L21, L22)
					rho = dens1(rs, E, L)
					rho[isnan(rho)] = 0
					rho_total += rho * rs**2
			rho = rho_total
			#print 2*E1, L1**2/rs**2, 2*pot(rs)
			#print "->", sqrt(2*E1 - L1**2/rs**2 - 2*pot(rs))
			print rho
			#print rs[mask]
			#y = rho[mask]*rs[mask]**2
			y = rho
			y /= sum(y)# * 1.2
			graph(rs, y)
			xlim(min(rs), max(rs))
			ylim(0, max(y))
			orbitnr = self.dfgrid.index_to_orbitnr(I1, I2, 0)
			print orbitnr
			print self.storage_3d.moments3d.shape
			rho = self.storage_3d.moments3d[orbitnr, 0]
			#rho /= sum(rho)
			graph(self.storage_3d.x, rho, color="red", linestyle="dash")
		else:
			logrs = arange(-3, 3, 0.01/4)
			rs = 10**logrs
			#rs = rs = self.storage_3d.x
			#rho_total = rs* 0
			global rho_total
			rho_totals = mab.utils.numpy.mmapzeros((self.dfgrid.n_I1, self.dfgrid.n_I2, len(rs)))
			weights = self.solution.findsolution()
			weights = weights.reshape((self.dfgrid.n_I2, self.dfgrid.n_I1))
			#indexedimage(weights)
			#draw()
			cores = 48
			@mab.parallelize.parallelize(cores=cores, info=info)
			def do(i, j):
					global rho_total
			#for i in range(self.dfgrid.n_I1):
			#	for j in range(self.dfgrid.n_I2):
					if 1:
						E1 = self.dfgrid.E_borders[i]
						E2 = self.dfgrid.E_borders[i+1]
						L11 = self.dfgrid.l_borders[j] * self.profile_model.Lmax_at_E(E1)
						L12 = self.dfgrid.l_borders[j+1]*  self.profile_model.Lmax_at_E(E1)
						L21 = self.dfgrid.l_borders[j] * self.profile_model.Lmax_at_E(E2)
						L22 = self.dfgrid.l_borders[j+1]*  self.profile_model.Lmax_at_E(E2)
						dE = E2-E1
						dL = L22-L21

						dither = self.dfgrid.dither
						for d1 in range(dither):
							for d2 in range(dither):
								E1 = self.dfgrid.subgrid.E_borders[i*dither+d1]
								E2 = self.dfgrid.subgrid.E_borders[(i)*dither+d1+1]
								L11 = self.dfgrid.subgrid.l_borders[j*dither+d2] * self.profile_model.Lmax_at_E(E1)
								L12 = self.dfgrid.subgrid.l_borders[(j)*dither+d2+1]*  self.profile_model.Lmax_at_E(E1)
								L21 = self.dfgrid.subgrid.l_borders[j*dither+d2] * self.profile_model.Lmax_at_E(E2)
								L22 = self.dfgrid.subgrid.l_borders[(j)*dither+d2+1]*  self.profile_model.Lmax_at_E(E2)
								dEd = E2-E1
								dLd = L22-L21
								rho = dens(rs, E1, E2, L11, L12, L21, L22)  #* dLd # *( dE * dL)
								#E = self.dfgrid.subgrid.Es[i*dither+d1]
								#L = self.dfgrid.subgrid.ls[j*dither+d2] * self.profile_model.Lmax_at_E(E)
								#rho = dens1(rs, E, L) # * 1 / (dE * dL / (dEd * dLd))
								rho[isnan(rho)] = 0
								rho = rho
								if sum(rho) > 0:
									rho /= sum(rho)
								rho_totals[i,j] += rho  * rs**1 * weights[j,i]
					else:
						E1 = self.dfgrid.E_borders[i]
						E2 = self.dfgrid.E_borders[i+1]
						L11 = self.dfgrid.l_borders[j] * self.profile_model.Lmax_at_E(E1)
						L12 = self.dfgrid.l_borders[j+1]*  self.profile_model.Lmax_at_E(E1)
						L21 = self.dfgrid.l_borders[j] * self.profile_model.Lmax_at_E(E2)
						L22 = self.dfgrid.l_borders[j+1]*  self.profile_model.Lmax_at_E(E2)
						#dE = E2-E1
						#dL = L22-L21
						#E = self.dfgrid.Es[i]
						#L = self.dfgrid.ls[j] * self.profile_model.Lmax_at_E(E)
						#rho = dens1(rs, E, L) * weights[j,i]# * dE * dL
						#d = quad(lambda r: dens(rs, E1, E2, L11, L12, L21, L22),
						d = dens(rs, E1, E2, L11, L12, L21, L22)
						rho = d * weights[j,i]# *( dE * dL)
						rho[isnan(rho)] = 0
						rho_totals[i,j] += rho * rs**1

			nE = self.dfgrid.n_I1
			nL = self.dfgrid.n_I2
			I1s = [i for i in range(nE) for j in range(nL)] #[:1]
			I2s = [j for i in range(nE) for j in range(nL)] #[:1]
			#print len(I1s), len(dfgrid.subgrid.E_borders), len(dfgrid.subgrid.Es)
			do(I1s, I2s)
			rho = sum(sum(rho_totals, axis=0), axis=0)
			print rho.shape
			rho = rho.reshape((len(rho)/8, 8))
			print rho.shape
			rho = sum(rho, axis=1)
			print rho.shape

			y = rho
			y /= sum(y)*8
			logdr = logrs[1] - logrs[0]
			graph(logrs[::8]+logdr*4, y)
			y = self.light_model.densityr(rs)*rs**3/1e4
			y /= sum(y)
			graph(logrs, y, color="red")
		#xlim(min(logrs), max(logrs))
		#select(1)

		draw()

max_stars_repo_name: maartenbreddels/mab

max_issues_repo_name: maartenbreddels/mab

max_forks_repo_name: maartenbreddels/mab

reason: from numpy


text:
def detect_cells(parameters):
    '''detect cells in volumes'''

    import os
    import h5py
    import time
    import numpy as np
    from types import SimpleNamespace
    from voluseg._steps.step4a import define_blocks
    from voluseg._steps.step4b import process_block_data
    from voluseg._steps.step4c import initialize_block_cells
    from voluseg._steps.step4d import nnmf_sparse
    from voluseg._tools.ball import ball
    from voluseg._tools.constants import hdf
    from voluseg._tools.evenly_parallelize import evenly_parallelize

    # set up spark
    from pyspark.sql.session import SparkSession
    spark = SparkSession.builder.getOrCreate()
    sc = spark.sparkContext

    p = SimpleNamespace(**parameters)

    ball_diam, ball_diam_xyz0 = ball(1.0 * p.diam_cell, p.affine_mat)

    # load timepoints
    fullname_timemean = os.path.join(p.dir_output, 'mean_timeseries')
    with h5py.File(fullname_timemean+hdf, 'r') as file_handle:
        timepoints = file_handle['timepoints'][()]

    # load plane filename
    for color_i in range(p.n_colors):
        fullname_cells = os.path.join(p.dir_output, 'cells%s_clean'%(color_i))
        if os.path.isfile(fullname_cells+hdf):
            continue

        dir_cell = os.path.join(p.dir_output, 'cells', str(color_i))
        os.makedirs(dir_cell, exist_ok=True)

        fullname_volmean = os.path.join(p.dir_output, 'volume%d'%(color_i))
        with h5py.File(fullname_volmean+hdf, 'r') as file_handle:
            volume_mean = file_handle['volume_mean'][()].T
            volume_mask = file_handle['volume_mask'][()].T
            volume_peak = file_handle['volume_peak'][()].T
            if 'n_blocks' in file_handle.keys():
                flag = 0
                n_voxels_cell = file_handle['n_voxels_cell'][()]
                n_blocks = file_handle['n_blocks'][()]
                block_valids = file_handle['block_valids'][()]
                xyz0 = file_handle['block_xyz0'][()]
                xyz1 = file_handle['block_xyz1'][()]
            else:
                flag = 1

        # broadcast volume peaks (for initialization) and volume_mean (for renormalization)
        bvolume_peak = sc.broadcast(volume_peak)
        bvolume_mean = sc.broadcast(volume_mean)

        # dimensions and resolution
        lxyz = volume_mean.shape
        rxyz = np.diag(p.affine_mat)[:3]

        # compute number of blocks (do only once)
        if flag:
            lx, ly, lz = lxyz
            rx, ry, rz = rxyz

            # get number of voxels in cell
            if (lz == 1) or (rz >= p.diam_cell):
                # area of a circle
                n_voxels_cell = np.pi * ((p.diam_cell / 2.0)**2) / (rx * ry)
            else:
                # volume of a cylinder (change to sphere later)
                n_voxels_cell = p.diam_cell * np.pi * ((p.diam_cell / 2.0)**2) / (rx * ry * rz)

            n_voxels_cell = np.round(n_voxels_cell).astype(int)

            # get number of voxels in each cell
            n_blocks, block_valids, xyz0, xyz1 = \
                define_blocks(lx, ly, lz, p.n_cells_block, n_voxels_cell, volume_mask)

            # save number and indices of blocks
            with h5py.File(fullname_volmean+hdf, 'r+') as file_handle:
                file_handle['n_voxels_cell'] = n_voxels_cell
                file_handle['n_blocks'] = n_blocks
                file_handle['block_valids'] = block_valids
                file_handle['block_xyz0'] = xyz0
                file_handle['block_xyz1'] = xyz1

        print('number of blocks, total: %d.'%(block_valids.sum()))

        for ii in np.where(block_valids)[0]:
            try:
                fullname_block = os.path.join(dir_cell, 'block%05d'%(ii))
                with h5py.File(fullname_block+hdf, 'r') as file_handle:
                    if ('completion' in file_handle.keys()) and file_handle['completion'][()]:
                        block_valids[ii] = 0
            except (NameError, OSError):
                pass

        print('number of blocks, remaining: %d.'%(block_valids.sum()))
        ix = np.where(block_valids)[0]
        block_ixyz01 = list(zip(ix, xyz0[ix], xyz1[ix]))

        # detect individual cells with sparse nnmf algorithm
        def detect_cells_block(tuple_i_xyz0_xyz1):
            os.environ['MKL_NUM_THREADS'] = '1'

            ii, xyz0, xyz1 = tuple_i_xyz0_xyz1[1]

            voxel_xyz, voxel_timeseries, peak_idx, voxel_similarity_peak = \
                process_block_data(xyz0, xyz1, parameters, color_i, lxyz, rxyz,
                                   ball_diam, bvolume_mean, bvolume_peak, timepoints)

            n_voxels_block = len(voxel_xyz)                        # number of voxels in block

            voxel_fraction_peak = np.argsort(((voxel_timeseries[peak_idx])**2).mean(1)) / len(peak_idx)
            for fraction in np.r_[1:0:-0.05]:
                try:
                    peak_valids = (voxel_fraction_peak >= (1 - fraction))   # valid voxel indices

                    n_cells = np.round(peak_valids.sum() / (0.5 * n_voxels_cell)).astype(int)
                    print((fraction, n_cells))

                    tic = time.time()
                    voxel_timeseries_valid, voxel_xyz_valid, cell_weight_init_valid, \
                        cell_neighborhood_valid, cell_sparseness = initialize_block_cells(
                            n_voxels_cell, n_voxels_block, n_cells, voxel_xyz, voxel_timeseries, peak_idx,
                            peak_valids, voxel_similarity_peak, lxyz, rxyz, ball_diam, ball_diam_xyz0)
                    print('cell initialization: %.1f minutes.\n'%((time.time() - tic) / 60))

                    tic = time.time()
                    cell_weights_valid, cell_timeseries_valid, d = nnmf_sparse(
                        voxel_timeseries_valid, voxel_xyz_valid, cell_weight_init_valid,
                        cell_neighborhood_valid, cell_sparseness, timepoints=timepoints,
                        miniter=10, maxiter=100, tolfun=1e-3)

                    success = 1
                    print('cell factorization: %.1f minutes.\n'%((time.time() - tic) / 60))
                    break
                except ValueError as msg:
                    print('retrying factorization of block %d: %s'%(ii, msg))
                    success = 0

            # get cell positions and timeseries, and save cell data
            fullname_block = os.path.join(dir_cell, 'block%05d'%(ii))
            with h5py.File(fullname_block+hdf, 'w') as file_handle:
                if success:
                    for ci in range(n_cells):
                        ix = cell_weights_valid[:, ci] > 0
                        xyzi = voxel_xyz_valid[ix]
                        wi = cell_weights_valid[ix, ci]
                        bi = np.sum(wi * bvolume_mean.value[tuple(zip(*xyzi))]) / np.sum(wi)
                        ti = bi * cell_timeseries_valid[ci] / np.mean(cell_timeseries_valid[ci])

                        file_handle['/cell/%05d/xyz'%(ci)] = xyzi
                        file_handle['/cell/%05d/weights'%(ci)] = wi
                        file_handle['/cell/%05d/timeseries'%(ci)] = ti

                file_handle['n_cells'] = n_cells
                file_handle['completion'] = 1

        if block_valids.any():
            evenly_parallelize(block_ixyz01).foreach(detect_cells_block)

        # collect_blocks(color_i, parameters, lxyz)

max_stars_repo_name: jingxlim/voluseg

max_issues_repo_name: jingxlim/voluseg

max_forks_repo_name: jingxlim/voluseg


text:
from PIL import Image
from torch.autograd import Variable
import numpy as np
import matplotlib.pyplot as plt
from examples.classification.cat_dog_net import CatDogNet
from pytorch_toolbox.transformations.image import Resize, Normalize, NumpyImage2Tensor
from pytorch_toolbox.transformations.to_float import ToFloat
from pytorch_toolbox.transformations.compose import Compose
from pytorch_toolbox.probe.activation import compare_activations

if __name__ == '__main__':

    backend = "cpu"     # cpu|cuda
    model_path = "cat_dog_model.pth.tar"

    #
    #   Instantiate model and load checkpoint
    #
    model = CatDogNet()
    model.load(model_path)
    model.eval()
    cat_path = "images/cat.jpg"

    # Load test images and prepare it for input in the loaded network
    cat_img_numpy = np.array(Image.open(cat_path).convert('RGB'))
    # Occluded version of the image
    cat_img_occluded_numpy = cat_img_numpy.copy()
    cat_img_occluded_numpy[121:150, 70:135, :] = 0
    cat_img_occluded_numpy[20:100, 150:210, :] = 0
    cat_img_occluded_numpy[150:240, 140:220, :] = 0

    imagenet_mean = [123, 116, 103]
    imagenet_std = [58, 57, 57]
    transformations = Compose([Resize((128, 128)),
                                ToFloat(),
                                NumpyImage2Tensor(),
                                Normalize(mean=imagenet_mean, std=imagenet_std)])
    cat_img = Variable(transformations(cat_img_numpy).unsqueeze(0))
    cat_occluded_img = Variable(transformations(cat_img_occluded_numpy).unsqueeze(0))

    """
    Show activation difference between two images, in this case we compare a cat and an occluded cat
    """
    # Show input images
    titles = ["Non occluded cat", "Occluded cat"]
    fig, axes = plt.subplots(1, 2)
    axes[0].imshow(cat_img_numpy)
    axes[0].set_title(titles[0])
    axes[1].imshow(cat_img_occluded_numpy)
    axes[1].set_title(titles[1])
    fig.tight_layout()
    plt.show()

    # Show activation difference
    prediction, prediction_occ = compare_activations(model, cat_img, cat_occluded_img,
                                                     error_func=lambda a, b: np.abs(a-b), cmin=0, cmax=2.5)
    plt.show()

    # Show predictions
    prediction = prediction.data.cpu().numpy()
    prediction_occ = prediction_occ.data.cpu().numpy()
    fig, axes = plt.subplots(2)
    predictions = [prediction, prediction_occ]
    for ax, prediction, title in zip(axes, predictions, titles):
        ax.bar(np.arange(2), prediction[0])
        ax.set_ylim([0, 1])
        ax.set_xticks([0, 1])
        ax.set_xticklabels(["cat", "dog"])
        ax.set_ylabel("probability")
        ax.set_title(title)
    plt.show()


text:
from sklearn.cluster import DBSCAN
import numpy as np
import firebase_admin
from firebase_admin import credentials
from firebase_admin import firestore

# Use the application default credentials
cred = credentials.ApplicationDefault()
firebase_admin.initialize_app(cred, {
  'projectId': "covid-cluster",
})

store = firestore.client()

def cluster(data, context):
    store = firestore.client()
    doc_ref = store.collection(u'Locations')
    points = []

    docs = doc_ref.get()
    for doc in docs:
        points.append(np.array(list(doc.to_dict().values())))

    points = np.array(points)

    db = DBSCAN(eps=3, min_samples=15).fit(points)

    core_samples_mask = np.zeros_like(db.labels_, dtype=bool)
    core_samples_mask[db.core_sample_indices_] = True
    labels = db.labels_

    # Number of clusters in labels, ignoring noise if present.
    n_clusters_ = len(set(labels)) - (1 if -1 in labels else 0)

    for i in range(n_clusters_):
        points_of_cluster = points[labels==i,:]
        c = np.mean(points_of_cluster, axis=0)
        store.collection(u'Parties').add({"Longtitude": c[0], "Latitude": c[1]})


text:
import os
from thor.utils.spice import getSPICEKernels
import pytest
import pandas as pd
from astropy import units as u
from astropy.time import Time

from ...testing import testOrbits
from ...utils import KERNELS_DE440
from ...utils import getSPICEKernels
from ...utils import setupSPICE
from ...utils import getMPCObservatoryCodes
from ..state import getObserverState

DATA_DIR = os.path.join(
    os.path.dirname(os.path.abspath(__file__)),
    "../../testing/data"
)

def test_getObserverState_heliocentric():
    """
    Read the test dataset for heliocentric state vectors of each observatory at each observation time.
    Use THOR to find heliocentric state vectors of each observatory at each observation time.
    Compare the resulting state vectors and test how well they agree with the ones pulled from Horizons.
    """
    getSPICEKernels(KERNELS_DE440)
    setupSPICE(KERNELS_DE440, force=True)

    # Make sure the latest version of the MPC observatory codes
    # has been downloaded
    getMPCObservatoryCodes()

    # Read observatory states from test data file
    observer_states_df = pd.read_csv(
        os.path.join(DATA_DIR, "observer_states.csv"),
        index_col=False
    )

    origin = "heliocenter"
    for observatory_code in observer_states_df["observatory_code"].unique():
        observatory_mask = observer_states_df["observatory_code"].isin([observatory_code])

        times = Time(
            observer_states_df[observatory_mask]["mjd_utc"].values,
            scale="utc",
            format="mjd"
        )
        observer_states = observer_states_df[observatory_mask][["x", "y", "z", "vx", "vy", "vz"]].values

        observer_states_thor = getObserverState(
            [observatory_code],
            times,
            origin=origin,
        )
        observer_states_thor = observer_states_thor[["obs_x", "obs_y", "obs_z", "obs_vx", "obs_vy", "obs_vz"]].values

        # Test that each state agrees with Horizons
        # to within the tolerances below
        testOrbits(
            observer_states_thor,
            observer_states,
            position_tol=(20*u.m),
            velocity_tol=(1*u.cm/u.s),
            magnitude=True
        )
    return

def test_getObserverState_barycentric():
    """
    Read the test dataset for barycentric state vectors of each observatory at each observation time.
    Use THOR to find barycentric state vectors of each observatory at each observation time.
    Compare the resulting state vectors and test how well they agree with the ones pulled from Horizons.
    """
    getSPICEKernels(KERNELS_DE440)
    setupSPICE(KERNELS_DE440, force=True)

    # Make sure the latest version of the MPC observatory codes
    # has been downloaded
    getMPCObservatoryCodes()

    # Read observatory states from test data file
    observer_states_df = pd.read_csv(
        os.path.join(DATA_DIR, "observer_states_barycentric.csv"),
        index_col=False
    )

    origin = "barycenter"
    for observatory_code in observer_states_df["observatory_code"].unique():
        observatory_mask = observer_states_df["observatory_code"].isin([observatory_code])

        times = Time(
            observer_states_df[observatory_mask]["mjd_utc"].values,
            scale="utc",
            format="mjd"
        )
        observer_states = observer_states_df[observatory_mask][["x", "y", "z", "vx", "vy", "vz"]].values

        observer_states_thor = getObserverState(
            [observatory_code],
            times,
            origin=origin,
        )
        observer_states_thor = observer_states_thor[["obs_x", "obs_y", "obs_z", "obs_vx", "obs_vy", "obs_vz"]].values

        # Test that each state agrees with Horizons
        # to within the tolerances below
        testOrbits(
            observer_states_thor,
            observer_states,
            position_tol=(20*u.m),
            velocity_tol=(1*u.cm/u.s),
            magnitude=True
        )
    return

def test_getObserverState_raises():

    times = Time([59000], scale="utc", format="mjd")

    with pytest.raises(ValueError):
        # Raise error for incorrect frame
        observer_states = getObserverState(["500"], times, frame="eccliptic")

        # Raise error for incorrect origin
        observer_states = getObserverState(["500"], times, origin="heeliocenter")

    with pytest.raises(TypeError):
        # Raise error for non-astropy time
        observer_states = getObserverState(["500"], times.tdb.mjd)

max_stars_repo_name: KatKiker/thor

max_issues_repo_name: KatKiker/thor

max_forks_repo_name: KatKiker/thor

reason: from astropy