Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

cuve_fit results are inconsistent over different CPU / Linux distro #20686

Closed
seongsujeong opened this issue May 10, 2024 · 3 comments
Closed
Labels
query A question or suggestion that requires further information scipy.optimize

Comments

@seongsujeong
Copy link

seongsujeong commented May 10, 2024

Describe your issue.

It has come to our attention that the curve_fit function in scipy.optimize produces different fitting results across platforms, even when running identical code and using the same data samples as shown in the sample code below.

The function in question models a trimodal distribution, which is the sum of three Gaussian functions with varying means, standard deviations, and amplitudes (or peak heights).

The fitting results differ across platforms, as detailed in the Error message section below.

Reproducing Code Example

'''
Example trimodal fitting code to be reported to SciPy

Code and sample data based on works from Junkyo Jung (JPL)
'''

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit

def map_param(param_name_list, param_array):
    '''
    A helper function to map the list of parameters to keyword parameters

    Parameters
    ----------
    param_name_list: list
        List of parameter names
    
    param_array: list or np.ndarray
        array of parameters

    returns
    -------
    param_dict: dict
        keyword parameter
    '''
    param_dict = {}

    for i_param, param_name in enumerate(param_name_list):
        param_dict[param_name] = param_array[i_param]
    
    return param_dict

# Define trimodal function
def gauss(array, mu, sigma, amplitude):
    """
    Generate gaussian distribution with given mean and std

    Parameters
    ----------

    array: np.ndarray
        Input value for the Gaussian function
    
    mu: float
        Mean value of the Gaussian function
    
    sigma: float
        Standard deviation of the Gaussian function

    amplitude: float
        Amplitude (or, "peak height") of the Gaussian function
    
    
    Returns
    -------
    _: np.ndarray
        Output of the Gaussian value
    """
    return amplitude * np.exp(-(array - mu)**2 / (2 * sigma**2))


def trimodal(array, mu1, sigma1, amplitude1,
             mu2, sigma2, amplitude2,
              mu3, sigma3, amplitude3):
    """
    Generate trimodal gaussian distribution with given means and stds

    Parameters
    ----------

    array: np.ndarray
        Input of the trimodal function

    mu1, sigma1, amplitude1: float
        mean value, standard deviation of the 1st Gaussian function
    
    mu2, sigma2, amplitude2: float
        mean value, standard deviation of the 2nd Gaussian function

    mu3, sigma3, amplitude3: float
        mean value, standard deviation of the 3rd Gaussian function
    
    Returns
    -------
    _ : np.ndarray
        Out of the trimodal function
    """
    return gauss(array, mu1, sigma1, amplitude1) + \
        gauss(array, mu2, sigma2, amplitude2) + \
        gauss(array, mu3, sigma3, amplitude3)


# Define the sample data, which is a histogram
bin_str = """
-35.  -34.9 -34.8 -34.7 -34.6 -34.5 -34.4 -34.3 -34.2 -34.1 -34.  -33.9
 -33.8 -33.7 -33.6 -33.5 -33.4 -33.3 -33.2 -33.1 -33.  -32.9 -32.8 -32.7
 -32.6 -32.5 -32.4 -32.3 -32.2 -32.1 -32.  -31.9 -31.8 -31.7 -31.6 -31.5
 -31.4 -31.3 -31.2 -31.1 -31.  -30.9 -30.8 -30.7 -30.6 -30.5 -30.4 -30.3
 -30.2 -30.1 -30.  -29.9 -29.8 -29.7 -29.6 -29.5 -29.4 -29.3 -29.2 -29.1
 -29.  -28.9 -28.8 -28.7 -28.6 -28.5 -28.4 -28.3 -28.2 -28.1 -28.  -27.9
 -27.8 -27.7 -27.6 -27.5 -27.4 -27.3 -27.2 -27.1 -27.  -26.9 -26.8 -26.7
 -26.6 -26.5 -26.4 -26.3 -26.2 -26.1 -26.  -25.9 -25.8 -25.7 -25.6 -25.5
 -25.4 -25.3 -25.2 -25.1 -25.  -24.9 -24.8 -24.7 -24.6 -24.5 -24.4 -24.3
 -24.2 -24.1 -24.  -23.9 -23.8 -23.7 -23.6 -23.5 -23.4 -23.3 -23.2 -23.1
 -23.  -22.9 -22.8 -22.7 -22.6 -22.5 -22.4 -22.3 -22.2 -22.1 -22.  -21.9
 -21.8 -21.7 -21.6 -21.5 -21.4 -21.3 -21.2 -21.1 -21.  -20.9 -20.8 -20.7
 -20.6 -20.5 -20.4 -20.3 -20.2 -20.1 -20.  -19.9 -19.8 -19.7 -19.6 -19.5
 -19.4 -19.3 -19.2 -19.1 -19.  -18.9 -18.8 -18.7 -18.6 -18.5 -18.4 -18.3
 -18.2 -18.1 -18.  -17.9 -17.8 -17.7 -17.6 -17.5 -17.4 -17.3 -17.2 -17.1
 -17.  -16.9 -16.8 -16.7 -16.6 -16.5 -16.4 -16.3 -16.2 -16.1 -16.  -15.9
 -15.8 -15.7 -15.6 -15.5 -15.4 -15.3 -15.2 -15.1 -15.  -14.9 -14.8 -14.7
 -14.6 -14.5 -14.4 -14.3 -14.2 -14.1 -14.  -13.9 -13.8 -13.7 -13.6 -13.5
 -13.4 -13.3 -13.2 -13.1 -13.  -12.9 -12.8 -12.7 -12.6 -12.5 -12.4 -12.3
 -12.2 -12.1 -12.  -11.9 -11.8 -11.7 -11.6 -11.5 -11.4 -11.3 -11.2 -11.1
 -11.  -10.9 -10.8 -10.7 -10.6 -10.5 -10.4 -10.3 -10.2 -10.1 -10.   -9.9
  -9.8  -9.7  -9.6  -9.5  -9.4  -9.3  -9.2  -9.1  -9.   -8.9  -8.8  -8.7
  -8.6  -8.5  -8.4  -8.3  -8.2  -8.1  -8.   -7.9  -7.8  -7.7  -7.6  -7.5
  -7.4  -7.3  -7.2  -7.1  -7.   -6.9  -6.8  -6.7  -6.6  -6.5  -6.4  -6.3
  -6.2  -6.1  -6.   -5.9  -5.8  -5.7  -5.6  -5.5  -5.4  -5.3  -5.2  -5.1
  -5.   -4.9  -4.8  -4.7  -4.6  -4.5  -4.4  -4.3  -4.2  -4.1  -4.   -3.9
  -3.8  -3.7  -3.6  -3.5  -3.4  -3.3  -3.2  -3.1  -3.   -2.9  -2.8  -2.7
  -2.6  -2.5  -2.4  -2.3  -2.2  -2.1  -2.   -1.9  -1.8  -1.7  -1.6  -1.5
  -1.4  -1.3  -1.2  -1.1  -1.   -0.9  -0.8  -0.7  -0.6  -0.5  -0.4  -0.3
  -0.2  -0.1   0.    0.1   0.2   0.3   0.4   0.5   0.6   0.7   0.8   0.9
   1.    1.1   1.2   1.3   1.4   1.5   1.6   1.7   1.8   1.9   2.    2.1
   2.2   2.3   2.4   2.5   2.6   2.7   2.8   2.9   3.    3.1   3.2   3.3
   3.4   3.5   3.6   3.7   3.8   3.9   4.    4.1   4.2   4.3   4.4   4.5
   4.6   4.7   4.8   4.9   5.    5.1   5.2   5.3   5.4   5.5   5.6   5.7
   5.8   5.9   6.    6.1   6.2   6.3   6.4   6.5   6.6   6.7   6.8   6.9
   7.    7.1   7.2   7.3   7.4   7.5   7.6   7.7   7.8   7.9   8.    8.1
   8.2   8.3   8.4   8.5   8.6   8.7   8.8   8.9   9.    9.1   9.2   9.3
   9.4   9.5   9.6   9.7   9.8   9.9
"""

count_str = """
0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
 0.    0.    0.    0.    0.    0.    0.    0.    0.004 0.004 0.    0.
 0.    0.004 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.004
 0.    0.    0.    0.    0.    0.    0.008 0.004 0.004 0.    0.004 0.008
 0.008 0.004 0.    0.    0.004 0.008 0.012 0.004 0.    0.016 0.016 0.016
 0.024 0.024 0.02  0.036 0.028 0.024 0.068 0.084 0.08  0.096 0.124 0.156
 0.176 0.136 0.116 0.112 0.1   0.088 0.056 0.052 0.048 0.032 0.024 0.056
 0.036 0.044 0.008 0.04  0.028 0.02  0.028 0.016 0.012 0.016 0.016 0.016
 0.024 0.012 0.016 0.016 0.02  0.016 0.012 0.02  0.02  0.028 0.012 0.012
 0.02  0.004 0.004 0.016 0.012 0.012 0.012 0.012 0.032 0.024 0.004 0.012
 0.004 0.02  0.012 0.024 0.012 0.024 0.016 0.016 0.012 0.016 0.016 0.024
 0.02  0.028 0.016 0.024 0.04  0.028 0.032 0.024 0.044 0.02  0.02  0.036
 0.048 0.024 0.02  0.028 0.032 0.032 0.02  0.016 0.044 0.02  0.036 0.024
 0.02  0.004 0.016 0.036 0.036 0.032 0.012 0.016 0.012 0.008 0.024 0.036
 0.06  0.052 0.116 0.116 0.072 0.088 0.072 0.116 0.092 0.092 0.072 0.06
 0.064 0.048 0.052 0.036 0.044 0.04  0.072 0.064 0.052 0.072 0.044 0.076
 0.092 0.092 0.088 0.108 0.1   0.108 0.104 0.072 0.124 0.112 0.128 0.16
 0.144 0.196 0.16  0.22  0.18  0.22  0.204 0.22  0.156 0.188 0.184 0.128
 0.12  0.152 0.132 0.124 0.076 0.068 0.036 0.032 0.028 0.06  0.06  0.032
 0.028 0.012 0.008 0.016 0.    0.    0.008 0.012 0.004 0.008 0.012 0.
 0.    0.004 0.004 0.004 0.004 0.008 0.008 0.    0.012 0.004 0.    0.004
 0.004 0.004 0.    0.008 0.    0.008 0.004 0.    0.    0.    0.004 0.
 0.    0.    0.    0.012 0.    0.    0.    0.    0.    0.    0.    0.
 0.    0.004 0.    0.004 0.    0.004 0.004 0.    0.    0.    0.004 0.
 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.004 0.
 0.    0.004 0.    0.    0.    0.    0.008 0.    0.    0.    0.    0.
 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.004
 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
 0.    0.    0.004 0.    0.    0.    0.    0.    0.    0.    0.    0.
 0.    0.    0.    0.    0.    0.004 0.    0.    0.    0.    0.    0.
 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
 0.    0.    0.    0.    0.    0.   
"""

count = np.array([float(number) for number in count_str.split()],
                 dtype=np.float64)
bins = np.array([float(number) for number in bin_str.split()],
                dtype=np.float64)

#initial value of the parameters
initvel = [-18.296757, 0.5, 0.1759999999999975,
            -11.750452, 0.5, 0.2200000000000047,
            -15.023604393005371, 0.5, 0.1]

# boundary to constrain the parameters
bounds = ([-35, 0, 0.01,
           -35, 0, 0.01,
           -35, 0, 0.01],
          [5, 5, 0.95,
           5, 5, 0.95,
           5, 5, 0.95])

# Run curve_fit
params, cov, infodict, mesg, ier = curve_fit(trimodal,
                                             bins,
                                             count,
                                             initvel,
                                             full_output=True,
                                             method='trf',
                                             bounds=bounds)

# print out the fitting result
print('SciPy curve_fit result:')
param_names = ['mu1', 'sigma1', 'amplitude1',
               'mu2', 'sigma2', 'amplitude2',
               'mu3', 'sigma3', 'amplitude3']
for i_param, param_name in enumerate(param_names):
    print(f'{param_name}: {params[i_param]:0.8f}', end='\t')
    if i_param % 3 == 2:
        print()

# Prepare to plot the fitting result and other reference data
y_fit_cf = trimodal(bins, **map_param(param_names, params))
y_fit_init = trimodal(bins, **map_param(param_names, initvel))

plt.plot(bins, count, 'r', alpha=0.5)
plt.plot(bins, y_fit_init, 'k--', alpha=0.5)
plt.plot(bins, y_fit_cf, 'g', alpha=0.5)

plt.legend(['Histogram', 'Initial parameters', 'curve_fit'])
plt.grid()
plt.show()

Error message

Intel Xeon Platinum 8259CL, Oracle Linux Server 8.9

SciPy curve_fit result:
mu1: -18.45068882	sigma1: 0.01106007	amplitude1: 0.02774448	
mu2: -9.40209468	sigma2: 1.04460988	amplitude2: 0.17496974	
mu3: -18.10321619	sigma3: 5.00000000	amplitude3: 0.03973587

AMD EPYC 7663, CentOS Linux 7

SciPy curve_fit result:
mu1: 4.48314458	sigma1: 0.01034696	amplitude1: 0.01507123	
mu2: -9.40203992	sigma2: 1.04450772	amplitude2: 0.17497326	
mu3: -18.10214019	sigma3: 5.00000000	amplitude3: 0.03973711	


### SciPy/NumPy/Python version and system information

```shell
For the Intel w/ Oracle linux:
1.12.0 1.26.4 sys.version_info(major=3, minor=12, micro=2, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: pkgconfig
    found: true
    include directory: $HOME/miniconda3/envs/isce3_0.20.0/include
    lib directory: $HOME/miniconda3/envs/isce3_0.20.0/lib
    name: blas
    openblas configuration: unknown
    pc file directory: $HOME/miniconda3/envs/isce3_0.20.0/lib/pkgconfig
    version: 3.9.0
  lapack:
    detection method: pkgconfig
    found: true
    include directory: $HOME/miniconda3/envs/isce3_0.20.0/include
    lib directory: $HOME/miniconda3/envs/isce3_0.20.0/lib
    name: lapack
    openblas configuration: unknown
    pc file directory: $HOME/miniconda3/envs/isce3_0.20.0/lib/pkgconfig
    version: 3.9.0
  pybind11:
    detection method: pkgconfig
    include directory: $HOME/miniconda3/envs/isce3_0.20.0/include
    name: pybind11
    version: 2.11.1
Compilers:
  c:
    args: -march=nocona, -mtune=haswell, -ftree-vectorize, -fPIC, -fstack-protector-strong,
      -fno-plt, -O2, -ffunction-sections, -pipe, -isystem, $HOME/miniconda3/envs/isce3_0.20.0/include,
      -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/work=/usr/local/src/conda/scipy-split-1.12.0,
      -fdebug-prefix-map=$HOME/miniconda3/envs/isce3_0.20.0=/usr/local/src/conda-prefix,
      -march=nocona, -mtune=haswell, -ftree-vectorize, -fPIC, -fstack-protector-strong,
      -fno-plt, -O2, -ffunction-sections, -pipe, -isystem, $HOME/miniconda3/envs/isce3_0.20.0/include,
      -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/work=/usr/local/src/conda/scipy-split-1.12.0,
      -fdebug-prefix-map=$HOME/miniconda3/envs/isce3_0.20.0=/usr/local/src/conda-prefix,
      -DNDEBUG, -D_FORTIFY_SOURCE=2, -O2, -isystem, $HOME/miniconda3/envs/isce3_0.20.0/include,
      -DNDEBUG, -D_FORTIFY_SOURCE=2, -O2, -isystem, $HOME/miniconda3/envs/isce3_0.20.0/include
    commands: /home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/_build_env/bin/x86_64-conda-linux-gnu-cc
    linker: ld.bfd
    linker args: -Wl,-O2, -Wl,--sort-common, -Wl,--as-needed, -Wl,-z,relro, -Wl,-z,now,
      -Wl,--disable-new-dtags, -Wl,--gc-sections, -Wl,--allow-shlib-undefined, -Wl,-rpath,$HOME/miniconda3/envs/isce3_0.20.0/lib,
      -Wl,-rpath-link,$HOME/miniconda3/envs/isce3_0.20.0/lib, -L$HOME/miniconda3/envs/isce3_0.20.0/lib,
      -Wl,-O2, -Wl,--sort-common, -Wl,--as-needed, -Wl,-z,relro, -Wl,-z,now, -Wl,--disable-new-dtags,
      -Wl,--gc-sections, -Wl,--allow-shlib-undefined, -Wl,-rpath,$HOME/miniconda3/envs/isce3_0.20.0/lib,
      -Wl,-rpath-link,$HOME/miniconda3/envs/isce3_0.20.0/lib, -L$HOME/miniconda3/envs/isce3_0.20.0/lib,
      -march=nocona, -mtune=haswell, -ftree-vectorize, -fPIC, -fstack-protector-strong,
      -fno-plt, -O2, -ffunction-sections, -pipe, -isystem, $HOME/miniconda3/envs/isce3_0.20.0/include,
      -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/work=/usr/local/src/conda/scipy-split-1.12.0,
      -fdebug-prefix-map=$HOME/miniconda3/envs/isce3_0.20.0=/usr/local/src/conda-prefix,
      -march=nocona, -mtune=haswell, -ftree-vectorize, -fPIC, -fstack-protector-strong,
      -fno-plt, -O2, -ffunction-sections, -pipe, -isystem, $HOME/miniconda3/envs/isce3_0.20.0/include,
      -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/work=/usr/local/src/conda/scipy-split-1.12.0,
      -fdebug-prefix-map=$HOME/miniconda3/envs/isce3_0.20.0=/usr/local/src/conda-prefix,
      -DNDEBUG, -D_FORTIFY_SOURCE=2, -O2, -isystem, $HOME/miniconda3/envs/isce3_0.20.0/include,
      -DNDEBUG, -D_FORTIFY_SOURCE=2, -O2, -isystem, $HOME/miniconda3/envs/isce3_0.20.0/include
    name: gcc
    version: 12.3.0
  c++:
    args: -fvisibility-inlines-hidden, -fmessage-length=0, -march=nocona, -mtune=haswell,
      -ftree-vectorize, -fPIC, -fstack-protector-strong, -fno-plt, -O2, -ffunction-sections,
      -pipe, -isystem, $HOME/miniconda3/envs/isce3_0.20.0/include, -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/work=/usr/local/src/conda/scipy-split-1.12.0,
      -fdebug-prefix-map=$HOME/miniconda3/envs/isce3_0.20.0=/usr/local/src/conda-prefix,
      -fvisibility-inlines-hidden, -fmessage-length=0, -march=nocona, -mtune=haswell,
      -ftree-vectorize, -fPIC, -fstack-protector-strong, -fno-plt, -O2, -ffunction-sections,
      -pipe, -isystem, $HOME/miniconda3/envs/isce3_0.20.0/include, -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/work=/usr/local/src/conda/scipy-split-1.12.0,
      -fdebug-prefix-map=$HOME/miniconda3/envs/isce3_0.20.0=/usr/local/src/conda-prefix,
      -DNDEBUG, -D_FORTIFY_SOURCE=2, -O2, -isystem, $HOME/miniconda3/envs/isce3_0.20.0/include,
      -DNDEBUG, -D_FORTIFY_SOURCE=2, -O2, -isystem, $HOME/miniconda3/envs/isce3_0.20.0/include
    commands: /home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/_build_env/bin/x86_64-conda-linux-gnu-c++
    linker: ld.bfd
    linker args: -Wl,-O2, -Wl,--sort-common, -Wl,--as-needed, -Wl,-z,relro, -Wl,-z,now,
      -Wl,--disable-new-dtags, -Wl,--gc-sections, -Wl,--allow-shlib-undefined, -Wl,-rpath,$HOME/miniconda3/envs/isce3_0.20.0/lib,
      -Wl,-rpath-link,$HOME/miniconda3/envs/isce3_0.20.0/lib, -L$HOME/miniconda3/envs/isce3_0.20.0/lib,
      -Wl,-O2, -Wl,--sort-common, -Wl,--as-needed, -Wl,-z,relro, -Wl,-z,now, -Wl,--disable-new-dtags,
      -Wl,--gc-sections, -Wl,--allow-shlib-undefined, -Wl,-rpath,$HOME/miniconda3/envs/isce3_0.20.0/lib,
      -Wl,-rpath-link,$HOME/miniconda3/envs/isce3_0.20.0/lib, -L$HOME/miniconda3/envs/isce3_0.20.0/lib,
      -fvisibility-inlines-hidden, -fmessage-length=0, -march=nocona, -mtune=haswell,
      -ftree-vectorize, -fPIC, -fstack-protector-strong, -fno-plt, -O2, -ffunction-sections,
      -pipe, -isystem, $HOME/miniconda3/envs/isce3_0.20.0/include, -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/work=/usr/local/src/conda/scipy-split-1.12.0,
      -fdebug-prefix-map=$HOME/miniconda3/envs/isce3_0.20.0=/usr/local/src/conda-prefix,
      -fvisibility-inlines-hidden, -fmessage-length=0, -march=nocona, -mtune=haswell,
      -ftree-vectorize, -fPIC, -fstack-protector-strong, -fno-plt, -O2, -ffunction-sections,
      -pipe, -isystem, $HOME/miniconda3/envs/isce3_0.20.0/include, -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/work=/usr/local/src/conda/scipy-split-1.12.0,
      -fdebug-prefix-map=$HOME/miniconda3/envs/isce3_0.20.0=/usr/local/src/conda-prefix,
      -DNDEBUG, -D_FORTIFY_SOURCE=2, -O2, -isystem, $HOME/miniconda3/envs/isce3_0.20.0/include,
      -DNDEBUG, -D_FORTIFY_SOURCE=2, -O2, -isystem, $HOME/miniconda3/envs/isce3_0.20.0/include
    name: gcc
    version: 12.3.0
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 3.0.7
  fortran:
    args: -march=nocona, -mtune=haswell, -ftree-vectorize, -fPIC, -fstack-protector-strong,
      -fno-plt, -O2, -ffunction-sections, -pipe, -isystem, $HOME/miniconda3/envs/isce3_0.20.0/include,
      -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/work=/usr/local/src/conda/scipy-split-1.12.0,
      -fdebug-prefix-map=$HOME/miniconda3/envs/isce3_0.20.0=/usr/local/src/conda-prefix
    commands: /home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/_build_env/bin/x86_64-conda-linux-gnu-gfortran
    linker: ld.bfd
    linker args: -Wl,-O2, -Wl,--sort-common, -Wl,--as-needed, -Wl,-z,relro, -Wl,-z,now,
      -Wl,--disable-new-dtags, -Wl,--gc-sections, -Wl,--allow-shlib-undefined, -Wl,-rpath,$HOME/miniconda3/envs/isce3_0.20.0/lib,
      -Wl,-rpath-link,$HOME/miniconda3/envs/isce3_0.20.0/lib, -L$HOME/miniconda3/envs/isce3_0.20.0/lib,
      -Wl,-O2, -Wl,--sort-common, -Wl,--as-needed, -Wl,-z,relro, -Wl,-z,now, -Wl,--disable-new-dtags,
      -Wl,--gc-sections, -Wl,--allow-shlib-undefined, -Wl,-rpath,$HOME/miniconda3/envs/isce3_0.20.0/lib,
      -Wl,-rpath-link,$HOME/miniconda3/envs/isce3_0.20.0/lib, -L$HOME/miniconda3/envs/isce3_0.20.0/lib,
      -march=nocona, -mtune=haswell, -ftree-vectorize, -fPIC, -fstack-protector-strong,
      -fno-plt, -O2, -ffunction-sections, -pipe, -isystem, $HOME/miniconda3/envs/isce3_0.20.0/include,
      -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/work=/usr/local/src/conda/scipy-split-1.12.0,
      -fdebug-prefix-map=$HOME/miniconda3/envs/isce3_0.20.0=/usr/local/src/conda-prefix
    name: gcc
    version: 12.3.0
  pythran:
    include directory: ../../_h_env_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold/lib/python3.12/site-packages/pythran
    version: 0.15.0
Machine Information:
  build:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
  cross-compiled: false
  host:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
Python Information:
  path: $HOME/miniconda3/envs/isce3_0.20.0/bin/python
  version: '3.12'



for AMD with CentOS:
1.12.0 1.26.4 sys.version_info(major=3, minor=12, micro=2, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: pkgconfig
    found: true
    include directory: $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/include
    lib directory: $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/lib
    name: blas
    openblas configuration: unknown
    pc file directory: $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/lib/pkgconfig
    version: 3.9.0
  lapack:
    detection method: pkgconfig
    found: true
    include directory: $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/include
    lib directory: $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/lib
    name: lapack
    openblas configuration: unknown
    pc file directory: $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/lib/pkgconfig
    version: 3.9.0
  pybind11:
    detection method: pkgconfig
    include directory: $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/include
    name: pybind11
    version: 2.11.1
Compilers:
  c:
    args: -march=nocona, -mtune=haswell, -ftree-vectorize, -fPIC, -fstack-protector-strong,
      -fno-plt, -O2, -ffunction-sections, -pipe, -isystem, $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/include,
      -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/work=/usr/local/src/conda/scipy-split-1.12.0,
      -fdebug-prefix-map=$HOME/tools/mambaforge/envs/isce3_0.19.2.conda=/usr/local/src/conda-prefix,
      -march=nocona, -mtune=haswell, -ftree-vectorize, -fPIC, -fstack-protector-strong,
      -fno-plt, -O2, -ffunction-sections, -pipe, -isystem, $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/include,
      -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/work=/usr/local/src/conda/scipy-split-1.12.0,
      -fdebug-prefix-map=$HOME/tools/mambaforge/envs/isce3_0.19.2.conda=/usr/local/src/conda-prefix,
      -DNDEBUG, -D_FORTIFY_SOURCE=2, -O2, -isystem, $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/include,
      -DNDEBUG, -D_FORTIFY_SOURCE=2, -O2, -isystem, $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/include
    commands: /home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/_build_env/bin/x86_64-conda-linux-gnu-cc
    linker: ld.bfd
    linker args: -Wl,-O2, -Wl,--sort-common, -Wl,--as-needed, -Wl,-z,relro, -Wl,-z,now,
      -Wl,--disable-new-dtags, -Wl,--gc-sections, -Wl,--allow-shlib-undefined, -Wl,-rpath,$HOME/tools/mambaforge/envs/isce3_0.19.2.conda/lib,
      -Wl,-rpath-link,$HOME/tools/mambaforge/envs/isce3_0.19.2.conda/lib,
      -L$HOME/tools/mambaforge/envs/isce3_0.19.2.conda/lib, -Wl,-O2,
      -Wl,--sort-common, -Wl,--as-needed, -Wl,-z,relro, -Wl,-z,now, -Wl,--disable-new-dtags,
      -Wl,--gc-sections, -Wl,--allow-shlib-undefined, -Wl,-rpath,$HOME/tools/mambaforge/envs/isce3_0.19.2.conda/lib,
      -Wl,-rpath-link,$HOME/tools/mambaforge/envs/isce3_0.19.2.conda/lib,
      -L$HOME/tools/mambaforge/envs/isce3_0.19.2.conda/lib, -march=nocona,
      -mtune=haswell, -ftree-vectorize, -fPIC, -fstack-protector-strong, -fno-plt,
      -O2, -ffunction-sections, -pipe, -isystem, $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/include,
      -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/work=/usr/local/src/conda/scipy-split-1.12.0,
      -fdebug-prefix-map=$HOME/tools/mambaforge/envs/isce3_0.19.2.conda=/usr/local/src/conda-prefix,
      -march=nocona, -mtune=haswell, -ftree-vectorize, -fPIC, -fstack-protector-strong,
      -fno-plt, -O2, -ffunction-sections, -pipe, -isystem, $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/include,
      -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/work=/usr/local/src/conda/scipy-split-1.12.0,
      -fdebug-prefix-map=$HOME/tools/mambaforge/envs/isce3_0.19.2.conda=/usr/local/src/conda-prefix,
      -DNDEBUG, -D_FORTIFY_SOURCE=2, -O2, -isystem, $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/include,
      -DNDEBUG, -D_FORTIFY_SOURCE=2, -O2, -isystem, $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/include
    name: gcc
    version: 12.3.0
  c++:
    args: -fvisibility-inlines-hidden, -fmessage-length=0, -march=nocona, -mtune=haswell,
      -ftree-vectorize, -fPIC, -fstack-protector-strong, -fno-plt, -O2, -ffunction-sections,
      -pipe, -isystem, $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/include,
      -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/work=/usr/local/src/conda/scipy-split-1.12.0,
      -fdebug-prefix-map=$HOME/tools/mambaforge/envs/isce3_0.19.2.conda=/usr/local/src/conda-prefix,
      -fvisibility-inlines-hidden, -fmessage-length=0, -march=nocona, -mtune=haswell,
      -ftree-vectorize, -fPIC, -fstack-protector-strong, -fno-plt, -O2, -ffunction-sections,
      -pipe, -isystem, $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/include,
      -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/work=/usr/local/src/conda/scipy-split-1.12.0,
      -fdebug-prefix-map=$HOME/tools/mambaforge/envs/isce3_0.19.2.conda=/usr/local/src/conda-prefix,
      -DNDEBUG, -D_FORTIFY_SOURCE=2, -O2, -isystem, $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/include,
      -DNDEBUG, -D_FORTIFY_SOURCE=2, -O2, -isystem, $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/include
    commands: /home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/_build_env/bin/x86_64-conda-linux-gnu-c++
    linker: ld.bfd
    linker args: -Wl,-O2, -Wl,--sort-common, -Wl,--as-needed, -Wl,-z,relro, -Wl,-z,now,
      -Wl,--disable-new-dtags, -Wl,--gc-sections, -Wl,--allow-shlib-undefined, -Wl,-rpath,$HOME/tools/mambaforge/envs/isce3_0.19.2.conda/lib,
      -Wl,-rpath-link,$HOME/tools/mambaforge/envs/isce3_0.19.2.conda/lib,
      -L$HOME/tools/mambaforge/envs/isce3_0.19.2.conda/lib, -Wl,-O2,
      -Wl,--sort-common, -Wl,--as-needed, -Wl,-z,relro, -Wl,-z,now, -Wl,--disable-new-dtags,
      -Wl,--gc-sections, -Wl,--allow-shlib-undefined, -Wl,-rpath,$HOME/tools/mambaforge/envs/isce3_0.19.2.conda/lib,
      -Wl,-rpath-link,$HOME/tools/mambaforge/envs/isce3_0.19.2.conda/lib,
      -L$HOME/tools/mambaforge/envs/isce3_0.19.2.conda/lib, -fvisibility-inlines-hidden,
      -fmessage-length=0, -march=nocona, -mtune=haswell, -ftree-vectorize, -fPIC,
      -fstack-protector-strong, -fno-plt, -O2, -ffunction-sections, -pipe, -isystem,
      $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/include, -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/work=/usr/local/src/conda/scipy-split-1.12.0,
      -fdebug-prefix-map=$HOME/tools/mambaforge/envs/isce3_0.19.2.conda=/usr/local/src/conda-prefix,
      -fvisibility-inlines-hidden, -fmessage-length=0, -march=nocona, -mtune=haswell,
      -ftree-vectorize, -fPIC, -fstack-protector-strong, -fno-plt, -O2, -ffunction-sections,
      -pipe, -isystem, $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/include,
      -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/work=/usr/local/src/conda/scipy-split-1.12.0,
      -fdebug-prefix-map=$HOME/tools/mambaforge/envs/isce3_0.19.2.conda=/usr/local/src/conda-prefix,
      -DNDEBUG, -D_FORTIFY_SOURCE=2, -O2, -isystem, $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/include,
      -DNDEBUG, -D_FORTIFY_SOURCE=2, -O2, -isystem, $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/include
    name: gcc
    version: 12.3.0
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 3.0.7
  fortran:
    args: -march=nocona, -mtune=haswell, -ftree-vectorize, -fPIC, -fstack-protector-strong,
      -fno-plt, -O2, -ffunction-sections, -pipe, -isystem, $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/include,
      -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/work=/usr/local/src/conda/scipy-split-1.12.0,
      -fdebug-prefix-map=$HOME/tools/mambaforge/envs/isce3_0.19.2.conda=/usr/local/src/conda-prefix
    commands: /home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/_build_env/bin/x86_64-conda-linux-gnu-gfortran
    linker: ld.bfd
    linker args: -Wl,-O2, -Wl,--sort-common, -Wl,--as-needed, -Wl,-z,relro, -Wl,-z,now,
      -Wl,--disable-new-dtags, -Wl,--gc-sections, -Wl,--allow-shlib-undefined, -Wl,-rpath,$HOME/tools/mambaforge/envs/isce3_0.19.2.conda/lib,
      -Wl,-rpath-link,$HOME/tools/mambaforge/envs/isce3_0.19.2.conda/lib,
      -L$HOME/tools/mambaforge/envs/isce3_0.19.2.conda/lib, -Wl,-O2,
      -Wl,--sort-common, -Wl,--as-needed, -Wl,-z,relro, -Wl,-z,now, -Wl,--disable-new-dtags,
      -Wl,--gc-sections, -Wl,--allow-shlib-undefined, -Wl,-rpath,$HOME/tools/mambaforge/envs/isce3_0.19.2.conda/lib,
      -Wl,-rpath-link,$HOME/tools/mambaforge/envs/isce3_0.19.2.conda/lib,
      -L$HOME/tools/mambaforge/envs/isce3_0.19.2.conda/lib, -march=nocona,
      -mtune=haswell, -ftree-vectorize, -fPIC, -fstack-protector-strong, -fno-plt,
      -O2, -ffunction-sections, -pipe, -isystem, $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/include,
      -fdebug-prefix-map=/home/conda/feedstock_root/build_artifacts/scipy-split_1706041474046/work=/usr/local/src/conda/scipy-split-1.12.0,
      -fdebug-prefix-map=$HOME/tools/mambaforge/envs/isce3_0.19.2.conda=/usr/local/src/conda-prefix
    name: gcc
    version: 12.3.0
  pythran:
    include directory: ../../_h_env_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold_placehold/lib/python3.12/site-packages/pythran
    version: 0.15.0
Machine Information:
  build:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
  cross-compiled: false
  host:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
Python Information:
  path: $HOME/tools/mambaforge/envs/isce3_0.19.2.conda/bin/python
  version: '3.12'
@seongsujeong seongsujeong added the defect A clear bug or issue that prevents SciPy from being installed or used as expected label May 10, 2024
@rkern
Copy link
Member

rkern commented May 10, 2024

Not entirely surprising. Floating point arithmetic can differ between CPUs. You have a degenerate model with a number of symmetries, so the local optimizer's path will end up having a number of branching points in the loss function that differing floating point calculations could send it down, creating large differences in the final result despite the small differences in any given floating point calculation.

You could try to help this by reparameterizing your model and applying constraints so that you don't have as many symmetries. For example, consider representing the means as mu1 (same as before), dmu2 = mu2 - mu1 and dmu3 = mu3 - mu2. Then optimize the parameter set mu1, dmu2, dmu3, ... etc. with dmu2 and dmu3 bounded from below by 0. That will break some of the symmetries, though the loss surface may still have a number of those branching points remaining.

I don't think there's much we can do about it on our end. CPU differences in floating point arithmetic are a thing, and we don't have much control over it.

@nickodell
Copy link
Member

nickodell commented May 10, 2024

This strikes me as a problem which curve_fit() is going to have a hard time solving. Fitting gaussians is a problem which has a lot of local optima. Any local optimizer is going to struggle. Since you already have bounds on reasonable values, why not use a global optimizer, such as dual_annealing()?

Example:

'''
Example trimodal fitting code to be reported to SciPy

Code and sample data based on works from Jun
'''

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit, minimize, Bounds, dual_annealing


def map_param(param_name_list, param_array):
    return dict(zip(param_name_list, param_array))


# Define trimodal function
def gauss(array, mu, sigma, amplitude):
    return amplitude * np.exp(-(array - mu)**2 / (2 * sigma**2))


def trimodal(array, mu1, sigma1, amplitude1,
             mu2, sigma2, amplitude2,
              mu3, sigma3, amplitude3):
    return gauss(array, mu1, sigma1, amplitude1) + \
        gauss(array, mu2, sigma2, amplitude2) + \
        gauss(array, mu3, sigma3, amplitude3)



# Define the sample data, which is a histogram
bin_str = """
-35.  -34.9 -34.8 -34.7 -34.6 -34.5 -34.4 -34.3 -34.2 -34.1 -34.  -33.9
 -33.8 -33.7 -33.6 -33.5 -33.4 -33.3 -33.2 -33.1 -33.  -32.9 -32.8 -32.7
 -32.6 -32.5 -32.4 -32.3 -32.2 -32.1 -32.  -31.9 -31.8 -31.7 -31.6 -31.5
 -31.4 -31.3 -31.2 -31.1 -31.  -30.9 -30.8 -30.7 -30.6 -30.5 -30.4 -30.3
 -30.2 -30.1 -30.  -29.9 -29.8 -29.7 -29.6 -29.5 -29.4 -29.3 -29.2 -29.1
 -29.  -28.9 -28.8 -28.7 -28.6 -28.5 -28.4 -28.3 -28.2 -28.1 -28.  -27.9
 -27.8 -27.7 -27.6 -27.5 -27.4 -27.3 -27.2 -27.1 -27.  -26.9 -26.8 -26.7
 -26.6 -26.5 -26.4 -26.3 -26.2 -26.1 -26.  -25.9 -25.8 -25.7 -25.6 -25.5
 -25.4 -25.3 -25.2 -25.1 -25.  -24.9 -24.8 -24.7 -24.6 -24.5 -24.4 -24.3
 -24.2 -24.1 -24.  -23.9 -23.8 -23.7 -23.6 -23.5 -23.4 -23.3 -23.2 -23.1
 -23.  -22.9 -22.8 -22.7 -22.6 -22.5 -22.4 -22.3 -22.2 -22.1 -22.  -21.9
 -21.8 -21.7 -21.6 -21.5 -21.4 -21.3 -21.2 -21.1 -21.  -20.9 -20.8 -20.7
 -20.6 -20.5 -20.4 -20.3 -20.2 -20.1 -20.  -19.9 -19.8 -19.7 -19.6 -19.5
 -19.4 -19.3 -19.2 -19.1 -19.  -18.9 -18.8 -18.7 -18.6 -18.5 -18.4 -18.3
 -18.2 -18.1 -18.  -17.9 -17.8 -17.7 -17.6 -17.5 -17.4 -17.3 -17.2 -17.1
 -17.  -16.9 -16.8 -16.7 -16.6 -16.5 -16.4 -16.3 -16.2 -16.1 -16.  -15.9
 -15.8 -15.7 -15.6 -15.5 -15.4 -15.3 -15.2 -15.1 -15.  -14.9 -14.8 -14.7
 -14.6 -14.5 -14.4 -14.3 -14.2 -14.1 -14.  -13.9 -13.8 -13.7 -13.6 -13.5
 -13.4 -13.3 -13.2 -13.1 -13.  -12.9 -12.8 -12.7 -12.6 -12.5 -12.4 -12.3
 -12.2 -12.1 -12.  -11.9 -11.8 -11.7 -11.6 -11.5 -11.4 -11.3 -11.2 -11.1
 -11.  -10.9 -10.8 -10.7 -10.6 -10.5 -10.4 -10.3 -10.2 -10.1 -10.   -9.9
  -9.8  -9.7  -9.6  -9.5  -9.4  -9.3  -9.2  -9.1  -9.   -8.9  -8.8  -8.7
  -8.6  -8.5  -8.4  -8.3  -8.2  -8.1  -8.   -7.9  -7.8  -7.7  -7.6  -7.5
  -7.4  -7.3  -7.2  -7.1  -7.   -6.9  -6.8  -6.7  -6.6  -6.5  -6.4  -6.3
  -6.2  -6.1  -6.   -5.9  -5.8  -5.7  -5.6  -5.5  -5.4  -5.3  -5.2  -5.1
  -5.   -4.9  -4.8  -4.7  -4.6  -4.5  -4.4  -4.3  -4.2  -4.1  -4.   -3.9
  -3.8  -3.7  -3.6  -3.5  -3.4  -3.3  -3.2  -3.1  -3.   -2.9  -2.8  -2.7
  -2.6  -2.5  -2.4  -2.3  -2.2  -2.1  -2.   -1.9  -1.8  -1.7  -1.6  -1.5
  -1.4  -1.3  -1.2  -1.1  -1.   -0.9  -0.8  -0.7  -0.6  -0.5  -0.4  -0.3
  -0.2  -0.1   0.    0.1   0.2   0.3   0.4   0.5   0.6   0.7   0.8   0.9
   1.    1.1   1.2   1.3   1.4   1.5   1.6   1.7   1.8   1.9   2.    2.1
   2.2   2.3   2.4   2.5   2.6   2.7   2.8   2.9   3.    3.1   3.2   3.3
   3.4   3.5   3.6   3.7   3.8   3.9   4.    4.1   4.2   4.3   4.4   4.5
   4.6   4.7   4.8   4.9   5.    5.1   5.2   5.3   5.4   5.5   5.6   5.7
   5.8   5.9   6.    6.1   6.2   6.3   6.4   6.5   6.6   6.7   6.8   6.9
   7.    7.1   7.2   7.3   7.4   7.5   7.6   7.7   7.8   7.9   8.    8.1
   8.2   8.3   8.4   8.5   8.6   8.7   8.8   8.9   9.    9.1   9.2   9.3
   9.4   9.5   9.6   9.7   9.8   9.9
"""

count_str = """
0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
 0.    0.    0.    0.    0.    0.    0.    0.    0.004 0.004 0.    0.
 0.    0.004 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.004
 0.    0.    0.    0.    0.    0.    0.008 0.004 0.004 0.    0.004 0.008
 0.008 0.004 0.    0.    0.004 0.008 0.012 0.004 0.    0.016 0.016 0.016
 0.024 0.024 0.02  0.036 0.028 0.024 0.068 0.084 0.08  0.096 0.124 0.156
 0.176 0.136 0.116 0.112 0.1   0.088 0.056 0.052 0.048 0.032 0.024 0.056
 0.036 0.044 0.008 0.04  0.028 0.02  0.028 0.016 0.012 0.016 0.016 0.016
 0.024 0.012 0.016 0.016 0.02  0.016 0.012 0.02  0.02  0.028 0.012 0.012
 0.02  0.004 0.004 0.016 0.012 0.012 0.012 0.012 0.032 0.024 0.004 0.012
 0.004 0.02  0.012 0.024 0.012 0.024 0.016 0.016 0.012 0.016 0.016 0.024
 0.02  0.028 0.016 0.024 0.04  0.028 0.032 0.024 0.044 0.02  0.02  0.036
 0.048 0.024 0.02  0.028 0.032 0.032 0.02  0.016 0.044 0.02  0.036 0.024
 0.02  0.004 0.016 0.036 0.036 0.032 0.012 0.016 0.012 0.008 0.024 0.036
 0.06  0.052 0.116 0.116 0.072 0.088 0.072 0.116 0.092 0.092 0.072 0.06
 0.064 0.048 0.052 0.036 0.044 0.04  0.072 0.064 0.052 0.072 0.044 0.076
 0.092 0.092 0.088 0.108 0.1   0.108 0.104 0.072 0.124 0.112 0.128 0.16
 0.144 0.196 0.16  0.22  0.18  0.22  0.204 0.22  0.156 0.188 0.184 0.128
 0.12  0.152 0.132 0.124 0.076 0.068 0.036 0.032 0.028 0.06  0.06  0.032
 0.028 0.012 0.008 0.016 0.    0.    0.008 0.012 0.004 0.008 0.012 0.
 0.    0.004 0.004 0.004 0.004 0.008 0.008 0.    0.012 0.004 0.    0.004
 0.004 0.004 0.    0.008 0.    0.008 0.004 0.    0.    0.    0.004 0.
 0.    0.    0.    0.012 0.    0.    0.    0.    0.    0.    0.    0.
 0.    0.004 0.    0.004 0.    0.004 0.004 0.    0.    0.    0.004 0.
 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.004 0.
 0.    0.004 0.    0.    0.    0.    0.008 0.    0.    0.    0.    0.
 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.004
 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
 0.    0.    0.004 0.    0.    0.    0.    0.    0.    0.    0.    0.
 0.    0.    0.    0.    0.    0.004 0.    0.    0.    0.    0.    0.
 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
 0.    0.    0.    0.    0.    0.   
"""

count = np.array([float(number) for number in count_str.split()],
                 dtype=np.float64)
bins = np.array([float(number) for number in bin_str.split()],
                dtype=np.float64)

param_names = ['mu1', 'sigma1', 'amplitude1',
               'mu2', 'sigma2', 'amplitude2',
               'mu3', 'sigma3', 'amplitude3']


#initial value of the parameters
initvel = [-18.296757, 0.5, 0.1759999999999975,
            -11.750452, 0.5, 0.2200000000000047,
            -15.023604393005371, 0.5, 0.1]

# boundary to constrain the parameters
# NJO: Note: changed bounds to forbid zero sigma
bounds = Bounds([-35, 0.01, 0.01,
           -35, 0.01, 0.01,
           -35, 0.01, 0.01],
          [5, 5, 0.95,
           5, 5, 0.95,
           5, 5, 0.95])


def loss(params, x, y_true):
    y_pred = trimodal(x, **map_param(param_names, params))
    return np.sqrt(np.mean((y_pred - y_true)**2))


result_da = dual_annealing(loss, x0=initvel, bounds=bounds, args=(bins, count))
result = None

print('SciPy dual_annealing result:')
for i_param, param_name in enumerate(param_names):
    print(f'{param_name}: {result_da.x[i_param]:0.8f}', end='\t')
    if i_param % 3 == 2:
        print()

params, cov, infodict, mesg, ier = curve_fit(trimodal,
                                             bins,
                                             count,
                                             initvel,
                                             full_output=True,
                                             method='trf',
                                             bounds=bounds)

print('SciPy curve_fit result:')
for i_param, param_name in enumerate(param_names):
    print(f'{param_name}: {params[i_param]:0.8f}', end='\t')
    if i_param % 3 == 2:
        print()

# Prepare to plot the fitting result and other reference data
y_fit_cf = trimodal(bins, **map_param(param_names, params))
y_fit_da = trimodal(bins, **map_param(param_names, result_da.x))
y_fit_init = trimodal(bins, **map_param(param_names, initvel))

plt.plot(bins, count, label='Histogram')
plt.plot(bins, y_fit_init, 'k--', alpha=0.5, label='Initial parameters')
plt.plot(bins, y_fit_cf, label='minimize')
plt.plot(bins, y_fit_da, label='dual_annealing')
plt.legend()

plt.grid()
plt.show()


Result:

Untitled

This looks better to me than either of the results that curve_fit() was giving - it can recognize the peak at -22, which neither curve_fit() result found.

@lucascolley lucascolley added query A question or suggestion that requires further information and removed defect A clear bug or issue that prevents SciPy from being installed or used as expected labels May 10, 2024
@seongsujeong
Copy link
Author

@rkern @nickodell Thank you for the comments and opinions. Maybe I have to admit that would not be possible to completely avoid the small discrepancy coming from CPU architecture.

I've had a chance to test dual_annealing based on the sample code @nickodell has posted here. The same code has run on AMD/Linux, Intel/Linux, and M1/MacOS, and all of the three results are consistent.
Screenshot 2024-05-10 at 19 56 44

I suspected numerical error when numerically computing the gradient, so I've come up with providing analytical Jacobian function into curve_fit. Below is the fitting result.
Screenshot 2024-05-10 at 19 52 43

@rkern rkern closed this as not planned Won't fix, can't repro, duplicate, stale May 11, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
query A question or suggestion that requires further information scipy.optimize
Projects
None yet
Development

No branches or pull requests

5 participants