Ein alter Klassiker mit neuer Perspektive

[ratwolf]

Der "Hopalong" Attractor (Martin's Map/Orbit) visualisiert mittels Density Heatmap und in 3D.

https://github.com/ratwolfzero/hopalong_python

https://github.com/ratwolfzero/hopalong ... 130_02.png
https://github.com/ratwolfzero/hopalong ... asic_1.png
https://github.com/ratwolfzero/hopalong ... c_3D_1.png

Code: Alles auswählen

import matplotlib.pyplot as plt
import numpy as np
from numba import njit
from math import copysign, sqrt, fabs
import time
import resource 


def validate_input(prompt, input_type=float, check_positive_non_zero=False, min_value=None):
    # Prompt for and return user input validated by type and positive/non-zero checks.
    while True:
        user_input = input(prompt)
        try:
            # Parse input as float first to handle scientific notation
            value = float(user_input)
            
            # Ensure the input is an integer, if expected
            if input_type == int:
                if not value.is_integer():
                    print('Invalid input. Please enter an integer.')
                    continue
                value = int(value)

            # Check if input is a positive non-zero value
            if check_positive_non_zero and value <= 0:
                print('Invalid input. The value must be a positive non-zero number.')
                continue

            # Then, check minimum value
            if min_value is not None and value < min_value:
                print(f'Invalid input. The value should be at least {min_value}.')
                continue

            return value
        except ValueError:
            print(f'Invalid input. Please enter a valid {input_type.__name__} value.')


def get_attractor_parameters():
    a = validate_input('Enter a float value for "a": ', float)
    b = validate_input('Enter a float value for "b": ', float)
    while True:
        c = validate_input('Enter a float value for "c": ', float)
        if (a == 0 and b == 0 and c == 0) or (a == 0 and c == 0):
            print('Invalid combination of parameters. The following combinations are not allowed:\n'
                  '- a = 0, b = 0, c = 0\n'
                  '- a = 0, b = any value, c = 0\n'
                  'Please enter different values.')
        else:
            break
    n = validate_input('Enter a positive integer value > 1000 for "n": ', int, check_positive_non_zero=True, min_value=1000)
    return {'a': a, 'b': b, 'c': c, 'n': n}


@njit #njit is an alias for nopython=True
def compute_trajectory_extents(a, b, c, n):
    # Dynamically compute and track the minimum and maximum extents of the trajectory over 'n' iterations.
    x = np.float64(0.0)
    y = np.float64(0.0)

    min_x = np.inf  # ensure that the initial minimum is determined correctly
    max_x = -np.inf # ensure that the initial maximum is determined correctly
    min_y = np.inf
    max_y = -np.inf

    for _ in range(n):
    # selective min/max update using direct comparisons avoiding min/max function
        if x < min_x:
            min_x = x
        if x > max_x:
            max_x = x
        if y < min_y:
            min_y = y
        if y > max_y:
            max_y = y
        # signum function respecting the behavior of floating point numbers according to IEEE 754 (signed zero)
        xx = y - copysign(1.0, x) * sqrt(fabs(b * x - c))
        yy = a-x
        x = xx
        y = yy
        
    return min_x, max_x, min_y, max_y

# Dummy call to ensure the function is pre-compiled by the JIT compiler before it's called by the interpreter.
_ = compute_trajectory_extents(1.0, 1.0, 1.0, 2)


@njit
def compute_trajectory_and_image(a, b, c, n, extents, image_size):
    # Compute the trajectory and populate the image with trajectory points
    image = np.zeros(image_size, dtype=np.uint64)
    
    # pre-compute image scale factors
    min_x, max_x, min_y, max_y = extents
    scale_x = (image_size[1] - 1) / (max_x - min_x) # column
    scale_y = (image_size[0] - 1) / (max_y - min_y) # row
    
    x = np.float64(0.0)
    y = np.float64(0.0)
    
    for _ in range(n):
        # map trajectory points to image pixel coordinates
        px = np.uint64((x - min_x) * scale_x)
        py = np.uint64((y - min_y) * scale_y)
        # populate the image arrayy "on the fly" with each computed point
        image[py, px] += 1  # respecting row/column convention, update # of hits

        # Update the trajectory "on the fly"
        xx = y - copysign(1.0, x) * sqrt(fabs(b * x - c))
        yy = a-x
        x = xx
        y = yy
        
    return image

# Dummy call to ensure the function is pre-compiled by the JIT compiler before it's called by the interpreter.
_ = compute_trajectory_and_image(1.0, 1.0, 1.0, 2, (-1, 0, 0, 1), (2, 2))


def render_trajectory_image(image, extents, params, color_map, mode='2D'):
    if mode == '2D':
        # Render in 2D using imshow
        fig = plt.figure(figsize=(8, 8), facecolor='gainsboro')
        ax = fig.add_subplot(1, 1, 1)
        
        img = ax.imshow(image, origin='lower', cmap=color_map, extent=extents, interpolation='none')
        ax.set_title(f'Hopalong Attractor@ratwolf@2024\nParams: a={params["a"]}, b={params["b"]}, c={params["c"]}, n={params["n"]:_}')
        ax.set_xlabel('X (Cartesian)')
        ax.set_ylabel('Y (Cartesian)')

        cbar = fig.colorbar(img, ax=ax, location='bottom')
        cbar.set_label('Pixel Density. (Scale = 0 - max)')
        max_hit_count = np.max(image)
        cbar.set_ticks(np.linspace(0, max_hit_count, num=8))
        plt.tight_layout()
        plt.show()
    elif mode == '3D':
        # Render in 3D using contourf3D
        x = np.linspace(extents[0], extents[1], image.shape[1])
        y = np.linspace(extents[2], extents[3], image.shape[0])
        x, y = np.meshgrid(x, y)
        z = image / np.max(image) if np.max(image) > 0 else image
        
        fig = plt.figure(figsize=(8, 8))
        ax = fig.add_subplot(111, projection='3d')
        ax.contourf3D(x, y, z, levels=100, cmap=color_map)
        
        ax.set_title(f'Hopalong Attractor - 3D Density (Z) Plot\nParams: a={params["a"]}, b={params["b"]}, c={params["c"]}, n={params["n"]:_}')
        ax.set_xlabel('X')
        ax.set_ylabel('Y')
        ax.set_zlabel('Z')
        ax.view_init(elev=75, azim=-95)
        plt.show()
    else:
        print("Invalid mode. Please choose '2D' or '3D'.")

def main(image_size=(1000, 1000), color_map='hot', mode='2D'):
    # Main execution process
    try:
        params = get_attractor_parameters()
        
        # Start the time measurement
        start_time = time.process_time()

        extents = compute_trajectory_extents(params['a'], params['b'], params['c'], params['n'])
        image = compute_trajectory_and_image(params['a'], params['b'], params['c'], params['n'], extents, image_size)
        render_trajectory_image(image, extents, params, color_map, mode=mode)

        # End the time measurement
        end_time = time.process_time()

        # Calculate the CPU user and system time
        cpu_sys_time_used = end_time - start_time

        # Calculate the memory resources used
        memMb = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024.0 / 1024.0
        
        print(f'CPU User&System time: {cpu_sys_time_used:.2f} seconds')
        print(f'Memory (RAM): {memMb:.2f} MByte used')
        
    except Exception as e:
        print(f'An error occurred: {e}')


# Main execution
if __name__ == '__main__':
    mode = input("Choose visualization mode (2D/3D): ").strip().upper()
    if mode not in ['2D', '3D']:
        print("Invalid choice. Defaulting to 2D.")
        mode = '2D'
    main(mode=mode)

[Sirius3]

Explizite Typangaben sind unnötig, wenn es sich eh um die defaults handelt. Ich bin eher ein Freund von built-in-Funktionen, wenn sie das selbe tun, wie das, was aus math extra importiert werden muß. Mit tuple-Assignment spart man sich Hilfsvariablen. Beim Umwandeln von float nach int sollte man runden um nummerische Fehler zu vermeiden. Bei `extends` ist es so wichtig, dass sie zu den Parametern passen, dass man diese innerhalb der Funktion ausrechnen sollte.

Code: Alles auswählen

@njit
def compute_trajectory_extents(a, b, c, n):
    """ Dynamically compute and track the minimum and maximum extents of the trajectory over 'n' iterations. """
    min_x = max_x = x = 0.0
    min_y = max_y = y = 0.0

    for _ in range(n):
        if x < min_x:
            min_x = x
        if x > max_x:
            max_x = x
        if y < min_y:
            min_y = y
        if y > max_y:
            max_y = y
        d = copysign(abs(b * x - c)**0.5, x)
        (x, y) = (y - d, a - x)
        
    return min_x, max_x, min_y, max_y

@njit
def compute_trajectory_and_image(a, b, c, n, image_size):
    """ Compute the trajectory and populate the image with trajectory points. """
    image = np.zeros(image_size, dtype=np.uint64)
    
    # pre-compute image scale factors
    min_x, max_x, min_y, max_y = compute_trajectory_extents(a, b, c, n)
    scale_x = (image_size[1] - 1) / (max_x - min_x) # column
    scale_y = (image_size[0] - 1) / (max_y - min_y) # row
    
    x = y = 0.0
    for _ in range(n):
        # map trajectory points to image pixel coordinates
        px = round((x - min_x) * scale_x)
        py = round((y - min_y) * scale_y)
        # populate the image array "on the fly" with each computed point
        image[py, px] += 1  # respecting row/column convention, update # of hits

        d = copysign(abs(b * x - c)**0.5, x)
        (x, y) = (y - d, a - x)
        
    return image, min_x, max_x, min_y, max_y

Mit `continue` sollte man möglichst sparsam sein, vor allem wenn else den selben Zweck erfüllt. Ein try-Block sollte möglichst kurz sein.

Code: Alles auswählen

def validate_input(prompt, input_type=float, check_positive_non_zero=False, min_value=None):
    """ Prompt for and return user input validated by type and positive/non-zero checks. """
    while True:
        user_input = input(prompt)
        try:
            value = float(user_input)
        except ValueError:
            print(f'Invalid input. Please enter a valid {input_type.__name__} value.')
        else:
            converted_value = input_type(value)
            if converted_value != value:
                print('Invalid input. Please enter an integer.')
            elif check_positive_non_zero and converted_value <= 0:
                print('Invalid input. The value must be a positive non-zero number.')
            elif min_value is not None and converted_value < min_value:
                print(f'Invalid input. The value should be at least {min_value}.')
            else:
                return converted_value

[snafu]

Man könnte die Logik zur Umwandlung auch in eine eigene Klasse stecken:

Code: Alles auswählen

#!/usr/bin/env python3
class Converter:
    def __init__(self, convert, interval=None):
        self.convert = convert
        self.interval = interval

    def __call__(self, value):
        converted = self.convert(value)
        if self.interval is None:
            return converted
        low, high = self.interval
        if low <= converted <= high:
            return converted
        raise ValueError(converted)


def get_input(prompt, convert, error_message):
    while True:
        user_input = input(prompt)
        try:
            return convert(user_input)
        except ValueError:
            print(error_message.format(user_input))


def main():
    number = get_input(
        "Please choose a number (1-100) ",
        Converter(int, (1, 100)),
        "{} is not a valid number"
    )
    print(f"Your choice was {number}")


if __name__ == "__main__":
    main()

Und anstelle des interval-Arguments wäre eine Validator-Klasse denkbar, um die Flexibilität zu erhöhen. Und schon hat man sich ein kleines Framework gebaut.

[Sirius3]

@snafu: sobald man wirklich verschiedene Converter braucht, könnte man sich auch überlegen, dafür extra mehrere Klassen zu definieren. Selbstzweck sollten aber Klassen nicht sein. Aber schön zu sehen, was alles geht.

[ratwolf]

Danke für die Rückmeldungen soweit.

Die explizite Deklaration von Datentypen ist zumindest Best Practice und in anderen Programmiersprachen notwendig.

Die Integer-Funktion heißt nicht so, weil sie rundet, sondern weil sie den Dezimalteil abschneidet und nur den Integer-Teil zurückgibt.
Will man runden, rundet man; will man den Integer-Teil, verwendet man die Integer-Funktion.

Streng genommen wird beim Runden aus einer Dezimalzahl kein Integer-Datentyp, sondern eine Dezimalzahl, die wie eine Integer-Zahl aussieht.

Das Modul „math“ wurde importiert, weil „copysign“ mit seiner spezifischen Definition zur Anwendung kommen soll, statt der üblichen Vorzeichenfunktion, die Python nicht „built-in“ hat.
„fabs“ bot sich hier an, da es immer einen Fließkommawert zurückgibt und jede weitere interne „int-float“-Konvertierung überflüssig macht.

Ob „sqrt“ oder „^0.5“ schneller ist, sei dahingestellt, aber meiner Erfahrung nach ist das nicht der Fall.

Hilfsvariablen zu vermeiden, ja, das ist eloquent, aber ich bevorzuge immer sauberen, geradlinigen Code. Bei der Verwendung von JIT verbessert dies die Leistung ohnehin nicht wirklich. Das ist vielleicht eher bei Standard Python der Fall.

Anspruchsvolle Routinen zur Validierung der Eingabedaten sind „nett zu haben“ …

Aber ich muss neidlos zugeben, dass das Runden auf die „nächste Ganzzahl“ in Bezug auf die Dichteapproximation im kontinuierlichen Raum tatsächlich vorteilhafter ist.

[__blackjack__]

@ratwolf: Ich würde das ja nicht als Deklaration bezeichnen wenn man ein Objekt von einem Datentyp erstellt. So etwas wie ``np.float64(0.0)`` ist in Python jedenfalls keine „best practice“. Deklarationen gibt es Python nur recht wenige (zum Beispiel ``global`` und ``nonlocal`` — keine von Datentypen), und es gibt Typannotationen, die hier aber nicht verwendet werden, und auf den Programmablauf auch keinen Einfluss haben.

[ratwolf]

Naja, eine dedizierte main()-Funktion als Einstiegspunkt gibt es in Python eigentlich auch nicht, wie in C.
Aber auch hier kann man sich vielleicht auf „good practice“ oder „common practice“ einigen – es ist nicht verkehrt, das Programm nach dem „One Responsibility Principle“ in sinnvolle Funktionen (Defs) aufzuteilen und dann als Einstiegspunkt eine main()-Funktion zu definieren, die die Programmausführung koordiniert. Natürlich kann ich in Python auch alles in einem Brei runter schreiben..

[Sirius3]

Was hat eine sinnvolle Aufteilung in Funktion mit überflüssigen Typumwandlungen zu tun?

[ratwolf]

frohes Fest

[__blackjack__]

@ratwolf: Boah, ey, erschreck mich doch nicht so! Weihnachten kommt immer so plötzlich und unerwartet, aber es sind ja noch 11 Tage bis ich die Geschenke besorgt haben muss.