Identische Spalten in numpy array suchen

[Tyrax]

Hallo,

ich versuche mich gerade an einem Programm, das n annihilierende random walker auf einem m-dimensionalen, endlichen lattice simuliert. Für den 2D-Fall plane ich zur Zeit, einen numpy array mit mit 3 Zeilen (zwei für die Koordinaten, eine für die Bezeichnung des jeweiligen random walkers) und n Spalten (je eine pro random walker).

Nach jeden Zeitschritt muss ich prüfen, ob sich random walker am gleichen Ort befinden, also ob Spalten meiner Matrix in den ersten zwei Zeilen identisch sind. Die Frage ist nun: Wie kann ich das effizient (und vielleicht sogar noch elegant) implementieren?

Mein bisheriger Plan sieht so aus:
- bewege jeden random walker (zur Zeit noch spaltenweise via for-Schleife...),
- sortiere die Matrix nach Spalten ( zB via sort(Matrix, axis=-1, kind='mergesort')) ),
- vergleiche aufeinander folgende Spalten in den ersten beiden Zeilen.

Meine Erfahrung mit numpy ist noch recht übersichtlich, daher würde ich mich über Tipps freuen.

Danke und beste Grüße, Tyrax

[BlackJack]

@Tyrax: Wo ist denn da der Vorteil von einem `numpy.array`? Welche Operationen werden da schnell(er) drauf ausgeführt im Gegensatz zu einer „normalen” Implementierung?

[Tyrax]

Hi BlackJack,

die Einträge in dem array sind ganzzahlig. Meinst Du mit der normalen Implementierung Listen? Ich war davon ausgegangen, dass ein numpy.array, den ich auf integer einstelle, effizienter ist als eine flexiblere Liste. Ich hatte zudem im Hinterkopf, ein paar Sachen über Matrixoperationen zu erledigen.

Mit welchen Mitteln würdest Du denn das beschriebene Projekt angehen?

Grüße, Tyrax

[BlackJack]

@Tyrax: Mit normal meine ich ein ordentliches, lesbares Programm, bei dem die Werte nicht auf so komisch auf Listen oder Arrays verteilt werden, sondern zu sinvollen Datentypen zusammengefasst werden. Man wird dann am Ende sicher irgendwo eine Liste mit Walker-Objekten haben. Oder aber ein Wörterbuch das Koordinaten auf Walker abbildet (oder auf Listen von Walkern).

`numpy`-Arrays mit Zahlen sind Listen gegenüber effizienter was den Speicherplatz angeht und in der Laufzeit wenn man die Operationen auf Arrays benutzen kann, die viele Elemente betreffen.

Wenn man anfängt heterogene Daten in Arrays zu speichern, wo man kommentieren muss was welche Zeile/Spalte bedeutet, dann ist Python vielleicht nicht die richtige Sprache für das Problem.

[Tyrax]

@BlackJack:
Vielleicht lerne ich bei der Gelegenheit ja mal was über sinnvolle Datentypen. Eine Matrix in der jede Spalte die Koordinaten eines Objekts stehen erschien mir sinnvoll. Die Matrix gibt mir den vollständigen Zustand des Systems zu diesem Zeitpunkt an. Warum ist das nicht ordentlich, wie geht es besser?

Was die heterogenen Daten angeht: Damit meinst Du vermutlich die dritte Zeile meiner Matrix. Die habe ich eingeführt um beim Umsortieren der Spalten die Identität der random walker verfolgen zu können. Was besseres ist mir noch nicht eingefallen.

Grüße, Tyrax

[BlackJack]

@Tyrax: Die Matrix gibt den Zustand wieder, aber der Code dafür ist halt nicht so schön, und nur mit dem Datenobjekt kann man ohne weitere Erklärung überhaupt nichts anfangen. Es ist halt einfach ein zweidimensionales Array mit ganzen Zahlen. Die Strukturierung ist Wissen, welches in den Funktionen versteckt ist, die damit arbeiten. Du schreibst das Programm um einen Datentyp herum und wählst nicht die Datentypen nach dem gestellten Problem.

Bei heterogen meinte ich alle Zeilen. Die erste Zeile bedeutet etwas anderes als die zweite, auch wenn beides Koordinaten sind.

Das so in Zeilen zu organisieren ist übrigens auch nicht so gut. Damit sieht man die Struktur in der ersten Dimension als eine Sammlung von X-Koordinaten, Y-Koordinaten, und Walker-IDs. Natürlicher wäre es zu sagen man hat eine Sammlung von Walkern, von denen jeder eine X- und eine Y-Koordinate und eine ID besitzt.

Was mir an dem Ansatz mit dem Array nicht gefällt ist, dass die Effizienzfrage vor allem anderen gestellt wird, und dadurch eine *deutlich* andere Lösung heraus kommt als man eigentlich in einer objektorientierten Programmiersprache entwickeln würde. Falls diese Lösung dann zu langsam oder speicherhungrig sein sollte, halte ich Python für die falsche Sprache für das Problem, weil vergleichbarer, objektorientierter Code in anderen, statisch kompilierten Sprachen letztendlich zu einer sehr ähnlichen interenen Repräsentation führt wie das `numpy.array` — aber eben bei deutlich anderem, IMHO besseren Quelltext. Welchem der die Lösung in einem Vokabular formuliert das dem Problem entspricht und nicht eines das sich auf eine interne Repräsentation der Daten bezieht.

[Tyrax]

@BlackJack:
Danke, jetzt vestehe ich, worum es Dir geht. Ich muss gestehen, dass ich außer Python keine Programmiersprache beherrsche (soweit ich das für Python überhaupt sagen kann ). Daher bleibe ich erstmal auch für solche Aufgaben bei Python.

Deinem zweiten Absatz kann ich noch nicht so ganz folgen. Mir erscheint es recht naheliegend, die Spalten einer Matrix als Objekte zu begreifen. Stehe ich in Deinen Augen besser da, wenn ich Zeilen als Objekte definiere (im Sinne der Repräsentation einer Matrix als Liste von Zeilen)?


@all:
Falls noch jemand ein paar Ideen für eine bessere Implementierung hat (gern auch mit eleganterer Struktur), immer her damit. Ich freue mich über Anregungen.

Tyrax

[BlackJack]

@Tyrax: Die erste Dimension sind aber nicht die Spalten sondern die Zeilen. Schreib mal eine ``for``-Schleife die nacheinander alle Werte für die einzelnen Walker ausgibt für jede der beiden Möglichkeiten und vergleiche welche einfacher ist. Und vor allem bei welcher man problemlos alle Werte die zu einem Walker gehören an einen Namen binden kann.

[Tyrax]

@BlackJack:
Hmm ja, verstehe. (Allerdings ist die Sortierung der Spalten einfacher, als die Sortierung der Zeilen. Man sucht ja entlang der Zeilen, um die Spalten zu sortieren.)

[BlackJack]

@Tyrax: Auch das Sortieren ist einfacher wenn man die Daten für die Objekte in Zeilen organisiert. Beziehungsweise ist das vom Quelltext für's sortieren letztendlich egal, denn man ruft ja einfach nur die Sortierfunktion mit der entsprechenden Axennummer als Argument auf, aber intern müssen dann nicht-zusammenhängende Daten vertauscht werden wenn man den Platz von zwei Walkern tauscht. Wo es dann den Unterschied macht, ist das iterieren über die Walker-„Objekte” in einer ``for``-Schleife was mit drei Spalten statt drei Zeilen deutlich einfacher geht.

Ich weiss nicht ob wir hier irgendwie aneinander vorbei reden, aber wenn man sich tatsächlich mal eine Schleife schreibt, die nacheinander etwas mit den Daten von allen Walkern macht, beispielsweise sie ausgeben, dann sehen die beiden Varianten so aus (ungetestet):

Code: Alles auswählen

# 
# Drei Spalten, also pro Zeile ein Walker:
# 
for walker in walkers:
    print 'x: {}, y: {}, id: {}'.format(*walker)

# 
# Drei Zeilen: X-Koordinaten, Y-Koordinaten, IDs:
# 
for i in xrange(len(walkers[0])):
    print 'x: {}, y: {}, id: {}'.format(
        *(walkers[j, i] for j in xrange(len(walkers)))
    )

Irgendwie ist die erste Variante mit den drei Spalten einfacher, oder? Welche Operation ist denn Deiner Meinung nach einfacher wenn man die Daten in drei Zeilen organisiert?

[Tyrax]

Hi BlackJack,

auf http://docs.scipy.org/doc/numpy/referen ... .sort.html steht unter anderem

All the sort algorithms make temporary copies of the data when sorting along any but the last axis. Consequently, sorting along the last axis is faster and uses less space than sorting along any other axis.

Ich gehe davon aus, dass das damit zusammenhängt, dass die Sortierung eines mehrdimensionalen arrays auf einer lexikographischen Ordnung beruht. Im Falle einer 2D-Matrix guckt sich der Algorithmus zur Sortierung der Spalten also zunächst die erste Zeile an.

Zugegeben, dass ist am Ende wieder ein Effizienz-Argument, und ich kann Deinen Ansatz das Programm eher am Problem als an den Datentypen zu orientieren durchaus verstehen.

Besten Dank jedenfalls für die Hilfe, Grüße, Tyrax

[Tyrax]

Hallo nochmal,

ich hab da offensichtlich was missverstanden. Die Sortierungs-Algos für numpy-arrays sortieren nicht lexikographisch sondern jede Zeile (oder Spalte) einzeln. Das bringt mich leider überhaupt nicht weiter, da ich ja die Spalten (oder auch Zeilen) untereinander vergleichen und danach sortieren will.

lexsort scheint genau das richtige zu sein, ich checke aber nicht, wie ich das im obigen Sinne anwenden kann.

Ich kann natürlich ohne numpy eine Liste von Listen schreiben und diese sortieren lassen, doch eigentlich würde ich gerne mit numpy arbeiten (Hilft mir an anderer Stelle). Kann man mit numpy lexikographisch sortieren? Alles, was ich gefunden habe, sah sehr nach workaround aus. Ich muss diese Sortierung sehr oft ausführen...

Danke und beste Grüße, Tyrax

[BlackJack]

Ich weiss nicht ob ich fachlich das umgesetzt habe was hier eigentlich modelliert werden soll, aber hier ist mal ein Random Walk wie man ihn in OOP umsetzen könnte. An den Rändern des Koordinatensystems verlässt ein Walker das System auf der einen Seite und erscheint wieder auf der gegenüber liegenden Seite. Und zu „annihilierend” habe ich auf die Schnelle nur Zusammenfassungen von Papers gefunden und einfach umgesetzt was ich mir darunter vorstelle: Wenn am Ende vom Bewegungsschritt Walker auf der gleichen Position stehen, werden alle alle Beteiligten gekillt.

Code: Alles auswählen

#!/usr/bin/env python
# coding: utf-8
"""Simulation of an annihilating random walk in a finite lattice.
"""
from collections import defaultdict
from itertools import izip
from operator import mul
from random import choice, randrange, sample


class CoordinateSystem(object):
    """N dimensional coordinate system.
    
    The object can be seen as a sequence of all coordinates within the
    coordinate system.
    
    :ivar shape: a non empty sequence of number of points per dimension.
    
    >>> coordinate_system = CoordinateSystem([10, 10, 10])
    >>> coordinate_system.dimension_count
    3
    >>> len(coordinate_system)
    1000
    >>> coordinate_system[0]
    (0, 0, 0)
    >>> coordinate_system[42]
    (2, 4, 0)
    >>> coordinate_system[999]
    (9, 9, 9)
    >>> (0, 0, 0) in coordinate_system
    True
    >>> (42, 23, 0) in coordinate_system
    False
    
    """
    def __init__(self, shape):
        """Initialises a coordinate system of given `shape`.
        
        .. seealso:: :attr:`~CoordinateSystem.shape`
        """
        if not shape or any(n < 0 for n in shape):
            raise ValueError(
                'shape must not be empty or contain negative values'
            )
        self.shape = shape
    
    def __len__(self):
        """Returns number of coordinates in the system."""
        return reduce(mul, self.shape, 1)
    
    def __getitem__(self, index):
        """Gets the coordinate at given `index`."""
        if index < 0:
            raise IndexError('negative index')
        result = list()
        for extend in self.shape:
            index, part = divmod(index, extend)
            result.append(part)
        if index != 0:
            raise IndexError('index too large (>{})'.format(len(self)))
        return tuple(result)
    
    def __contains__(self, coordinate):
        """Tests if given coordinate is within the coordinate system."""
        return (
            len(coordinate) == self.dimension_count
            and all(0 <= n < s for n, s in izip(coordinate, self.shape))
        )
    
    @property
    def dimension_count(self):
        """Number of dimensions of the coordinate system."""
        return len(self.shape)
    
    def get_sample(self, count):
        """Gets a sample of `count` coordinates from the coordinate system."""
        return sample(self, count)


class Walker(object):
    """A random walker in a lattice.
    
    :ivar coordinate_system: coordinate system in which the walker operates.
    :ivar coordinate: current coordinate of the walker.
    :ivar id: ID of the walker.
    """
    def __init__(self, coordinate_system, coordinate, id_):
        if coordinate not in coordinate_system:
            raise ValueError('coordinate system and coordinate do not match')
        self.coordinate_system = coordinate_system
        self.coordinate = coordinate
        self.id = id_
    
    def move(self, delta):
        """Moves the walker by given relative coordinate `delta`.
        
        If the new coordinate would be outside the coordinate system of the
        walker the coordinates ”wrap around”, i.e. if the walker leaves a
        2D grid on the left it reappears on the right.
        
        >>> walker = Walker(CoordinateSystem([10, 10]), (0, 5), 42)
        >>> walker.coordinate
        (0, 5)
        >>> walker.move((0, 2))
        >>> walker.coordinate
        (0, 7)
        >>> walker.move((-1, 0))
        >>> walker.coordinate
        (9, 7)
        >>> walker.move((5, 0))
        >>> walker.coordinate
        (4, 7)
        
        """
        self.coordinate = tuple(
            (i + d) % j
            for i, d, j
            in izip(self.coordinate, delta, self.coordinate_system.shape)
        )
    
    def random_move(self):
        """Makes a step in a random dimension into a random direction.
        
        .. seealso:: :meth:`move`
        """
        dimension_count = len(self.coordinate)
        delta = [0] * dimension_count
        delta[randrange(dimension_count)] = choice([-1, 1])
        self.move(delta)


class Simulation(object):
    """A simulation of the random walkers.
    
    The object has a :meth:`step` method to calculate the next step of the
    simulation and is an iterable over the :class:`Walker` objects living
    in the simulation.
    """
    def __init__(self, walkers):
        self.walkers = walkers
    
    def __len__(self):
        return len(self.walkers)
    
    def __iter__(self):
        return iter(self.walkers)
    
    def step(self):
        """Moves every :class:`Walker` one step in a random direction.
        
        Walkers ending up in the same spot at the end of the step are
        annihilated.
        """
        coordinate2walkers = defaultdict(list)
        for walker in self:
            walker.random_move()
            coordinate2walkers[walker.coordinate].append(walker)
        self.walkers = [
            ws[0] for ws in coordinate2walkers.itervalues() if len(ws) == 1
        ]
    
    @classmethod
    def with_randomly_placed_walkers(cls, coordinate_system, count):
        """Creates a simulation with `count` randomly placed
        :class:`Walker`\s within the given `coordinate_system`.
        """
        return cls(
            [
                Walker(coordinate_system, coordinate, i)
                for i, coordinate
                in enumerate(coordinate_system.get_sample(count))
            ]
        )


def main():
    """Runs a little test simulation."""
    simulation = Simulation.with_randomly_placed_walkers(
        CoordinateSystem([100, 100]), 500
    )
    for i in xrange(100):
        simulation.step()
        print '{0} walkers left in step {1}.'.format(len(simulation), i)
        if len(simulation) <= 1:
            break


if __name__ == '__main__':
    main()

Stimmt die Richtung davon wenigstens grob mit dem überein was Du machen möchtest? Und von welchen Dimensionen und Walkeranzahlen reden wir hier eigentlich?

[Tyrax]

Hi BlackJack,

ich werde einen Moment brauchen, um genau nachzuvollziehen, was Du da programmiert hast. Vielen Dank jedenfalls für die Mühe. Zumindest kann ich anhand des Codes mal meine OOP Fähigkeiten aufbessern und die Klasse selbst an mein Modell anpassen. Ist sowieso höchste Zeit.

Ich würde natürlich schon gerne noch wissen, ob lexikographische Spaltensortierung in einer numpy-Matrix zB mit lexsort möglich ist -- schon aus Prinzip.

Beste Grüße, Tyrax

[Sirius3]

ich konnte der Versuchung nicht wiederstehen, Blackjacks Programm mit numpy umzuschreiben

Code: Alles auswählen

#!/usr/bin/env python
# coding: utf-8
"""Simulation of an annihilating random walk in a finite lattice.
"""
from collections import Counter
from operator import mul
from random import sample
import numpy

class CoordinateSystem(object):
    """N dimensional coordinate system.
   
    The object can be seen as a sequence of all coordinates within the
    coordinate system.
   
    :ivar shape: a non empty sequence of number of points per dimension.
   
    >>> coordinate_system = CoordinateSystem([10, 10, 10])
    >>> coordinate_system.dimension_count
    3
    >>> len(coordinate_system)
    1000
    >>> coordinate_system[0]
    (0, 0, 0)
    >>> coordinate_system[42]
    (2, 4, 0)
    >>> coordinate_system[999]
    (9, 9, 9)
    >>> (0, 0, 0) in coordinate_system
    True
    >>> (42, 23, 0) in coordinate_system
    False
   
    """
    def __init__(self, shape):
        """Initialises a coordinate system of given `shape`.
       
        .. seealso:: :attr:`~CoordinateSystem.shape`
        """
        if not shape or any(n < 0 for n in shape):
            raise ValueError(
                'shape must not be empty or contain negative values'
            )
        self.shape = shape
   
    def __len__(self):
        """Returns number of coordinates in the system."""
        return reduce(mul, self.shape, 1)
   
    def __getitem__(self, index):
        """Gets the coordinate at given `index`."""
        if index < 0:
            raise IndexError('negative index')
        result = list()
        for extend in self.shape:
            index, part = divmod(index, extend)
            result.append(part)
        if index != 0:
            raise IndexError('index too large (>{})'.format(len(self)))
        return tuple(result)
   
    def __contains__(self, coordinate):
        """Tests if given coordinate is within the coordinate system."""
        return (
            len(coordinate) == self.dimension_count
            and all(0 <= n < s for n, s in izip(coordinate, self.shape))
        )
   
    @property
    def dimension_count(self):
        """Number of dimensions of the coordinate system."""
        return len(self.shape)
   
    def get_sample(self, count):
        """Gets a sample of `count` coordinates from the coordinate system."""
	idxprod = numpy.cumprod([1]+self.shape[:-1])
	return numpy.array(sample(xrange(len(self)),count),ndmin=2).T/idxprod%self.shape


class WalkerBulk(object):
    """Many random walkers in a lattice.
   
    :ivar coordinate_system: coordinate system in which the walkers operates.
    :ivar coordinates_and_id: current coordinates of the walkers.
      Last column is the id.
    """
    def __init__(self, coordinate_system, coordinates_and_id):
	assert coordinate_system.dimension_count == coordinates_and_id.shape[1]-1
        self.coordinate_system = coordinate_system
	self.coordinates_and_id = coordinates_and_id

	dim = self.coordinate_system.dimension_count
	self._steps = numpy.vstack([numpy.eye(dim,dtype=int),-numpy.eye(dim,dtype=int)])
	self._idxprod = numpy.cumprod([1]+self.coordinate_system.shape[:-1])


    def move(self, delta):
        """Moves the walkers by given relative coordinates `delta`.
       
        If the new coordinate would be outside the coordinate system of the
        walker the coordinates ”wrap around”, i.e. if the walker leaves a
        2D grid on the left it reappears on the right.
	"""
        self.coordinates_and_id[:,:-1] += delta
	self.coordinates_and_id[:,:-1] %= self.coordinate_system.shape
   
    def random_move(self):
        """Makes a step in a random dimension into a random direction.
       
        .. seealso:: :meth:`move`
        """
	pm = numpy.random.randint(len(self._steps),size=len(self.coordinates_and_id))
	self.move(self._steps[pm])

    def __len__(self):
        return len(self.coordinates_and_id)
   
    def __iter__(self):
        return iter(self.coordinates_and_id)
	
    def step(self):
        """Moves every Walker one step in a random direction.
       
        Walkers ending up in the same spot at the end of the step are
        annihilated.
        """
	self.random_move()
        coords = (self._idxprod*self.coordinates_and_id[:,:-1]).sum(1)
	counts = Counter(coords)
	keep = [counts[co] == 1 for co in coords]
        self.coordinates_and_id = self.coordinates_and_id.compress(keep,0)
   
    @classmethod
    def with_randomly_placed_walkers(cls, coordinate_system, count):
        """Creates a simulation with `count` randomly placed
        Walkers within the given `coordinate_system`.
        """
        return cls(coordinate_system,
	    numpy.hstack([coordinate_system.get_sample(count),
	        numpy.arange(count).reshape(count,1)]))


def main():
    """Runs a little test simulation."""
    simulation = WalkerBulk.with_randomly_placed_walkers(
        CoordinateSystem([100, 100]), 500
    )
    for i in xrange(100):
        simulation.step()
        print '{0} walkers left in step {1}.'.format(len(simulation), i)
        if len(simulation) <= 1:
            break

if __name__ == '__main__':
    main()

Kleiner Geschwindigkeitsvergleich bei mir, mit 100 Steps:
Reine Python-Version: 1.66s
Numpy-Version: 0.36s

Edit: Wie BlackJack richtig bemerkt hat, waren in der vorhergehenden Version die Startpositionen nicht
duplikatfrei. Jetzt korrigiert.

[BlackJack]

@Sirius3: Mit welchen Grössen sind denn diese Zeiten entstanden?

Wenn ich das richtig sehe stellst Du nicht sicher das alle Walker alleine auf einer Koordinate stehen‽ Je dichter man das Koordinatensystem besetzt um so weiter dürfte die Startsituation der beiden Implementierungen abweichen. Und ist ein Schritt eines Waklers bei Dir auch nur in einer Dimension?

[Sirius3]

@Blackjack: Du hast Recht. Die Startbedingungen waren nicht optimal. Ich hab's korrigiert.

Die möglichen Schritt-Richtungen sind in der Instanzvariable _steps gespeichert, ±1 in jede Richtung,
und werden in random_move zufällig ausgewählt. Damit braucht man pro Walker und Step nur jeweils
eine statt zwei Zufallszahlen, was die Geschwindigkeit um ~15% steigert.

Rein aus Interesse, wie sähenoch ein Geschwindigkeitsvergleich mit einer einfachen c-Version aus?

[BlackJack]

@Sirius3: Ich bekomme folgende Laufzeiten heraus:

Code: Alles auswählen

Python: 0.141s
+Numpy: 0.100s

Ist mein Rechner so schnell oder Deiner so langsam? Kann ich fast nicht glauben. Ich habe hier ja fast den Eindruck dass die beiden bei mir deutlich näher beieinander sind, liegt einfach am konstanten Aufwand der hier dominiert und dadurch der Vorteil von `numpy` nicht so sichtbar wird.

Bei einer einfachen C-Version bekommt man das Problem, dass man eine Hashtable selber implementieren oder aus einer Bibliothek verwenden muss, wenn man eine äquivalente Implementierung liefern möchte.

Ausserdem müsste man sich `random.sample()` mal anschauen um das auch effizient zu implementieren.

@Tyrax: Kannst Du mal ungefähre Grössenordnungen verraten mit denen man das testen kann um Deinem Szenario nahe zu kommen.

[EyDu]

Bei einer einfachen C-Version bekommt man das Problem, dass man eine Hashtable selber implementieren oder aus einer Bibliothek verwenden muss, wenn man eine äquivalente Implementierung liefern möchte.

Ausserdem müsste man sich `random.sample()` mal anschauen um das auch effizient zu implementieren.

In C++ sollte das schnell gemacht sein. Dictionaries benötigen als hash-Funktion für Indizes nur den <-Operator, welcher bei Koordinaten ja trivial zu implementieren ist. `random.sample()` kenne ich unter C++ nicht, auch nicht bei boost, wenn ich mich richtig erinnere, steht aber irgendwo bei Knuth etwas dazu.

Edit: je nach Datenlage könnte man random_shuffle verwenden. Und hier noch ein Stackoverflow-Thread, sogar mit dem Hinweis auf Knuth.

[BlackJack]

@EyDu: Das was den ``<``-Operator braucht ist keine Hashtable und damit hat das auch nichts mit einer Hashfunktion zu tun. Der Datentyp aus der STL ist als Baum implementiert.