#   Title: z_b_CalibrationPEST_iTHydroPlus.py
# 	Description: This program is designed to prepare PEST inputs for HydroPlus.
#   Output: PEST files with extensions .pst, .tpl, and .ins files and a HydroPlusConfig.xml file
#   Input: HydroPlus_PEST_Commands.csv provides several rows of commands, each row progressing with a variable keyword, variable value, and perhaps other variables
#   Example Input:
#       row 1:HydroPlus_exe,C:\iTree\HydroPlus\HydroPlus\Release\HydroPlus.exe
#       row 2:HydroPlus_input,C:\iTree\HydroPlus\TestingFilesAndScript\TestCases\Calibration\Bioretention\input
#       row 3:SVDMODE,1
#       row 4:FLAG_LOG_MODE,0
#       row 5:PEST_date_YYYYMMDDHHMM,SimulationNumericalParams/StartDate_YYYYMMDD,201205081200
#       row 6:PEST_date_YYYYMMDDHHMM,SimulationNumericalParams/StopDate_YYYYMMDD,201205082355
#       row 7:PEST_parameter,DataOrganizer/DataDrawer/DataFolder/Parameter_m_KsatExpDecay,0.2,0.1,0.5,1
#       row 8:PEST_parameter,DataOrganizer/DataDrawer/DataFolder/Routing_PerviousArea_Alpha_h,0.04,0.01,1.0
#       row 9:PEST_observation_variable,Qobs
#       row 10:PEST_observation_file,C:\iTree\HydroPlus\TestingFilesAndScript\TestCases\Calibration\Bioretention\input\Runoff_Observed_for_Calibration.csv
#       row 11:HydroPlus_prediction,C:\iTree\HydroPlus\TestingFilesAndScript\TestCases\Calibration\Bioretention\output\Runoff_Predicted_for_Calibration.csv
#       row 12:Flag_PEST_observation_CSV,1
#       row 13:Flag_PEST_observation_TimeSeries,0
#       row 14:Flag_PEST_observation_ColumnID,1
#       row 15:Flag_PEST_prediction_ColumnID,2
#       row 16:Flag_PEST_If_ObservationIsZero_Then_WeightIsZero,1
#       row 17:Flag_PEST_observation_MapSeries,0
#       row 18:nrows,78
#       row 19:ncols,68
#   Run: Prepare PEST with Python 3+ at command line interface %>python z_b_CalibrationPEST_iTHydroPlus.py [path_to_HydroPlus_PEST_Commands.csv]
#   Run: Launch PEST with command line interface in directory input %>pest HydroPlus_PEST_ControlFile.pst or %>i64pest HydroPlus_PEST_ControlFile.pst 
#   Notes: 
#       0. PEST_observation_file data must be defined with flags below: CSV separated, TimeSeries, ColumnID
#       0a. If PEST_observation_file data are not a time series, a unique file is needed for each variable, named after the variable
#       0b. If PEST_observation_file data are time series, it is expected the 1st and 2nd column contain date and time YYYYMMDD,HH:MM:SS
#       1. If changing from use of Runoff_Predicted_for_Calibration.csv, then change structure in code with command PEST_01_ins_file.write
#       2. Calibration at weekly or monthly timestep requires preparing PEST_observation_file for that timestep;
#       3   . Refactor to: a) move within HydroPlus.sln to simplify management with changing calibration timestep; 
#           ... b) enable calibration of DataFolder variables outside of 1st DataFolder;
#           ... c) ensure multiple output files can be used in calibration;
#           ... d) 
#		4. Keyword PEST_date_YYYYMMDDHHMM is used to write observation data and represent StartDate_YYYYMMDD and StopDate_YYYYMMDD.
#		5. Keyword PEST_parameter is used to identify variables to adjust, and its initial value, minimum value, maximum value, and instance (1st or 2nd) in HydroPlusCong 
#          Note: HydroPlus_PEST_Commands.csv optionally indicates if variable is 1st or 2nd instance in HydroPlusConfig.xml
#          Note: HydroPlus_PEST_Commands.csv lists the variable twice if both BulkArea and GI are to be adjusted, giving instance 1 to BulkArea, instance 2 to GI
#       6. Keyword PEST_observation_variable will truncate all variable names to the first four (4) characters
#           ... a) It determines the number of observation groups and the number of PEST_observation_file and HydroPlus_prediction that are read into the .pst
#       7. Keyword Flag_PEST_observation_CSV, set to 1 for CSV files, 0 otherwise [but script is not ready for this case of 0]
#       8. Keyword Flag_PEST_observation_TimeSeries, set to 1 for time series, 0 for single observations, & overridden if Flag_PEST_observation_MapSeries = 1
#       9. Keyword Flag_PEST_observation_ColumnID, set to column with observations starting with ID=0 for first column
#       10. Keyword Flag_PEST_prediction_ColumnID, set to column with predictions starting with ID=0 for first column
#       11. Keyword Flag_PEST_If_ObservationIsZero_Then_WeightIsZero, set to 1 to give 0 weight to observations with 0 value
#       12. Keyword Flag_PEST_observation_MapSeries, set to 1 for map series, 0 for time series
#       13. Keywords nrows and ncols, set to map extent when Flag_PEST_observation_MapSeries = 1
#       14. Python script below presumes only one model output file is read for comparison with observations, using NINSFLE = 1; This could be changed to NINSFLE = NOBSGP 
#
#   Test with PEST tools >inschek HydroPlus_PEST_Output_0.ins C:\iTree\HydroPlus\TestingFilesAndScript\TestCases\Calibration\WaterTable2Obs\input
#
#   Created: T. Endreny, 2022, te@esf.edu
#
# ========================================================================================================================================================

#import Python libraries
import subprocess
import os
import pandas as pd
import time
from random import uniform
import numpy as np
import csv 
import shutil
import glob
import sys
from filecmp import dircmp
import lxml.etree as LXML_File
from datetime import datetime
from itertools import islice

arguments_count = len(sys.argv)
if arguments_count < 2:
    print("Warning: Command line argument inputs = {}, and appear to be missing the path to inputs.".format(arguments_count))
    print("Exiting: This warning triggers the program to exit the script.")
    print("Correction: To properly run this program, provide two arguments after the python command ...")
    print("%>python z_b_CalibrationPEST_iTHydroPlus.py [path_to_HydroPlus_PEST_Commands.csv]")
    sys.exit()

print("\n")
print("Notice: The last line of data in HydroPlus_PEST_Commands.csv is read only if it is followed by an empty line.")

#Define PEST_input_file_str as input HydroPlus_PEST_Commands.csv
PEST_input_file_str = "HydroPlus_PEST_Commands.csv"
#Define HydroPlus.exe path
HydroPlus_exe = "HydroPlus_exe"

#Define Project input folder path
HydroPlus_input = "HydroPlus_input"
#Define HydroPlusConfigXML_input_str as HydroPlusConfig.xml
HydroPlusConfigXML_input_str = "HydroPlusConfig.xml"
#Define PEST date string keyword used to identify calibration dates within PEST_input_file_str
PEST_date_YYYYMMDDHHMM = "PEST_date_YYYYMMDDHHMM"
#Define XML_StartDate  using HydroPlusConfig.xml element key
XML_StartDate = "StartDate_YYYYMMDD"
#Define XML_StopDate  using HydroPlusConfig.xml element key
XML_StopDate = "StopDate_YYYYMMDD"
#Define PEST parameter string keyword PEST_parameter used to identify calibration input parameters
PEST_parameter = "PEST_parameter"
#Define PEST PEST_observation_variable string keyword PEST_observation_variable used to identify calibration observed data
PEST_observation_variable = "PEST_observation_variable"
#Define PEST PEST_observation_file string keyword PEST_observation_file used to identify calibration observed file path
PEST_observation_file = "PEST_observation_file"
#Define HydroPlus_prediction string keyword HydroPlus_prediction used to identify calibration predicted data
HydroPlus_prediction = "HydroPlus_prediction"
#Define calibration_symbol for use in .tpl file
PEST_key_symbol = "#"
#Define PEST_HydroPlus_tpl_str as HydroPlusConfigXML_input_str.tpl
PEST_tpl_str = HydroPlusConfigXML_input_str + ".tpl"
#Define PEST_01_ins_str as PEST_Output_01.ins
PEST_01_ins_str = "HydroPlus_PEST_Output_01.ins"
#Define PEST_pst_str as HydroPlus_PEST_ControlFile.pst
PEST_pst_str = "HydroPlus_PEST_ControlFile.pst"

#Set path to PEST management folder
#project_path defined as 1st command line argument, sys.argv[0]
project_path = sys.argv[1]
#determine if last character of project_path has backslash
last_character_argv00 = project_path[-1]
#If last_character_argv00 is not \ then
#Note: "\\" is used due to backslash also serving as an escape character
if (last_character_argv00 != "\\"):
    #add backslash to path to working directory
    project_path = project_path + "\\"

#PEST_inputs_path is the project_path + PEST_input_file_str
#Note: PEST_input_file_str is the HydroPlus_PEST_Commands.csv
PEST_inputs_path = project_path + PEST_input_file_str

#Define variables
#Variables from HydroPlusConfigXML_input_str
variable_names = []
variable_values = []
calibration_date = []
variable_initial = []
variable_min = []
variable_max = []
variable_instance = []
#Variables from Observed_data_01_str
Observed_data_path = []
HydroPlus_prediction_path = []
Observed_variable_full_name = []
NOBSPGP = []
OBSNME = []
PARNME = []
observation_date = []
observation_time = []
observation_value = []
#Define XML_Key_Elements_String as array
XML_Key_Elements_String = []
XML_Key_Elements_Date = []
XML_Key_Elements_PEST_Parameter = []
#Define XML_Key_Elements_Replacement_String as array
XML_Key_Elements_Replacement_String = []
XML_Key_Elements_Replacement_Date = []
XML_Key_Elements_Replacement_PEST_Parameter = []
#Define XML_File_Tree_Element_Vector to hold multiple instances
XML_File_Tree_Element_Vector = []

#Flag_PEST_observation_CSV and similar flags set to default values of zero
Flag_PEST_observation_CSV = 0
Flag_PEST_observation_TimeSeries = 0
Flag_PEST_observation_ColumnID = 0
Flag_PEST_prediction_ColumnID = 0
Flag_PEST_observation_MapSeries = 0
Flag_PEST_If_ObservationIsZero_Then_WeightIsZero = 0
#pst_observation_start_index initiated to -1 to serve as flag if start date is not found in data
pst_observation_start_index = -1
#pst_observation_stop_index initiated to -1 to serve as flag if stop date is not found in data
pst_observation_stop_index = -1

###################################################################
#Algorithm to read in data from PEST_input_file
###################################################################
#Initialize NPAR, number of parameters to 0
NPAR = 0
#Initialize NOBSGP, number of observation group data sets to 0
NOBSGP = 0
#open PEST_inputs_path as input file as variable PEST_input_file
#Note: PEST_inputs_path is the HydroPlus_PEST_Commands.csv
with open(PEST_inputs_path) as PEST_inputs_file:
    #get each rows of input, separate rows vertically by split of line break
    #Note: If HydroPlus_PEST_Commands.csv does not have an empty line at end of file, last line is not read
    lines = PEST_inputs_file.read().split('\n')
    #initiate row_counter
    row_counter = 0  
    #for each of the rows of input, starting at 1st line at position 0
    for line in lines[0:-1]:
        #line is defined without white space using line.strip()
        line = line.strip()
        #if line is a conditional checking if there are any elements in the line
        #Note: this helps avoid reading in empty lines at the end of a file
        if line:
            #assign input to input_data, separate the row horizontally by split of comma
            input_data = line.split(",")
            #Append to variable_names the 1st element in row
            variable_names.append(input_data[0])
            #Append to variable_values the 2nd element in row 
            variable_values.append(input_data[1])
            #If substring of variable_names[row_counter] equals HydroPlus_exe then
            if variable_names[row_counter] == HydroPlus_exe:
                #HydroPlus_exe_path equals variable_values[row_counter], which contains path
                HydroPlus_exe_path = variable_values[row_counter]
            #If substring of variable_names[row_counter] equals HydroPlus_input then
            if variable_names[row_counter] == HydroPlus_input:
                #HydroPlus_exe_path equals variable_values[row_counter], which contains path
                HydroPlus_input_path = variable_values[row_counter]
                #determine if last character of HydroPlus_input_path has backslash
                last_character_argv00 = HydroPlus_input_path[-1]
                #If last_character_argv00 is not \ then
                #Note: "\\" is used due to backslash also serving as an escape character
                if (last_character_argv00 != "\\"):
                    #add backslash to path to working directory
                    HydroPlus_input_path = HydroPlus_input_path + "\\"
            #If substring of variable_names[row_counter] equals FLAG_LOG_MODE then
            if variable_names[row_counter] == "FLAG_LOG_MODE":
                #FLAG_LOG_MODE set equal to variable_values[row_counter]
                #FLAG_LOG_MODE is choice on variable transform mode from interface, 1 = log mode, 0 = no transform
                FLAG_LOG_MODE = variable_values[row_counter]
            #If substring of variable_names[row_counter] equals SVDMODE then
            if variable_names[row_counter] == "SVDMODE":
                #SVDMODE set equal to variable_values[row_counter]
                #SVDMODE flag for use of singular value decomposition, 0 is inactive, 1 is active and assured stability, 2 is active but less stable
                SVDMODE = variable_values[row_counter]
            #If substring of variable_names[row_counter] equals PEST_date_YYYYMMDDHHMM then
            if variable_names[row_counter] == PEST_date_YYYYMMDDHHMM:
                #Append to calibration_date the 3rd element in row 
                calibration_date.append(input_data[2])
            #If substring of variable_names[row_counter] equals PEST_parameter then
            #Note: PEST_parameter should be followed by: 1) variable path in HydroPlusConfig.xml, ...
            #Note: ... 2) initial value, 3) minimum value, 4) maximum value, and optionally ... 
            #Note: ... 5) variable variable, 1 = 1st instance, e.g., Bulk Area or 2 = 2nd instance, e.g., GI
            if variable_names[row_counter] == PEST_parameter:
                #NPAR is number of parameters
                #NPAR is increased by 1
                NPAR = int(NPAR) + 1
                #Append to variable_initial the 3rd element in row 
                variable_initial.append(input_data[2])
                #Append to variable_min the 4th element in row 
                variable_min.append(input_data[3])
                #Append to variable_max the 5th element in row 
                variable_max.append(input_data[4])
                #Try to find and append optional data on instance
                try:
                    #If len(input_data) > 5 then the 6th and optional element was provided
                    if len(input_data) > 5:
                        # Check if the length of input_data[5] is greater than 0
                        #If len(input_data[5] > 0 then optional instance identified and append instance 
                        if len(input_data[5]) > 0:
                            #If input_data[5] value > 2 then change to 2
                            if (int(input_data[5]) > 2): 
                                input_data[5] = 2
                            #If input_data[5] value < 1 then change to 1
                            if (int(input_data[5]) < 1): 
                                input_data[5] = 1
                            #variable_instance appends the value provided, either 1 or 2
                            variable_instance.append(input_data[5])
                    #Else len(input_data) < 5 then the 6th and optional element was not provided
                    else:
                        #variable_instance appends the value of 1, indicating 1st instance, which is assumed
                        variable_instance.append(1)
                #Except IndexError when the input_data doesn't have a valid index, such as fewer than 6 elements
                except IndexError:
                    #variable_instance appends the value of 1, indicating 1st instance, which is assumed
                    variable_instance.append(1)
            #If substring of variable_names[row_counter] equals PEST_observation_variable then
            if variable_names[row_counter] == PEST_observation_variable:
                #NOBSGP is increased by 1
                NOBSGP = int(NOBSGP) + 1
                #Observed_data_str set equal to variable_values[row_counter]
                Observed_variable_AllCharacter_str = variable_values[row_counter]
                Observed_variable_4Character_str = variable_values[row_counter][0:4]
                OBSNME.append(Observed_variable_4Character_str)
                #Observed_variable_full_name list is appended with Observed_variable_AllCharacter_str, which has full variable name
                #Note: Observed_variable_full_name only needed when searching primary marker in Summary_Budget_Water_mm.csv
                Observed_variable_full_name.append(Observed_variable_AllCharacter_str)
            #If substring of variable_names[row_counter] equals PEST_observation_file then
            if variable_names[row_counter] == PEST_observation_file:
                #Observed_data_file_path set equal to variable_values[row_counter]
                Observed_data_file_path = variable_values[row_counter]
                #Observed_data_path appends HydroPlus_input_path 
                Observed_data_path.append(Observed_data_file_path) 
            #If substring of variable_names[row_counter] equals HydroPlus_prediction then
            if variable_names[row_counter] == HydroPlus_prediction:
                #Predicted_data_file_path set equal to variable_values[row_counter]
                Predicted_data_file_path = variable_values[row_counter]
                #HydroPlus_prediction_path appends Predicted_data_file_path
                HydroPlus_prediction_path.append(Predicted_data_file_path) 
            #If substring of variable_names[row_counter] equals Flag_PEST_observation_CSV then
            if variable_names[row_counter] == "Flag_PEST_observation_CSV":
                #PEST_observation_file flags: CSV separated, TimeSeries, ColumnID
                #Change: Flag_PEST_observation_CSV (1 = CSV) identifies the line separation format for the data used in calibration
                Flag_PEST_observation_CSV = int(variable_values[row_counter])
            #If substring of variable_names[row_counter] equals Flag_PEST_observation_TimeSeries then
            if variable_names[row_counter] == "Flag_PEST_observation_TimeSeries":
                #Change: Flag_PEST_observation_TimeSeries (0 = no TimeSeries; 1 = TimeSeries) identifies whether data used in calibration are time series
                Flag_PEST_observation_TimeSeries = int(variable_values[row_counter])
                if Flag_PEST_observation_TimeSeries > 0:
                    #Flag_PEST_observation_MapSeries set to default 0 in case not with inputs
                    Flag_PEST_observation_MapSeries = 0
            #If substring of variable_names[row_counter] equals Flag_PEST_observation_ColumnID then
            if variable_names[row_counter] == "Flag_PEST_observation_ColumnID":
                #Change: Flag_PEST_observation_ColumnID (starts at 0) identifies the column with the data used in calibration
                Flag_PEST_observation_ColumnID = int(variable_values[row_counter])
            #If substring of variable_names[row_counter] equals Flag_PEST_prediction_ColumnID then
            if variable_names[row_counter] == "Flag_PEST_prediction_ColumnID":
                #Change: Flag_PEST_prediction_ColumnID (starts at 0) identifies the column with the data used in calibration
                Flag_PEST_prediction_ColumnID = int(variable_values[row_counter])
            #If substring of variable_names[row_counter] equals Flag_PEST_If_ObservationIsZero_Then_WeightIsZero then
            if variable_names[row_counter] == "Flag_PEST_If_ObservationIsZero_Then_WeightIsZero":
                #Change: Flag_PEST_If_ObservationIsZero_Then_WeightIsZero = 1 assigns weight of 0 to observations of 0
                Flag_PEST_If_ObservationIsZero_Then_WeightIsZero = int(variable_values[row_counter])
            #If substring of variable_names[row_counter] equals Flag_PEST_observation_MapSeries then
            if variable_names[row_counter] == "Flag_PEST_observation_MapSeries":
                #Change: Flag_PEST_observation_MapSeries (0 = no MapSeries; 1 = MapSeries) identifies whether data used in calibration are map series
                Flag_PEST_observation_MapSeries = int(variable_values[row_counter])
                if Flag_PEST_observation_MapSeries > 0:
                    #Flag_PEST_observation_MapSeries set to default 0 in case not with inputs
                    Flag_PEST_observation_TimeSeries = 0

            #advance row_counter
            row_counter = row_counter + 1

#Set path to XML_Input_and_Template_path as HydroPlus_input_path + HydroPlusConfigXML_input_str
XML_Input_and_Template_path = HydroPlus_input_path + HydroPlusConfigXML_input_str
#Create XML_File_Tree to hold all HydroPlusConfig.xml elements by using XML_File.parse with XML_Input_and_Template_path
XML_File_Tree = LXML_File.parse(XML_Input_and_Template_path)
#Set path to PEST_tpl_file_path as HydroPlus_input_path + PEST_tpl_str
PEST_tpl_file_path = HydroPlus_input_path + PEST_tpl_str
#Set path to PEST_ins_file_path as HydroPlus_input_path + PEST_01_ins_str
PEST_ins_file_path = HydroPlus_input_path + PEST_01_ins_str
#Set path to PEST_pst_file_path as HydroPlus_input_path + PEST_pst_str
PEST_pst_file_path = HydroPlus_input_path + PEST_pst_str

###################################################################
#Algorithm to extract calibration data from and place parameters in XML_File_Tree
###################################################################

# Function to generate unique variable name
def generate_unique_variable_name(instance_index, PEST_parameter_counter):

    return f"variable_name_{instance_index}_{PEST_parameter_counter}"
    
#initialize row_counter
row_counter = 0
#initialize PEST_parameter_counter
PEST_parameter_counter = 0
#initialize PEST_date_counter
PEST_date_counter = 0
#for each row of data in variable_names list
#Note: variable_names list contains all of the variables that will be calibrated
for PEST_input_row in variable_names:
    #Algorithm to extract PEST_date_YYYYMMDDHHMM from XML file
    #if PEST_input_row equals PEST_date_YYYYMMDDHHMM then
    if PEST_input_row == PEST_date_YYYYMMDDHHMM:
        #location_variable_len defined using int and rfind functions with variable_values[row_counter] and character /, finding length to last instance
        location_variable_len = int(variable_values[row_counter].rfind("/"))
        #root_len defined using len function with variable_values[row_counter]
        root_len = len(variable_values[row_counter])
        #date_name_str is defined as substring of variable_values[row_counter], to extract variable name after character /
        date_name_str = variable_values[row_counter][int(location_variable_len+1):root_len]
        #XML_Key_Elements_Date appends variable_values[row_counter], containing 
        XML_Key_Elements_Date.append(variable_values[row_counter])
        #XML_Key_Elements_Replacement_Date appends calibration_date[PEST_date_counter][0:10] to get YYYYMMDDHH
        XML_Key_Elements_Replacement_Date.append(calibration_date[PEST_date_counter][0:10])
        #Define XML_File_Tree_Element the value XML_Key_Elements_Date[PEST_date_counter]
        #Note: XML_File_Tree.findall function an option, and then search vector of output for all instances, or last instance, etc. 
        #Note: XML_File_Tree_Element is known to XML_File_Tree as a member, and can be altered and known by XML_File_Tree
        XML_File_Tree_Element = XML_File_Tree.find(XML_Key_Elements_Date[PEST_date_counter])
        #If date_name_str equals XML_StartDate then get value
        if (date_name_str == XML_StartDate):
            #XML_StartDate_GD equals XML_Key_Elements_Replacement_Date[PEST_date_counter], value going into HydroPlusConfig.xml element key
            XML_StartDate_GD = calibration_date[PEST_date_counter]
        #If date_name_str equals XML_StopDate then get value
        if (date_name_str == XML_StopDate):
            #XML_StopDate_GD equals XML_Key_Elements_Replacement_Date[PEST_date_counter], value going into HydroPlusConfig.xml element key
            XML_StopDate_GD = calibration_date[PEST_date_counter]
        #Redefine XML_File_Tree_Element to new value, XML_Key_Elements_Replacement_Date, within HydroPlusConfg.xml file
        #Note: XML_File_Tree_Element is known to XML_File_Tree as a member, and can be altered and known by XML_File_Tree
        XML_File_Tree_Element.text = str(XML_Key_Elements_Replacement_Date[PEST_date_counter])
        #advance PEST_date_counter
        PEST_date_counter = PEST_date_counter + 1

    #Algorithm to place PEST_parameter variable names into XML file
    #if PEST_input_row equals PEST_parameter then
    if PEST_input_row == PEST_parameter:
        #rfind to locate last instance of character / in string variable_values
        location_variable_len = int(variable_values[row_counter].rfind("/"))
        root_len = len(variable_values[row_counter])
        #variable_name_str is defined as substring of variable_values[row_counter], to extract variable name after character /
        variable_name_str = variable_values[row_counter][int(location_variable_len+1):root_len]

        #Algorithm to limit PEST parameter variable names to 12 characters within the .tpl file
        #var_spacebar_max is set to 12, accounting for 2 more spaces allocated to bookends with PEST_key_symbol
        var_spacebar_max = 12
        #If len(variable_name_str) > var_spacebar_max then
        if (len(variable_name_str) > int(var_spacebar_max)):
            #If PEST_parameter_counter < 10 then
            if PEST_parameter_counter < 10:
                #variable_name_str equals substring of variable_name_str[0:12] appended to str(PEST_parameter_counter)
                variable_name_str = variable_name_str[0:11] + str(PEST_parameter_counter)
            #If PEST_parameter_counter < 10 then
            elif PEST_parameter_counter < 100:
                #variable_name_str equals substring of variable_name_str[0:11] appended to str(PEST_parameter_counter)
                variable_name_str = variable_name_str[0:10] + str(PEST_parameter_counter)
        #var_spacebar has default length of zero
        var_spacebar = ""
        #var_spacebvar_spacebarar_loop equals var_spacebar_max minus len(variable_name_str)
        var_spacebar_loop = var_spacebar_max - len(variable_name_str)
        #If var_spacebar_loop > 0 then
        if int(var_spacebar_loop) > 0:
            #for i_space in range from 0 to var_spacebar_loop
            for i_space in range(0,var_spacebar_loop):
                #var_spacebar increased by a space
                var_spacebar = var_spacebar + " "
                
        #PARNME appends variable_name_str; this name is written in .tpl and .pst files
        PARNME.append(variable_name_str)
        #XML_Key_Elements_PEST_Parameter appends variable_values[row_counter]
        XML_Key_Elements_PEST_Parameter.append(variable_values[row_counter])
        #XML_Key_Elements_Replacement_PEST_Parameter appends PEST_key_symbol + variable_name_str + var_spacebar + PEST_key_symbol
        XML_Key_Elements_Replacement_PEST_Parameter.append(PEST_key_symbol + variable_name_str + var_spacebar + PEST_key_symbol)

        #Define XML_File_Tree_Element the value XML_Key_Elements_PEST_Parameter[PEST_parameter_counter]
        #Note: XML_File_Tree.findall function when >1 instance of element in HydroPlusConfig.xml due to BulkArea and GI, creating a vector of similar names
        #Note: XML_File_Tree.find function when only 1 instance or only want 1st instance of element in HydroPlusConfig.xml
        #Note: XML_File_Tree_Element_Vector[], is known to XML_File_Tree as a member, and can be altered and known by XML_File_Tree
        XML_File_Tree_Element_Vector = XML_File_Tree.findall(XML_Key_Elements_PEST_Parameter[PEST_parameter_counter])

        #If int(variable_instance[PEST_parameter_counter]) equals 1 then
        #Note: To modify 1st and 2nd instance of a variable, just list it twice in HydroPlus_PEST_Commands.csv ...
        #Note: ... and for east listing, provide either 1 or 2 in the 6th position, after the maximum value
        if (int(variable_instance[PEST_parameter_counter]) == 1):
            #XML_File_Tree_Element equals XML_File_Tree_Element_Vector[0], where 0 indicates first instance
            XML_File_Tree_Element = XML_File_Tree_Element_Vector[0]
            #Define variable_name_value using XML_File_Tree_Element.text 
            variable_name_value = XML_File_Tree_Element.text 
            #Redefine XML_File_Tree_Element to new value, XML_Key_Elements_Replacement_PEST_Parameter, within HydroPlusConfg.xml file
            #Note: XML_File_Tree_Element is known to XML_File_Tree as a member, and can be altered and known by XML_File_Tree
            XML_File_Tree_Element.text = str(XML_Key_Elements_Replacement_PEST_Parameter[PEST_parameter_counter])
        #Else If int(variable_instance[PEST_parameter_counter]) == 2 and XML_File_Tree_Element_Vector contains 2 members then
        elif (int(variable_instance[PEST_parameter_counter]) == 2 and len(XML_File_Tree_Element_Vector) > 1):
            #XML_File_Tree_Element equals XML_File_Tree_Element_Vector[1], where 1 indicates second or GI instance
            XML_File_Tree_Element = XML_File_Tree_Element_Vector[1]
            #Define variable_name_value using XML_File_Tree_Element.text 
            variable_name_value = XML_File_Tree_Element.text 
            #Redefine XML_File_Tree_Element to new value, XML_Key_Elements_Replacement_PEST_Parameter, within HydroPlusConfg.xml file
            #Note: XML_File_Tree_Element is known to XML_File_Tree as a member, and can be altered and known by XML_File_Tree
            XML_File_Tree_Element.text = str(XML_Key_Elements_Replacement_PEST_Parameter[PEST_parameter_counter])
        #Else int(variable_instance[PEST_parameter_counter]) == 2 but len(XML_File_Tree_Element_Vector) < 1
        else:
            #XML_File_Tree_Element equals XML_File_Tree_Element_Vector[0], where 0 indicates first instance
            XML_File_Tree_Element = XML_File_Tree_Element_Vector[0]
            #Define variable_name_value using XML_File_Tree_Element.text 
            variable_name_value = XML_File_Tree_Element.text 
            #Redefine XML_File_Tree_Element to new value, XML_Key_Elements_Replacement_PEST_Parameter, within HydroPlusConfg.xml file
            #Note: XML_File_Tree_Element is known to XML_File_Tree as a member, and can be altered and known by XML_File_Tree
            XML_File_Tree_Element.text = str(XML_Key_Elements_Replacement_PEST_Parameter[PEST_parameter_counter])

        #advance PEST_parameter_counter
        PEST_parameter_counter = PEST_parameter_counter + 1
    #advance row_counter
    row_counter = row_counter + 1

#XML_StartDate_YYYYMMDD is defined as substring of XML_StartDate_GD[0,8], taking first 8 characters
XML_StartDate_YYYYMMDD = XML_StartDate_GD[0:8]
#XML_StopDate_YYYYMMDD is defined as substring of XML_StopDate_GD[0,8], taking first 8 characters
XML_StopDate_YYYYMMDD = XML_StopDate_GD[0:8]
#XML_StartDate_HHMM is defined as substring of XML_StartDate_GD[8,12], taking last 4 characters
XML_StartDate_HHMM = XML_StartDate_GD[8:12]
#XML_StopDate_HHMM is defined as substring of XML_StopDate_GD[8,12], taking last 4 characters
XML_StopDate_HHMM = XML_StopDate_GD[8:12]

#If XML_StartDate_HHMM is not entered, set it to zeros
if (len(XML_StartDate_HHMM) <=0):
    XML_StartDate_HHMM = "0000"
#If XML_StopDate_HHMM is not entered, set it to zeros
if (len(XML_StopDate_HHMM) <=0):
    XML_StopDate_HHMM = "0000"
#If XML_StartDate_HHMM is 2 digits long, add MM zeros
if (len(XML_StartDate_HHMM) <3):
    XML_StartDate_HHMM = XML_StartDate_HHMM + "00"
#If XML_StopDate_HHMM is 2 digits long, add MM zeros
if (len(XML_StopDate_HHMM) <3):
    XML_StopDate_HHMM = XML_StopDate_HHMM + "00"

#Obtain YYYYMMDD values from XML_StartDate_GD and XML_StopDate_GD
print("XML_StartDate_GD = {}; XML_StopDate_GD = {}".format(XML_StartDate_GD, XML_StopDate_GD))
print("XML_StartDate_YYYYMMDD = {}; XML_StopDate_YYYYMMDD = {}".format(XML_StartDate_YYYYMMDD, XML_StopDate_YYYYMMDD))
print("XML_StartDate_HHMM = {}; XML_StopDate_HHMM = {}".format(XML_StartDate_HHMM, XML_StopDate_HHMM))

###################################################################
#Algorithm to extract observed data from Observed_data_file
###################################################################
print("Flag_PEST_observation_CSV = {}".format(Flag_PEST_observation_CSV))
print("Flag_PEST_observation_TimeSeries = {}".format(Flag_PEST_observation_TimeSeries))
print("Flag_PEST_observation_ColumnID = {}".format(Flag_PEST_observation_ColumnID))
print("Flag_PEST_prediction_ColumnID = {}".format(Flag_PEST_prediction_ColumnID))
print("Flag_PEST_observation_MapSeries = {}".format(Flag_PEST_observation_MapSeries))

#def format_time_to_HHMM function takes time_str and uses datetime.strptime function to get the H M S objects
def format_time_to_HHMM(time_str):
    #time_obj is the datetime object result of parsing the time_str with the datetime.strptime function
    time_obj = datetime.strptime(time_str, "%H:%M:%S")
    #time_HHMM_str is the formatted datetime object in the desired string format "HHMM"
    time_HHMM_str = time_obj.strftime("%H%M")
    #return time_HHMM_str
    return time_HHMM_str

#for loop with OBSGP in range from 0 to NOBSGP
for OBSGP in range(0, NOBSGP):
    #Initialize Flag_Advance_counter to 0
    Flag_Advance_counter = 0
    #If Flag_PEST_observation_MapSeries is 1, then
    if (Flag_PEST_observation_MapSeries == 1):
        #row_start is 0, reading from 1st line of header
        row_start = 0
    else:
        #row_start is 0, skipping 1st line which is header
        row_start = 1
    #open observed data file as variable PEST_inputs_file
    with open(Observed_data_path[OBSGP]) as Observed_data_file:

        #Initialize Observation_counter to zero
        Observation_counter = 0
        #get each rows of input, separate rows vertically by split of line break
        lines = Observed_data_file.read().split('\n')

        #initiate row_counter
        row_counter = 0
        #for each of the rows of input, starting at 1st line at position 1, skipping header
        #Note: lines[-1:] will read to end of file, but use [1:-1] to read until penultimate item in file
        for line in lines[row_start:]:
            #Note: This allows files to have extra line after all data fields
            #Note: lines[-1:] will read to end of file, but use [1:-1] to read until penultimate item in file
            if line.strip():  
                #assign input to observed_01_data, separate the row horizontally by split of comma
                #Note: Use Flag_PEST_observation_CSV == 1 to use line.split(",")
                if (Flag_PEST_observation_CSV == 1):
                    observed_01_data = line.split(",")
                #Note: Use Flag_PEST_observation_CSV == 0 to use line.split()
                if (Flag_PEST_observation_CSV != 1):
                    observed_01_data = line.split()
                    
                #If conditional based on flags is Flag_PEST_observation_TimeSeries == 1 and Flag_PEST_observation_ColumnID == 2 and Flag_PEST_observation_MapSeries != 1
                #Note: Flag_PEST_observation_ColumnID starts at 0, so 2 is third column
                #Note: If len(observed_01_data) <= 0 then line in lines[1:] contained an empty row
                if (Flag_PEST_observation_TimeSeries == 1 and Flag_PEST_observation_ColumnID > 1 and len(observed_01_data) > 0 and Flag_PEST_observation_MapSeries != 1):
                    #Append to variable_names the 1st element in row 
                    #Note: Values are in YYYYMMDD
                    observation_date.append(observed_01_data[0])
                    #observation_time appends values returned from the format_time_to_HHMM function, sending observed_01_data[1]
                    observation_time.append(format_time_to_HHMM(observed_01_data[1]))
                    #Check that observed_01_data elements equal or exceed Flag_PEST_observation_ColumnID
                    num_columns = len(observed_01_data)
                    num_commas = num_columns - 1
                    if (Flag_PEST_observation_ColumnID > num_commas):
                        print("The observed data file has too few columns and commas compared with the Flag_PEST_observation_ColumnID")
                        # Terminate the script
                        sys.exit()                       
                    #Append to variable_values the Flag_PEST_observation_ColumnID element in row
                    #Note: Flag_PEST_observation_ColumnID is presumed to be separated by commas
                    observation_value.append(observed_01_data[Flag_PEST_observation_ColumnID])
                    #If observation_date equals XML_StartDate_YYYYMMDD and observation_time equals XML_StartDate_HHMM then 
                    if (observation_date[row_counter] == XML_StartDate_YYYYMMDD and observation_time[row_counter] == XML_StartDate_HHMM):
                        print("observervation_start_date = XML_StartDate_YYYYMMDD, {} and observation_time = XML_StartDate_HHMM, {}".format(XML_StartDate_YYYYMMDD,XML_StartDate_HHMM))
                        #pst_observation_start_index = row_counter when XML_StartDate_YYYYMMDD and XML_StartDate_HHMM identified
                        pst_observation_start_index = row_counter
                        #Flag_Advance_counter changed to 1
                        Flag_Advance_counter = 1
                    #If Flag_Advance_counter not zero then advance counter
                    if (Flag_Advance_counter > 0):
                        #Observation_counter increased by 1
                        Observation_counter = int(Observation_counter) + 1
                    #If observation_date equals XML_StopDate_YYYYMMDD and observation_time equals XML_StopDate_HHMM then 
                    if (observation_date[row_counter] == XML_StopDate_YYYYMMDD and observation_time[row_counter] == XML_StopDate_HHMM):
                        print("observervation_stop_date = XML_StopDate_YYYYMMDD, {} and observation_time = XML_StopDate_HHMM, {}".format(XML_StopDate_YYYYMMDD,XML_StopDate_HHMM))
                        #pst_observation_stop_index = row_counter when XML_StartDate_YYYYMMDD and XML_StartDate_HHMM identified
                        pst_observation_stop_index = row_counter
                        #Flag_Advance_counter changed to 0
                        Flag_Advance_counter = 0

                    #advance row_counter
                    row_counter = row_counter + 1
                    
                #If conditional based on flags is Flag_PEST_observation_TimeSeries == 0 and Flag_PEST_observation_ColumnID == 1 and Flag_PEST_observation_MapSeries != 1
                #Note: Flag_PEST_observation_ColumnID starts at 0, so 2 is third column
                #Note: If len(observed_01_data) > 2 then line in lines[1:] has two or more commas and has left the header
                elif (Flag_PEST_observation_TimeSeries == 0 and Flag_PEST_observation_ColumnID == 1 and len(observed_01_data) >= 3 and Flag_PEST_observation_MapSeries != 1):
                    #Append to variable_values the 3rd element in row
                    #Note: Values are in m/h
                    observation_value.append(observed_01_data[1])
                    #pst_observation_start_index = row_counter, which should be one given Flag_PEST_observation_TimeSeries = 0
                    pst_observation_start_index = row_counter
                    #pst_observation_stop_index = row_counter, which should be one given Flag_PEST_observation_TimeSeries = 0
                    pst_observation_stop_index = row_counter
                    #advance row_counter
                    row_counter = row_counter + 1
                    #Observation_counter increased by 1
                    Observation_counter = 1
                #If conditional based on flags is Flag_PEST_observation_MapSeries == 1 and Flag_PEST_observation_TimeSeries != 1
                #Note: Flag_PEST_observation_MapSeries == 1 presumes reading ArcGIS ASCII map with 6 header rows
                elif (Flag_PEST_observation_MapSeries == 1 and Flag_PEST_observation_TimeSeries != 1):
                    #If row_counter equals 0, ncols in header
                    if (row_counter == 0):
                        #read number of cols as 2nd element in 2nd line of header
                        ncols = int(observed_01_data[1])
                        print(f"OBSGP={OBSGP}: ncols={ncols}")
                    #If row_counter equals 1, nrows of header
                    if (row_counter == 1):
                        #read number of rows as 2nd element in 1st line of header
                        nrows = int(observed_01_data[1])
                        print(f"OBSGP={OBSGP}: nrows={nrows}")
                    #If row_counter > 5 then read data below header rows; 0-5 are header rows
                    if (row_counter > 5):
                        #observation_value = [list(map(float, line.split())) for line in islice(Observed_data_file, nrows)]
                        observation_value.extend(observed_01_data)
                        #Observation_counter = int(Observation_counter) + len(observed_01_data) ensures that Observation_counter increase by length of each line
                        Observation_counter = int(Observation_counter) + len(observed_01_data)
                    #advance row_counter
                    row_counter = row_counter + 1
                #Else, print error about invalid input
                else: print("Invalid input. Flag_PEST_observation_TimeSeries={} and Flag_PEST_observation_MapSeries={} cannot both be 1 or 0".format(Flag_PEST_observation_TimeSeries,Flag_PEST_observation_MapSeries))
                
    #If pst_observation_start_index = -1 and Flag_PEST_observation_TimeSeries == 1 and Flag_PEST_observation_ColumnID > 1
    #Note: If pst_observation_start_index is empty, then there is date disagreement between HydroPlusConfig.xml and observed time series
    if (pst_observation_start_index == -1 and Flag_PEST_observation_TimeSeries == 1 and Flag_PEST_observation_ColumnID > 1):
        #Explain the problem
        print("\n")
        print("Warning: XML_StartDate_YYYYMMDD is not contained within dates of PEST_observation_file.")
        print("Solution: Modify the XML_StartDate_YYYYMMDD variable, or the observation date and time.")
        print("PEST_observation_file: ... {}".format(Observed_data_path[OBSGP]))
        #Terminate the script
        os._exit(0)

    #If pst_observation_stop_index = -1 and Flag_PEST_observation_TimeSeries == 1 and Flag_PEST_observation_ColumnID > 1
    #Note: If pst_observation_stop_index is empty, then there is date disagreement between HydroPlusConfig.xml and observed time series
    if (pst_observation_stop_index == -1 and Flag_PEST_observation_TimeSeries == 1 and Flag_PEST_observation_ColumnID > 1):
        #Explain the problem
        print("\n")
        print("Warning: XML_StopDate_YYYYMMDD is not contained within dates of PEST_observation_file.")
        print("Solution: Modify the XML_StopDate_YYYYMMDD variable, or the observation date and time.")
        print("PEST_observation_file: ... {}".format(Observed_data_path[OBSGP]))
        #Terminate the script
        os._exit(0)
    
    #If (Flag_PEST_observation_TimeSeries == 1 or Flag_PEST_observation_MapSeries == 1) then counter is updated
    if (Flag_PEST_observation_TimeSeries == 1 or Flag_PEST_observation_MapSeries == 1):
        #Observation_counter is length of observation_value for observation group
        print(f"OBSGP={OBSGP}: Observation_counter={Observation_counter}")

    #NOBSPGP appends Observation_counter for this observation group
    NOBSPGP.append(Observation_counter)
    #Close Observed_data_file
    Observed_data_file.close()


##############################################################
#Algorithm to organize PEST .pst file variables
##############################################################
#define .pst variables, as described in PEST manual chapter 4
#control data 1st row
PSTHDRCD = "* control data"

#control data 2nd row
RSTFLE = "restart" #option to maintain restart data
PESTMODE = "estimation" #mode of PEST operation

#control data 3rd row
#number of parameters
NPAR = NPAR
#total number of observations per group, initiated to zero, updated below using NOBSPGP[]
NOBS = 0
#number of parameter groups, ranges from 1 to NPAR, choice is to use 1
NPARGP = 1 
#number of prior data, typically 0
NPRIOR = 0 
#number of observation groups
NOBSGP = NOBSGP

#control data 4th row
#number of input files with parameters
NTPLFLE = 1 
#number of instruction files for reading output; typically equal to NOBSGP or 1
NINSFLE = NOBSGP 
#precision of PEST parameters, with 12 characters maintained
PRECIS = "single" 
#decimal point is maintained
DPOINT = "point" 
#value associated with DERCOM, default is 1
NUMCOM = 1 
#optional parameters from 4th row
JACFILE = 0 
#optional parameters from 4th row
MESSFILE = 0 

#control data 5th row
#initial value of Marquardt lambda used in calculating fit, default is 10
RLAMBDA1 = 10#  
#value to adjust Marquardt lambda value, default is 2
RLAMFAC = 2
#phi ratio sufficient value, default is 0.3
PHIRATSUF = 0.3 
#alternative criterion to PHIRATSUF, default is 0.03
PHIREDLAM = 0.03 
#upper limit on number of lambdas PEST will test during iteration
NUMLAM = 10 
#JACUPDATE, LAMFORGIVE, DERFORGIVE #leaving off optional parameters from 5th row

#control data 6th row
#parameter change factor defined with parameter, for relative adjustment, default is 10
RELPARMAX = 3#  
#parameter change factor defined with parameter, for factor adjustment, default is 10 (does not allow for sign change)
FACPARMAX = 3#  
#parameter lower bound factor used in product with parameter original value, default is 0.001
FACORIG = 0.001 
#leaving off optional parameter in 6th row, parameter change factor defined with parameter, for absolute adjustment
#ABSPARMAX 
#iteration at which parameter sticks to its upper or lower bound, default is 0 which is not sticking
IBOUNDSTICK = 0 
#UPVECBEND #leaving off optional parameter in 6th row

#control data 7th row
#criteria of change in objective function for switching to higher order derivatives, used when FORCEN = switch, default is 0.1
PHIREDSWH = 0.1 
#NOPTSWITCH, SPLITSWH, DOAUI, DOSENREUSE, BOUNDSCALE #leaving off optional parameters in 7th row, related to switching to higher order derivatives

#control data 8th row
#maximum number of iterations for PEST per estimation run, default is 30, high value is 50
NOPTMAX = 30 
#criteria of change in objective function, default is 0.005
PHIREDSTP = 0.005 
#regulates number of iterations when below PHIREDSTP, default is 4
NPHISTP = 4 
#number of iterations before PEST terminates if objective function not lowering, default is 3
NPHINORED = 3 
#criteria of relative change in objective function, default is 0.005
RELPARSTP = 0.005 
#regulates number of iterations when below RELPARSTP, default is 4
NRELPAR = 4 
#PHISTOPTHRESH, LASTRUN and PHIABANDON #leaving off optional parameters in 8th row

#control data 9th row
#flag to write covariance coefficient matrix based on Jacobian matrix, default is 1
ICOV = 1 
#flag to write correlation coefficient matrix based on Jacobian matrix, default is 1
ICOR = 1 
#flag to write eignevectors of posterior covariance matrix based on Jacobian matrix, default is 1
IEIG = 1 
#IRES, JCOSAVE, JCOSAVEITN, VERBOSEREC, REISAVEITN, PARSAVEITN, PARSAVERUN #leaving off optional parameters in 9th row

#SVD 1st row; SVD-Assist is complementary section; if SVD is too slow, try LSQR as alternative
PSTHDRSVD = "* singular value decomposition"
#SVD 2nd row
#SVDMODE read from PEST_inputs_file 
#SVD 3rd row
#number of singular values before truncation, default is >= number of parameters
MAXSING = NPAR 
#ratio of lowest to highest eigenvalue of matrix at which singular truncation, default is 5E-y
EIGTHRESH = 0.0000005 
#SVD 4th row
#flag to write only singular value of eignvector to case.svd file
EIGWRITE = 1 

#parameter group 1st row
PSTHDRPARGP = "* parameter groups"
#parameter group 2nd row, written for each parameter value
#parameter group name, choice is to use 1 group
PARGPNME = "prgp01" 
#increment type for derivative of parameter, default is "relative", "absolute" is option
INCTYP = "relative"
#increment amount for INCTYP, default is 0.01
DERINC = 0.01 
#lower bound on increment amount, default is 0.0
DERINCLB = 0#  
#determines if derivatives are always forward difference, default is switch, which switches between forward difference and three-point methods
FORCEN = "switch" 
#multiplier by which DERINC is multiplied to change increment amount, default is 2
DERINCMUL = 2# 
#method of derivatives calculation, default is "parabolic" when FORCEN = "switch"
DERMTHD = "parabolic" 
#SPLITHRESH, SPLITRELDIFF and SPLITACTION #leaving off optional parameters in 2nd row

#parameter data 1st row
PSTHDRPARDATA = "* parameter data"
#parameter data 2nd row
#parameter name
#if int(FLAG_LOG_MODE) < 1:
#Note: get choice on FLAG_LOG_MODE variable transform mode from interface, 1 = log mode, 0 = no transform
if (int(FLAG_LOG_MODE) < 1):
    #PARNME = PARNME  
    #default value is "none", options include "log", "fixed", and "tied"; if using "tied" another line of parameter data is required.
    PARTRANS = "none" 
else:
    PARTRANS = "log"

#default value is "relative", but when PARTRANS = log, PARCHGLIM = "factor"
if (PARTRANS == "log"):
    PARCHGLIM = "factor"
else:
    PARCHGLIM = "relative"  
#initial value of parameter, read from user supplied inputs
PARVALI = 0 
#lower bound of parameter, read from user supplied inputs
PARLBND = 0 
#upper bound of parameter, read from user supplied inputs
PARUBND = 0 
#name of parameter group, could have unique name for each parameter controlled differently, default using 1 group
PARGP = PARGPNME 
#before writing parameter to input file, multiply the parameter by this value (conceal real value from inversion), default is 1
PSCALE = 1
#before writing parameter to input file, offset the parameter by this value (conceal real value from inversion), default is 0
POFFSET = 0
#derivative commands to run the model, default is 1
DERCOM = 1 

#observation group 1st row
PSTHDROBSGP = "* observation groups"
#observation group 2nd row; repeat for NOBSGP
#create OBGNME(), string with name of observation group
OBGNME = []
#For loop i_group in range 0 to NOBSGP
for i_group in range(0, NOBSGP):
    #observation group for observation, need a unique name for each NOBSGP
    OBGNME.append("obsgp0" + str(i_group))
    #NOBS is sum of int(NOBS) + NOBSPGP[i_group]
    NOBS = int(NOBS) + NOBSPGP[i_group]    

#observation data 1st row
PSTHDROBSDTA = "* observation data"
#observation group 2nd row; repeat for NOBS
#observation name, unique for each observation, assigned via loop below
#OBSNME = OBSNME
#observation value, provided by user, read in via loop below
OBSVAL = 0
#observation weight, default is 1
OBSWEIGHT = 1
##observation group for observation, need a unique name for each NOBSGP, defined above as vector
#OBGNME 

#model command line section
MODELCLIHDR = "* model command line"
#Path to model executable and command line argument as string
MODELCLI = HydroPlus_exe_path + " " + HydroPlus_input_path

#model input output section
MODELINOUTHDR = "* model input/output"
#PEST tpl file and Model parameter file as string
PESTTPL_MODELXML = PEST_tpl_str + " " + HydroPlusConfigXML_input_str

#PEST ins file and associated Model observations for 1st output
#cte 2025 not used?: PESTINS_MODELOBS = PEST_01_ins_str + " " + HydroPlus_prediction_path

#####################################################
#Algorithm to write PEST .pst file
#####################################################
#open PEST_pst_file_path file for writing output as PEST_pst_file
with open (PEST_pst_file_path, 'w') as PEST_pst_file:
    #PEST_pst_file write function places keyword pcf and line break
    #Note: pcf is PEST control file information
    PEST_pst_file.write("pcf " + "\n")
    
    #cd or control data section 
    #1st row
    #Note: PSTHDRCD is "* control data"
    PEST_pst_file.write(PSTHDRCD + "\n")
    #2nd row
    #Note: RSTFLE is option to maintain restart data; PESTMODE is mode of PEST operation
    PEST_pst_file.write("{: <15} {: <13} \n".format(RSTFLE, PESTMODE))
    #3rd row
    #Note: NPAR is # of parameters; NOBS is # of observations total; NPARGP is # of parameter groups; 
    #Note: NPRIOR is number of prior data; NOBSGP is # of observation groups
    PEST_pst_file.write("{: <15} {: <13} {: <13} {: <13} {: <13} \n".format(NPAR, NOBS, NPARGP, NPRIOR, NOBSGP))
    #4th row
    #Note: NTPLFLE is # of input files with parameters; NINSFLE is # of instruction files; PRECIS is precision of PEST parameters
    #Note: DPOINT is decimal point is maintained; NUMCOM is value associated with DERCOM or derivative commands to run the model; 
    #Note: JACFILE is optional parameters; MESSFILE is optional parameters
    PEST_pst_file.write("{: <15} {: <13} {: <13} {: <13} {: <13} {: <13} {: <13} \n".format(NTPLFLE, NINSFLE, PRECIS, DPOINT, NUMCOM, JACFILE, MESSFILE))
    #5th row
    #Note: RLAMBDA1 is value of Marquardt lambda used in calculating fit; RLAMFAC is value to adjust Marquardt lambda value
    #Note: PHIRATSUF is phi ratio sufficient value; PHIREDLAM is alternative criterion to PHIRATSUF
    #Note: NUMLAM is upper limit on number of lambdas PEST will test during iteration
    PEST_pst_file.write("{: <15} {: <13} {: <13} {: <13} {: <13} \n".format(RLAMBDA1, RLAMFAC, PHIRATSUF, PHIREDLAM, NUMLAM))
    #6th row
    #Note: RELPARMAX is parameter change factor defined with parameter; FACPARMAX is parameter change factor defined with parameter
    #Note: FACORIG is parameter lower bound factor used in product with parameter original value
    #Note: IBOUNDSTICK is iteration at which parameter sticks to its upper or lower bound
    PEST_pst_file.write("{: <15} {: <13} {: <13} {: <13} \n".format(RELPARMAX, FACPARMAX, FACORIG, IBOUNDSTICK))
    #7th row
    #Note: PHIREDSWH is criteria of change in objective function for switching to higher order derivatives, used when FORCEN = switch
    PEST_pst_file.write("{: <15} \n".format(PHIREDSWH))
    #8th row
    #Note: NOPTMAX is maximum number of iterations for PEST per estimation run; PHIREDSTP is criteria of change in objective function; 
    #Note: NPHISTP regulates number of iterations when below PHIREDSTP; NPHINORED is number of iterations before PEST terminates if objective function not lowering; 
    #Note: RELPARSTP is criteria of relative change in objective function; NRELPAR is regulates number of iterations when below RELPARSTP
    PEST_pst_file.write("{: <15} {: <13} {: <13} {: <13} {: <13} {: <13} \n".format(NOPTMAX, PHIREDSTP, NPHISTP, NPHINORED, RELPARSTP, NRELPAR))
    #9th row
    #Note: ICOV is flag to write covariance coefficient matrix based on Jacobian matrix; 
    #Note: ICOR is flag to write correlation coefficient matrix based on Jacobian matrix; IEIG is flag to write eignevectors of posterior covariance matrix based on Jacobian matrix
    PEST_pst_file.write("{: <15} {: <13} {: <13} \n".format(ICOV, ICOR, IEIG))
    
    #If int(SVDMODE) > 0 then write parameters
    if (int(SVDMODE) > 0):
        #svd or singular value decomposition section
        #1st row
        #Note: PSTHDRSVD is string "* singular value decomposition"
        PEST_pst_file.write(PSTHDRSVD + "\n")
        #2nd row
        #Note: SVDMODE is flag for use of singular value decomposition
        PEST_pst_file.write("{: <15} \n".format(SVDMODE))
        #3rd row
        #Note: MAXSING is number of singular values before truncation, default is >= number of parameters
        #Note: EIGTHRESH is ratio of lowest to highest eigenvalue of matrix at which singular truncation, default is 5E-y
        PEST_pst_file.write("{: <15} {: <13} \n".format(MAXSING, EIGTHRESH))
        #4th row
        #Note: EIGWRITE is flag to write only singular value of eignvector to case.svd file
        PEST_pst_file.write("{: <15} \n".format(EIGWRITE))
    
    #pg or parameter groups section
    #1st row
    PEST_pst_file.write(PSTHDRPARGP + "\n")
    #2nd row
    PEST_pst_file.write("{: <15} {: <13} {: <13} {: <13} {: <13} {: <13} {: <13} \n".format(PARGPNME, INCTYP, DERINC, DERINCLB, FORCEN, DERINCMUL, DERMTHD))

    #pd or parameter data section
    #1st row
    PEST_pst_file.write(PSTHDRPARDATA + "\n")
   
    #For loop i_par in range 0 to len(PARNME)
    for i_par in range(0, len(PARNME)):
        #PARVALI equals variable_initial[i_par]
        PARVALI = variable_initial[i_par]
        #PARLBND equals variable_min[i_par]
        PARLBND = variable_min[i_par]
        #PARUBND equals variable_max[i_par]
        PARUBND = variable_max[i_par]
        #next row(s)
        PEST_pst_file.write("{: <15} {: <13} {: <13} {: <13} {: <13} {: <13} {: <13} {: <13} {: <13} {: <13} \n".format(PARNME[i_par], PARTRANS, PARCHGLIM, PARVALI, PARLBND, PARUBND, PARGP, PSCALE, POFFSET, DERCOM))

    #og or observation groups section
    #1st row
    PEST_pst_file.write(PSTHDROBSGP + "\n")
    for i_group in range(0, NOBSGP):
        #next row(s)
        #Note: OBGNME[i_group] is observation group name
        PEST_pst_file.write("{: <15} \n".format(OBGNME[i_group]))

    #od or observation data section
    #1st row
    #Note: = PSTHDROBSDTA is defined as "* observation data"
    PEST_pst_file.write(PSTHDROBSDTA + "\n")
    
    #obs_cursor runs from 0 to len(observation_value)-1
    obs_cursor = 0
    #For loop i_group in range 0 to NOBSGP
    for i_group in range(0, NOBSGP):
        #For loop i_obs in range 0 to NOBSPGP[i_group]
        for i_obs in range(0, NOBSPGP[i_group]):
            #next row(s)
            #If Flag_PEST_observation_TimeSeries equals 1 then it is a time series and use i_obs will increment
            if (Flag_PEST_observation_TimeSeries == 1 and Flag_PEST_observation_MapSeries == 0):
                OBSVAL = observation_value[i_obs]
                #If Flag_PEST_If_ObservationIsZero_Then_WeightIsZero equals 1 and OBSVAl equals zero then set weight to zero
                if (Flag_PEST_If_ObservationIsZero_Then_WeightIsZero == 1 and OBSVAL != "0"):
                    OBSWEIGHT = 0
                #Else If Flag_PEST_If_ObservationIsZero_Then_WeightIsZero equals 0 then set weight to one
                elif (Flag_PEST_If_ObservationIsZero_Then_WeightIsZero != 1):
                    OBSWEIGHT = 1
            #ElIf Flag_PEST_observation_TimeSeries equals 0 then it is not a time series and use i_group to increment
            #Note: i_obs remains zero for all groups
            elif (Flag_PEST_observation_TimeSeries == 0 and Flag_PEST_observation_MapSeries == 0):
                OBSVAL = observation_value[i_group]
            #ElIf Flag_PEST_observation_MapSeries equals 1 then it is a map series and use i_obs will increment
            #Note: i_obs remains zero for all groups
            elif (Flag_PEST_observation_MapSeries == 1):
                #OBSVAL = observation_value[obs_cursor]
                #Note: obs_cursor extends across all observation groups, beyond i_obs within any one group
                OBSVAL = observation_value[obs_cursor]
                #obs_cursor += 1 to advance within observation_value vector
                obs_cursor += 1
                if (OBSVAL == "-9999"): 
                    OBSWEIGHT = 0
                else: 
                    OBSWEIGHT = 1
            
            #Note: OBSNME[i_group] + str(i_obs) is observation group + number i_obs, followed by OBSVAL w/ observation value, and OBSWEIGHT w/ observation weight, 
            PEST_pst_file.write("{: <15} {: <13} {: <13} {: <13} \n".format(OBSNME[i_group] + str(i_obs), OBSVAL, OBSWEIGHT, OBGNME[i_group]))

    #mlc or model command line section
    #1st row
    PEST_pst_file.write(MODELCLIHDR + "\n")
    #2nd row
    PEST_pst_file.write(MODELCLI + "\n")

    #mio or model input output section
    #1st row
    PEST_pst_file.write(MODELINOUTHDR + "\n")
    #2nd row
    PEST_pst_file.write(PESTTPL_MODELXML + "\n")

    #HydroPlus_prediction_path configuration: Normalize predictions to a list of length NOBSGP
    if isinstance(HydroPlus_prediction_path, str):
        HydroPlus_prediction_paths = [HydroPlus_prediction_path]
    else:
        # e.g., list/tuple -> list
        HydroPlus_prediction_paths = list(HydroPlus_prediction_path)  

    #If only one path is provided, replicate it across all groups
    if len(HydroPlus_prediction_paths) == 1 and NOBSGP > 1:
        HydroPlus_prediction_paths = HydroPlus_prediction_paths * NOBSGP


    #For loop i_group in range 0 to NOBSGP
    for i_group in range(0,NOBSGP):

        #Define PEST_01_ins_str as PEST_Output_01.ins
        PEST_01_ins_str = f"HydroPlus_PEST_Output_{i_group}.ins" 
        
        #Set path to PEST_ins_file_path as HydroPlus_input_path + PEST_01_ins_str
        PEST_ins_file_path = os.path.join(HydroPlus_input_path, PEST_01_ins_str)
        
        #PESTINS_MODELOBS f"{PEST_01_ins_str} {HydroPlus_prediction_paths[i_group]}"; contains PEST ins file and associated model observation for this output
        #Note: PEST expects ins file and associated Model observations for output
        PESTINS_MODELOBS = f"{PEST_01_ins_str} {HydroPlus_prediction_paths[i_group]}"
        #cte 2025
        print(f"HydroPlus_prediction_paths={HydroPlus_prediction_paths[i_group]}")
        
        #3rd row to PEST_pst_file needs to be within loop
        PEST_pst_file.write(PESTINS_MODELOBS + "\n")

        #####################################################
        #Algorithm to write PEST .ins files
        #####################################################
        #Note: For the .ins file, each instruction line must begin with either a primary marker or a line advance item. 
        #Note: The primary marker is a string of characters, bracketed at each end by a marker delimiter. 
        #Note: ... If a marker is the first item on an instruction line, then it is a primary marker; if it occurs later in the line, following other instruction items, it is a secondary marker.
        #Note: ... PEST reads the model output file, line by line, searching for the string between the marker delimiter characters.
        #Note: ... A primary marker may be the only item on a PEST instruction line, or it may precede a number of other items directing further processing of the line containing the marker. 
        #Note: ... In the former case the purpose of the primary marker is simply to establish a reference point for further downward movement within the model output file as set out in subsequent instruction lines.
        #Note: The line advance item is “ln” where n is the number of lines to advance; note that “l” is “el”, the twelfth letter of the alphabet, not “one”.
        #Note: ... The line advance will more rapidly arrive at a line than the primary marker, so if the output file has a fixed line with the primary marker, use line advance for efficiency.
        #Note: ... Use a secondary marker to isolate a variable in the instruction file, as shown next: 
        #Note: ... Output file line: SOIL WATER CONTENT (NO CORRECTION)=21.345634%
        #Note: ... Ins file line: l5 *=* !sws! *%*
        #Note: ... Use whitespace to isolate a variable in the instruction file, as shown next: 
        #Note: ... Output file line: 4.33 -20.3 23.392093 3.394382
        #Note: ... Ins file line: l10 w w w !obs1!
        #open PEST_ins_file_path file for writing output as PEST_01_ins_file
        #Note: Prepare to write .ins instruction file
        with open (PEST_ins_file_path, 'w') as PEST_01_ins_file:
            #PEST_01_ins_file write function places keyword pif, PEST_key_symbol, and line break
            #Note: pif is PEST instruction file
            PEST_01_ins_file.write("pif " + PEST_key_symbol + "\n")

            #Create string with secondary marker
            #def repeat_symbols function takes the argument Flag_PEST_prediction_ColumnID
            def repeat_symbols(Flag_PEST_prediction_ColumnID):
                #marker_symbol_str is the string of symbols used as a marker to find the output variable, in this case comma
                marker_symbol_str = "%s,%s" %(PEST_key_symbol,PEST_key_symbol)
                #result is output from whitespace " " and join, with for loop w/ _ as throwaway variable in range of Flag_PEST_prediction_ColumnID
                #Note: range(Flag_PEST_prediction_ColumnID) creates a sequence of integers from 0 to Flag_PEST_prediction_ColumnID - 1
                #Note: [marker_symbol_str for _ in range(Flag_PEST_prediction_ColumnID)] is a list comprehension that creates a list of num_times copies of the string symbols. 
                #Note: " ".join(...) is used to join the elements of the list created in step 2 into a single string, separated by a space " "
                result = " ".join([marker_symbol_str for _ in range(Flag_PEST_prediction_ColumnID)])
                #return result to call below
                return result
            #try function
            try:
                #If Flag_PEST_prediction_ColumnID < 0 then warn user
                if Flag_PEST_prediction_ColumnID < 0:
                    print("Please enter a positive integer equal or greater than 0 for Flag_PEST_prediction_ColumnID.")
                #Else Flag_PEST_prediction_ColumnID >=0 and proceed
                else:
                    #marker_symbol_combined_str is the combined markers for use in the .ins file, called from function repeat_symbols
                    marker_symbol_combined_str = repeat_symbols(Flag_PEST_prediction_ColumnID)
                    #print(marker_symbol_combined_str)
            except ValueError:
                print("Invalid input. Please enter a valid integer for Flag_PEST_prediction_ColumnID.")

            #If Flag_PEST_observation_TimeSeries equals 1, then time series instructions are simple, advance each line
            #Note: Section works for time series observation using chronological pst_observation_start_index & pst_observation_stop_index
            if (Flag_PEST_observation_TimeSeries == 1 and Flag_PEST_observation_MapSeries != 1):
                #PEST_01_ins_file.write command for skipping first line of data ("l1 " + "\n")
                PEST_01_ins_file.write("l1 " + "\n")

                #For loop irow in range 0 to int(pst_observation_stop_index - pst_observation_start_index) + 1
                for irow in range(0,int(pst_observation_stop_index - pst_observation_start_index) + 1):
                    #PEST_01_ins_file.write the commands for reading predicted output orl1 " must + 1OBSNME[i_group] + str(irow) + "!" + "\n")
                    #Note: Instructions are simply "l1", ell one, which proceeds 1 line and read; 1st row of file is skipped, header containing number of rows
                    #Note: Instructions can easily handle more complex instructions, e.g., identify a keyword in file to find and then begin reading next line, ...
                    #Note: ... placing keyword characters within PEST_key_symbol, e.g., #YYYYMMDD#, and then for each row number of commas to pass, e.g., l1 #,# #,# #,# !Cl1!
                    PEST_01_ins_file.write("l1 " + marker_symbol_combined_str + " !" + OBSNME[i_group] + str(irow) + "!" + "\n")
            elif (Flag_PEST_observation_TimeSeries == 1 and Flag_PEST_observation_MapSeries == 1): 
                print("Invalid input. Flag_PEST_observation_TimeSeries and Flag_PEST_observation_MapSeries cannot both be 1")
            #If Flag_PEST_observation_TimeSeries not equal to 1, then write each instruction as one line to unique files
            if (Flag_PEST_observation_TimeSeries != 1 and Flag_PEST_observation_MapSeries != 1):
                #Target output is: #Primary Marker,# !Variable! #SecondaryMarker#"
                #Note: # is the key to surround primary and secondary markers, and ',' must be at end of primary, start of secondary
                PEST_01_ins_file.write("#{},# !{}0! #,#\n".format(Observed_variable_full_name[i_group],OBSNME[i_group]))
            elif (Flag_PEST_observation_TimeSeries == 0 and Flag_PEST_observation_MapSeries == 0): 
                print("Invalid input. Flag_PEST_observation_TimeSeries or Flag_PEST_observation_MapSeries must be 1")

            #If Flag_PEST_observation_MapSeries equals 1, then map instructions are to write each row and column
            #Note: Section works for map series observation using spatial row and column indicators
            if (Flag_PEST_observation_MapSeries == 1 and Flag_PEST_observation_TimeSeries != 1):
                #PEST_01_ins_file.write command for primary marker, end of header ("#NODATA_value#" + "\n")
                PEST_01_ins_file.write("#NODATA_value#" + "\n")
                #For loop through nrows
                for irow in range(nrows):
                    #line_ins starts with l1, ell one, and space
                    line_ins = "l1 "
                    #For loop through ncols
                    for icol in range(1, ncols + 1):
                        #if not last ncols+1
                        if (icol < ncols):
                            #line_ins has added OBSNME[i_group] + str(irow * ncols + icol - 1) + w, for whitespace
                            line_ins += ("!" + OBSNME[i_group] + str(irow * ncols + icol - 1) + "! w ")
                        else:
                            #line_ins has added OBSNME[i_group] + str(irow * ncols + icol - 1) + \n for line return
                            line_ins += ("!" + OBSNME[i_group] + str(irow * ncols + icol - 1) + "!" + "\n")
                    
                    #PEST_01_ins_file has line written, with the trailing space removed at end of ncol loop
                    PEST_01_ins_file.write(line_ins.rstrip() + '\n')

        #Close PEST_01_ins_file
        PEST_01_ins_file.close()

        #print output to monitor
        print("PEST calibration output for i_group {}, NOBSPGP = {}, NOBS = {}".format(i_group, NOBSPGP[i_group], NOBS))

#####################################################
#Algorithm to write PEST .tpl files
#####################################################
#Write XML_File_Tree as new HydroPlusConfig.xml file to folder defined by path_to_TestCase_HydroPlusConfigXML
#Note: XML_File_Tree is created above using LMXL_File.parse(XML_Input_and_Template_path)
#Note: LMXL_File is from import lmxl.etree as LMXL_File
XML_File_Tree.write(PEST_tpl_file_path)

#open PEST_tpl_file_path file for input as project_XML_output_file
with open(PEST_tpl_file_path) as PEST_tpl_file_file:
    #XML_file_lines equals all XML data i PEST_tpl_file_file
    XML_file_lines = PEST_tpl_file_file.read()

#open PEST_tpl_file_path file for writing output as PEST_tpl_file_file
with open (PEST_tpl_file_path, 'w') as PEST_tpl_file_file:
    #PEST_tpl_file_file write function places keyword ptf, PEST_key_symbol, and line break
    PEST_tpl_file_file.write("ptf " + PEST_key_symbol + "\n")
    #PEST_tpl_file_file write function places all XML_file_lines 
    PEST_tpl_file_file.write(XML_file_lines)

#####################################################
#Close files and write comments
#####################################################
#Close PEST_inputs_file
PEST_inputs_file.close()
#Close PEST_tpl_file_file
PEST_tpl_file_file.close()
#Close PEST_pst_file
PEST_pst_file.close()

print("\nPEST preparation complete. Follow these steps to use the PEST files:")
print("1. Open XTerm window in directory containing the file {}".format(PEST_pst_str))
print("2a. With 64-bit PEST, enter in that XTerm command line >i64pest {}".format(PEST_pst_str))
print("2b. With 32-bit PEST, enter in that XTerm command line >pest {}".format(PEST_pst_str))
print("\nMay PEST Calibration Help Us Improve Our World with i-Tree Tools.")
