distributed_representation_of_APF/fig_kappa_effect_on_RF.py at main · songfei-psy/distributed_representation_of_APF · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
import numpy as np
import os
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
from scipy.special import iv  # Import the modified Bessel function of the first kind

def encode_angle(phi, theta, kappa):
    return np.exp(kappa * np.cos(phi - theta)) / np.exp(kappa)

# NRF
def decode_angle(N, kappa, phig, phio):
    theta = np.linspace(0, 2 * np.pi, N, endpoint=False)
    o = encode_angle(phio, theta, kappa)
    g = encode_angle(phig, theta, 1.5)
    s = g + o
    weighted_sum = np.sum(s * np.exp(1j * theta))
    decoded_phi = np.angle(weighted_sum)
    return decoded_phi


if __name__ == '__main__':
    # folder
    project_folder = os.path.dirname(os.path.abspath(__file__))
    data_folder = os.path.join(project_folder, 'data', 'method')
    if not os.path.exists(data_folder):
        os.makedirs(data_folder)

    # -------------------------- Define extended ranges of kappa and scale ratio ----------------------------
    # Define the range of kappa values
    kappa = np.linspace(0.001, 50, 500)
    # Calculate the exact ratio I1(kappa)/exp(kappa)
    scaling_factor = iv(1, kappa) / np.exp(kappa)
    # Approximation for small kappa: kappa / 2
    approx_small_kappa = kappa / 2
    # Asymptotic approximation for large kappa: 1 / sqrt(2 * pi * kappa)
    approx_large_kappa = 1 / np.sqrt(2 * np.pi * kappa)

    # Find the maximum value of the exact ratio and its corresponding kappa
    max_idx = np.argmax(scaling_factor)
    max_kappa = kappa[max_idx]
    max_value = scaling_factor[max_idx]
    # Plot the results
    plt.figure(figsize=(10, 8))
    fontsize = 16
    font = FontProperties()
    font.set_size(fontsize-2)
    font.set_weight('bold')
    plt.plot(kappa, scaling_factor, label=r"Exact $\frac{I_1(\kappa)}{e^\kappa}$", linewidth=2,
             color=plt.cm.viridis(2 / 6))
    plt.plot(kappa, approx_small_kappa, label=r"Approximation $\frac{\kappa}{2}$", linestyle='--', linewidth=2,
             color='orange')
    plt.plot(kappa, approx_large_kappa, label=r"Asymptotic $\frac{1}{\sqrt{2\pi\kappa}}$", linestyle='--', linewidth=2,
             color='green')

    # Highlight the maximum point
    plt.scatter([max_kappa], [max_value], color='red', label=f"Maximum at $\kappa = {max_kappa:.2f}$", zorder=5)

    # Add axis labels and title
    plt.xlabel(r'$\kappa$', fontsize=fontsize, fontweight='bold')
    plt.ylabel(r'Amplitude Ratio', fontsize=fontsize, fontweight='bold')
    plt.ylim(0, 0.3)
    plt.title(r'Amplitude Ratio for Different Values of $\kappa$', fontsize=fontsize, fontweight='bold')
    # Add grid and legend
    plt.grid(alpha=0.2)
    plt.xticks(fontsize=fontsize-1, fontweight='bold')
    plt.yticks(fontsize=fontsize-1, fontweight='bold')
    plt.legend(loc='upper right', prop=font)
    # save
    fts = {'.png', '.eps'}
    for ft in fts:
        filename = 'kappa_influence_ratio' + ft
        plt.savefig(os.path.join(data_folder, filename), dpi=600)

    # -------------------------------- Define extended ranges of kappa and N --------------------------------
    N = 100
    kappa_group = [[1.5, 5, 10, 20, 50, 250],[1.5, 0.54, 0.34, 0.22, 0.13, 0.056]]
    # kappa_values = [1.5, 0.5, 0.12, 0.087, 0.056, 0.025]

    # Plot decoded angle for different kappas
    phig = 0  # goal angle
    phio = np.linspace(0, np.pi, 100, endpoint=False) # obstacle angle
    phio = phio[1:] # remove the first element to avoid division by zero
    sub_names = ['a', 'b']
    for kappa_values, sub_name in zip(kappa_group, sub_names):
        fig = plt.figure(figsize=(10, 8))
        plt.fill_between(np.rad2deg(phio), np.rad2deg((phig)*np.ones((len(phio),))), np.rad2deg((phio + phig) / 2),
                         color='blue',
                         alpha=0.1,
                         label="Region between \n $\phi_g$ and $\\frac{\phi_g+\phi_o}{2}$")
        plt.plot(np.rad2deg(phio), np.rad2deg((phig)*np.ones((len(phio),))), 'r--', label="$\phi_g$")
        for k, kappa in enumerate(kappa_values):
            decoded = np.zeros((len(phio),))
            color_value = plt.cm.viridis(k / len(kappa_values))
            for i, phi in enumerate(phio):
                arg = decode_angle(N, kappa, phig, phi)
                decoded[i] = np.rad2deg(arg)
            plt.plot(np.rad2deg(phio), decoded, marker='o', color=color_value, label=f'$\kappa$ = {kappa}')

        plt.xlabel(r"$\phi_o$ (degrees)",fontsize=fontsize, fontweight='bold')
        plt.ylabel('Decoded Angle (degrees)', fontsize=fontsize, fontweight='bold')
        plt.xticks(fontsize=fontsize-1, fontweight='bold')
        plt.yticks(fontsize=fontsize-1, fontweight='bold')
        plt.legend(loc='upper left', prop=font)
        # save
        fts = {'.png', '.eps'}
        for ft in fts:
            filename = 'RF_decoded_angle_' + sub_name + ft
            plt.savefig(os.path.join(data_folder, filename), dpi=600)

    plt.show()