# Copyright (c) 2018 Kelvin Say
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

library(lubridate)

source('preprocess_read_profiles.R')
source('process_technical.R')
source('process_financial.R')
source('visualise.R')
source('utility_functions.R')

# ---------------------------------------------------------------------------------------------------------------------------
# Configuration Options for the Research Scenario
# ---------------------------------------------------------------------------------------------------------------------------
g.force_update <- FALSE
g.save_detailed <- FALSE

g.scenario_name <-
  'P1.xxx'

g.scenario_simulation <- list(
  'Time Step (s)' = 900,
  'Starting Year' = 2018,
  'Number of Years' = 16,
  'PV Start (Wp)' = 0,
  'PV End (Wp)' = 10000,
  'PV Delta (Wp)' = 500,
  'Battery Start (Wh)' = 0,
  'Battery End (Wh)' = 20000,
  'Battery Delta (Wh)' = 1000,
  'NPV Skim Top (%)' = 0.05)

g.scenario_load <- list(
  'Profile' = 'SunWiz: Double Peak',
  'Annual Consumption (MWh)' = 5.2,
  'Filename' = 'pattern_double_peak.csv',
  'Filepath' = file.path(getwd(), '_Scenarios/ConsumptionPatterns_SunWiz'))

g.scenario_insolation <- list(
  'Profile' = 'PVWatts: Perth',
  'Type' = 'Fixed',
  'Tilt' = 20,
  'Azimuth' = 0,
  'SysLoss (%)' = 0.14, # xxx - from PVWatts, need to remove it and get the PV info to apply it, along with tilt/azimuth...
  'Latitude' = 31.95,
  'Longitude' = 115.77,
  'Filename' = 'pvwatts_perth_hourly.csv',
  'Filepath' = file.path(getwd(), '_Scenarios/Insolation_PVWatts'))

g.scenario_pv <- list(
  'Panel Type' = 'Generic',
  'System Losses (%)' = 0.14,
  '25 Year Capacity (%)' = 0.8)

g.scenario_battery <- list(
  'Type' = "Tesla Powerwall 2",
  'Input Loss (%)' = 1 - sqrt(0.89),
  'Output Loss (%)' = 1 - sqrt(0.89),
  'Input Limit (W)' = 5000,
  'Output Limit (W)' = 5000,
  '10 Year Capacity (%)' = 0.7)

g.scenario_tariff <- list(
  'Import ($/kWh)' = 0.27, # 0.26474
  'Export ($/kWh)' = 0.07, # 0.07135
  'Import Inflation (%)' = 0.05,
  'Export Inflation (%)' = 0.05,
  'FiT System Limit (Wp)' = 5000)

g.scenario_costs <- list(
  'PV Cost ($/kWp)' = 1400,
  'Battery Cost ($/kWh)' = 900,
  'PV Cost Change (%)' = -0.059,
  'Battery Cost Change (%)' = -0.08)

g.scenario_finance <- list(
  'Analysis Years' = 10,
  'Discount Rate (%)' = 0.06)

# ---------------------------------------------------------------------------------------------------------------------------
# Scenario Execution
# ---------------------------------------------------------------------------------------------------------------------------
run <- function() {
  graphics.off()

  print(g.scenario_name)

  start_msg <- ""
  new_technical_simulation_required <- TRUE
  if (!g.force_update) {
    if (exists("technical_data") & exists("technical_summary")) {
      if ((all.equal(g.scenario_simulation, technical_data$Inputs$Scenario$Simulation)[1] == TRUE) &
          (all.equal(g.scenario_load, technical_data$Inputs$Scenario$Load)[1] == TRUE) &
          (all.equal(g.scenario_insolation, technical_data$Inputs$Scenario$Insolation)[1] == TRUE) &
          (all.equal(g.scenario_pv, technical_data$Inputs$Scenario$PV)[1] == TRUE) &
          (all.equal(g.scenario_battery, technical_data$Inputs$Scenario$Battery)[1] == TRUE) &
          g.scenario_finance$'Analysis Years' == technical_data$Inputs$Scenario$Financial$'Analysis Years')
      {
        new_technical_simulation_required <- FALSE
        technical_data$Inputs$Scenario$Name <- g.scenario_name
        technical_data$Inputs$Scenario$Financial$'Discount Rate (%)' <- g.scenario_finance$'Discount Rate (%)'
        technical_summary$Inputs$Scenario$Name <- g.scenario_name
        start_msg <- "No technical changes, reusing technical data"
      } else {
        start_msg <- "Technical changes detected, re-generating technical data"
      }
    } else {
      start_msg <- "Generating technical data"
    }
  } else {
    start_msg <- "Force updating"
  }
  print(start_msg)

  if (new_technical_simulation_required) {
    technical_data <- process_Scenario_Technical(g.scenario_simulation, g.scenario_finance, g.scenario_load, g.scenario_insolation, g.scenario_pv, g.scenario_battery)
    assign("technical_data", technical_data, envir = .GlobalEnv)
    technical_summary <- process_Summarise_Technical(technical_data)
    assign("technical_summary", technical_summary, envir = .GlobalEnv)
  }

  economic_summary <- process_Scenario_Economic(technical_data, technical_summary, g.scenario_simulation, g.scenario_finance, g.scenario_tariff, g.scenario_costs)
  assign("economic_summary", economic_summary, envir = .GlobalEnv)

  # Save the outcome data
  path_scenario_results <- file.path(getwd(), "_Outputs", g.scenario_name)
  if (!g.save_detailed) {
    saveRDS(technical_data, file.path(path_scenario_results, "data_technical_data.rds"))
  }
  saveRDS(technical_summary, file.path(path_scenario_results, "data_technical_summary.rds"))
  saveRDS(economic_summary, file.path(path_scenario_results, "data_economic_summary.rds"))

  # Print out all input elements for this scenario
  sink(file.path(path_scenario_results, "data_inputs.txt"))
  print(economic_summary$Scenario)
  sink()

  print("Scenario analysis complete")
}

# ---------------------------------------------------------------------------------------------------------------------------
#
# ---------------------------------------------------------------------------------------------------------------------------
process_Scenario_Technical <- function(simulation_info, financial_info, load_info, insolation_info, pv_info, battery_info) {
  start_time <- now()
  time_a <- start_time

  # Store the raw scenario inputs (for later comparison)
  scenario_inputs <- list(Name = g.scenario_name,
                          Simulation = simulation_info,
                          Financial = financial_info,
                          Load = load_info,
                          Insolation = insolation_info,
                          PV = pv_info,
                          Battery = battery_info)

  # Create the time line vector (assuming a 365 day year of 2017)
  num_steps <- 365 * 3600 * 24 / simulation_info$'Time Step (s)'
  v_time_series <- ymd_hms("2017-01-01 00:00:00") + seconds(simulation_info$'Time Step (s)' * (0:(num_steps-1)))

  # Obtain the reference insolation_info and baseload profiles # xxx - obtain from a specific file/folder specified above
  preprocess_Read_Scenario(i.file_name_load = load_info$'Filename',
                           i.file_path_load = load_info$'Filepath',
                           i.file_name_insolation = insolation_info$'Filename',
                           i.file_path_insolation = insolation_info$'Filepath',
                           i.sample_period_sec = simulation_info$'Time Step (s)',
                           i.annual_consumption_kWh = (load_info$'Annual Consumption (MWh)' * 1000))
  insolation_data <- fread(file.path(getwd(), "_Inputs", "UnitInsolation.csv"))
  v_solar_insolation_W <- insolation_data$'Power' / 1000
  baseload_profile_data <- fread(file.path(getwd(), "_Inputs", "HouseholdLoad.csv"))
  v_baseload_W <- baseload_profile_data$'NetPower'

  # Setup Technical Data
  # --------------------
  # Prepare the technical process' input objects

  # Time information
  technical_timeline <- list(Info = list('Operational Years' = financial_info$'Analysis Years',
                                         'Time Step (s)' = simulation_info$'Time Step (s)'),
                             Data = list('Time' = v_time_series))

  # Solar and PV information
  technical_solar <- list(Info = list('Nominal Capacity (Wp)' = NULL,
                                      'Original Sample Period (s)' = NULL), # xxx - should come from insolation_data field
                          Data = list('Generation (W)' = NULL,
                                      'Insolation Reference (W/Wp)' = v_solar_insolation_W))
  technical_solar$Info <- append(technical_solar$Info, insolation_info)
  technical_solar$Info <- append(technical_solar$Info, pv_info)

  # Baseload information
  technical_load <- list(Info = list('Original Sample Period (s)' = NULL), # xxx - should come from baseload_profile_data field
                         Data = list('Load (W)' = v_baseload_W))
  technical_load$Info <- append(technical_load$Info, load_info)

  # Battery information
  technical_battery <- list(Info = list('Rated Capacity (Wh)' = NULL),
                            Data = list('Capacity (Wh)' = NULL),
                            Initial = list('SoC (Wh)' = NULL))
  technical_battery$Info <- append(technical_battery$Info, battery_info)

  # Store the input data away for future reference
  technical_inputs <- list(Timeline = technical_timeline,
                           Solar = technical_solar,
                           Baseload = technical_load,
                           Battery = technical_battery,
                           Scenario = scenario_inputs)

  # Model the technical operation of various combinations for this scenario
  # -----------------------------------------------------------------------
  # Run through the range of PV and Battery combinations
  pv_num_elements <- (simulation_info$'PV End (Wp)' - simulation_info$'PV Start (Wp)') / simulation_info$'PV Delta (Wp)' + 1
  battery_num_elements <- (simulation_info$'Battery End (Wh)' - simulation_info$'Battery Start (Wh)') / simulation_info$'Battery Delta (Wh)' + 1

  technical_results <- matrix(list(list()), nrow = pv_num_elements, ncol = battery_num_elements)

  print(paste0("--> Initialised: ", round(now() - time_a,2), " seconds"))
  time_a <- now()

  pv_size_Wp <- simulation_info$'PV Start (Wp)'
  for (i in 1:pv_num_elements) {
    technical_solar$Info$'Nominal Capacity (Wp)' <- pv_size_Wp

    battery_size_Wh <- simulation_info$'Battery Start (Wh)'
    for (j in 1:battery_num_elements) {
      technical_battery$Info$'Rated Capacity (Wh)' <- battery_size_Wh

      results <- process_MeterLoad_Years(simulation_info, technical_timeline, technical_load, technical_solar, technical_battery)
      technical_results[i,j][[1]]$'PV Capacity (Wp)' <- pv_size_Wp
      technical_results[i,j][[1]]$'Battery Capacity (Wh)' <- battery_size_Wh
      technical_results[i,j][[1]]$'Daily' <- results$'Daily'
      if (g.save_detailed) {
        technical_results[i,j][[1]]$'Detailed' <- results$'Operation' # This ONLY contains the observation data for the final year
      }
      technical_results[i,j][[1]]$'Summary' <- results$'Summary'

      print(paste0("--> Processed [", pv_size_Wp, ",", battery_size_Wh, "]: ", round(now() - time_a,2), " seconds"))
      time_a <- now()

      battery_size_Wh <- battery_size_Wh + simulation_info$'Battery Delta (Wh)'
    }
    pv_size_Wp <- pv_size_Wp + simulation_info$'PV Delta (Wp)'
  }

  print(paste0("====> COMPLETED <====="))
  print(round(now() - start_time))

  return(list('Inputs' = technical_inputs,
              'Operational' = technical_results))
}


# ---------------------------------------------------------------------------------------------------------------------------
#
# ---------------------------------------------------------------------------------------------------------------------------
process_Scenario_Economic <- function(technical_data, technical_summary, simulation_info, finance_info, tariff_info, costs_info) {

  year <- simulation_info$'Starting Year'
  num_years <- simulation_info$'Number of Years'
  init_import_tariff <- tariff_info$'Import ($/kWh)'
  init_export_tariff <- tariff_info$'Export ($/kWh)'
  init_pv_cost <- costs_info$'PV Cost ($/kWp)'
  init_battery_cost <- costs_info$'Battery Cost ($/kWh)'
  baseload <- technical_summary$Total$'Load (kWh)'
  num_pv <- length(technical_summary$'PV Capacity (Wp)')
  npv_skim <- simulation_info$'NPV Skim Top (%)'
  path_scenario_results <- file.path(getwd(), "_Outputs", g.scenario_name)
  init_tariff_info <- tariff_info
  init_costs_info <- costs_info

  dir.create(path_scenario_results, showWarnings = FALSE)
  do.call(file.remove, list(list.files(path_scenario_results, full.names = TRUE)))

  output <- data.table(year = simulation_info$'Starting Year':(simulation_info$'Starting Year' + simulation_info$'Number of Years' - 1),
                       cfgtype = rep('No system', num_years),
                       descr = rep('None', num_years),
                       npv = rep(NaN, num_years),
                       capital = rep(NaN, num_years),
                       dependence = rep(100, num_years),
                       dependence.min = rep(100, num_years),
                       dependence.max = rep(100, num_years),
                       gen_dispatch = rep(0, num_years),
                       gen_dispatch.min = rep(0, num_years),
                       gen_dispatch.max = rep(0, num_years),
                       gen_nondispatch = rep(0, num_years),
                       gen_nondispatch.min = rep(0, num_years),
                       gen_nondispatch.max = rep(0, num_years),
                       payback = rep(NaN, num_years),
                       payback.min = rep(NaN, num_years),
                       payback.max = rep(NaN, num_years),
                       npv_max = rep(NaN, num_years),
                       npv_max.min = rep(NaN, num_years),
                       npv_max.max = rep(NaN, num_years),
                       exported = rep(0, num_years))

  input <- data.table(year = simulation_info$'Starting Year':(simulation_info$'Starting Year' + simulation_info$'Number of Years' - 1),
                      growth_tariff_import = rep(0, num_years),
                      growth_tariff_export = rep(0, num_years),
                      growth_pv_cost = rep(0, num_years),
                      growth_batt_cost = rep(0, num_years))

  for (i in 1:num_years) {

    financial_data <- process_Financial_Viability(technical_summary, finance_info, tariff_info, costs_info, FALSE, FALSE)

    # Collate the results for export later
    if (i == 1) {
      detailed_financials <- list(financial_data)
      names(detailed_financials) <- year
    } else {
      name_list <- c(names(detailed_financials), year)
      detailed_financials[[i]] <- financial_data
      names(detailed_financials) <- name_list
    }

    # Determine the PV + battery combination that has the highest NPV
    maxnpv <- max(financial_data$'System.NPV')
    if (maxnpv > 0) {
      v_max <- which(financial_data$'System.NPV'==maxnpv)
      max_battery_index <- ceiling(v_max[1]/num_pv)
      max_pv_index <- (v_max[1]-1) %% num_pv + 1

      # Record the system results from the maxima
      if (technical_summary$'Battery Capacity (Wh)'[max_battery_index] == 0) {
        cfgtype <- "PV"
        descr <- paste0(technical_summary$'PV Capacity (Wp)'[max_pv_index]/1000, "kWp")
      } else {
        cfgtype <- "PV-Battery"
        descr <- paste0(technical_summary$'PV Capacity (Wp)'[max_pv_index]/1000, "kWp + ",technical_summary$'Battery Capacity (Wh)'[max_battery_index]/1000, "kWh")
      }
      if (length(v_max) > 1) {
        descr <- paste0(descr, "*")
      }
      capital <- financial_data$'System.Costs'[max_pv_index, max_battery_index]
      dependence <- technical_summary$Total$'Dependence (%)'[max_pv_index, max_battery_index]
      battery_discharged <- technical_summary$Total$'Battery Discharged (kWh)'[max_pv_index, max_battery_index]
      exported <- technical_summary$Total$'Exported (kWh)'[max_pv_index, max_battery_index]
      payback <- financial_data$'System.Payback'[max_pv_index, max_battery_index]
      npv_max <- maxnpv

      output$descr[i] <- descr
      output$npv[i] <- maxnpv
      output$npv_max[i] <- npv_max
      output$payback[i] <- payback
      output$cfgtype[i] <- cfgtype
      output$capital[i] <- round(capital,2)
      output$dependence[i] <- round(dependence*100,2)
      output$gen_dispatch[i] <- round(battery_discharged/baseload*100,2)
      output$gen_nondispatch[i] <- round(exported/baseload*100,2)
      output$exported[i] <- round(exported,0)
    }

    # Record the input financial values
    growth_tariff_import <- (tariff_info$'Import ($/kWh)'-init_import_tariff)/init_import_tariff
    if (init_export_tariff > 0) {
      growth_tariff_export <- (tariff_info$'Export ($/kWh)'-init_export_tariff)/init_export_tariff
    } else {
      growth_tariff_export <- 0
    }
    growth_pv_cost <- (costs_info$'PV Cost ($/kWp)'-init_pv_cost)/init_pv_cost
    growth_batt_cost <- (costs_info$'Battery Cost ($/kWh)'-init_battery_cost)/init_battery_cost
    input$year[i] <- year
    input$growth_tariff_import[i] <- round(growth_tariff_import*100,2)
    input$growth_tariff_export[i] <- round(growth_tariff_export*100,2)
    input$growth_pv_cost[i] <- round(growth_pv_cost*100,2)
    input$growth_batt_cost[i] <- round(growth_batt_cost*100,2)

    # Record the values of the region around the maximum NPV system (defined by the "skim %" parameter)
    v_region <- which( (financial_data$'System.NPV' >= (maxnpv * (1-npv_skim))) &  (financial_data$'System.NPV' > 0) )
    if (length(v_region) > 0) {
      output$dependence.min[i] <- Inf
      output$dependence.max[i] <- -Inf
      output$gen_dispatch.min[i] <- Inf
      output$gen_dispatch.max[i] <- -Inf
      output$gen_nondispatch.min[i] <- Inf
      output$gen_nondispatch.max[i] <- -Inf
      output$payback.min[i] <- Inf
      output$payback.max[i] <- -Inf
      output$npv_max.min[i] <- Inf
      output$npv_max.max[i] <- -Inf

      for (j in 1:length(v_region)) {
        battery_index <- ceiling(v_region[j]/num_pv)
        pv_index <- (v_region[j]-1) %% num_pv + 1

        dependence <- round(technical_summary$Total$'Dependence (%)'[pv_index, battery_index]*100, 2)
        output$dependence.min[i] <- min(output$dependence.min[i], dependence)
        output$dependence.max[i] <- max(output$dependence.max[i], dependence)

        gen_dispatch <- round(technical_summary$Total$'Battery Discharged (kWh)'[pv_index, battery_index]/baseload*100, 2)
        output$gen_dispatch.min[i] <- min(output$gen_dispatch.min[i], gen_dispatch)
        output$gen_dispatch.max[i] <- max(output$gen_dispatch.max[i], gen_dispatch)

        gen_nondispatch <- round(technical_summary$Total$'Exported (kWh)'[pv_index, battery_index]/baseload*100, 2)
        output$gen_nondispatch.min[i] <- min(output$gen_nondispatch.min[i], gen_nondispatch)
        output$gen_nondispatch.max[i] <- max(output$gen_nondispatch.max[i], gen_nondispatch)

        payback <- round(financial_data$'System.Payback'[pv_index, battery_index], 2)
        output$payback.min[i] <- min(output$payback.min[i], payback)
        output$payback.max[i] <- max(output$payback.max[i], payback)

        npv <- round(financial_data$'System.NPV'[pv_index, battery_index], 2)
        output$npv_max.min[i] <- min(output$npv_max.min[i], npv)
        output$npv_max.max[i] <- max(output$npv_max.max[i], npv)
      }

    }

    mypath <- file.path(path_scenario_results, paste(year, ".jpg", sep = ""))
    jpeg(filename=mypath, width=909, height=432, pointsize=12, quality=100, res=96)
    plot_npv <- vis_Annual_NPV(simulation_info, technical_summary, financial_data)
    plot_optimal <- vis_Annual_Optimal(simulation_info, technical_summary, financial_data, v_region, output$year[i])
    multiplot(plot_npv, plot_optimal, cols=2)
    dev.off()

    tariff_info$'Import ($/kWh)' <- tariff_info$'Import ($/kWh)' * (1 + tariff_info$'Import Inflation (%)')
    tariff_info$'Export ($/kWh)' <- tariff_info$'Export ($/kWh)' * (1 + tariff_info$'Export Inflation (%)')
    costs_info$'PV Cost ($/kWp)' <- costs_info$'PV Cost ($/kWp)' * (1 + costs_info$'PV Cost Change (%)')
    costs_info$'Battery Cost ($/kWh)' <- costs_info$'Battery Cost ($/kWh)' * (1 + costs_info$'Battery Cost Change (%)')
    year <- year + 1
  }

  # Save the economic scenario data for future reference
  economic_summary <- list('Financial Assumptions' = input,
                           'Economic Outcomes' = output,
                           'Detailed Financials' = detailed_financials,
                           'Scenario' = technical_data$Inputs$'Scenario')
  economic_summary$'Scenario' <- append(economic_summary$'Scenario', list('System Costs' = init_costs_info))
  economic_summary$'Scenario' <- append(economic_summary$'Scenario', list('Electricity Tariffs' = init_tariff_info))

  # Plot the change in energy dependence and available battery capacity
  mypath <- file.path(path_scenario_results, "Results_EnergyDependence.jpg")
  jpeg(filename=mypath, width=909, height=432, pointsize=12, quality=100, res=96)
  plot(vis_Scenario_EnergyDependence(economic_summary))
  dev.off()

  # Plot the changing NPV amounts
  mypath <- file.path(path_scenario_results, "Results_NPV.jpg")
  jpeg(filename=mypath, width=909, height=216, pointsize=12, quality=100, res=96)
  plot(vis_scenario_ResultsNPV(output))
  dev.off()

  #Plot the changing payback periods
  mypath <- file.path(path_scenario_results, "Results_Payback.jpg")
  jpeg(filename=mypath, width=909, height=216, pointsize=12, quality=100, res=96)
  plot(vis_Scenario_ResultsPayback(output))
  dev.off()

  # Plot the input financial assumptions
  mypath <- file.path(path_scenario_results, "Financial_Assumptions.jpg")
  jpeg(filename=mypath, width=909, height=216, pointsize=12, quality=100, res=96)
  plot(vis_Scenario_FinancialAssumptions(input))
  dev.off()

  # Plot the energy profile across a single day without a battery or PV installed
  mypath <- file.path(path_scenario_results, "Energy_Profile.jpg")
  jpeg(filename=mypath, width=909, height=432, pointsize=12, quality=100, res=96)
  plot(vis_Day_EnergyProfile(technical_data, 1, 2, 2000, 6000))
  dev.off()

  # Backup the insolation and load files into the scenario directory
  flist <- list.files(file.path(getwd(), "_Inputs"), full.names = TRUE)
  file.copy(flist, file.path(path_scenario_results))

  # Display the changing energy depedence
  plot(vis_Scenario_EnergyDependence(economic_summary))

  return(economic_summary)
}