text:
"""
    multifrequency_opf(
        folder::String,
        obj::String;
        gen_areas=[],
        area_interface=[],
        gen_zones=[],
        zone_interface=[],
        print_results::Bool=false,
        override_param::Dict{Any}=Dict(),
        fix_f_override::Bool=false,
        direct_pq::Bool=true,
        master_subnet::Int64=1,
        suffix::String="",
        start_vals=Dict{String, Dict}("sn"=>Dict()),
        no_converter_loss::Bool=false,
        uniform_gen_scaling::Bool=false,
        unbounded_pg::Bool=false
    )

Models and solves the OPF for a single network with data contained in `folder`.

# Arguments
- `folder::String`: the directory containing all subnetwork data, *subnetworks.csv*, and *interfaces.csv*
- `obj::String`: the objective function to use, from the following:
   - "mincost": minimize generation cost
   - "areagen": minimize generation in the areas specified in `gen_areas`
   - "zonegen": minimize generation in the zones specified in `gen_zones`
   - "minredispatch": minimize the change in generator dispatch from the initial values defined in the network data
- `gen_areas`: integer array of all areas in which generation should be minimized if `obj=="areagen"`
- `area_interface`: two areas with power transfer between them that should be saved and plotted. Results for P, Q, S, and loss are saved for power transfer between the two areas. Must have exactly two elements.
- `gen_zones`: integer array of all zones in which generation should be minimized if `obj=="zonegen"`
- `zone_interface`: two zones with power transfer between them that should be saved and plotted. Results for P, Q, S, and loss are saved for power transfer between the two zones. Must have exactly two elements.
- `print_results::Bool`: if true, print the DataFrames containing the output values for buses, branches, generators, and interfaces. These values are always saved to the output *.csv* files whether true or false.
- `override_param::Dict{Any}`: values to override in the network data defined in `folder`. Must follow the same structure as the full network data dictionary, beginning with key "sn". Default empty Dict.
- `fix_f_override::Bool`: if true, fix the frequency in every subnetwork to the base value, overriding the `variable_f` parameter to `variable_f=false` for every subnetwork. Default false.
- `direct_pq::Bool`: If direct_pq is false, then the interface is treated as a single node and power flow respects Kirchoff Laws, by constraining the voltage magnitude and angle on each side to be equal and enforcing reactive power balance. Default true.
- `master_subnet::Int64`: if `direct_pq==false`, the angle reference must be defined for exactly one subnetwork, since the other subnetwork angles are coupled through the interfaces. Value of `master_subnet` defines which subnetwork provides this reference. Default 1.
- `suffix::String`: suffix to add to the output directory when saving results. Default empty string.
- `start_vals`: Nested dictionary populated with values to be used as a starting point in the optimization model. Applies to bus `vm` and `va`, gen `pg` and `qg`, branch `pt`, `pf`, `qt` and `qf` and subnet `f`. Any of these values which are present in the dictionary will be applied; other values will be ignored. A full network data dictionary can be used. Default `Dict{String, Dict}("sn"=>Dict())`.
"""
function multifrequency_opf(
   folder::String,
   obj::String;
   gen_areas=[],
   area_interface=[],
   gen_zones=[],
   zone_interface=[],
   print_results::Bool=false,
   override_param::Dict=Dict(),
   fix_f_override::Bool=false,
   direct_pq::Bool=true,
   master_subnet::Int64=1,
   suffix::String="",
   start_vals=Dict{String, Dict}("sn"=>Dict()),
   no_converter_loss::Bool=false,
   uniform_gen_scaling::Bool=false,
   unbounded_pg::Bool=false,
   output_to_files::Bool=true,
   output_location_base::String="",
   output_top_folder::String="",
   regularize_f::Float64=0.0,
   ipopt_max_iter::Int64=10000
   )

   println("read_sn_data($folder)")
   mn_data = read_sn_data(folder, no_converter_loss=no_converter_loss)

   folder = abspath(folder)
   if length(output_location_base) > 0
      output_folder = output_location_base
   else
      folder_split = splitpath(folder)
      toplevels = folder_split[1:end-3]
      output_folder = joinpath(toplevels...,"results/$(folder_split[end-1])/$(folder_split[end])$suffix")
   end
   output_folder = joinpath(output_folder, output_top_folder)
   if !isdir(output_folder)
      mkpath(output_folder)
   end

   # toplevel = split(folder,"/")[1]
   # output_folder = join([(i==toplevel ? "results/" : "$i/") for i in split(folder,"/")])*suffix
   println("output folder: $output_folder")
   if !isdir(output_folder)
      mkpath(output_folder)
   end

   (output_dict, res_summary, solution_pm, binding_cnstr_dict) = multifrequency_opf(
      mn_data,
      output_folder,
      obj,
      gen_areas,
      area_interface,
      gen_zones,
      zone_interface,
      print_results,
      override_param,
      fix_f_override,
      direct_pq,
      master_subnet,
      start_vals,
      uniform_gen_scaling,
      unbounded_pg,
      output_to_files,
      regularize_f=regularize_f,
      ipopt_max_iter=ipopt_max_iter
      )

      return (output_dict, res_summary, solution_pm, binding_cnstr_dict)
   end

   function multifrequency_opf(
      folder::String,
      obj::String,
      gen_areas,
      area_interface,
      gen_zones,
      zone_interface,
      print_results::Bool,
      override_param::Dict=Dict(),
      fix_f_override::Bool=false,
      direct_pq::Bool=true,
      master_subnet::Int64=1,
      suffix::String="",
      start_vals=Dict{String, Dict}("sn"=>Dict());
      regularize_f::Float64=0.0,
      ipopt_max_iter::Int64=10000
      )

      (output_dict, res_summary, solution_pm, binding_cnstr_dict) = multifrequency_opf(
         folder,
         obj;
         gen_areas=gen_areas,
         area_interface=area_interface,
         gen_zones=gen_zones,
         zone_interface=zone_interface,
         print_results=print_results,
         override_param=override_param,
         fix_f_override=fix_f_override,
         direct_pq=direct_pq,
         master_subnet=master_subnet,
         suffix=suffix,
         start_vals=start_vals,
         regularize_f=regularize_f,
         ipopt_max_iter=ipopt_max_iter
      )

      return (output_dict, res_summary, solution_pm, binding_cnstr_dict)
   end

   function multifrequency_opf(
      mn_data::Dict,
      output_folder::String,
      obj::String,
      gen_areas=[],
      area_interface=[],
      gen_zones=[],
      zone_interface=[],
      print_results::Bool=false,
      override_param::Dict=Dict(),
      fix_f_override::Bool=false,
      direct_pq::Bool=true,
      master_subnet::Int64=1,
      start_vals=Dict{String, Dict}("sn"=>Dict()),
      uniform_gen_scaling::Bool=false,
      unbounded_pg::Bool=false,
      output_to_files::Bool=true;
      regularize_f::Float64=0.0,
      ipopt_max_iter::Int64=10000
      )

   # If direct_pq is false, then interface flow respects Kirchoff
   # Laws, constraining the voltage magnitude and angle on each side to be equal
   # and enforcing reactive power balance.
   # If false, the angle reference must be defined for only
   # one subnetwork, since the other subnetwork angles are coupled through the
   # interfaces. This "master subnet" is defined by input master_subnet

   # If fix_f_override is true, all subnetwork frequencies will be fixed,
   # regardless of what's specified in the subnetworks file

   # Instantiate a Solver
   #---------------------
   # pass options to IPOPT
   # make sure IPOPT logs to file so we can grep time, residuals and number of iterations
   ipopt_print_level = 0
   ipopt_log_to_file = true
   ipopt_file_log_level = 3
   ipopt_log_file = "$output_folder/ipopt_log"
   if ipopt_log_to_file
      0 < ipopt_file_log_level < 3 && @warn("`file_print_level` should be 0 or ≥ 3 for IPOPT to report elapsed time, final residuals and number of iterations")
   end
   # Callback
   # callback === nothing || setIntermediateCallback(model, callback)
   nlp_optimizer = optimizer_with_attributes(
      Ipopt.Optimizer,
      "print_level" =>  ipopt_print_level,
      "output_file"=> ipopt_log_file,
      "file_print_level" => ipopt_file_log_level,
      "max_iter" => ipopt_max_iter,
      # "derivative_test" => "first-order",
      # "print_info_string" => "yes"
   )
   # nlp_optimizer = optimizer_with_attributes(
   #    () -> AmplNLWriter.Optimizer(Ipopt_jll.amplexe),
   #    "print_level" =>  ipopt_print_level,
   #    "output_file"=> ipopt_log_file,
   #    "file_print_level" => ipopt_file_log_level,
   #    "max_iter" => ipopt_max_iter,
   #    # "halt_on_ampl_error" => "yes",
   #    # "derivative_test" => "first-order",
   #    # "print_info_string" => "yes"
   # )
   PowerModels.logger_config!("error")
   # nlp_optimizer = with_optimizer(Ipopt.Optimizer)
   # note: print_level changes the amount of solver information printed to the terminal

      # Add fix_f_override to override_params if true
   if fix_f_override
      for (subnet_idx, subnet) in mn_data["sn"]
         set_nested!(override_param, ["sn","$subnet_idx","variable_f"], false)
      end
   end

   # if any override parameter subnetwork index has a value of -1, apply it to all variable frequency networks
   if "sn" in keys(override_param)
      if "-1" in keys(override_param["sn"])
         subdict = pop!(override_param["sn"],"-1")
         for (subnet_idx, subnet) in mn_data["sn"]
            if subnet["variable_f"]
               override_param["sn"]["$subnet_idx"] = deepcopy(subdict)
               # println("Set override of subnet $subnet_idx to $subdict")
            end
         end
      end
   end

   # if any override parameters exist in the override_param dictionary, apply them
   combine_nested!(mn_data, override_param)

   if output_to_files
      stringnet = JSON.json(mn_data)
      open("$output_folder/network_data.json", "w") do f
         write(f, stringnet)
      end
   end
   # use build_ref to filter out inactive components
   ref = Dict{Int64,Any}()
   for (subnet_idx,subnet) in mn_data["sn"]
      subnet_idx_int = parse(Int64,subnet_idx)
      ref_temp = PowerModels.build_ref(subnet)
      if :it in keys(ref_temp)
         ref[subnet_idx_int] = ref_temp[:it][:pm][:nw][0]
      else
         ref[subnet_idx_int] = ref_temp[:nw][0]
      end

   end

   # note: ref contains all the relevant system parameters needed to build the OPF model
   # and it is a Dict with the ref data for each subnetwork
   # When we introduce constraints and variable bounds below, we use the parameters in ref.

   # Export ref to file for debugging
   if output_to_files
      # println("Set ref data. Saving JSON.")
      stringref = JSON.json(ref)
      open("$output_folder/ref_data.json", "w") do f
         write(f, stringref)
      end
   end

   # Use this to print any generators which have intial values outside their limits
   print_init_gen_viol = false
   if print_init_gen_viol
      for (subnet_idx, ref_subnet) in ref
         for (i,gen) in ref_subnet[:gen]
            if gen["pg"] > gen["pmax"]
               println("gen $(gen["index"]) at bus $(gen["gen_bus"]): pg = $(gen["pg"]);\t pmax = $(gen["pmax"])\t (status = $(gen["gen_status"]))")
            end
         end
      end
   end

   # print_summary(mn_data["sn"]["1"])
   ###############################################################################
   # 1. Building the Optimal Power Flow Model
   ###############################################################################

   # Initialize a JuMP Optimization Model
   #-------------------------------------
   t_start = time_ns()
   model = Model(nlp_optimizer)

   va = Dict{Int64,Any}()
   vm = Dict{Int64,Any}()
   pg = Dict{Int64,Any}()
   qg = Dict{Int64,Any}()
   p = Dict{Int64,Any}()
   q = Dict{Int64,Any}()
   p_dc = Dict{Int64,Any}()
   q_dc = Dict{Int64,Any}()
   g = Dict{Int64,Any}()
   b = Dict{Int64,Any}()
   b_fr = Dict{Int64,Any}()
   b_to = Dict{Int64,Any}()
   g_fr = Dict{Int64,Any}()
   g_to = Dict{Int64,Any}()
   l_series = Dict{Int64,Any}()
   c_inv_series = Dict{Int64,Any}()
   f = Dict{Int64,Any}()
   b_shunt = Dict{Int64,Any}()
   cost_pg = Dict{Int64,Any}()
   cost_dcline = Dict{Int64,Any}()
   area_pg = Dict{Int64,Any}()
   zone_pg = Dict{Int64,Any}()
   redispatch_pg = Dict{Int64,Any}()
   redispatch_qg = Dict{Int64,Any}()
   redispatch_vm = Dict{Int64,Any}()
   redispatch_va = Dict{Int64,Any}()
   p_i = Dict{Int64,Any}()
   q_i = Dict{Int64,Any}()
   s_i = Dict{Int64,Any}()
   i_rms_sq = Dict{Int64,Any}()
   i_mabs = Dict{Int64,Any}()
   i_dc = Dict{Int64,Any}()
   s_M = Dict{Int64,Any}()
   conv_loss = Dict{Int64,Any}()

   # constraint references
   constraints = Dict{Symbol,Any}(
      :vm_ulim     => Dict{Int64,Any}(),
      :vm_llim     => Dict{Int64,Any}(),
      :pg_ulim     => Dict{Int64,Any}(),
      :pg_llim     => Dict{Int64,Any}(),
      :qg_ulim     => Dict{Int64,Any}(),
      :qg_llim     => Dict{Int64,Any}(),
      :s_ulim      => Dict{Int64,Any}(),
      :s_llim      => Dict{Int64,Any}(),
      :theta_llim  => Dict{Int64,Any}(),
      :theta_ulim  => Dict{Int64,Any}(),
      :theta_ref  => Dict{Int64,Any}()
   )

   for (subnet_idx, ref_subnet) in ref
      add_vars!(
         mn_data,
         ref_subnet,
         subnet_idx,
         model,
         va, vm, pg, qg, p, q, p_dc, q_dc, g, b, b_fr, b_to, g_fr, g_to, l_series,
         c_inv_series, f, b_shunt, cost_pg, cost_dcline, area_pg, zone_pg, redispatch_pg,
         redispatch_qg, redispatch_vm, redispatch_va, p_i, q_i,
         constraints,
         )
   end
   # Add converter interface variables and constraints
   for (conv_idx, conv_params) in mn_data["converter"]
      p_i[conv_idx] = @variable(
         model,
         [(subnet, bus) in conv_params["converter_buses"]],
         base_name = "p_i_$(conv_idx)",
      )

      q_i[conv_idx] = @variable(
         model,
         [(subnet, bus) in conv_params["converter_buses"]],
         base_name = "q_i_$(conv_idx)",
      )
      for (subnet, bus) in conv_params["converter_buses"]
         smax = ref[subnet][:bus][bus]["converter_vmax"][conv_idx]*ref[subnet][:bus][bus]["converter_imax"][conv_idx]
         set_lower_bound(p_i[conv_idx][(subnet, bus)], -smax)
         set_upper_bound(p_i[conv_idx][(subnet, bus)], smax)
         set_lower_bound(q_i[conv_idx][(subnet, bus)], -smax)
         set_upper_bound(q_i[conv_idx][(subnet, bus)], smax)
      end

      # Keep track of DC subnets connected to the converter in conv_dc_subnets
      # so that they can be excluded from the loss model (the loss parameters
      # for an AC-DC converter should be applied only on the AC side)
      conv_dc_subnets = []
      for (subnet, bus) in conv_params["converter_buses"]
         dc_subnet = ((!ref[subnet][:variable_f]) && (ref[subnet][:f_base] == 0))
         if (:f_max in keys(ref[subnet]))
            # Also DC if max frequency is 0
            dc_subnet = dc_subnet || (ref[subnet][:f_max] == 0)
         end
         if dc_subnet
            push!(conv_dc_subnets, subnet)
         end
      end

      i_rms_sq[conv_idx] = @variable(
         model,
         [(subnet, bus) in conv_params["converter_buses"]],
         base_name = "i_rms_sq_$(conv_idx)",
         lower_bound = 0
      )
      for (subnet, bus) in conv_params["converter_buses"]
         if !(subnet in conv_dc_subnets)
            # if (ref[subnet][:bus][bus]["converter_c1"][conv_idx] == 0) &&
            #       (ref[subnet][:bus][bus]["converter_sw1"][conv_idx] == 0)
            #    @constraint(
            #       model,
            #       i_rms_sq[conv_idx][(subnet, bus)] == 0
            #    )
            # else
               M = ref[subnet][:bus][bus]["converter_M"][conv_idx]
               p_im = p_i[conv_idx][(subnet, bus)]
               q_im = q_i[conv_idx][(subnet, bus)]
               v_i = vm[subnet][bus]
               @NLconstraint(
                  model,
                  i_rms_sq[conv_idx][(subnet, bus)] == (
                     M^2*p_im^2 / (18*v_i^2)
                     + (p_im^2 + q_im^2) / (36*v_i^2)
                     )
               )
            # end
         end
      end

      # Interface power balance, including converter loss as a linear constraint if all c1, c2, c3, sw1, sw2, sw3 are not all zero
      if all(all(values(ref[subnet][:bus][bus]["converter_c1"][conv_idx]) .== 0) for (subnet, bus) in conv_params["converter_buses"]) &&
            all(all(values(ref[subnet][:bus][bus]["converter_c2"][conv_idx]) .== 0) for (subnet, bus) in conv_params["converter_buses"]) &&
            all(all(values(ref[subnet][:bus][bus]["converter_c3"][conv_idx]) .== 0) for (subnet, bus) in conv_params["converter_buses"]) &&
            all(all(values(ref[subnet][:bus][bus]["converter_sw1"][conv_idx]) .== 0) for (subnet, bus) in conv_params["converter_buses"]) &&
            all(all(values(ref[subnet][:bus][bus]["converter_sw2"][conv_idx]) .== 0) for (subnet, bus) in conv_params["converter_buses"]) &&
            all(all(values(ref[subnet][:bus][bus]["converter_sw3"][conv_idx]) .== 0) for (subnet, bus) in conv_params["converter_buses"]) &&
            false
         @constraint(
            model,
            sum(p_i[conv_idx][(subnet, bus)]
                for (subnet, bus) in conv_params["converter_buses"] if !(subnet in conv_dc_subnets)) == 0
         )
      else
         i_mabs[conv_idx] = @variable(
            model,
            [(subnet, bus) in conv_params["converter_buses"]],
            base_name = "i_mabs_$(conv_idx)",
            lower_bound = 0
         )

         if !(conv_idx in keys(s_M))
            s_M[conv_idx] = Dict{Any,Any}()
         end
         if !(conv_idx in keys(s_i))
            s_i[conv_idx] = Dict{Any,Any}()
         end
         for (subnet, bus) in conv_params["converter_buses"]
            if !(subnet in conv_dc_subnets)
               if (ref[subnet][:bus][bus]["converter_c2"][conv_idx] == 0) &&
                     (ref[subnet][:bus][bus]["converter_sw2"][conv_idx] == 0)
                  @constraint(
                     model,
                     i_mabs[conv_idx][(subnet, bus)] == 0
                  )
               else
                  s_i[conv_idx][(subnet, bus)] = @variable(
                        model,
                        base_name = "s_i_$(conv_idx)_$(subnet)_$(bus)",
                        lower_bound = 0,
                     )
                  @constraint(
                     model,
                     p_i[conv_idx][(subnet, bus)]^2 + q_i[conv_idx][(subnet, bus)]^2 == s_i[conv_idx][(subnet, bus)]^2
                  )
                  M = ref[subnet][:bus][bus]["converter_M"][conv_idx]
                  p_im = p_i[conv_idx][(subnet, bus)]
                  q_im = q_i[conv_idx][(subnet, bus)]
                  v_i = vm[subnet][bus]
                  if M == 1
                     @NLconstraint(
                        model,
                        i_mabs[conv_idx][(subnet, bus)] ==
                        2/(3*v_i*sqrt(2))*(
                           p_im^2/s_i[conv_idx][(subnet, bus)] + abs(q_im)
                           )
                     )
                  else
                     s_M[conv_idx][(subnet, bus)] = @variable(
                           model,
                           base_name = "s_M_$(conv_idx)_$(subnet)_$(bus)",
                           lower_bound = 0,
                           upper_bound = ref[subnet][:bus][bus]["converter_vmax"][conv_idx]*ref[subnet][:bus][bus]["converter_imax"][conv_idx],
                        )
                     @constraint(
                        model,
                        p_i[conv_idx][(subnet, bus)]^2*(1-ref[subnet][:bus][bus]["converter_M"][conv_idx]^2) + q_i[conv_idx][(subnet, bus)]^2 == s_M[conv_idx][(subnet, bus)]^2
                     )
                     @NLconstraint(
                        model,
                        s_i[conv_idx][(subnet, bus)]*i_mabs[conv_idx][(subnet, bus)] ==
                        2/(3*v_i*sqrt(2))*(
                           M*p_im*(M*p_im) + s_i[conv_idx][(subnet, bus)]*s_M[conv_idx][(subnet, bus)]
                           )
                     )
                  end
               end
            end
         end

         i_dc[conv_idx] = @variable(
            model,
            [(subnet, bus) in conv_params["converter_buses"]],
            base_name = "i_dc$(conv_idx)",
            lower_bound = 0
         )
         for (subnet, bus) in conv_params["converter_buses"]
            if !(subnet in conv_dc_subnets)
               if (ref[subnet][:bus][bus]["converter_c3"][conv_idx] == 0) &&
                     (ref[subnet][:bus][bus]["converter_sw3"][conv_idx] == 0)
                  @constraint(
                     model,
                     i_dc[conv_idx][(subnet, bus)] == 0
                  )
               else
                  M = ref[subnet][:bus][bus]["converter_M"][conv_idx]
                  p_im = p_i[conv_idx][(subnet, bus)]
                  v_i = vm[subnet][bus]
                  @NLconstraint(
                     model,
                     i_dc[conv_idx][(subnet, bus)] == M*abs(p_im) / (sqrt(2) * v_i)
                  )
               end
            end
         end

         conv_loss[conv_idx] = @variable(
            model,
            base_name = "conv_loss_$(conv_idx)",
            lower_bound = 0
         )

         # Interface power balance
         @constraint(
            model,
            conv_loss[conv_idx] == 6*(
               sum(ref[subnet][:bus][bus]["converter_c1"][conv_idx]*i_rms_sq[conv_idx][(subnet, bus)] for (subnet, bus) in conv_params["converter_buses"])
               + sum(ref[subnet][:bus][bus]["converter_c2"][conv_idx]*i_mabs[conv_idx][(subnet, bus)] for (subnet, bus) in conv_params["converter_buses"])
               + sum(ref[subnet][:bus][bus]["converter_c3"][conv_idx]*i_dc[conv_idx][(subnet, bus)] for (subnet, bus) in conv_params["converter_buses"])
               + sum(ref[subnet][:bus][bus]["converter_sw1"][conv_idx]*i_rms_sq[conv_idx][(subnet, bus)] for (subnet, bus) in conv_params["converter_buses"])
               + sum(ref[subnet][:bus][bus]["converter_sw2"][conv_idx]*i_mabs[conv_idx][(subnet, bus)] for (subnet, bus) in conv_params["converter_buses"])
               + sum(ref[subnet][:bus][bus]["converter_sw3"][conv_idx] for (subnet, bus) in conv_params["converter_buses"])
            )
         )

         @constraint(
            model,
            sum(p_i[conv_idx][(subnet, bus)] for (subnet, bus) in conv_params["converter_buses"]) + conv_loss[conv_idx] == 0
         )
      end
      # Converter current limits
      for (subnet, bus) in conv_params["converter_buses"]
         i_lim = ref[subnet][:bus][bus]["converter_imax"][conv_idx]
         if i_lim < Inf
            @constraint(
               model,
               i_rms_sq[conv_idx][(subnet, bus)] - i_lim^2*vm[subnet][bus]^2 <= 0
            )
         end
      end

      # Interface coupling of vm and va to model control of frequency alone (no direct PQ control)
      ############################################################################################
      if !direct_pq
         println("No direct PQ control")
         @constraint(
            model,
            sum(q_i[conv_idx[1]][(subnet, bus)]
                for (subnet, bus) in conv_params["converter_buses"]) == 0
         )
         for ((subnet1, bus1), (subnet2, bus2)) in combinations(conv_params["converter_buses"], 2)

            @constraint(
               model,
               vm[subnet1][bus1] == vm[subnet2][bus2]
            )
            @constraint(
               model,
               va[subnet1][bus1] == va[subnet2][bus2]
            )
         end
      end
      ############################################################################################

   end

   # Set apparent power limit for a 0 Hz system based on DC current,
   # incorporating constants k_ins and k_cond if true
   # If false, apply the same apparent power limit to the active power
   dc_current_limit = false

   for (subnet_idx, ref_subnet) in ref

      add_constraints!(
         mn_data,
         ref_subnet,
         subnet_idx,
         model,
         va, vm, pg, qg, p, q, p_dc, q_dc, g, b, b_fr, b_to, g_fr, g_to, l_series,
         c_inv_series, f, b_shunt, cost_pg, cost_dcline, area_pg, zone_pg, redispatch_pg,
         redispatch_qg, redispatch_vm, redispatch_va, p_i, q_i,
         constraints,
         obj, gen_areas, area_interface,
         gen_zones, zone_interface,
         direct_pq,
         dc_current_limit,
         master_subnet
         )

   end

   if uniform_gen_scaling
      println("APPLYING UNIFORM GENERATOR SCALING")
      (alpha_upstream,
      alpha_downstream) = add_gen_zone_scaling_constraint!(
         ref,
         model,
         pg,
         constraints,
         unbounded_pg,
         gen_zones
         )
   end

   # Set starting values if given
   # ----------------------------
   set_startvals!(
      start_vals,
      va, vm, pg, qg, p, q, f
      )

   # Add Objective Function
   # ----------------------
   if regularize_f > 0.0
      f_regularizer = @variable(
         model,
         base_name = "f_regularizer",
         lower_bound = 0.0
      )
      @constraint(
         model,
         f_regularizer >= sum(regularize_f/ref_subnet[:f_base]^2*(ref_subnet[:f_base]-f[subnet_idx])^2 for (subnet_idx, ref_subnet) in ref if ref_subnet[:variable_f])
      )
   end
   if obj=="mincost"
      # Minimize the cost of active power generation and cost of HVDC line usage
      if regularize_f > 0.0
         @objective(model, Min,
            sum(cost_pg[subnet_idx] for (subnet_idx, ref_subnet) in ref) +
            sum(cost_dcline[subnet_idx] for (subnet_idx, ref_subnet) in ref) +
            f_regularizer
         )
      else
         @objective(model, Min,
            sum(cost_pg[subnet_idx] for (subnet_idx, ref_subnet) in ref) +
            sum(cost_dcline[subnet_idx] for (subnet_idx, ref_subnet) in ref)
         )
      end
   elseif obj=="minredispatch"
      # Minimize the deviations from the base case dispatch
      if regularize_f > 0
      @objective(model, Min,
         sum(redispatch_pg[subnet_idx] for (subnet_idx, ref_subnet) in ref) +
         # sum(redispatch_qg[subnet_idx] for (subnet_idx, ref_subnet) in ref) +
         1e9*sum(redispatch_vm[subnet_idx] for (subnet_idx, ref_subnet) in ref) +
         # sum(redispatch_va[subnet_idx] for (subnet_idx, ref_subnet) in ref)
         f_regularizer
      )
      else
         @objective(model, Min,
            sum(redispatch_pg[subnet_idx] for (subnet_idx, ref_subnet) in ref) +
            # sum(redispatch_qg[subnet_idx] for (subnet_idx, ref_subnet) in ref) +
            1e9*sum(redispatch_vm[subnet_idx] for (subnet_idx, ref_subnet) in ref)
            # sum(redispatch_va[subnet_idx] for (subnet_idx, ref_subnet) in ref)
         )
      end
   elseif obj=="areagen"
      # Minimize the generation in the specified areas
      if regularize_f > 0.0
         @objective(model, Min,
            sum(area_pg[subnet_idx] for (subnet_idx, ref_subnet) in ref) +
            f_regularizer
         )
      else
         @objective(model, Min,
            sum(area_pg[subnet_idx] for (subnet_idx, ref_subnet) in ref)
         )
      end
   elseif obj=="zonegen"
      # Minimize the generation in the specified zones
      if regularize_f > 0.0
         @objective(model, Min,
            sum(zone_pg[subnet_idx] for (subnet_idx, ref_subnet) in ref) +
            f_regularizer
         )
      else
         @objective(model, Min,
            sum(zone_pg[subnet_idx] for (subnet_idx, ref_subnet) in ref)
         )
      end
   else
      println("The specified objective has not been implemented. Using minumum cost as the objective.")
      obj = "mincost"
      # Minimize the cost of active power generation and cost of HVDC line usage
      if regularize_f > 0.0
         @objective(model, Min,
            sum(cost_pg[subnet_idx] for (subnet_idx, ref_subnet) in ref) +
            sum(cost_dcline[subnet_idx] for (subnet_idx, ref_subnet) in ref) +
            f_regularizer
         )
      else
         @objective(model, Min,
            sum(cost_pg[subnet_idx] for (subnet_idx, ref_subnet) in ref) +
            sum(cost_dcline[subnet_idx] for (subnet_idx, ref_subnet) in ref)
         )
      end
   end

   #########################################
   # 2. Solve the Optimal Power Flow Model #
   #########################################
   # Solve the optimization problem
   # io = open("$output_folder/model.nl","w")
   # mywrite_nl_file(io,model)
   # close(io)
   open("$output_folder/model.lp", "w") do f
      print(f, model)
   end

   # open("$output_folder/variables.txt", "w") do f
   #    all_indices = MOI.get(model, MOI.ListOfVariableIndices())::Vector{MOI.VariableIndex}
   #    for idx in all_indices
   #       v = VariableRef(model, idx)
   #       print(f, "$idx: $v\n")
   #    end
   # end
   #
   # println()
   #
   # open("$output_folder/constraints.txt", "w") do f
   #    for (f_type,set_type) in JuMP.list_of_constraint_types(model)
   #       f_type = moi_function_type(f_type)
   #       # set_type = MathOptInterface.LessThan{Float64}
   #       for idx in MathOptInterface.get(model, MathOptInterface.ListOfConstraintIndices{f_type,set_type}())
   #          cnstr = JuMP.constraint_ref_with_index(model, idx)
   #          print(f, "$idx: $cnstr\n")
   #       end
   #    end
   # end

   optimize!(model)

   t_end = time_ns()
   t_build_solve = (t_end - t_start)/1e9
   # Check that the solver terminated without an error
   println("The solver termination status is $(termination_status(model))")

   ipopt_output = readlines(ipopt_log_file)

   cpu_time = 0.0
   dual_feas = primal_feas = Inf
   iter = -1
   for line in ipopt_output
      if occursin("CPU secs", line)
         cpu_time += Meta.parse(split(line, "=")[2])
      elseif occursin("Dual infeasibility", line)
         dual_feas = Meta.parse(split(line)[4])
      elseif occursin("Constraint violation", line)
         primal_feas = Meta.parse(split(line)[4])
      elseif occursin("Number of Iterations....", line)
         iter = Meta.parse(split(line)[4])
      end
   end

   #########################
   # 3. Review the Results #
   #########################

   # Save the resulting setpoints in PowerModels format
   solution_pm = deepcopy(mn_data)
   for (subnet_idx, ref_subnet) in ref
      dc_subnet = ((!ref_subnet[:variable_f]) && (ref_subnet[:f_base] == 0))
      if (:f_max in keys(ref_subnet))
         # Also DC if max frequency is 0
         dc_subnet = dc_subnet || (ref_subnet[:f_max] == 0)
      end

      for (i,gen) in ref_subnet[:gen]
         solution_pm["sn"]["$subnet_idx"]["gen"]["$i"]["pg"] = value(pg[subnet_idx][gen["index"]])
         if !dc_subnet
            solution_pm["sn"]["$subnet_idx"]["gen"]["$i"]["qg"] = value(qg[subnet_idx][gen["index"]])
         end
      end
      for (i,bus) in ref_subnet[:bus]
         solution_pm["sn"]["$subnet_idx"]["bus"]["$i"]["vm"] = value(vm[subnet_idx][bus["bus_i"]])
         if !dc_subnet
            solution_pm["sn"]["$subnet_idx"]["bus"]["$i"]["va"] = value(va[subnet_idx][bus["bus_i"]])
         end
      end
      for (i,branch) in ref_subnet[:branch]
         f_idx = (i, branch["f_bus"], branch["t_bus"])
         t_idx = (i, branch["t_bus"], branch["f_bus"])
         solution_pm["sn"]["$subnet_idx"]["branch"]["$i"]["pf"] = value(p[subnet_idx][f_idx])
         solution_pm["sn"]["$subnet_idx"]["branch"]["$i"]["pt"] = value(p[subnet_idx][t_idx])
         if !dc_subnet
            solution_pm["sn"]["$subnet_idx"]["branch"]["$i"]["qf"] = value(q[subnet_idx][f_idx])
            solution_pm["sn"]["$subnet_idx"]["branch"]["$i"]["qt"] = value(q[subnet_idx][t_idx])
            va_fr = va[subnet_idx][branch["f_bus"]]
            va_to = va[subnet_idx][branch["t_bus"]]
            solution_pm["sn"]["$subnet_idx"]["branch"]["$i"]["angle"] = value(va_fr) - value(va_to)
         end
      end
      for (i,dcline) in ref_subnet[:dcline]
         f_idx = (i, dcline["f_bus"], dcline["t_bus"])
         t_idx = (i, dcline["t_bus"], dcline["f_bus"])
         solution_pm["sn"]["$subnet_idx"]["dcline"]["$i"]["pf"] = value(p_dc[subnet_idx][f_idx])
         solution_pm["sn"]["$subnet_idx"]["dcline"]["$i"]["pt"] = value(p_dc[subnet_idx][t_idx])
         solution_pm["sn"]["$subnet_idx"]["dcline"]["$i"]["qf"] = value(q_dc[subnet_idx][f_idx])
         solution_pm["sn"]["$subnet_idx"]["dcline"]["$i"]["qt"] = value(q_dc[subnet_idx][t_idx])
      end
   end
   # Export network to file for debugging
   if output_to_files
      # println("Wrote solution to PowerModels network. Saving JSON.")
      stringnet = JSON.json(solution_pm)
      open("$output_folder/network_solved.json", "w") do f
         write(f, stringnet)
      end
   end

   # Check the value of the objective function
   objval = objective_value(model)
   println("The objective value is $(objval).")

   println("CPU time (s): $cpu_time")
   println("Dual infeasibility: $dual_feas")
   println("Constraint violation: $primal_feas")
   println("Iterations: $iter")

   # cost = (sum(value(cost_pg[subnet_idx]) for (subnet_idx, ref_subnet) in ref) +
   #         sum(value(cost_dcline[subnet_idx]) for (subnet_idx, ref_subnet) in ref)
   #         )

   # total_load = sum(sum(ref_subnet[:load][l] for l in ref_subnet[:bus_loads][i]) for (i,bus) in ref_subnet[:bus])
   res_summary = Dict{Int64, DataFrame}()

   # Define arrays for summary table
   n_subnets = length(ref)
   subnets_arr=Array{Int64}(undef, n_subnets)
   f_arr=Array{Float64}(undef, n_subnets)
   fmin_arr=Array{Float64}(undef, n_subnets)
   fmax_arr=Array{Float64}(undef, n_subnets)
   min_v_arr=Array{Float64}(undef, n_subnets)
   max_v_arr=Array{Float64}(undef, n_subnets)
   avg_v_arr=Array{Float64}(undef, n_subnets)
   ngen_Qll_arr=Array{Int64}(undef, n_subnets)
   ngen_Qul_arr=Array{Int64}(undef, n_subnets)
   ngen_arr=Array{Int64}(undef, n_subnets)
   minQloss_arr=Array{Float64}(undef, n_subnets)
   maxQloss_arr=Array{Float64}(undef, n_subnets)
   n_negQloss_arr=Array{Int64}(undef, n_subnets)
   n_s_bound_arr=Array{Int64}(undef, n_subnets)
   max_ang_diff_arr=Array{Float64}(undef, n_subnets)
   min_ang_diff_arr=Array{Float64}(undef, n_subnets)
   nbranch_arr=Array{Int64}(undef, n_subnets)
   loss_arr=Array{Float64}(undef, n_subnets)
   cost_arr=Array{Float64}(undef, n_subnets)
   obj_arr=Array{Float64}(undef, n_subnets)
   time_s_arr=Array{Float64}(undef, n_subnets)
   status_arr=Array{TerminationStatusCode}(undef, n_subnets)
   subnet_load=Array{Float64}(undef, n_subnets)
   subnet_generation=Array{Float64}(undef, n_subnets)
   subnet_line_loss=Array{Float64}(undef, n_subnets)
   subnet_total_loss=Array{Float64}(undef, n_subnets)
   converter_loss=Array{Float64}(undef, n_subnets)

   summary_idx = 1
   for (subnet_idx, ref_subnet) in ref
      dc_subnet = ((!ref_subnet[:variable_f]) && (ref_subnet[:f_base] == 0))
      if (:f_max in keys(ref_subnet))
         # Also DC if max frequency is 0
         dc_subnet = dc_subnet || (ref_subnet[:f_max] == 0)
      end

      print("\nsubnetwork $(subnet_idx)\n================\n")
      if :f_fixed in keys(ref_subnet)
         freq_res = ref_subnet[:f_fixed]
         freq_min = ref_subnet[:f_min]
         freq_max = ref_subnet[:f_max]
      elseif ref_subnet[:variable_f]
         freq_res = value(f[subnet_idx])
         freq_min = ref_subnet[:f_min]
         freq_max = ref_subnet[:f_max]
      else
         freq_res = ref_subnet[:f_base]
         freq_min = ref_subnet[:f_base]
         freq_max = ref_subnet[:f_base]
      end
      println("frequency: $(freq_res) Hz\n================\n")
      # Iterate over buses
      bus_load = Dict{Int64,Any}()
      # gen_string = "\ngen output\ngen,\tbus,\tp,\tq\n"
      # bus_string = "\nbus voltage\nbus,\tvm,\tva\n"

      # Define arrays for bus table
      n_bus = length(ref_subnet[:bus])
      bus_arr = Array{Int64}(undef, n_bus)
      vm_arr = Array{Float64}(undef, n_bus)
      va_arr = Array{Float64}(undef, n_bus)
      vmin_arr = Array{Float64}(undef, n_bus)
      vmax_arr = Array{Float64}(undef, n_bus)

      # Define arrays for gen table
      n_gen = length(ref_subnet[:gen])
      gen_arr = Array{Int64}(undef, n_gen)
      gen_bus_arr = Array{Int64}(undef, n_gen)
      p_arr = Array{Float64}(undef, n_gen)
      pmin_arr = Array{Float64}(undef, n_gen)
      pmax_arr = Array{Float64}(undef, n_gen)
      q_arr = Array{Float64}(undef, n_gen)
      qmin_arr = Array{Float64}(undef, n_gen)
      qmax_arr = Array{Float64}(undef, n_gen)
      pbase_arr = Array{Float64}(undef, n_gen)
      qbase_arr = Array{Float64}(undef, n_gen)

      bus_idx = 1
      gen_i = 1
      for (i,bus) in ref_subnet[:bus]
         # println("bus_idx: $bus_idx")
         bus_arr[bus_idx] = bus["index"]
         # println("bus[\"index\"]: $(bus["index"])")
         vm_arr[bus_idx] = value(vm[subnet_idx][i])
         va_arr[bus_idx] = dc_subnet ? NaN : value(va[subnet_idx][i]) * 180/pi
         vmin_arr[bus_idx] = bus["vmin"]
         vmax_arr[bus_idx] = bus["vmax"]
         # bus_string *= "$i,\t$(vm_tmp),\t$(va_tmp)\n"
         busgens = ref_subnet[:bus_gens][i]
         for gen_idx in busgens
            # println("gen_i: $gen_i")
            # println("gen_idx: $gen_idx")
            gen_arr[gen_i] = gen_idx
            # println("bus[\"index\"]: $(bus["index"])")
            gen_bus_arr[gen_i] = bus["index"]
            p_arr[gen_i] = value(pg[subnet_idx][gen_idx])
            q_arr[gen_i] = dc_subnet ? NaN : value(qg[subnet_idx][gen_idx])
            pmin_arr[gen_i] = ref_subnet[:gen][gen_idx]["pmin"]
            pmax_arr[gen_i] = ref_subnet[:gen][gen_idx]["pmax"]
            qmin_arr[gen_i] = ref_subnet[:gen][gen_idx]["qmin"]
            qmax_arr[gen_i] = ref_subnet[:gen][gen_idx]["qmax"]
            pbase_arr[gen_i] = ref_subnet[:gen][gen_idx]["pg"]
            qbase_arr[gen_i] = ref_subnet[:gen][gen_idx]["qg"]
            # gen_string *= "$(gen_idx),\t$i,\t$(pg_tmp),\t$(qg_tmp)\n"
            gen_i += 1
         end
         busloads = ref_subnet[:bus_loads][i]
         if !isempty(busloads)
            bus_load[i] = sum([ref_subnet[:load][l]["pd"] for l in busloads])
         end
         bus_idx += 1
      end
      bus_df = DataFrame(
         bus = bus_arr, vm = vm_arr, va = va_arr, vmin = vmin_arr, vmax = vmax_arr
      )
      gen_df = DataFrame(
         gen = gen_arr, bus = gen_bus_arr, pg = p_arr, qg = q_arr,
         pgmin = pmin_arr, pgmax = pmax_arr, qgmin = qmin_arr, qgmax = qmax_arr,
         pbase = pbase_arr, qbase = qbase_arr
      )

      # Iterate over branches
      n_branch = length(ref_subnet[:branch])

      # Define arrays for branch table
      branch_arr = Array{Int64}(undef, n_branch)
      from_bus_arr = Array{Int64}(undef, n_branch)
      to_bus_arr = Array{Int64}(undef, n_branch)
      B_arr = Array{Float64}(undef, n_branch)
      G_arr = Array{Float64}(undef, n_branch)
      Xseries_arr = Array{Float64}(undef, n_branch)
      Bfr_arr = Array{Float64}(undef, n_branch)
      Bto_arr = Array{Float64}(undef, n_branch)
      Gfr_arr = Array{Float64}(undef, n_branch)
      Gto_arr = Array{Float64}(undef, n_branch)
      rate_a_arr = Array{Float64}(undef, n_branch)
      p_from_arr = Array{Float64}(undef, n_branch)
      p_to_arr = Array{Float64}(undef, n_branch)
      p_loss_arr = Array{Float64}(undef, n_branch)
      q_from_arr = Array{Float64}(undef, n_branch)
      q_to_arr = Array{Float64}(undef, n_branch)
      q_loss_arr = Array{Float64}(undef, n_branch)
      s_from_arr = Array{Float64}(undef, n_branch)
      s_to_arr = Array{Float64}(undef, n_branch)
      angle_diff_arr = Array{Float64}(undef, n_branch)
      angle_min_arr = Array{Float64}(undef, n_branch)
      angle_max_arr = Array{Float64}(undef, n_branch)
      for (i,(idx,branch)) in enumerate(ref_subnet[:branch])
         branch_arr[i] = idx
         from_bus_arr[i] = branch["f_bus"]
         to_bus_arr[i] = branch["t_bus"]
         rate_a_arr[i] = "rate_a" in keys(branch) ? branch["rate_a"] : Inf
         B_arr[i] = typeof(b[subnet_idx][idx])==VariableRef ? value(b[subnet_idx][idx]) : b[subnet_idx][idx]
         G_arr[i] = typeof(g[subnet_idx][idx])==VariableRef ? value(g[subnet_idx][idx]) : g[subnet_idx][idx]
         Xseries_arr[i] = ref_subnet[:variable_f] ? (value(f[subnet_idx])*2pi*l_series[subnet_idx][idx] - c_inv_series[subnet_idx][idx]/(value(f[subnet_idx])*2pi)) : branch["br_x"]

         Bfr_arr[i] = typeof(b_fr[subnet_idx][idx])==VariableRef ? value(b_fr[subnet_idx][idx]) : b_fr[subnet_idx][idx]
         Bto_arr[i] = typeof(b_to[subnet_idx][idx])==VariableRef ? value(b_to[subnet_idx][idx]) : b_to[subnet_idx][idx]
         Gfr_arr[i] = typeof(g_fr[subnet_idx][idx])==VariableRef ? value(g_fr[subnet_idx][idx]) : g_fr[subnet_idx][idx]
         Gto_arr[i] = typeof(g_to[subnet_idx][idx])==VariableRef ? value(g_to[subnet_idx][idx]) : g_to[subnet_idx][idx]
         f_idx = (idx, branch["f_bus"], branch["t_bus"])
         t_idx = (idx, branch["t_bus"], branch["f_bus"])
         p_from_arr[i] = value(p[subnet_idx][f_idx])
         q_from_arr[i] = dc_subnet ? 0 : value(q[subnet_idx][f_idx])
         p_to_arr[i] = value(p[subnet_idx][t_idx])
         q_to_arr[i] = dc_subnet ? 0 : value(q[subnet_idx][t_idx])
         s_from_arr[i] = sqrt(p_from_arr[i]^2 + q_from_arr[i]^2)
         s_to_arr[i] = sqrt(p_to_arr[i]^2 + q_to_arr[i]^2)
         p_loss_arr[i] = abs(p_from_arr[i] + p_to_arr[i])
         q_loss_arr[i] = q_from_arr[i] + q_to_arr[i]
         angle_diff_arr[i] = dc_subnet ? NaN : (value(va[subnet_idx][branch["f_bus"]]) - value(va[subnet_idx][branch["t_bus"]])) * 180/pi
         angle_min_arr[i] = branch["angmin"] * 180/pi
         angle_max_arr[i] = branch["angmax"] * 180/pi
      end
      br_df = DataFrame(
         branch = branch_arr, from_bus = from_bus_arr, to_bus = to_bus_arr,
         X = Xseries_arr, B = B_arr, G = G_arr, Bfr = Bfr_arr, Bto = Bto_arr,
         # Gfr = Gfr_arr, Gto = Gto_arr,
         p_from = p_from_arr,
         p_to = p_to_arr, q_from = q_from_arr, q_to = q_to_arr,
         p_loss = p_loss_arr, q_loss = q_loss_arr, s_from = s_from_arr,
         s_to = s_to_arr, rate_a = rate_a_arr,
         ang_deg = angle_diff_arr, angmin = angle_min_arr, angmax = angle_max_arr,
      )

      # Sum of generation
      subnet_load[summary_idx] = sum(Float64[load for (i,load) in bus_load])
      subnet_generation[summary_idx] = sum(Float64[value(pg[subnet_idx][i]) for i in keys(ref_subnet[:gen])])
      # full_gen += total_generation
      subnet_line_loss[summary_idx] = sum(Float64[abs(value(p[subnet_idx][(i,br["f_bus"],br["t_bus"])]) + value(p[subnet_idx][(i,br["t_bus"],br["f_bus"])])) for (i,br) in ref_subnet[:branch]])
      subnet_total_loss[summary_idx] = subnet_generation[summary_idx] - subnet_load[summary_idx]

      converter_loss[summary_idx] = sum(Float64[value(conv_loss[conv_idx]) for (conv_idx, conv_params) in mn_data["converter"]])

      if print_results
         print("\nbus voltages\n$(bus_df)\n")
         print("\ngeneration\n$(gen_df)\n")
         print("\nbranch values\n$(br_df)\n")

         println("\ntotal load: $(subnet_load)")

         println("total subnet generation: $(subnet_generation)")
         println("total subnet losses: $(subnet_loss)")
      end

      if output_to_files
         CSV.write("$output_folder/res_bus_$subnet_idx.csv", bus_df)
         CSV.write("$output_folder/res_gen_$subnet_idx.csv", gen_df)
         CSV.write("$output_folder/res_branch_$subnet_idx.csv", br_df)
      end

      subnets_arr[summary_idx] = subnet_idx
      f_arr[summary_idx] = freq_res
      fmin_arr[summary_idx] = freq_min
      fmax_arr[summary_idx] = freq_max
      min_v_arr[summary_idx]=min(vm_arr...)
      max_v_arr[summary_idx]=max(vm_arr...)
      avg_v_arr[summary_idx]=mean(vm_arr)
      ngen_Qll_arr[summary_idx]=sum(q_arr .≈ qmin_arr)
      ngen_Qul_arr[summary_idx]=sum(q_arr .≈ qmax_arr)
      ngen_arr[summary_idx]=length(ref_subnet[:gen])
      minQloss_arr[summary_idx]=min(q_loss_arr...)
      ngen_arr[summary_idx]=length(ref_subnet[:gen])
      minQloss_arr[summary_idx]=min(q_loss_arr...)
      maxQloss_arr[summary_idx]=max(q_loss_arr...)
      n_negQloss_arr[summary_idx]=sum(q_loss_arr.<0.0)
      n_s_bound_arr[summary_idx]=sum((s_to_arr .≈ rate_a_arr) .| (s_from_arr .≈ rate_a_arr))
      max_ang_diff_arr[summary_idx]=max(angle_diff_arr...)
      min_ang_diff_arr[summary_idx]=min(angle_diff_arr...)
      nbranch_arr[summary_idx]=length(ref_subnet[:branch])
      cost_arr[summary_idx]=(value(cost_pg[subnet_idx]) + value(cost_dcline[subnet_idx]))
      obj_arr[summary_idx]=objval
      time_s_arr[summary_idx]=t_build_solve
      status_arr[summary_idx]=termination_status(model)
      summary_idx += 1
   end
   res_summary = DataFrame(
      subnet=subnets_arr,
      f=f_arr, fmin=fmin_arr, fmax=fmax_arr,
      min_v=min_v_arr, max_v=max_v_arr, avg_v=avg_v_arr,
      ngen_Qll=ngen_Qll_arr, ngen_Qul=ngen_Qul_arr,
      ngen=ngen_arr, minQloss=minQloss_arr,
      maxQloss=maxQloss_arr, n_negQloss=n_negQloss_arr,
      n_s_bound=n_s_bound_arr,
      max_ang_diff=max_ang_diff_arr, min_ang_diff=min_ang_diff_arr,
      nbranch=nbranch_arr,
      total_loss=subnet_total_loss,
      converter_loss=converter_loss,
      cost=cost_arr,
      objective=obj_arr,
      time_s=time_s_arr,
      status=status_arr
   )
   if output_to_files
      CSV.write("$output_folder/res_summary.csv", res_summary)
   end

   tot_load = 0.0
   tot_gen = 0.0
   tot_shunt_p = 0.0
   tot_branch_loss = 0.0
   tot_interface_p = 0.0
   tot_dc_loss = 0.0
   tot_bus_flow = 0.0
   for (subnet_idx, ref_subnet) in ref
      sn_load = 0.0
      sn_gen = 0.0
      sn_shunt_p = 0.0
      sn_branch_loss = 0.0
      sn_interface_p = 0.0
      sn_dc_loss = 0.0
      sn_bus_flow = 0.0
      for (i,bus) in ref_subnet[:bus]
         if haskey(bus, "converter_index")
            bus_interfaces = bus["converter_index"]
         else
            bus_interfaces = Array{Int64,1}()
         end
         bus_flows = sum(Float64[value(p[subnet_idx][a]) for a in ref_subnet[:bus_arcs][i]])
         bus_loads = [ref_subnet[:load][l] for l in ref_subnet[:bus_loads][i]]
         bus_shunts = [ref_subnet[:shunt][s] for s in ref_subnet[:bus_shunts][i]]
         bus_i_gen = sum(Float64[value(pg[subnet_idx][g]) for g in ref_subnet[:bus_gens][i]])                  # sum of active power generation at bus i
         bus_i_load = sum(Float64[load["pd"] for load in bus_loads])                                    # sum of active load consumption at bus i
         bus_i_shunt_p = sum(Float64[shunt["gs"] for shunt in bus_shunts])*value(vm[subnet_idx][i])^2   # sum of active shunt element injections at bus i
         bus_i_interface_p = sum(Float64[value(p_i[conv_idx][(subnet_idx, i)]) for (conv_idx,busi,busj) in bus_interfaces])
         bus_i_dc_p = sum(Float64[value(p_dc[subnet_idx][a_dc]) for a_dc in ref_subnet[:bus_arcs_dc][i]])
         sn_load += bus_i_load
         sn_gen += bus_i_gen
         sn_shunt_p += bus_i_shunt_p
         sn_interface_p += bus_i_interface_p
         sn_dc_loss += bus_i_dc_p
         sn_bus_flow += bus_flows
         bus_bal = -bus_flows+bus_i_gen+bus_i_interface_p-bus_i_shunt_p-bus_i_load-bus_i_dc_p
         if bus_bal > 1e-6
            println("Subnet $subnet_idx bus $i power balance: $(bus_bal)")
         end
      end
      subnet_br_loss = sum(Float64[abs(value(p[subnet_idx][(i,br["f_bus"],br["t_bus"])]) + value(p[subnet_idx][(i,br["t_bus"],br["f_bus"])])) for (i,br) in ref_subnet[:branch]])
      sn_branch_loss += subnet_br_loss
      if output_to_files
         open("$output_folder/balance_$(subnet_idx).csv", "w") do io
            write(io, "total load, $(sn_load)\ntotal generation, $(sn_gen)\ntotal shunt losses, $(sn_shunt_p)\ntotal branch losses, $(sn_branch_loss)\ntotal interface losses, $(sn_interface_p)\ntotal dc line losses, $(sn_dc_loss)\ntotal flow from buses, $(sn_bus_flow)")
         end
      end
      tot_load += sn_load
      tot_gen += sn_gen
      tot_shunt_p += sn_shunt_p
      tot_interface_p += sn_interface_p
      tot_dc_loss += sn_dc_loss
      tot_bus_flow += sn_bus_flow
   end
   if output_to_files
      open("$output_folder/balance_tot.csv", "w") do io
         write(io, "total load, $(tot_load)\ntotal generation, $(tot_gen)\ntotal shunt losses, $(tot_shunt_p)\ntotal branch losses, $(tot_branch_loss)\ntotal interface losses, $(tot_interface_p)\ntotal dc line losses, $(tot_dc_loss)\ntotal flow from buses, $(tot_bus_flow)")
      end
   end

   # Define arrays for interface table
   n_iface_bus = sum(Int64[length(interface["converter_buses"]) for (i,interface) in mn_data["converter"]])
   iface_arr = Array{Int64}(undef, n_iface_bus)
   subnet_arr = Array{Int64}(undef, n_iface_bus)
   bus_arr = Array{Int64}(undef, n_iface_bus)
   p_i_arr = Array{Float64}(undef, n_iface_bus)
   q_i_arr = Array{Float64}(undef, n_iface_bus)
   s_i_arr = Array{Float64}(undef, n_iface_bus)
   # Iterate over interfaces
   i = 1
   for (iface_idx,interface) in mn_data["converter"]
      for conv_bus in interface["converter_buses"]
         iface_arr[i] = iface_idx[1]
         subnet_arr[i] = conv_bus[1]
         bus_arr[i] = conv_bus[2]
         p_i_arr[i] = value(p_i[iface_idx][conv_bus])
         q_i_arr[i] = value(q_i[iface_idx][conv_bus])
         s_i_arr[i] = sqrt(p_i_arr[i]^2 + q_i_arr[i]^2)
         i += 1
      end
   end
   iface_df = DataFrame(
      interface = iface_arr, subnetwork = subnet_arr, bus = bus_arr,
      p_int = p_i_arr, q_int = q_i_arr, s_int = s_i_arr
   )
   if print_results
      print("\ninterface transfers\n$(iface_df)\n")

      # println("total system generation: $(full_gen)")
      # println("total system losses: $(full_loss)")
   end
   if output_to_files
      CSV.write("$output_folder/res_interfaces.csv", iface_df)
   end
   # Apply the PV bus changes and run the power flow
   # for (subnet_idx, ref_subnet) in ref
   #    for (i,bus) in ref_subnet[:bus]
   #       # println("bus $i:")
   #       # println("base value: $(mn_data["sn"]["$subnet_idx"]["bus"]["$(bus["index"])"]["vm"])")
   #       mn_data["sn"]["$subnet_idx"]["bus"]["$(bus["index"])"]["vm"] = value(vm[subnet_idx][i])
   #       # println("changed to: $(value(vm[subnet_idx][i]))")
   #    end
   #    for (i,gen) in ref_subnet[:gen]
   #       # println("gen $i:")
   #       # println("base value: $(mn_data["sn"]["$subnet_idx"]["gen"]["$(gen["index"])"]["pg"])")
   #       mn_data["sn"]["$subnet_idx"]["gen"]["$(gen["index"])"]["pg"] = value(pg[subnet_idx][i])
   #       # println("changed to: $(value(pg[subnet_idx][i]))")
   #    end
   # end
   # pm_1 = build_model(mn_data["sn"]["1"], ACPPowerModel, PowerModels.post_pf, setting = Dict("output" => Dict("branch_flows" => true)))
   # result = optimize_model!(pm_1, with_optimizer(Ipopt.Optimizer))
   # # result=run_pf(ref, ACPPowerModel, with_optimizer(Ipopt.Optimizer))
   # print_summary(result["solution"])
   # key = collect(keys(res_summary))[1]

   # Define output dictionary
   output_dict = Dict{String, Any}(
      "cost"=>cost_arr,
      "objective"=>obj_arr,
      "frequency (Hz)"=>f_arr,
      "subnet"=>[subnet_idx for (subnet_idx, ref_subnet) in ref],
      "status"=>termination_status(model),
      "total time (s)"=>t_build_solve,
      "CPU time (s)"=>cpu_time,
      "iterations"=>iter,
      "generation"=>subnet_generation,
      "total loss"=>subnet_total_loss,
      "line loss"=>subnet_line_loss,
      "converter loss"=>converter_loss
      )

   if obj=="areagen"
      areagen = [value(area_pg[subnet_idx]) for (subnet_idx, ref_subnet) in ref]
      n_areas = length(gen_areas)
      if n_areas > 1
         gen_areas_sorted = sort(gen_areas)
         gen_label = "generation in areas"
         for i in 1:(n_areas-2)
            gen_label *= " $(gen_areas_sorted[i]),"
         end
         gen_label *=  " $(gen_areas_sorted[end-1]) and $(gen_areas_sorted[end]) (p.u.)"
      else
         gen_label = "generation in areas $(gen_areas[1]) (p.u.)"
      end
      output_dict[gen_label] = areagen
   elseif length(gen_areas) > 0
      n_areas = length(gen_areas)
      # areagens = keys(filter(p->(ref_subnet[:bus][last(p)["gen_bus"]]["area"] in gen_areas), ref_subnet[:gen]))
      # areagen = []
      # for (subnet_idx, ref_subnet) in ref
      #    println("subnet: $subnet_idx")
      #    areagens = keys(filter(p->(ref_subnet[:bus][last(p)["gen_bus"]]["area"] in gen_areas), ref_subnet[:gen]))
      #    println("areagens: $areagens")
      #    println("sum: $(Float64[value(pg[subnet_idx][i]) for i in areagens])")
      #    println("sum: $(sum(Float64[value(pg[subnet_idx][i]) for i in areagens]))")
      #    append!(areagen,sum(Float64[value(pg[subnet_idx][i]) for i in areagens]))
      # end
      areagen = [sum(Float64[value(pg[subnet_idx][i]) for i in keys(filter(p->(ref_subnet[:bus][last(p)["gen_bus"]]["area"] in gen_areas), ref_subnet[:gen]))]) for (subnet_idx, ref_subnet) in ref]
      if n_areas > 1
         gen_areas_sorted = sort(gen_areas)
         gen_label = "generation in areas"
         for i in 1:(n_areas-2)
            gen_label *= " $(gen_areas_sorted[i]),"
         end
         gen_label *=  " $(gen_areas_sorted[end-1]) and $(gen_areas_sorted[end]) (p.u.)"
      else
         gen_label = "generation in areas $(gen_areas[1]) (p.u.)"
      end
      output_dict[gen_label] = areagen
   end

   if obj=="zonegen"
      zonegen = [value(zone_pg[subnet_idx]) for (subnet_idx, ref_subnet) in ref]
      n_zones = length(gen_zones)
      if n_zones > 1
         gen_zones_sorted = sort(gen_zones)
         gen_label = "generation in zones"
         for i in 1:(n_zones-2)
            gen_label *= " $(gen_zones_sorted[i]),"
         end
         gen_label *=  " $(gen_zones_sorted[end-1]) and $(gen_zones_sorted[end]) (p.u.)"
      else
         gen_label = "generation in zones $(gen_zones[1]) (p.u.)"
      end
      output_dict[gen_label] = zonegen
   end

   if length(area_interface) == 2
      boundary_p = zeros(Float64,length(ref))
      boundary_q = zeros(Float64,length(ref))
      boundary_s = Array{Float64}(undef, length(ref))
      boundary_s_rev = Array{Float64}(undef, length(ref))
      boundary_slack = Array{Float64}(undef, length(ref))
      boundary_p_rev = zeros(Float64,length(ref))
      boundary_q_rev = zeros(Float64,length(ref))
      boundary_capacity = zeros(Float64,length(ref))
      boundary_loss = Array{Float64}(undef, length(ref))
      for (i,(subnet_idx, ref_subnet)) in enumerate(ref)
         dc_subnet = ((!ref_subnet[:variable_f]) && (ref_subnet[:f_base] == 0))
         if (:f_max in keys(ref_subnet))
            # Also DC if max frequency is 0
            dc_subnet = dc_subnet || (ref_subnet[:f_max] == 0)
         end
         boundary_arcs = Array{Tuple{Int64,Int64,Int64},1}()
         boundary_arcs_reverse = Array{Tuple{Int64,Int64,Int64},1}()
         f_area_buses = filter(p->(last(p)["area"]==area_interface[1]),ref_subnet[:bus])
         # println("f_area_buses: $f_area_buses")
         # bus_arcs = getindex.(Ref(ref_subnet[:bus_arcs]), f_area_buses)
         for (idx, bus) in f_area_buses
            bus_arcs = ref_subnet[:bus_arcs][idx]
            for arc in bus_arcs
               if ref_subnet[:bus][arc[3]]["area"] == area_interface[2]
                  push!(boundary_arcs, arc)
                  push!(boundary_arcs_reverse, (arc[1],arc[3],arc[2]))
                  boundary_p[i] += value(p[subnet_idx][arc])
                  boundary_q[i] += dc_subnet ? 0 : value(q[subnet_idx][arc])
                  boundary_p_rev[i] += value(p[subnet_idx][(arc[1],arc[3],arc[2])])
                  boundary_q_rev[i] += dc_subnet ? 0 : value(q[subnet_idx][(arc[1],arc[3],arc[2])])
                  boundary_capacity[i] += ref[subnet_idx][:branch][arc[1]]["rate_a"] * ref[subnet_idx][:branch][arc[1]]["br_status"]
               end
            end
         end
         # println("boundary_arcs: $boundary_arcs")
         # println("boundary_p: $boundary_p")
         boundary_s[i] = sqrt(boundary_p[i]^2 + boundary_q[i]^2)
         boundary_s_rev[i] = sqrt(boundary_p_rev[i]^2 + boundary_q_rev[i]^2)
         boundary_slack[i] = boundary_capacity[i] - max(abs(boundary_s[i]), abs(boundary_s_rev[i]))
         boundary_loss[i] = abs(boundary_p[i] + boundary_p_rev[i])
      end
      output_dict["P flow from area $(area_interface[1]) to $(area_interface[2]) (p.u.)"] = boundary_p
      output_dict["P flow from area $(area_interface[2]) to $(area_interface[1]) (p.u.)"] = boundary_p_rev
      output_dict["Q flow from area $(area_interface[1]) to $(area_interface[2]) (p.u.)"] = boundary_q
      output_dict["Q flow from area $(area_interface[2]) to $(area_interface[1]) (p.u.)"] = boundary_q_rev
      output_dict["S flow from area $(area_interface[1]) to $(area_interface[2]) (p.u.)"] = boundary_s
      output_dict["S flow from area $(area_interface[2]) to $(area_interface[1]) (p.u.)"] = boundary_s_rev
      output_dict["S slack from area $(area_interface[1]) to $(area_interface[2]) (p.u.)"] = boundary_slack
      output_dict["S capacity from area $(area_interface[1]) to $(area_interface[2]) (p.u.)"] = boundary_capacity
      output_dict["loss between area $(area_interface[1]) and $(area_interface[2]) (p.u.)"] = boundary_loss
   elseif length(area_interface) > 2
      println("Warning: area_interface has $(length(area_interface)) elements. Results are saved for power transfer between areas only when exactly two areas are given in area_interface.")
   end

   if length(zone_interface) == 2
      unspecified_to_zones = false
      if length(zone_interface[2]) == 0
         unspecified_to_zones = true
      end
      boundary_p = zeros(Float64,length(ref))
      boundary_q = zeros(Float64,length(ref))
      boundary_s = Array{Float64}(undef, length(ref))
      boundary_s_rev = Array{Float64}(undef, length(ref))
      boundary_slack = Array{Float64}(undef, length(ref))
      boundary_p_rev = zeros(Float64,length(ref))
      boundary_q_rev = zeros(Float64,length(ref))
      boundary_capacity = zeros(Float64,length(ref))
      boundary_loss = Array{Float64}(undef, length(ref))
      for (i,(subnet_idx, ref_subnet)) in enumerate(ref)
         boundary_arcs = Array{Tuple{Int64,Int64,Int64},1}()
         boundary_arcs_reverse = Array{Tuple{Int64,Int64,Int64},1}()
         f_area_buses = filter(p->(last(p)["zone"] in zone_interface[1]),ref_subnet[:bus])
         # println("f_area_buses: $f_area_buses")
         # bus_arcs = getindex.(Ref(ref_subnet[:bus_arcs]), f_area_buses)
         for (idx, bus) in f_area_buses
            bus_arcs = ref_subnet[:bus_arcs][idx]
            for arc in bus_arcs
               if (ref_subnet[:bus][arc[3]]["zone"] in zone_interface[2]) || (unspecified_to_zones && !(ref_subnet[:bus][arc[3]]["zone"] in zone_interface[2]))
                  push!(boundary_arcs, arc)
                  push!(boundary_arcs_reverse, (arc[1],arc[3],arc[2]))
                  boundary_p[i] += value(p[subnet_idx][arc])
                  boundary_q[i] += dc_subnet ? 0 : value(q[subnet_idx][arc])
                  boundary_p_rev[i] += value(p[subnet_idx][(arc[1],arc[3],arc[2])])
                  boundary_q_rev[i] += dc_subnet ? 0 : value(q[subnet_idx][(arc[1],arc[3],arc[2])])
                  br_rate_a = Inf
                  if "rate_a" in keys(ref[subnet_idx][:branch][arc[1]])
                     br_rate_a = ref[subnet_idx][:branch][arc[1]]["rate_a"]
                  end
                  boundary_capacity[i] += br_rate_a * ref[subnet_idx][:branch][arc[1]]["br_status"]
               end
            end
         end
         # println("boundary_arcs: $boundary_arcs")
         # println("boundary_p: $boundary_p")
         boundary_s[i] = sqrt(boundary_p[i]^2 + boundary_q[i]^2)
         boundary_s_rev[i] = sqrt(boundary_p_rev[i]^2 + boundary_q_rev[i]^2)
         boundary_slack[i] = boundary_capacity[i] - max(abs(boundary_s[i]), abs(boundary_s_rev[i]))
         boundary_loss[i] = abs(boundary_p[i] + boundary_p_rev[i])
      end
      output_dict["P flow from zones $(zone_interface[1]) to $(zone_interface[2]) (p.u.)"] = boundary_p
      output_dict["P flow from zones $(zone_interface[2]) to $(zone_interface[1]) (p.u.)"] = boundary_p_rev
      output_dict["Q flow from zones $(zone_interface[1]) to $(zone_interface[2]) (p.u.)"] = boundary_q
      output_dict["Q flow from zones $(zone_interface[2]) to $(zone_interface[1]) (p.u.)"] = boundary_q_rev
      output_dict["S flow from zones $(zone_interface[1]) to $(zone_interface[2]) (p.u.)"] = boundary_s
      output_dict["S flow from zones $(zone_interface[2]) to $(zone_interface[1]) (p.u.)"] = boundary_s_rev
      output_dict["S slack from zones $(zone_interface[1]) to $(zone_interface[2]) (p.u.)"] = boundary_slack
      output_dict["S capacity from zones $(zone_interface[1]) to $(zone_interface[2]) (p.u.)"] = boundary_capacity
      output_dict["loss between zones $(zone_interface[1]) and $(zone_interface[2]) (p.u.)"] = boundary_loss
   elseif length(zone_interface) > 2
      println("Warning: zone_interface has $(length(zone_interface)) elements. Results are saved for power transfer between zones only when exactly two arrays of zones are given in zone_interface.")
   end
   if uniform_gen_scaling
      output_dict["alpha_upstream"] = value(alpha_upstream)
      output_dict["alpha_downstream"] = value(alpha_downstream)
   end

   outdict_string = JSON.json(output_dict)
   if output_to_files
      open("$output_folder/output_values.json", "w") do f
          write(f, outdict_string)
      end
   end

   binding_cnstr_dict = Dict{Symbol, Any}()
   # Get duals for all constraints in the constraints reference dict
   if termination_status(model) in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, OPTIMAL, ALMOST_OPTIMAL]
      if has_duals(model)
         for (class_i,class) in constraints
            # println("$class_i constraints")
            for (subnet_i,subnet) in class
               # println("subnet: $subnet")
               # println("typeof(subnet): $(typeof(subnet))")
               if !isa(subnet, Dict)
                  subnet_dict = Dict(1=>subnet)
               else
                  subnet_dict = subnet
               end
               # println("subnet_dict: $subnet_dict")
               for (cnstr_i,cnstr) in subnet_dict
                  if abs(dual(cnstr)) >= 1e-6
                     if !(class_i in keys(binding_cnstr_dict))
                        binding_cnstr_dict[class_i] = Dict{String, Any}()
                     end
                     # if "ctg_label" in keys(mn_data["sn"]["$ctg_i"])
                     #    lbl = mn_data["sn"]["$ctg_i"]["ctg_label"]
                     # else
                     #    lbl = "no"
                     # end
                     if !("subnet $subnet_i" in keys(binding_cnstr_dict[class_i]))
                        binding_cnstr_dict[class_i]["subnet $subnet_i"] = Array{String,1}()
                     end
                     # println("dual: $(dual(cnstr))")
                     push!(binding_cnstr_dict[class_i]["subnet $subnet_i"], "$cnstr")
                  end
               end
            end
         end
      else
         println("model has no duals")
      end
   end
   binding_cnstr_dict_string = JSON.json(binding_cnstr_dict)
   if output_to_files
      open("$output_folder/binding_constraints.json", "w") do f
         write(f, binding_cnstr_dict_string)
      end
   end

   return (output_dict, res_summary, solution_pm, binding_cnstr_dict)
end

# function mywrite_nl_file(f::IO, m::JuMP.Model)
#     # JuMP.build(mj)
#     #
#     # m = internalmodel(mj)
#     # m.status = :NotSolved
#     # m.solve_exitcode = -1
#     # m.solve_result_num = -1
#     # m.solve_result = "?"
#     # m.solve_message = ""
#
#     # There is no non-linear binary type, only non-linear discrete, so make
#     # sure binary vars have bounds in [0, 1]
#     # for i in 1:m.nvar
#     #     if m.vartypes[i] == :Bin
#     #         if m.x_l[i] < 0
#     #             m.x_l[i] = 0
#     #         end
#     #         if m.x_u[i] > 1
#     #             m.x_u[i] = 1
#     #         end
#     #     end
#     # end
#     AmplNLWriter.make_var_index!(m)
#     AmplNLWriter.make_con_index!(m)
#
#     AmplNLWriter.write_nl_file(f, m.inner)
# end


text:
# need to load LineProperties first
include("structures/LineProperties.jl")
include("structures/SolarData.jl")
include("structures/SpecParams.jl")
include("structures/SynthWorkspace.jl")
include("structures/DiskParams.jl")

ext: jl

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.0
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# # PyPlot.jl

# ```
# This module provides a Julia interface to the Matplotlib plotting library from Python, and specifically to the matplotlib.pyplot module.
#
# PyPlot uses the Julia PyCall package to call Matplotlib directly from Julia with little or no overhead (arrays are passed without making a copy).
# ```

using Images
using PyPlot
using TestImages

# > Only the currently documented matplotlib.pyplot API is exported. To use other functions in the module, you can also call matplotlib.pyplot.foo(...) as plt.foo(...). For example, plt.plot(x, y) also works.

c = testimage("cameraman") |> channelview |> rawview
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5)) # you can call without `plt` e.g `subplots(1, 2)`
ax1.set_title("cameraman")
ax1.imshow(c, cmap=:gray)
ax2.set_title("hist")
ax2.hist(vec(c), density=true, bins=255);

m = testimage("lena") |> channelview |> rawview
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5))
ax1.set_title("lena")
ax1.imshow(permutedims(m, (2,3,1)), cmap=:gray)
ax2.set_title("hist")
ax2.hist(vec.([m[1,:,:], m[2,:,:], m[3,:,:]]), density=true, bins=25, color=[:red, :green, :blue]);


text:
function barzanLinInvoke(model, features)

	predictions = [features ones(size(features,1),1)] * model;
	return predictions

end

max_stars_repo_name: barzan/dbseer

max_issues_repo_name: barzan/dbseer

max_forks_repo_name: barzan/dbseer


text:
apply_laser_over_atoms(laser, atoms) = _apply_laser(laser, atoms)
apply_laser_over_sensors(laser, sensors) = _apply_laser(laser, sensors)

apply_laser_over_oneAtom(laser, oneAtom) = _apply_laser_over_oneCoordinate(laser, oneAtom)
apply_laser_over_oneSensor(laser, oneSensor) =
    _apply_laser_over_oneCoordinate(laser, oneSensor)

function _apply_laser(laser, spatial_coordinates)
    N = _get_number_elements(spatial_coordinates)
    _core_LaserFunction = _get_core_LaserFunction(laser)
    _spatial_matrix = _get_spatial_matrix(spatial_coordinates)
    E₀ = estimate_E₀(laser)

    E = zeros(ComplexF64, N)
    matrix_view = view(_spatial_matrix, :, :)

    Threads.@threads for n = 1:N
        oneCoordinate = matrix_view[:, n]
        @inbounds E[n] = _core_LaserFunction(oneCoordinate, E₀, laser)
    end

    return E
end
function _apply_laser_over_oneCoordinate(laser, oneCoordinate)
    _core_LaserFunction = _get_core_LaserFunction(laser)
    E₀ = estimate_E₀(laser)
    E = _core_LaserFunction(oneCoordinate, E₀, laser)

    return E
end

_get_number_elements(r::AbstractMatrix) = size(r, 2)
_get_number_elements(atoms::Atom{T}) where {T} = atoms.N

_get_spatial_matrix(r::AbstractMatrix) = r
_get_spatial_matrix(atoms::Atom{T}) where {T} = atoms.r

_get_core_LaserFunction(laser::Laser{PlaneWave3D}) = _core_PlaneWave3D
_get_core_LaserFunction(laser::Laser{Gaussian3D}) = _core_Gaussian3D

function _core_PlaneWave3D(oneAtom, E₀, laser)
    direction = view(laser.pump.direction, :)
    Eᵢ =
        E₀ * cis(
            +k₀ * (
                oneAtom[1] * direction[1] +
                oneAtom[2] * direction[2] +
                oneAtom[3] * direction[3]
            ),
        )
    return Eᵢ
end

function _core_Gaussian3D(oneAtom, E₀, laser)
    w₀ = laser.pump.w₀
    x, y, z = oneAtom[1], oneAtom[2], oneAtom[3]

    ## This formula is stable for z==0
    # Ref: Eq 3.11 from "CHAPTER 3. PROPAGATION AND FOCUSING OF OPTICAL FIELDS"
    denominator_factor = 1 + 2im * z / (k₀ * w₀^2)
    Eᵢ = E₀ * cis(+k₀ * z)
    Eᵢ = Eᵢ * exp(-(x^2 + y^2) / (denominator_factor * w₀^2))
    Eᵢ = Eᵢ / denominator_factor

    return Eᵢ
end


text:
# helper functions to check that the given system can be applied to the states/inputs
@inline _is_conformable_state(system, x) = statedim(system) == length(x)
@inline _is_conformable_input(system, u) = inputdim(system) == length(u)
@inline _is_conformable_noise(system, w) = noisedim(system) == length(w)
@inline _in_stateset(system, x) = x ∈ stateset(system)
@inline _in_inputset(system, u) = u ∈ inputset(system)
@inline _in_noiseset(system, w) = w ∈ noiseset(system)

@inline function _argument_error(sym, set=:none)
    if sym == :x
        vector = "state"
    elseif sym == :u
        vector = "input"
    elseif sym == :w
        vector = "noise"
    end
    if set == :none
        throw(ArgumentError("the $vector vector has the wrong dimensions"))
    else
        throw(ArgumentError("the $vector vector is not contained in the $(vector) set"))
    end
end

# ====================
# Generic _instantiate method
# ====================

"""
    _instantiate(system::AbstractSystem, x::AbstractVector;
                [check_constraints]=true)

Return the result of instantiating an `AbstractSystem` at the current state.

### Input

- `system`            -- `AbstractSystem`
- `x`                 -- state (it should be any vector type)
- `check_constraints` -- (optional, default: `true`) check if the state belongs to
                         the state set

### Output

The result of applying the system to state `x`.

### Notes

The `_instantiate` method generalizes the `successor` of an `AbstractDiscreteSystem`
and the `vector_field` of an `AbstractContinuousSystem` into a single method.
"""
function _instantiate(system::AbstractSystem, x::AbstractVector;
                     check_constraints::Bool=true)
    !_is_conformable_state(system, x) && _argument_error(:x)
    if isconstrained(system) && check_constraints
        !_in_stateset(system, x) && _argument_error(:x,:X)
    end

    if islinear(system)
        return state_matrix(system) * x

    elseif isaffine(system)
        return state_matrix(system) * x + affine_term(system)

    elseif ispolynomial(system) || isblackbox(system)
        return mapping(system)(x)

    else
        throw(ArgumentError("_instantiate not defined for type `$(typename(sys))`"))
    end
end

"""
    _instantiate(system::AbstractSystem, x::AbstractVector, u::AbstractVector;
                [check_constraints]=true)

Return the result of instantiating an `AbstractSystem` at the current state and
applying one input.

### Input

- `system`            -- `AbstractSystem`
- `x`                 -- state (it should be any vector type)
- `u`                 -- input (it should be any vector type) or noise, if `system`
                         is not controlled
- `check_constraints` -- (optional, default: `true`) check if the state belongs to
                         the state set

### Output

The result of applying the system to state `x` and input `u`.

### Notes

If the system is not controlled but noisy, the input `u` is interpreted as noise.

The `_instantiate` method generalizes the `successor` of an `AbstractDiscreteSystem`
and the `vector_field` of an `AbstractContinuousSystem` into a single method.
"""
function _instantiate(system::AbstractSystem, x::AbstractVector, u::AbstractVector;
                     check_constraints::Bool=true)
    # Figure out if `system` is a controlled or noisy system and set the function
    # `matrix` to either `input_matix` or `noise_matrix`
    if iscontrolled(system) && !isnoisy(system)
        input_var = :u; input_set = :U; matrix = input_matrix
        _is_conformable = _is_conformable_input
        _in_set = _in_inputset
    elseif isnoisy(system) && !iscontrolled(system)
        input_var = :w; input_set = :W; matrix = noise_matrix
        _is_conformable = _is_conformable_noise
        _in_set = _in_noiseset
    else
        throw(ArgumentError("successor function for $(typeof(system)) does not have 2 arguments"))
    end

    !_is_conformable_state(system, x) && _argument_error(:x)
    !_is_conformable(system, u) && _argument_error(input_var)

    if isconstrained(system) && check_constraints
        !_in_stateset(system, x) && _argument_error(:x,:X)
        !_in_set(system, u) && _argument_error(input_var, input_set)
    end

    if islinear(system)
        return state_matrix(system) * x + matrix(system) * u

    elseif isaffine(system)
        return state_matrix(system) * x + affine_term(system) + matrix(system) * u

    elseif ispolynomial(system) || isblackbox(system)
        return mapping(system)(x, u)

    else
        throw(ArgumentError("_instantiate not defined for type `$(typename(system))`"))
    end
end

"""
    _instantiate(system::AbstractSystem,
                x::AbstractVector, u::AbstractVector, w::AbstractVector; [check_constraints]=true)

Return the result of instantiating an `AbstractSystem` at the current state and
applying two inputs to an `AbstractSystem`.

### Input

- `system`            -- `AbstractSystem`
- `x`                 -- state (it should be any vector type)
- `u`                 -- input (it should be any vector type)
- `w`                 -- noise (it should be any vector type)
- `check_constraints` -- (optional, default: `true`) check if the state belongs to
                         the state set

### Output

The result of applying the system to state `x`, input `u` and noise `w`.

### Notes

The `_instantiate` method generalizes the `successor` of an `AbstractDiscreteSystem`
and the `vector_field` of an `AbstractContinuousSystem` into a single method.
"""
function _instantiate(system::AbstractSystem, x::AbstractVector, u::AbstractVector, w::AbstractVector;
                     check_constraints::Bool=true)
    !_is_conformable_state(system, x) && _argument_error(:x)
    !_is_conformable_input(system, u) && _argument_error(:u)
    !_is_conformable_noise(system, w) && _argument_error(:w)

    if isconstrained(system) && check_constraints
        !_in_stateset(system, x) && _argument_error(:x,:X)
        !_in_inputset(system, u) && _argument_error(:u,:U)
        !_in_noiseset(system, w) && _argument_error(:w,:W)
    end

    if islinear(system)
        return state_matrix(system) * x + input_matrix(system) * u + noise_matrix(system) * w

    elseif isaffine(system)
        return state_matrix(system) * x + affine_term(system) + input_matrix(system) * u + noise_matrix(system) * w

    elseif ispolynomial(system) || isblackbox(system)
        return mapping(system)(x, u, w)

    else
        throw(ArgumentError("_instantiate not defined for type `$(typename(system))`"))

    end
end


text:
using Documenter
using WorldMapData

# Generate documents
makedocs(
    modules   = [WorldMapData],
    doctest   = false,
    clean     = true,
    linkcheck = true,
    format    = Documenter.HTML(),
    sitename  = "WorldMapData.jl",
    authors   = "Duncan Eddy",
    pages     = Any[
        "Home" => "index.md",
    ]
)

# Generate plots
# makeplots()

deploydocs(
    repo = "github.com/sisl/WorldMapData.jl",
    devbranch = "master",
    devurl = "latest",
)


text:
using Documenter, KnetLayers

makedocs(

    modules = [KnetLayers],
    clean = false,              # do we clean build dir
    format = :html,
    sitename = "KnetLayers.jl",
    authors = "Ekin Akyurek and contributors.",
    doctest = true,
    pages = Any[ # Compat: `Any` for 0.4 compat
        "Home" => "index.md",
        "Function Documentation" => Any[
            "reference.md",
        ],
    ],
#    analytics = "UA-89508993-1",
#    linkcheck = !("skiplinks" in ARGS),
)

deploydocs(
    repo = "github.com/ekinakyurek/KnetLayers.jl.git",
    julia = "1.0",
    osname = "linux",
    target = "build",
    make = nothing,
    deps = nothing,
    #deps   = Deps.pip("mkdocs", "python-markdown-math"),
)


text:
function POMDPs.reward(mdp::PedCarMDP, s::PedCarMDPState, a::PedCarMDPAction, sp::PedCarMDPState)
    r = mdp.action_cost
    ego = sp.ego
    if sp.crash
        r += mdp.collision_cost
    end
    if ego.posF.s >= get_end(mdp.env.roadway[mdp.ego_goal]) &&
       get_lane(mdp.env.roadway, ego).tag == mdp.ego_goal
        r += mdp.goal_reward
    end
    return r
end


text:
module Data

using DrWatson
using Random, Statistics
using DynamicalSystems, OrdinaryDiffEq, LinearAlgebra
using DelimitedFiles, DrWatson
using StochasticDiffEq

# default constants used if no other is given
const default_N = 10000
const diffeq = (alg = Vern9(), reltol = 1e-12, abstol = 1e-12, maxiters = typemax(Int))
const diffeq_lowacc = (alg = Tsit5(), reltol = 1e-9, abstol = 1e-9, maxiters = typemax(Int))

function kaplanyorke_map(; N = default_N, λ = 0.2, kwargs...)
    f(u, λ, t) = SVector((2u[1]) % 1.0 + 1e-15rand(Random.GLOBAL_RNG), λ*u[2] + cospi(4u[1]))
    ds = DiscreteDynamicalSystem(f, SVector(0.15, 0.2), λ)
    tr = trajectory(ds, N; Ttr = 1000)
    standardize(tr)
end

function roessler_periodic(; N = default_N, Δt = 0.5, kwargs...)
    ds = Systems.roessler([0.1, -0.2, 0.1]; c = 3.0, a = 0.2, b = 0.2)
    tr = trajectory(ds, N*Δt; Δt, Ttr = 1000, diffeq...)
    standardize(tr)
end

function roessler_chaotic(; η=0, correlated=false, N = default_N, Δt = 0.5, kwargs...)
    ds = Systems.roessler([0.1, -0.2, 0.1]; c = 5.7, a = 0.2, b = 0.2)
    tr = trajectory(ds, N*Δt; Δt, Ttr = 1000, diffeq...)
    tr = standardize(tr)
    if η ≠ 0 # add static noise
        if !correlated
            seed = Random.MersenneTwister(5653435)
            for i in 1:length(tr)
                tr[i] += η*randn(seed, SVector{3})
            end
        else correlated
            a = arma_noise(N+1)
            for i in 1:length(tr)
                tr[i] = tr[i] .+ η .* a[i]
            end
        end
    end
    return tr
end

function roessler_embedding(; d = 3, kwargs...)
    x = roessler_chaotic(kwargs...)[:, 1]
    x = standardize(x)
    τ = estimate_delay(x, "mi_min")
    return embed(x, d, τ)
end

function roessler_rounded(; digits = 2, kwargs...)
    X = roessler_chaotic(; kwargs...)
    Z = [round.(x; digits = digits) for x in X]
    return Dataset(Z)
end

function roessler_dynoise(du,u, p, t)
    a, b, c, η = p
    du[1] = -u[2]-u[3]
    du[2] = u[1] + a*u[2]
    du[3] = b + u[3]*(u[1] - c)
end

# SDEProblem wants a specific function for the noise
# (This assumes diagonal noise, with the only noise acting on the
# 2nd variable)
function σ_roessler_dynoise(du, u, p, t)
    a, b, c, η = p
    du[1] = 0
    du[2] = η
    du[3] = 0
end

function roessler_dynamic(; N = default_N, Δt = 0.2, η=0.1, kwargs...)
    # using the SDE solve from DifferentialEquations here
    tspan = (0.0, N*Δt + 100)
    prob_sde_roessler = SDEProblem(
        roessler_dynoise, σ_roessler_dynoise, [0.1, -0.2, 0.1],
        tspan, [0.2, 0.2, 5.7, η]
    )
    saveat = 100.0:Δt:tspan[end]
    # I really don't have much experience with which solver to use
    sol = solve(prob_sde_roessler, SRA(); saveat, maxiters = typemax(Int))
    X = Dataset(hcat(sol.u...)')
    return standardize(X)
end

function roessler_lorenz(; N = default_N, kwargs...)
    tr1 = roessler_periodic(N = N÷2)
    tr2 = lorenz96_chaotic(; D = 4, N = N÷2)[:, 1:3]
    tr2 = trajectory(ds, (N/2)*Δt; Δt, Ttr = 1000, diffeq...)
    append!(standardize(tr1), standardize(tr2))
end

function roessler_sphere(; N = default_N, Δt = 0.5, kwargs...)
    ds = Systems.roessler([0.1, -0.2, 0.1]; c = 3.0, a = 0.2, b = 0.2)
    tr1 = standardize(trajectory(ds, (N/2)*Δt; Δt, Ttr = 1000, diffeq...))
    tr2 = standardize(uniform_sphere(; N = N÷2))
    return append!(tr1, tr2)
end

function roessler_nonstationary_eom(u, p, t)
    @inbounds begin
        a, b, dc = p
        du1 = -u[2]-u[3]
        du2 = u[1] + a*u[2]
        du3 = b + u[3]*(u[1] - u[4])
        du4 = dc
        return SVector{4}(du1, du2, du3, du4)
    end
end

function roessler_nonstationary(; N = default_N, Δt = 0.1, kwargs...)
    u0 = [-5.0, -5.0, 0.0, 3.0] # u0[4] is initial `c` value
    p0 = [0.2, 0.2, (20 - 3)/(N*Δt)]
    ds = ContinuousDynamicalSystem(roessler_nonstationary_eom, u0, p0)
    tr = trajectory(ds, N*Δt; Δt, Ttr = 0, diffeq...)
end

function torus2(; N = default_N, Δt = 0.1, ω = sqrt(3), R = 2.0, r = 1.0, kwargs...)
    function torus(u)
        θ, φ = u
        x = (R + r*cos(θ))*cos(φ)
        y = (R + r*cos(θ))*sin(φ)
        z = r*sin(θ)
        return SVector(x, y, z)
    end
    θs = range(0; step = 0.1, length = N)
    φs = ω .* θs
    return Dataset([torus(u) for u in zip(θs, φs)])
end

function koch(; maxk = 7, kwargs...)
    flakepoints = SVector{2}.([[0.0; 0.0], [0.5; sqrt(3)/2], [1; 0.0], [0.0; 0.0]])
    function innerkoch(points, maxk, α = sqrt(3)/2)
        Q = @SMatrix [0 -1; 1 0]
        for k = 1:maxk
            n = length(points)
            new_points = eltype(points)[]
            for i = 1:n-1
                p1, p2 = points[i], points[i+1]
                v = (p2 - p1) / 3
                q1 = p1 + v
                q2 = p1 + 1.5v + α * Q * v
                q3 = q1 + v
                push!(new_points, p1, q1, q2, q3)
            end
            push!(new_points, points[end])
            points = new_points
        end
        return points
    end
    kochpoints = innerkoch(flakepoints, maxk)
    kochdata = Dataset(unique!(kochpoints))
end

function henonheiles_quasi(; N = default_N, Δt = 0.2, kwargs...)
    u0 = [0.0, 0.1, 0.5, 0.0]
    ds = Systems.henonheiles(u0)
    tr = trajectory(ds, N*Δt; Δt, Ttr = 100, diffeq...)
    standardize(tr)
end

function henonheiles_chaotic(; N = default_N, Δt = 0.5, kwargs...)
    u0 = [0.0, -0.25, 0.42081, 0.0]
    ds = Systems.henonheiles(u0)
    tr = trajectory(ds, N*Δt; Δt, Ttr = 100, diffeq...)
    standardize(tr)
end

function uniform_sphere(; N = default_N, kwargs...)
    A = SVector{3, Float64}[]
    i = 0
    while i < N
        x = rand(SVector{3, Float64}) .- 0.5
        if norm(x) ≤ 1/2√3
            push!(A, x)
            i += 1
        end
    end
    return Dataset(A)
end

function standardmap_chaotic(; N = default_N, k = 64.0, kwargs...)
    ds = Systems.standardmap([0.08152, 0.122717]; k)
    tr = trajectory(ds, N; Ttr = 100)
end

function henon_chaotic(; N = default_N, a = 1.4, b = 0.3, kwargs...)
    ds = Systems.henon([0.08152, 0.122717]; a, b)
    tr = trajectory(ds, N; Ttr = 100)
    standardize(tr)
end

function towel_chaotic(; N = default_N, η = 0, correlated = false, kwargs...)
    ds = Systems.towel([0.085, -0.121, 0.075])
    tr = trajectory(ds, N; Ttr = 100)
    tr = standardize(tr)
    if η ≠ 0 # add static noise
        if !correlated
            seed = Random.MersenneTwister(5653435)
            for i in 1:length(tr)
                tr[i] += η*randn(seed, SVector{3})
            end
        else correlated
            a = arma_noise(N+1)
            for i in 1:length(tr)
                tr[i] = tr[i] .+ η .* a[i]
            end
        end
    end
    return tr
end

function arma_noise(N)
    η = randn(N)
    s = ones(N)
    for n in 4:N
        s[n] = 1.625s[n-1] - 0.284s[n-2] - 0.355s[n-3] + η[n] - 0.96η[n-1]
    end
    s ./= std(s)
    return s
end

function rule_towel_noisy(x, p, n)
    @inbounds x1, x2, x3 = x[1], x[2], x[3]
    η, seed = p
    SVector(
        3.8*x1*(1-x1) - 0.05*(x2+0.35)*(1-2*x3),
        0.1*( (x2+0.35)*(1-2*x3) - 1 + η*randn(seed) )*(1 - 1.9*x1),
        3.78*x3*(1-x3)+0.2*x2
    )
end
function towel_dynamic_noise(; N = default_N, η = 0.1, kwargs...)
    p = [η, MersenneTwister(919009)]
    ds = DiscreteDynamicalSystem(rule_towel_noisy, [0.085, -0.121, 0.075], p)
    tr = trajectory(ds, N; Ttr = 100)
    tr = standardize(tr)
end

function lorenz96_chaotic(; N = default_N, Δt = 0.2, D = 4, F = 24.0, kwargs...)
    ds = Systems.lorenz96(D, range(0; length = D, step = 0.1); F)
    tr = trajectory(ds, N*Δt; Δt, Ttr = 1000.0, diffeq...)
    standardize(tr)
end

function lorenz96_rounded(; digits = 2, kwargs...)
    X = lorenz96_chaotic(; kwargs...)
    Z = [round.(x; digits) for x in X]
    return Dataset(Z)
end

function lorenz96_nonstationary_f!(dx, x, p, t)
    dF = p[1]
    N = length(x)-1
    F = x[end]

    # 3 edge cases
    @inbounds dx[1] = (x[2] - x[N - 1]) * x[N] - x[1] + F
    @inbounds dx[2] = (x[3] - x[N]) * x[1] - x[2] + F
    @inbounds dx[N] = (x[1] - x[N - 2]) * x[N - 1] - x[N] + F
    # then the general case
    for n in 3:(N - 1)
      @inbounds dx[n] = (x[n + 1] - x[n - 2]) * x[n - 1] - x[n] + F
    end
    dx[end] = dF
    return nothing
end

function lorenz96_nonstationary(; N = default_N, Δt = 0.2, D = 6, F2=24.0, F1=1.0, kwargs...)
    u0 = range(0.1; length = D, step = 0.1)
    dF = (F2-F1)/(N*Δt)
    u0 = [u0..., F1]
    ds = ContinuousDynamicalSystem(lorenz96_nonstationary_f!, u0, [dF])
    tr = trajectory(ds, N*Δt; Δt, Ttr = 0,  diffeq...)
    # standardize(tr)
end

function lorenz96_nonstationary_2(; N = default_N, Δt = 0.2, D = 6, F2=40.0, F1=10.0, kwargs...)
    ds = Systems.lorenz96(D, range(0; length = D, step = 0.1); F = F1)
    tr1 = trajectory(ds, (N/2)*Δt; Δt, Ttr = 1000.0, diffeq...)

    ds = Systems.lorenz96(D, range(0; length = D, step = 0.1); F = F2)
    tr2 = trajectory(ds, (N/2)*Δt; Δt, Ttr = 1000.0, diffeq...)

    return standardize(append!(tr1, tr2))
end

logistic(u) = 4*u*(1-u)
function logistics!(un, u, p, t)
    L = length(u)
    k, r, η = p # η is dynamic noise
    for i in 2:L-1
        un[i] = logistic(u[i] + k*(u[i-1] - 2u[i] + u[i+1]) + η*randn())
    end
    # circling around
    un[1] = logistic(u[1] + k*(u[L] - 2u[1] + u[2]) + η*randn())
    un[L] = logistic(u[L] + k*(u[L-1] - 2u[L] + u[1]) + η*randn())
    return
end

function coupled_logistics(; N = default_N, D = 8, k = 0.1, r = 4.0, η = 0, kwargs...)
    u0 = range(0.1, 0.9; length = D) .+ 0.0001
    ds = DiscreteDynamicalSystem(logistics!, u0, [k, r, η])
    tr = trajectory(ds, N; Ttr = 1000)
end

function embed_system(; system, d = 4, kwargs...)
    data_producing_function = getfield(Data, system)
    x = data_producing_function(; kwargs...)[:, 1]
    x = standardize(x)
    τ = estimate_delay(x, "mi_min")
    return embed(x, d, τ)
end

using DelimitedFiles
function experimental_data(; N = nothing, name, kwargs...)
    file = datadir("experimental", name*".txt")
    if name == "electrochemical1"
        τ = (0, 26, 13, 5, 20)
        x = vec(readdlm(file))
        A = genembed(x, τ)
    elseif name == "electrochemical2"
        τ = (0, 25, 16, 148, 138, 87, 60, 105)
        x = vec(readdlm(file))
        A = genembed(x, τ)
    elseif name == "shinriki1"
        τ = (0, 19, 38, 57)
        # js = ones(length(τ)) # perhaps use multivariate here
        x = readdlm(file)[:, 1]
        A = genembed(x, τ)
    elseif name == "shinriki2"
        τ = (0, 20, 40, 60)
        # js = ones(length(τ)) # perhaps use multivariate here
        x = readdlm(file)[:, 1]
        A = genembed(x, τ)
    elseif name == "nifty50"
        τ = [j*43 for j in 0:5]
        x = vec(readdlm(file))
        A = genembed(x, τ)
    elseif name == "vostok"
        τ = [j*50 for j in 0:7]
        x = vec(readdlm(file))
        A = genembed(x, τ)
    elseif name == "roessler_embed"
        τ = [0, 6, 3, 14]      # from PECUZAL
        τ = [j*7 for j in 0:6] # from TRADITIONAL
        x = readdlm(file)[:, 1]
        A = genembed(x, τ)
    elseif name == "roessler_all"
        M = readdlm(file)
        A = Dataset(M) # just take all 28 oscillators
    elseif name == "doublependulum"
        M = readdlm(file)
        τx = [0, 51, 25, 39, 12]
        τy = [0, 60, 30, 46, 14]
        jx = fill(1, length(τx))
        jy = fill(2, length(τy))
        τd = [0, 61, 31, 47, 31, 0]
        jd = [2, 2, 2, 2, 1, 1]
        A = Dataset(M)
        genembed(A, τx, jx)
    else
        error("Experimental dataset $(name) is uknown")
    end
    A = standardize(A)
    if isnothing(N)
        return A
    else
        return Dataset(A[1:min(length(A), N)])
    end
end

function electronic_roessler_mean_field(;R, Y = 1, d = 4, kwargs...)
    tr = mean(readdlm(datadir("exp_pro","R$(R)_ST_10_$Y.dat")); dims = 2)
    x = standardize(vec(tr))
    τ = estimate_delay(x, "mi_min")
    return embed(x, d, τ)
end

function electronic_roessler(; R=1, Y = 1, dms = 1:28, kwargs...)
    tr = readdlm(datadir("exp_pro","R$(R)_ST_10_$Y.dat"))[:,dms]
    standardize(Dataset(tr))
end

end

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
export
    SWS_FAST_BILINEAR,
    SWS_BILINEAR,
    SWS_BICUBIC,
    SWS_X,
    SWS_POINT,
    SWS_AREA,
    SWS_BICUBLIN,
    SWS_GAUSS,
    SWS_SINC,
    SWS_LANCZOS,
    SWS_SPLINE,
    SWS_SRC_V_CHR_DROP_MASK,
    SWS_SRC_V_CHR_DROP_SHIFT,
    SWS_PARAM_DEFAULT,
    SWS_PRINT_INFO,
    SWS_FULL_CHR_H_INT,
    SWS_FULL_CHR_H_INP,
    SWS_DIRECT_BGR,
    SWS_ACCURATE_RND,
    SWS_BITEXACT,
    SWS_ERROR_DIFFUSION,
    SWS_CPU_CAPS_MMX,
    SWS_CPU_CAPS_MMXEXT,
    SWS_CPU_CAPS_MMX2,
    SWS_CPU_CAPS_3DNOW,
    SWS_CPU_CAPS_ALTIVEC,
    SWS_CPU_CAPS_BFIN,
    SWS_CPU_CAPS_SSE2,
    SWS_MAX_REDUCE_CUTOFF,
    SWS_CS_ITU709,
    SWS_CS_FCC,
    SWS_CS_ITU601,
    SWS_CS_ITU624,
    SWS_CS_SMPTE170M,
    SWS_CS_SMPTE240M,
    SWS_CS_DEFAULT,
    SwsVector,
    SwsFilter,
    SwsContext

const SWS_FAST_BILINEAR = 1
const SWS_BILINEAR = 2
const SWS_BICUBIC = 4
const SWS_X = 8
const SWS_POINT = 0x10
const SWS_AREA = 0x20
const SWS_BICUBLIN = 0x40
const SWS_GAUSS = 0x80
const SWS_SINC = 0x0100
const SWS_LANCZOS = 0x0200
const SWS_SPLINE = 0x0400
const SWS_SRC_V_CHR_DROP_MASK = 0x00030000
const SWS_SRC_V_CHR_DROP_SHIFT = 16
const SWS_PARAM_DEFAULT = 123456
const SWS_PRINT_INFO = 0x1000
const SWS_FULL_CHR_H_INT = 0x2000
const SWS_FULL_CHR_H_INP = 0x4000
const SWS_DIRECT_BGR = 0x8000
const SWS_ACCURATE_RND = 0x00040000
const SWS_BITEXACT = 0x00080000
const SWS_ERROR_DIFFUSION = 0x00800000
const SWS_CPU_CAPS_MMX = 0x80000000
const SWS_CPU_CAPS_MMXEXT = 0x20000000
const SWS_CPU_CAPS_MMX2 = 0x20000000
const SWS_CPU_CAPS_3DNOW = 0x40000000
const SWS_CPU_CAPS_ALTIVEC = 0x10000000
const SWS_CPU_CAPS_BFIN = 0x01000000
const SWS_CPU_CAPS_SSE2 = 0x02000000
const SWS_MAX_REDUCE_CUTOFF = 0.002
const SWS_CS_ITU709 = 1
const SWS_CS_FCC = 4
const SWS_CS_ITU601 = 5
const SWS_CS_ITU624 = 5
const SWS_CS_SMPTE170M = 5
const SWS_CS_SMPTE240M = 7
const SWS_CS_DEFAULT = 5

immutable SwsVector
    coeff::Ptr{Cdouble}
    length::Cint
end

immutable SwsFilter
    lumH::Ptr{SwsVector}
    lumV::Ptr{SwsVector}
    chrH::Ptr{SwsVector}
    chrV::Ptr{SwsVector}
end

typealias SwsContext Void


text:
qloss()=QuadLoss()
relu()=Relu()
sigm()=Sigm()
tanh()=Tanh()


text:
function construct_network!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    ::Type{CopperPlatePowerModel},
    template::OperationsProblemTemplate,
)
    buses = PSY.get_components(PSY.Bus, sys)
    bus_count = length(buses)

    if get_balance_slack_variables(optimization_container.settings)
        add_slacks!(optimization_container, CopperPlatePowerModel)
    end
    copper_plate(optimization_container, :nodal_balance_active, bus_count)
    return
end

function construct_network!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    ::Type{AreaBalancePowerModel},
    template::OperationsProblemTemplate,
)
    area_mapping = PSY.get_aggregation_topology_mapping(PSY.Area, sys)
    branches = get_available_components(PSY.Branch, sys)
    if get_balance_slack_variables(optimization_container.settings)
        throw(
            IS.ConflictingInputsError(
                "Slack Variables are not compatible with AreaBalancePowerModel",
            ),
        )
    end

    area_balance(optimization_container, :nodal_balance_active, area_mapping, branches)
    return
end

function construct_network!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    ::Type{StandardPTDFModel},
    template::OperationsProblemTemplate,
)
    buses = PSY.get_components(PSY.Bus, sys)
    ptdf = get_PTDF(optimization_container)

    if ptdf === nothing
        throw(ArgumentError("no PTDF matrix supplied"))
    end

    if get_balance_slack_variables(optimization_container.settings)
        add_slacks!(optimization_container, StandardPTDFModel)
    end

    copper_plate(optimization_container, :nodal_balance_active, length(buses))
    return
end

function construct_network!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    ::Type{T},
    template::OperationsProblemTemplate,
) where {T <: PTDFPowerModel}
    construct_network!(
        optimization_container,
        sys,
        T,
        template;
        instantiate_model = instantiate_nip_ptdf_expr_model,
    )

    add_pm_expr_refs!(optimization_container, T, sys)
    copper_plate(
        optimization_container,
        :nodal_balance_active,
        length(PSY.get_components(PSY.Bus, sys)),
    )

    return
end

function construct_network!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    ::Type{T},
    template::OperationsProblemTemplate;
    instantiate_model = instantiate_nip_expr_model,
) where {T <: PM.AbstractPowerModel}
    if T in UNSUPPORTED_POWERMODELS
        throw(
            ArgumentError(
                "$(T) formulation is not currently supported in PowerSimulations",
            ),
        )
    end

    if get_balance_slack_variables(optimization_container.settings)
        add_slacks!(optimization_container, T)
    end

    @debug "Building the $T network with $instantiate_model method"
    powermodels_network!(optimization_container, T, sys, template, instantiate_model)
    add_pm_var_refs!(optimization_container, T, sys)
    add_pm_con_refs!(optimization_container, T, sys)

    return
end

function construct_network!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    ::Type{T},
    template::OperationsProblemTemplate;
    instantiate_model = instantiate_bfp_expr_model,
) where {T <: PM.AbstractBFModel}
    if T in UNSUPPORTED_POWERMODELS
        throw(
            ArgumentError(
                "$(T) formulation is not currently supported in PowerSimulations",
            ),
        )
    end

    get_balance_slack_variables(optimization_container.settings) &&
        add_slacks!(optimization_container, T)

    @debug "Building the $T network with $instantiate_model method"
    powermodels_network!(optimization_container, T, sys, template, instantiate_model)
    add_pm_var_refs!(optimization_container, T, sys)
    add_pm_con_refs!(optimization_container, T, sys)
    return
end

function construct_network!(
    optimization_container::OptimizationContainer,
    sys::PSY.System,
    ::Type{T},
    template::OperationsProblemTemplate;
    instantiate_model = instantiate_vip_expr_model,
) where {T <: PM.AbstractIVRModel}
    if T in UNSUPPORTED_POWERMODELS
        throw(
            ArgumentError(
                "$(T) formulation is not currently supported in PowerSimulations",
            ),
        )
    end

    if get_balance_slack_variables(optimization_container.settings)
        add_slacks!(optimization_container, T)
    end

    @debug "Building the $T network with $instantiate_model method"
    powermodels_network!(optimization_container, T, sys, template, instantiate_model)
    add_pm_var_refs!(optimization_container, T, sys)
    add_pm_con_refs!(optimization_container, T, sys)
    return
end


text:
# Modelo de turbina sem sangria
type Turbina
	Turbina()=begin
		PP2=outers.PP2
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"Steam tables"
			)),
			Entalpia(),
			Eficiencia(Dict{Symbol,Any}(
				:Brief=>"Eficiencia da turbina"
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in}",
				:PosX=>0,
				:PosY=>0.25
			)),
			Potencia(Dict{Symbol,Any}(
				:Brief=>"Potencia da turbina",
				:PosX=>1,
				:PosY=>0.5
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{out}",
				:PosX=>1,
				:PosY=>1
			)),
			[
				:(H_IS = PP2.propPS(Fout.P,Fin.S)),
				:(Fout.H = (H_IS - Fin.H) * EF_T + Fin.H),
				:([Fout.S,Fout.T] = PP2.propPH(Fout.P,Fout.H)),
				:(Fin.F * (Fin.H - Fout.H) = POT_TURB),
				:(Fout.F = Fin.F),
			],
			[
				"","","","","",
			],
			[:PP2,],
			[:H_IS,:EF_T,:Fin,:POT_TURB,:Fout,]
		)
	end
	PP2::DanaPlugin
	H_IS::Entalpia
	EF_T::Eficiencia
	Fin::Corrente
	POT_TURB::Potencia
	Fout::Corrente
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Turbina
function setEquationFlow(in::Turbina)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
end
function atributes(in::Turbina,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/turbina"
	drive!(fields,_)
	return fields
end
Turbina(_::Dict{Symbol,Any})=begin
	newModel=Turbina()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
"""
    flatten(x::AbstractArray)

Reshape arbitrarly-shaped input into a matrix-shaped output,
preserving the size of the last dimension.

See also [`unsqueeze`](@ref).

# Examples
```jldoctest
julia> rand(3,4,5) |> Flux.flatten |> size
(12, 5)

julia> xs = rand(Float32, 10,10,3,7);

julia> m = Chain(Conv((3,3), 3 => 4, pad=1), Flux.flatten, Dense(400 => 33));

julia> xs |> m[1] |> size
(10, 10, 4, 7)

julia> xs |> m |> size
(33, 7)
```
"""
function flatten(x::AbstractArray)
  return reshape(x, :, size(x)[end])
end

"""
    normalise(x; dims=ndims(x), ϵ=1e-5)

Normalise `x` to mean 0 and standard deviation 1 across the dimension(s) given by `dims`.
Per default, `dims` is the last dimension.
`ϵ` is a small additive factor added to the denominator for numerical stability.
"""
@inline function normalise(x::AbstractArray; dims=ndims(x), ϵ=ofeltype(x, 1e-5))
  μ = mean(x, dims=dims)
  σ = std(x, dims=dims, mean=μ, corrected=false)
  return @. (x - μ) / (σ + ϵ)
end