ap1approx6.py

#
# Approximate AP1 cusp with 6 cos/sin coefficients
#
import numpy as np
from scipy.optimize import minimize
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# 100 nits AP1 data (ACES2 DRT model space M correlate)

# All x data is in radians
# All y data is M

x_corners = np.array([
    0.43200641091143477,
    1.8514061932331212,
    2.5711039444948653,
    3.3551054834647194,
    4.690868255638757,
    5.962357982293879
])

y_corners = np.array([
    118.92053782235328,
    76.75733362819672,
    110.89122468941676,
    65.37220909498073,
    58.17770736440632,
    100.49203357671745
])

x_ap1 = np.array([
    0.010,
    0.022,
    0.034,
    0.046,
    0.059,
    0.072,
    0.086,
    0.101,
    0.117,
    0.133,
    0.151,
    0.170,
    0.191,
    0.213,
    0.239,
    0.269,
    0.305,
    0.353,
    0.432,
    0.455,
    0.479,
    0.503,
    0.528,
    0.553,
    0.579,
    0.605,
    0.632,
    0.659,
    0.686,
    0.714,
    0.742,
    0.771,
    0.800,
    0.829,
    0.858,
    0.888,
    0.917,
    0.947,
    0.977,
    1.006,
    1.036,
    1.066,
    1.095,
    1.124,
    1.153,
    1.182,
    1.210,
    1.238,
    1.266,
    1.293,
    1.319,
    1.346,
    1.371,
    1.396,
    1.421,
    1.445,
    1.468,
    1.491,
    1.514,
    1.535,
    1.557,
    1.577,
    1.597,
    1.617,
    1.636,
    1.654,
    1.672,
    1.690,
    1.707,
    1.723,
    1.739,
    1.755,
    1.770,
    1.784,
    1.799,
    1.812,
    1.826,
    1.839,
    1.851,
    1.864,
    1.876,
    1.889,
    1.902,
    1.914,
    1.927,
    1.940,
    1.953,
    1.965,
    1.978,
    1.991,
    2.004,
    2.017,
    2.030,
    2.043,
    2.056,
    2.068,
    2.081,
    2.094,
    2.107,
    2.120,
    2.133,
    2.145,
    2.158,
    2.171,
    2.184,
    2.196,
    2.209,
    2.221,
    2.234,
    2.246,
    2.259,
    2.271,
    2.283,
    2.295,
    2.307,
    2.319,
    2.331,
    2.343,
    2.355,
    2.367,
    2.378,
    2.390,
    2.401,
    2.413,
    2.424,
    2.435,
    2.446,
    2.457,
    2.468,
    2.479,
    2.489,
    2.500,
    2.510,
    2.521,
    2.531,
    2.541,
    2.551,
    2.561,
    2.571,
    2.597,
    2.620,
    2.641,
    2.661,
    2.679,
    2.697,
    2.714,
    2.731,
    2.747,
    2.763,
    2.778,
    2.793,
    2.808,
    2.822,
    2.837,
    2.851,
    2.865,
    2.878,
    2.892,
    2.905,
    2.918,
    2.932,
    2.945,
    2.957,
    2.970,
    2.983,
    2.995,
    3.008,
    3.020,
    3.032,
    3.044,
    3.056,
    3.068,
    3.080,
    3.092,
    3.103,
    3.115,
    3.126,
    3.137,
    3.149,
    3.160,
    3.171,
    3.182,
    3.193,
    3.203,
    3.214,
    3.225,
    3.235,
    3.246,
    3.256,
    3.266,
    3.277,
    3.287,
    3.297,
    3.307,
    3.317,
    3.326,
    3.336,
    3.346,
    3.355,
    3.365,
    3.375,
    3.385,
    3.395,
    3.406,
    3.416,
    3.428,
    3.439,
    3.451,
    3.463,
    3.475,
    3.488,
    3.501,
    3.515,
    3.529,
    3.543,
    3.557,
    3.572,
    3.588,
    3.604,
    3.620,
    3.637,
    3.654,
    3.672,
    3.690,
    3.709,
    3.729,
    3.749,
    3.769,
    3.790,
    3.812,
    3.834,
    3.857,
    3.880,
    3.904,
    3.929,
    3.954,
    3.980,
    4.007,
    4.034,
    4.062,
    4.090,
    4.120,
    4.149,
    4.180,
    4.210,
    4.242,
    4.274,
    4.306,
    4.339,
    4.373,
    4.406,
    4.441,
    4.475,
    4.510,
    4.546,
    4.581,
    4.618,
    4.654,
    4.691,
    4.744,
    4.796,
    4.847,
    4.897,
    4.945,
    4.991,
    5.036,
    5.079,
    5.121,
    5.160,
    5.198,
    5.235,
    5.269,
    5.303,
    5.334,
    5.364,
    5.393,
    5.421,
    5.447,
    5.472,
    5.496,
    5.519,
    5.540,
    5.561,
    5.581,
    5.600,
    5.619,
    5.636,
    5.653,
    5.669,
    5.684,
    5.699,
    5.714,
    5.727,
    5.741,
    5.753,
    5.766,
    5.778,
    5.789,
    5.800,
    5.811,
    5.821,
    5.831,
    5.841,
    5.850,
    5.859,
    5.868,
    5.877,
    5.885,
    5.893,
    5.901,
    5.908,
    5.916,
    5.923,
    5.930,
    5.937,
    5.943,
    5.950,
    5.956,
    5.962,
    5.968,
    5.975,
    5.981,
    5.987,
    5.993,
    6.000,
    6.006,
    6.013,
    6.019,
    6.026,
    6.033,
    6.039,
    6.046,
    6.053,
    6.060,
    6.067,
    6.074,
    6.081,
    6.089,
    6.096,
    6.104,
    6.111,
    6.119,
    6.127,
    6.135,
    6.143,
    6.151,
    6.159,
    6.167,
    6.176,
    6.185,
    6.194,
    6.203,
    6.212,
    6.221,
    6.231,
    6.241,
    6.251,
    6.261,
    6.271,
    6.282
])

y_ap1 = np.array([
    110.201,
    110.560,
    110.927,
    111.302,
    111.684,
    112.075,
    112.473,
    112.880,
    113.295,
    113.720,
    114.155,
    114.601,
    115.060,
    115.536,
    116.033,
    116.562,
    117.144,
    117.839,
    118.921,
    114.550,
    110.511,
    106.770,
    103.298,
    100.072,
    97.072,
    94.280,
    91.682,
    89.264,
    87.015,
    84.925,
    82.985,
    81.186,
    79.522,
    77.986,
    76.571,
    75.271,
    74.082,
    72.998,
    72.014,
    71.126,
    70.329,
    69.620,
    68.994,
    68.447,
    67.975,
    67.575,
    67.244,
    66.977,
    66.771,
    66.624,
    66.532,
    66.491,
    66.500,
    66.554,
    66.652,
    66.791,
    66.968,
    67.181,
    67.427,
    67.705,
    68.013,
    68.347,
    68.707,
    69.091,
    69.497,
    69.923,
    70.369,
    70.832,
    71.311,
    71.805,
    72.313,
    72.834,
    73.366,
    73.909,
    74.463,
    75.025,
    75.595,
    76.173,
    76.757,
    76.835,
    76.925,
    77.029,
    77.146,
    77.277,
    77.422,
    77.580,
    77.753,
    77.941,
    78.143,
    78.360,
    78.592,
    78.839,
    79.101,
    79.380,
    79.674,
    79.984,
    80.310,
    80.652,
    81.011,
    81.387,
    81.780,
    82.189,
    82.616,
    83.060,
    83.522,
    84.002,
    84.499,
    85.014,
    85.548,
    86.100,
    86.670,
    87.259,
    87.867,
    88.494,
    89.140,
    89.805,
    90.490,
    91.195,
    91.919,
    92.664,
    93.429,
    94.214,
    95.019,
    95.846,
    96.693,
    97.562,
    98.452,
    99.364,
    100.297,
    101.253,
    102.231,
    103.231,
    104.254,
    105.301,
    106.371,
    107.465,
    108.582,
    109.724,
    110.891,
    106.969,
    103.787,
    101.084,
    98.723,
    96.622,
    94.727,
    93.001,
    91.415,
    89.949,
    88.586,
    87.315,
    86.125,
    85.007,
    83.954,
    82.960,
    82.020,
    81.130,
    80.286,
    79.484,
    78.721,
    77.995,
    77.304,
    76.645,
    76.016,
    75.416,
    74.843,
    74.296,
    73.773,
    73.274,
    72.796,
    72.340,
    71.904,
    71.486,
    71.088,
    70.707,
    70.343,
    69.995,
    69.663,
    69.346,
    69.043,
    68.754,
    68.479,
    68.217,
    67.967,
    67.729,
    67.503,
    67.288,
    67.084,
    66.891,
    66.708,
    66.535,
    66.371,
    66.216,
    66.071,
    65.934,
    65.806,
    65.686,
    65.574,
    65.469,
    65.372,
    64.848,
    64.324,
    63.802,
    63.281,
    62.762,
    62.245,
    61.730,
    61.218,
    60.709,
    60.202,
    59.699,
    59.200,
    58.705,
    58.215,
    57.729,
    57.249,
    56.775,
    56.307,
    55.847,
    55.394,
    54.949,
    54.513,
    54.086,
    53.671,
    53.266,
    52.873,
    52.493,
    52.128,
    51.777,
    51.442,
    51.125,
    50.826,
    50.547,
    50.289,
    50.054,
    49.843,
    49.658,
    49.501,
    49.373,
    49.277,
    49.214,
    49.187,
    49.197,
    49.248,
    49.341,
    49.479,
    49.666,
    49.904,
    50.196,
    50.545,
    50.955,
    51.431,
    51.975,
    52.593,
    53.290,
    54.071,
    54.943,
    55.912,
    56.987,
    58.178,
    57.533,
    57.061,
    56.746,
    56.575,
    56.534,
    56.611,
    56.795,
    57.074,
    57.440,
    57.882,
    58.393,
    58.964,
    59.589,
    60.261,
    60.974,
    61.724,
    62.505,
    63.314,
    64.146,
    64.999,
    65.868,
    66.753,
    67.649,
    68.555,
    69.470,
    70.391,
    71.317,
    72.247,
    73.179,
    74.112,
    75.046,
    75.980,
    76.912,
    77.843,
    78.772,
    79.698,
    80.620,
    81.539,
    82.455,
    83.366,
    84.272,
    85.174,
    86.071,
    86.963,
    87.850,
    88.731,
    89.607,
    90.478,
    91.344,
    92.203,
    93.057,
    93.906,
    94.749,
    95.586,
    96.418,
    97.244,
    98.064,
    98.879,
    99.688,
    100.492,
    100.625,
    100.761,
    100.902,
    101.046,
    101.194,
    101.347,
    101.503,
    101.664,
    101.829,
    101.999,
    102.172,
    102.351,
    102.534,
    102.721,
    102.914,
    103.111,
    103.313,
    103.520,
    103.732,
    103.949,
    104.171,
    104.399,
    104.632,
    104.871,
    105.115,
    105.364,
    105.620,
    105.881,
    106.148,
    106.421,
    106.700,
    106.986,
    107.278,
    107.576,
    107.880,
    108.191,
    108.509,
    108.833,
    109.165,
    109.503,
    109.848
])

def f(x, ax, bx, cx, ay, by, cy, off):
    hr = x
    hr2 = 2.0 * hr
    hr3 = 3.0 * hr
    M = (ax * np.cos(hr) +
         bx * np.cos(hr2) +
         cx * np.cos(hr3) +
         ay * np.sin(hr) +
         by * np.sin(hr2) +
         cy * np.sin(hr3) +
         off)
    return M

def fit_M(x_data, coeffs):
    ax, bx, cx, ay, by, cy, off = coeffs
    M_values = []
    for x in x_data:
        M_values.append(f(x, ax, bx, cx, ay, by, cy, off))
    return np.array(M_values)

# Define the error function to minimize
def error_function(coeffs, x_data, reference_data):
    M_values = fit_M(x_data, coeffs)
    return np.sum((M_values - reference_data) ** 2)

# Define the weighted error function to minimize
def weighted_error_function(coeffs, x_data, reference_data, weights):
    M_values = fit_M(x_data, coeffs)
    return np.sum(weights * (M_values - reference_data) ** 2)

# Define the error function with L2 regularization
def regularized_error_function(coeffs, x_data, reference_data, lambda_reg):
    M_values = fit_M(x_data, coeffs)
    error = np.sum((M_values - reference_data) ** 2)
    regularization_term = lambda_reg * np.sum(np.square(coeffs))
    return error + regularization_term

# Define the error function with penalty term
def penalized_error_function(coeffs, x_data, reference_data, penalty_factor):
    M_values = fit_M(x_data, coeffs)
    error = np.sum((M_values - reference_data) ** 2)
    penalty = np.sum(np.maximum(0, reference_data - M_values) ** 2)
    return error + penalty_factor * penalty

# Initial guess for the coefficients
initial_guess = [10, 10, 8.0, 8.0, -10, 8.0, 60.0]

lambda_reg = 0.01  # Regularization parameter
penalty_factor = 5.0  # Penalty factor

# Minimize the error function
#result = minimize(regularized_error_function, initial_guess, args=(x_ap1, y_ap1, lambda_reg)).x
#result = minimize(error_function, initial_guess, args=(x_ap1, y_ap1)).x
result = minimize(penalized_error_function, initial_guess, args=(x_ap1, y_ap1, penalty_factor)).x
#result, pcov = curve_fit(f, x_ap1, y_ap1, p0=initial_guess)

# Extract the optimized coefficients
optimized_coeffs = result
ax, bx, cx, ay, by, cy, off = optimized_coeffs

print("Optimized coefficients:")
print(f"ax: {ax}, bx: {bx}, cx: {cx}")
print(f"ay: {ay}, by: {by}, cy: {cy}")
print(f"off: {off}")

M_values = []
for x in x_ap1:
    M_values.append(f(x, ax, bx, cx, ay, by, cy, off))

plt.figure(figsize=(8, 8))
px = plt.subplot(111, projection='polar')
px.plot(x_ap1, y_ap1, label='AP1 cusp 100 nits')
px.plot(x_ap1, M_values, label='AP1 cusp approx 6 coeffs')
px.scatter(x_corners, y_corners, color='red', label='AP1 cusp corners')
px.legend(loc='upper right', bbox_to_anchor=(1.1, 1.1))
px.set_theta_zero_location('N')
px.set_theta_direction(-1)
plt.show()

fig, px = plt.subplots(figsize=(10, 6))
px.plot(x_ap1, y_ap1, label='AP1 cusp 100 nits')
px.plot(x_ap1, M_values, label='AP1 cusp approx 6 coeffs')
px.scatter(x_corners, y_corners, color='red', label='AP1 cusp corners')
plt.xlabel('x')
plt.ylabel('y')
px.legend()
px.grid(True)
plt.show()