<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 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]])