docs: def find_goodness_of_fit(rapid_qout_file, reach_id_file, observed_file,
                           out_analysis_file, daily=False):
      """
      Finds the goodness of fit comparing observed streamflow in a rapid Qout
      file with simulated flows in a csv file.
  
      Parameters
      ----------
      rapid_qout_file: str
          Path to the RAPID Qout file.
      reach_id_file: str
          ath to file with river reach ID's associate with the RAPID Qout file.
          It is in the format of the RAPID observed flows reach ID file.
      observed_file: str
          Path to input csv with with observed flows corresponding to the
          RAPID Qout. It is in the format of the RAPID observed flows file.
      out_analysis_file: str
          Path to the analysis output csv file.
      daily: bool, optional
          If True and the file is CF-Compliant, it will compare the
          *observed_file* with daily average flow from Qout. Default is False.
  
  
      Example with CF-Compliant RAPID Qout file:
  
      .. code:: python
  
          import os
          from RAPIDpy.postprocess import find_goodness_of_fit
  
          INPUT_DATA_PATH = '/path/to/data'
          reach_id_file = os.path.join(INPUT_DATA_PATH, 'obs_reach_id.csv')
          observed_file = os.path.join(INPUT_DATA_PATH, 'obs_flow.csv')
  
          cf_input_qout_file = os.path.join(COMPARE_DATA_PATH,
                                           'Qout_nasa_lis_3hr_20020830_CF.nc')
          cf_out_analysis_file = \
              os.path.join(OUTPUT_DATA_PATH,
                          'cf_goodness_of_fit_results-daily.csv')
          find_goodness_of_fit(cf_input_qout_file,
                               reach_id_file,
                               observed_file,
                               cf_out_analysis_file,
                               daily=True)
  
      """
      reach_id_list = np.loadtxt(reach_id_file,
                                 delimiter=",", usecols=(0,),
                                 ndmin=1, dtype=np.int32)
  
      data_nc = RAPIDDataset(rapid_qout_file)
  
      # analyze and write
      observed_table = np.loadtxt(observed_file,
                                  ndmin=2, delimiter=",",
                                  usecols=tuple(range(reach_id_list.size)))
      with open(out_analysis_file, 'w') as outcsv:
          writer = csvwriter(outcsv)
          writer.writerow(["reach_id",
                           "percent_bias",
                           "abs_percent_bias",
                           "rmse",
                           "mae",
                           "bias",
                           "NSE",
                           "likelihood",
                           "correlation_coeff",
                           "index_agreement",
                           "KGE"])
  
          for index, reach_id in enumerate(reach_id_list):
              observed_array = observed_table[:, index]
              simulated_array = data_nc.get_qout(reach_id, daily=daily)
              # make sure they are the same length
              simulated_array = simulated_array[:len(observed_array)]
              observed_array = observed_array[:len(simulated_array)]
              simulated_array, observed_array = \
                  filter_nan(simulated_array, observed_array)
              writer.writerow([reach_id,
                               pc_bias(simulated_array, observed_array),
                               apb(simulated_array, observed_array),
                               rmse(simulated_array, observed_array),
                               mae(simulated_array, observed_array),
                               bias(simulated_array, observed_array),
                               NS(simulated_array, observed_array),
                               L(simulated_array, observed_array),
                               correlation(simulated_array, observed_array),
                               index_agreement(simulated_array, observed_array),
                               KGE(simulated_array, observed_array)[0]])