PyPDEVS-and-JADEVS-DEVStone-benchmarks/devstone/plotNumberOfDCP.py

import argparse
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from model_properties import model_properties

# ==========================================
#               Model Formulas
# ==========================================

def model_dLI(width, depth, n):
    """
    dLI Model Formula:
    DCP = ((n + 1)(d - 1) + 1) * floor((w - 1) / n) + remaining
    remaining = (n_2 + 1)(d - 1) + 1  if n_2 > 0 else 0
    where n_2 = (w - 1) % n
    """
    w = np.atleast_1d(width)
    d = np.atleast_1d(depth)
    n_val = np.atleast_1d(n)

    # Broadcast scalar n if w/d are arrays
    if len(n_val) == 1 and len(w) > 1:
        n_val = np.full_like(w, n_val[0])

    # Ensure n >= 1
    n_val = np.maximum(n_val, 1)

    # Main Term
    floor_term = np.floor((w - 1) / n_val)
    term_main = ((n_val + 1) * (d - 1) + 1) * floor_term

    # Remaining Term
    n_2 = (w - 1) % n_val
    term_rem = (n_2 + 1) * (d - 1) + 1
    term_rem = np.where(n_2 > 0, term_rem, 0)

    return term_main + term_rem

def model_dHI(width, depth, n):
    """
    dHI Model Formula:
    DCP = c1 + ci * floor((w - n - 1) / n) + remaining
    c1 = (2n + 1)(d - 1) + 1
    ci = (2n + 2)(d - 1) + 1
    remaining = (2n_2 + 1)(d - 1) + 1  if n_2 > 0 else 0
    where n_2 = (w - n - 1) % n
    """
    w = np.atleast_1d(width)
    d = np.atleast_1d(depth)
    n_val = np.atleast_1d(n)

    if len(n_val) == 1 and len(w) > 1:
        n_val = np.full_like(w, n_val[0])

    # Safety: clamp n to width-1 to preserve physical meaning (chunk size cannot exceed width)
    # This prevents negative floor values when w < n+1
    n_eff = np.where(w > 1, np.minimum(n_val, w - 1), 1)

    # First term (c1)
    term_c1 = (2 * n_eff + 1) * (d - 1) + 1

    # Intermediate term (ci)
    val_for_floor = w - n_eff - 1
    floor_term = np.floor(val_for_floor / n_eff)
    term_ci = ((2 * n_eff + 2) * (d - 1) + 1) * floor_term

    # Remaining term
    n_2 = val_for_floor % n_eff
    term_rem = (2 * n_2 + 1) * (d - 1) + 1
    term_rem = np.where(n_2 > 0, term_rem, 0)

    return term_c1 + term_ci + term_rem


MODELS = {
    'dLI': model_dLI,
    'dHI': model_dHI,
}


def get_model_func(model_name):
    return MODELS.get(model_name)


# ==========================================
#               Execution Modes
# ==========================================

def run_vary_n(args):
    """Mode 1: Fix W/D, Vary n"""
    print(f"Running 'Vary N' | Model: {args.model}, Width: {args.width}, Depth: {args.depth}")

    func = get_model_func(args.model)
    if not func:
        print(f"Error: Model '{args.model}' not found.")
        return

    n_vals = np.arange(1, args.limit + 1)
    Y = func(args.width, args.depth, n_vals)

    plt.figure(figsize=(10, 6))
    plt.plot(n_vals, Y, marker='o', label=f'{args.model}: W={args.width}, D={args.depth}')
    plt.yscale('log')
    plt.xlabel('n (Models per cycle)')
    plt.ylabel('Number of DCP calls')
    plt.title(f'#DCP Calls vs n for {args.model}')
    plt.legend()
    plt.grid(True, which="both", ls="--")
    plt.tight_layout()
    plt.show()


def run_vary_dim(args):
    """Mode 2: Fix n, Vary Width/Depth (The 'Old Code' adaptation)"""
    print(f"Running 'Vary Dimensions' | Model: {args.model}, n: {args.n}")

    func = get_model_func(args.model)
    if not func:
        print(f"Error: Model '{args.model}' not found.")
        return

    # Generate grid of W and D
    w_vals = np.arange(2, args.limit + 1)
    d_vals = np.arange(2, args.limit + 1)
    width, depth = np.meshgrid(w_vals, d_vals)

    # Calculate X axis (Total Atomic Models)
    X = (width - 1) * (depth - 1) + 1

    # Calculate Y axis (DCP Calls)
    Y = func(width, depth, args.n)

    plt.figure(figsize=(10, 6))
    plt.scatter(X.flatten(), Y.flatten(), s=5, alpha=0.5, label=f'{args.model}, n={args.n}')
    plt.yscale('log')
    # plt.xscale('log') # Optional: log scale for X if desired
    plt.xlabel('Number of atomic models')
    plt.ylabel('Number of DCP calls')
    plt.title(f'#DCP Calls vs Atomic Models for {args.model}')
    plt.legend()
    plt.grid(True, which="both", ls="--")
    plt.tight_layout()
    plt.show()

def clean_data(df):
    df["total_time"] = (
        df["simulation_time"]
        + df["creation_time"]
        + df["setup_time"]
    )

    possible_group_cols = ["depth", "width", "number", "total_time"]
    grouping_cols = [c for c in possible_group_cols if c in df.columns]

    def clean(group):
        Q1 = group["total_time"].quantile(0.25)
        Q3 = group["total_time"].quantile(0.75)
        IQR = Q3 - Q1
        lower = Q1 - 1.5 * IQR
        upper = Q3 + 1.5 * IQR
        return group[(group["total_time"] >= lower) & (group["total_time"] <= upper)]

    cleaned = (
        df.groupby(grouping_cols, group_keys=False)[grouping_cols]
          .apply(clean)
    )

    return cleaned

def run_csv_mode(args):
    """Mode 3: Read CSV, calculate theoretical DCP, plot vs Runtime"""
    print(f"Running 'CSV Mode' | File: {args.file}")

    try:
        df = pd.read_csv(args.file)
    except FileNotFoundError:
        print(f"Error: File {args.file} not found.")
        return
    df = clean_data(df)

    # Create 'dcp_calls' column
    df['dcp_calls'] = 0.0

    # Determine model per row or globally
    if 'model_type' not in df.columns:
        df['model_type'] = args.model

    # Check for unknown models in CSV
    unique_models = df['model_type'].unique()
    if len(unique_models) > 1:
        print(f"Error: More than one model type found. Found: {unique_models}")
        return
    model = unique_models[0]
    if model not in MODELS:
        print(f"Error: CSV contains unknown model '{model}'.")
        return

    df['atomic_models'] = (df['width'] - 1) * (df['depth'] - 1) + 1

    # Calculate DCP calls by group (efficient for different models)
    for m_name, group in df.groupby('model_type'):
        func = get_model_func(m_name)
        # Fallback to default if model not found in dictionary
        if not func:
            func = get_model_func(args.model)

        # Use 'n' from CSV if available, otherwise CLI default
        n_vals = group['number'].values if 'number' in group.columns else args.n

        df.loc[group.index, 'dcp_calls'] = func(group['width'].values, group['depth'].values, n_vals)

    plt.figure(figsize=(10, 6))
    scatter = plt.scatter(df['dcp_calls'], df['total_time'],
                          c=df['atomic_models'], cmap='viridis', alpha=0.7, edgecolors='k', s=50)

    plt.colorbar(scatter, label='# Atomic Models')
    plt.xlabel('Theoretical DCP Calls')
    plt.ylabel('Total Simulation Time (s)')
    plt.title(f'Run Time vs DCP Calls for {model}')
    plt.grid(True, linestyle='--', alpha=0.6)
    plt.tight_layout()
    plt.show()


# ==========================================
#               CLI
# ==========================================

def main():
    parser = argparse.ArgumentParser(description="DCP Calculation Tool")
    subparsers = parser.add_subparsers(dest='mode', required=True, help='Operation Mode')

    # Common arguments
    parent_parser = argparse.ArgumentParser(add_help=False)
    parent_parser.add_argument('--model', type=str, default='dHI', choices=MODELS.keys(),
                               help='Model formula to use (default: dHI)')

    # 1. Vary N
    parser_n = subparsers.add_parser('created', parents=[parent_parser], help='Fix Width and Depth, Vary number of models created per cycle')
    parser_n.add_argument('--width', '-w', type=int, default=100)
    parser_n.add_argument('--depth', '-d', type=int, default=100)
    parser_n.add_argument('--limit', '-l', type=int, default=100, help='Max n')

    # 2. Vary Dimensions
    parser_dim = subparsers.add_parser('dimensions', parents=[parent_parser], help='Fix n, Vary W/D')
    parser_dim.add_argument('--n', type=int, default=2)
    parser_dim.add_argument('--limit', '-l', type=int, default=100, help='Max width/depth')

    # 3. CSV
    parser_csv = subparsers.add_parser('time', parents=[parent_parser], help='Read CSV and plot DCP calls against simulation time')
    parser_csv.add_argument('--file', '-f', type=str)
    parser_csv.add_argument('--n', type=int, default=1, help='Default n if not in CSV')

    args = parser.parse_args()

    if args.mode == 'created':
        run_vary_n(args)
    elif args.mode == 'dimensions':
        run_vary_dim(args)
    elif args.mode == 'time':
        run_csv_mode(args)

if __name__ == "__main__":
    main()