text:
type DBError
  # The field contents are ERROR, FATAL, or PANIC (in an error message),
  # or WARNING, NOTICE, DEBUG, INFO, or LOG (in a notice message),
  # or a localized translation of one of these.
  severity :: String

  # The SQLSTATE code for the error
  code :: Symbol

  # The primary human-readable error message.
  # This should be accurate but terse (typically one line).
  message :: String

  # An optional secondary error message carrying
  # more detail about the problem. Might run to multiple lines.
  detail :: Nullable{String}

  # An optional suggestion what to do about the problem.
  # This is intended to differ from Detail in that it offers advice
  # (potentially inappropriate) rather than hard facts.
  # Might run to multiple lines.
  hint :: Nullable{String}

  # Either indicating an error cursor position as an index into the original
  # query string or refering to a cursor position of an internally generated
  # command rather than the one submitted by the client. The first character
  # has index 1, and positions are measured in characters not bytes.
  position :: Nullable{Int}

  # The text of a failed internally-generated command.
  # This could be, for example, a SQL query issued by a PL/pgSQL function
  query :: Nullable{String}

  # An indication of the context in which the error occurred.
  # Presently this includes a call stack traceback of active
  # procedural language functions and internally-generated queries.
  # The trace is one entry per line, most recent first.
  where :: Nullable{String}

  # If the error was associated with a specific database object,
  # the name of the schema containing that object, if any.
  scheme :: Nullable{String}

  # If the error was associated with a specific table, the name of the table.
  # (Refer to the schema name field for the name of the table's schema.)
  table :: Nullable{String}

  # If the error was associated with a specific table column, the name of the
  # column. (Refer to the schema and table name fields to identify the table.)
  column :: Nullable{String}

  # If the error was associated with a specific data type,
  # the name of the data type.
  # (Refer to the schema name field for the name of the data type's schema.)
  datatype :: Nullable{String}

  # If the error was associated with a specific constraint,
  # the name of the constraint. Refer to fields listed above for the
  # associated table or domain. (For this purpose, indexes are treated
  # as constraints, even if they weren't created with constraint syntax.)
  constraint :: Nullable{String}

  # The file name of the source-code location where the error was reported.
  file :: Nullable{String}

  # The line number of the source-code location where the error was reported.
  line :: Nullable{String}

  # The name of the source-code routine reporting the error.
  routine :: Nullable{String}

  function DBError(msg :: MSG)
    dbe = new()
    fields = Dict{Char, ByteString}
    while length(msg) > 0
      field = readuntil(msg, '\0')[1:end-1]
      if length(field) > 1
        fields[field[1]] = field[2:end]
      end
    end

    # Check for query position and if not set for internal query
    p = get_nullable(fields, 'P')
    p = isnull(p) ? p = get_nullable(fields, 'p') : p
    pos = isnull(p) ? Nullable{Int}() : Nullable(parseint(get(p)))

    dbe.severity    = fields['S']
    dbe.code        = ERROR_CODES[fields['C']]
    dbe.message     = fields['M']

    dbe.detail      = get_nullable(fields, 'D')
    dbe.hint        = get_nullable(fields, 'H')

    dbe.position    = pos

    dbe.query       = get_nullable(fields, 'q')
    dbe.where       = get_nullable(fields, 'W')
    dbe.schema      = get_nullable(fields, 's')
    dbe.table       = get_nullable(fields, 't')
    dbe.column      = get_nullable(fields, 'c')
    dbe.datatype    = get_nullable(fields, 'd')
    dbe.constraint  = get_nullable(fields, 'n')
    dbe.file        = get_nullable(fields, 'F')
    dbe.line        = get_nullable(fields, 'L')
    dbe.routine     = get_nullable(fields, 'R')
  end
end

# taken from http://www.postgresql.org/docs/9.4/static/errcodes-appendix.html
const ERROR_CODES = @compat Dict{String, Symbol}(
  # Class 00 — Successful Completion
  "00000" => :successful_completion,

  # Class 01 — Warning
  "01000" => :warning,
  "0100C" => :dynamic_result_sets_returned,
  "01008" => :implicit_zero_bit_padding,
  "01003" => :null_value_eliminated_in_set_function,
  "01007" => :privilege_not_granted,
  "01006" => :privilege_not_revoked,
  "01004" => :string_data_right_truncation,
  "01P01" => :deprecated_feature,

  # Class 02 — No Data (this is also a warning class per the SQL standard)
  "02000" => :no_data,
  "02001" => :no_additional_dynamic_result_sets_returned,

  # Class 03 — SQL Statement Not Yet Complete
  "03000" => :sql_statement_not_yet_complete,

  # Class 08 — Connection Exception
  "08000" => :connection_exception,
  "08003" => :connection_does_not_exist,
  "08006" => :connection_failure,
  "08001" => :sqlclient_unable_to_establish_sqlconnection,
  "08004" => :sqlserver_rejected_establishment_of_sqlconnection,
  "08007" => :transaction_resolution_unknown,
  "08P01" => :protocol_violation,

  # Class 09 — Triggered Action Exception
  "09000" => :triggered_action_exception,

  # Class 0A — Feature Not Supported
  "0A000" => :feature_not_supported,

  # Class 0B — Invalid Transaction Initiation
  "0B000" => :invalid_transaction_initiation,

  # Class 0F — Locator Exception
  "0F000" => :locator_exception,
  "0F001" => :invalid_locator_specification,

  # Class 0L — Invalid Grantor
  "0L000" => :invalid_grantor,
  "0LP01" => :invalid_grant_operation,

  # Class 0P — Invalid Role Specification
  "0P000" => :invalid_role_specification,

  # Class 0Z — Diagnostics Exception
  "0Z000" => :diagnostics_exception,
  "0Z002" => :stacked_diagnostics_accessed_without_active_handler,

  # Class 20 — Case Not Found
  "20000" => :case_not_found,

  # Class 21 — Cardinality Violation
  "21000" => :cardinality_violation,

  # Class 22 — Data Exception
  "22000" => :data_exception,
  "2202E" => :array_subscript_error,
  "22021" => :character_not_in_repertoire,
  "22008" => :datetime_field_overflow,
  "22012" => :division_by_zero,
  "22005" => :error_in_assignment,
  "2200B" => :escape_character_conflict,
  "22022" => :indicator_overflow,
  "22015" => :interval_field_overflow,
  "2201E" => :invalid_argument_for_logarithm,
  "22014" => :invalid_argument_for_ntile_function,
  "22016" => :invalid_argument_for_nth_value_function,
  "2201F" => :invalid_argument_for_power_function,
  "2201G" => :invalid_argument_for_width_bucket_function,
  "22018" => :invalid_character_value_for_cast,
  "22007" => :invalid_datetime_format,
  "22019" => :invalid_escape_character,
  "2200D" => :invalid_escape_octet,
  "22025" => :invalid_escape_sequence,
  "22P06" => :nonstandard_use_of_escape_character,
  "22010" => :invalid_indicator_parameter_value,
  "22023" => :invalid_parameter_value,
  "2201B" => :invalid_regular_expression,
  "2201W" => :invalid_row_count_in_limit_clause,
  "2201X" => :invalid_row_count_in_result_offset_clause,
  "22009" => :invalid_time_zone_displacement_value,
  "2200C" => :invalid_use_of_escape_character,
  "2200G" => :most_specific_type_mismatch,
  "22004" => :null_value_not_allowed,
  "22002" => :null_value_no_indicator_parameter,
  "22003" => :numeric_value_out_of_range,
  "22026" => :string_data_length_mismatch,
  "22001" => :string_data_right_truncation,
  "22011" => :substring_error,
  "22027" => :trim_error,
  "22024" => :unterminated_c_string,
  "2200F" => :zero_length_character_string,
  "22P01" => :floating_point_exception,
  "22P02" => :invalid_text_representation,
  "22P03" => :invalid_binary_representation,
  "22P04" => :bad_copy_file_format,
  "22P05" => :untranslatable_character,
  "2200L" => :not_an_xml_document,
  "2200M" => :invalid_xml_document,
  "2200N" => :invalid_xml_content,
  "2200S" => :invalid_xml_comment,
  "2200T" => :invalid_xml_processing_instruction,

  # Class 23 — Integrity Constraint Violation
  "23000" => :integrity_constraint_violation,
  "23001" => :restrict_violation,
  "23502" => :not_null_violation,
  "23503" => :foreign_key_violation,
  "23505" => :unique_violation,
  "23514" => :check_violation,
  "23P01" => :exclusion_violation,

  # Class 24 — Invalid Cursor State
  "24000" => :invalid_cursor_state,

  # Class 25 — Invalid Transaction State
  "25000" => :invalid_transaction_state,
  "25001" => :active_sql_transaction,
  "25002" => :branch_transaction_already_active,
  "25008" => :held_cursor_requires_same_isolation_level,
  "25003" => :inappropriate_access_mode_for_branch_transaction,
  "25004" => :inappropriate_isolation_level_for_branch_transaction,
  "25005" => :no_active_sql_transaction_for_branch_transaction,
  "25006" => :read_only_sql_transaction,
  "25007" => :schema_and_data_statement_mixing_not_supported,
  "25P01" => :no_active_sql_transaction,
  "25P02" => :in_failed_sql_transaction,

  # Class 26 — Invalid SQL Statement Name
  "26000" => :invalid_sql_statement_name,

  # Class 27 — Triggered Data Change Violation
  "27000" => :triggered_data_change_violation,

  # Class 28 — Invalid Authorization Specification
  "28000" => :invalid_authorization_specification,
  "28P01" => :invalid_password,

  # Class 2B — Dependent Privilege Descriptors Still Exist
  "2B000" => :dependent_privilege_descriptors_still_exist,
  "2BP01" => :dependent_objects_still_exist,

  # Class 2D — Invalid Transaction Termination
  "2D000" => :invalid_transaction_termination,

  # Class 2F — SQL Routine Exception
  "2F000" => :sql_routine_exception,
  "2F005" => :function_executed_no_return_statement,
  "2F002" => :modifying_sql_data_not_permitted,
  "2F003" => :prohibited_sql_statement_attempted,
  "2F004" => :reading_sql_data_not_permitted,

  # Class 34 — Invalid Cursor Name
  "34000" => :invalid_cursor_name,

  # Class 38 — External Routine Exception
  "38000" => :external_routine_exception,
  "38001" => :containing_sql_not_permitted,
  "38002" => :modifying_sql_data_not_permitted,
  "38003" => :prohibited_sql_statement_attempted,
  "38004" => :reading_sql_data_not_permitted,

  # Class 39 — External Routine Invocation Exception
  "39000" => :external_routine_invocation_exception,
  "39001" => :invalid_sqlstate_returned,
  "39004" => :null_value_not_allowed,
  "39P01" => :trigger_protocol_violated,
  "39P02" => :srf_protocol_violated,

  # Class 3B — Savepoint Exception
  "3B000" => :savepoint_exception,
  "3B001" => :invalid_savepoint_specification,

  # Class 3D — Invalid Catalog Name
  "3D000" => :invalid_catalog_name,

  # Class 3F — Invalid Schema Name
  "3F000" => :invalid_schema_name,

  # Class 40 — Transaction Rollback
  "40000" => :transaction_rollback,
  "40002" => :transaction_integrity_constraint_violation,
  "40001" => :serialization_failure,
  "40003" => :statement_completion_unknown,
  "40P01" => :deadlock_detected,

  # Class 42 — Syntax Error or Access Rule Violation
  "42000" => :syntax_error_or_access_rule_violation,
  "42601" => :syntax_error,
  "42501" => :insufficient_privilege,
  "42846" => :cannot_coerce,
  "42803" => :grouping_error,
  "42P20" => :windowing_error,
  "42P19" => :invalid_recursion,
  "42830" => :invalid_foreign_key,
  "42602" => :invalid_name,
  "42622" => :name_too_long,
  "42939" => :reserved_name,
  "42804" => :datatype_mismatch,
  "42P18" => :indeterminate_datatype,
  "42P21" => :collation_mismatch,
  "42P22" => :indeterminate_collation,
  "42809" => :wrong_object_type,
  "42703" => :undefined_column,
  "42883" => :undefined_function,
  "42P01" => :undefined_table,
  "42P02" => :undefined_parameter,
  "42704" => :undefined_object,
  "42701" => :duplicate_column,
  "42P03" => :duplicate_cursor,
  "42P04" => :duplicate_database,
  "42723" => :duplicate_function,
  "42P05" => :duplicate_prepared_statement,
  "42P06" => :duplicate_schema,
  "42P07" => :duplicate_table,
  "42712" => :duplicate_alias,
  "42710" => :duplicate_object,
  "42702" => :ambiguous_column,
  "42725" => :ambiguous_function,
  "42P08" => :ambiguous_parameter,
  "42P09" => :ambiguous_alias,
  "42P10" => :invalid_column_reference,
  "42611" => :invalid_column_definition,
  "42P11" => :invalid_cursor_definition,
  "42P12" => :invalid_database_definition,
  "42P13" => :invalid_function_definition,
  "42P14" => :invalid_prepared_statement_definition,
  "42P15" => :invalid_schema_definition,
  "42P16" => :invalid_table_definition,
  "42P17" => :invalid_object_definition,

  # Class 44 — WITH CHECK OPTION Violation
  "44000" => :with_check_option_violation,

  # Class 53 — Insufficient Resources
  "53000" => :insufficient_resources,
  "53100" => :disk_full,
  "53200" => :out_of_memory,
  "53300" => :too_many_connections,
  "53400" => :configuration_limit_exceeded,

  # Class 54 — Program Limit Exceeded
  "54000" => :program_limit_exceeded,
  "54001" => :statement_too_complex,
  "54011" => :too_many_columns,
  "54023" => :too_many_arguments,

  # Class 55 — Object Not In Prerequisite State
  "55000" => :object_not_in_prerequisite_state,
  "55006" => :object_in_use,
  "55P02" => :cant_change_runtime_param,
  "55P03" => :lock_not_available,

  # Class 57 — Operator Intervention
  "57000" => :operator_intervention,
  "57014" => :query_canceled,
  "57P01" => :admin_shutdown,
  "57P02" => :crash_shutdown,
  "57P03" => :cannot_connect_now,
  "57P04" => :database_dropped,

  # Class 58 — System Error (errors external to PostgreSQL itself)
  "58000" => :system_error,
  "58030" => :io_error,
  "58P01" => :undefined_file,
  "58P02" => :duplicate_file,

  # Class F0 — Configuration File Error
  "F0000" => :config_file_error,
  "F0001" => :lock_file_exists,

  # Class HV — Foreign Data Wrapper Error (SQL/MED)
  "HV000" => :fdw_error,
  "HV005" => :fdw_column_name_not_found,
  "HV002" => :fdw_dynamic_parameter_value_needed,
  "HV010" => :fdw_function_sequence_error,
  "HV021" => :fdw_inconsistent_descriptor_information,
  "HV024" => :fdw_invalid_attribute_value,
  "HV007" => :fdw_invalid_column_name,
  "HV008" => :fdw_invalid_column_number,
  "HV004" => :fdw_invalid_data_type,
  "HV006" => :fdw_invalid_data_type_descriptors,
  "HV091" => :fdw_invalid_descriptor_field_identifier,
  "HV00B" => :fdw_invalid_handle,
  "HV00C" => :fdw_invalid_option_index,
  "HV00D" => :fdw_invalid_option_name,
  "HV090" => :fdw_invalid_string_length_or_buffer_length,
  "HV00A" => :fdw_invalid_string_format,
  "HV009" => :fdw_invalid_use_of_null_pointer,
  "HV014" => :fdw_too_many_handles,
  "HV001" => :fdw_out_of_memory,
  "HV00P" => :fdw_no_schemas,
  "HV00J" => :fdw_option_name_not_found,
  "HV00K" => :fdw_reply_handle,
  "HV00Q" => :fdw_schema_not_found,
  "HV00R" => :fdw_table_not_found,
  "HV00L" => :fdw_unable_to_create_execution,
  "HV00M" => :fdw_unable_to_create_reply,
  "HV00N" => :fdw_unable_to_establish_connection,

  # Class P0 — PL/pgSQL Error
  "P0000" => :plpgsql_error,
  "P0001" => :raise_exception,
  "P0002" => :no_data_found,
  "P0003" => :too_many_rows,

  # Class XX — Internal Error
  "XX000" => :internal_error,
  "XX001" => :data_corrupted,
  "XX002" => :index_corrupted,
)


text:
# # Basic Ploting Tutorial
# This tutorial lists the examples done in the basic Plots tutorial
# **Note**: it is possible that the functionality available in plots might not be available yet in Makie.
# The orignal tutorial can be found [here](http://docs.juliaplots.org/latest/tutorial/)

# Importing required packages
using Makie
using DataFrames
using StatsMakie

# ## Basic Plotting: Line Plots
scene = Scene()
lines!(scene, rand(10), color=:blue)
lines!(scene, rand(10), color=:orange)
lines!(scene, rand(10), color=:green)

display(scene)

# The `Scene` object holds everything in a plot. We can think of it like a container for everything.
# ### OUTPUT:
# ()[https://github.com/Akshat-mehrotra/codein/blob/master/JULIA/Makie%20Images/plot1.png]
#
# ## Plot Attributes
# In Makie.jl, the modifiers to plots are called attributes.
# These are documented at the [attributes page](http://makie.juliaplots.org/dev/plot-attributes.html)
# Here, we are using the attribute `linewidth` which contols the thickness of the line.
# We are also adding `legend` and `axisnames` to the scene.
# legend is not a ! function i.e. it needs to be concatinated to the scene using hbox/vbox
# We can edit the axis name and other attributes by just passing `Axis` as index to scene.

scene = Scene()
lines!(scene, rand(10), linewidth=3.0, color=:red)
lines!(scene, rand(10), linewidth=3.0, color=:green)
lines!(scene, rand(10), linewidth=3.0, color=:orange)

leg = legend([scene[i] for i in 2:length(scene)], ["line $(jk-1)" for jk in 2:length(scene)])

axis = scene[Axis]
axis[:names, :axisnames] = ("x label", "", "")

title_plot = title(scene,"Three Lines")
s = vbox(title_plot, leg)

display(s)
#
# ### OUTPUT:
# ()[https://github.com/Akshat-mehrotra/codein/blob/master/JULIA/Makie%20Images/plot2.png]
#
# ## Changing the Plotting Type
# Here we are introducting a new way of plotting data, a scatter plot
# We can easily generate scatter plots by using `scatter` or `scatter!`, if we want to add it to an existing scene
# In Makie a rule of thumb is that whenever a plotting function ends with a bang, a scene can be passed into that to plot on it.

scene = Scene()
scatter!(scene,rand(10), color=:red)
scatter!(scene, rand(10), color=:blue)

leg = legend([scene[i] for i in 2:length(scene)], ["point $(jk-1)" for jk in 2:length(scene)])

tit = title(scene,"My Scatter Plot")

s = hbox(leg, tit)
display(s)

# ### OUTPUT:
# ()[https://github.com/Akshat-mehrotra/codein/blob/master/JULIA/Makie%20Images/plot3.png]

# ## Combining Multiple Plots as Subplots
# We can combine multiple plots together as subplots using hbox/vbox. There are many methods for doing this, and we will show two simple methods for generating simple layouts
# hbox places the plots on top of each other and vbox places them side by side
# Example 1:
l1 = lines(rand(10), linewidth=3.0, color=:red)
l2 = lines(rand(10), linewidth=3.0, color=:green)
l3 = lines(rand(10), linewidth=3.0, color=:orange)

s = hbox(l1,l2,l3)
display(s)

# ### OUTPUT:
# ()[https://github.com/Akshat-mehrotra/codein/blob/master/JULIA/Makie%20Images/plot4.png]
# In the next way, we will create 4 diffrent scenes and place them in a 4x4 box
scene1 = Scene()
lines!(scene1, rand(10), color=:red)
lines!(scene1, rand(10), color=:blue)
lines!(scene1, rand(10), color=:green)
t = title(scene1,"subtitle")

scene2 = Scene()
lines!(scene2, rand(10), linewidth=3.0, color=:red)
lines!(scene2, rand(10), linewidth=3.0, color=:blue)
lines!(scene2, rand(10), linewidth=3.0, color=:green)
text!(scene2,"This one is labelled",textsize = 1)

scene3 = Scene()
scatter!(scene3, rand(10), color=:red)
scatter!(scene3, rand(10), color=:blue)

scene4 = Scene()
plot!(scene4,histogram, rand(10), color=:blue)
plot!(scene4,histogram, rand(10), color=:red)
plot!(scene4,histogram, rand(10), color=:green)

s = vbox(hbox(scene2,t),hbox(scene4,scene3))
display(s)

# ### OUTPUT:
# ()[https://github.com/Akshat-mehrotra/codein/blob/master/JULIA/Makie%20Images/plot5.png]
#
# ## Reading data from DataFrame
# We can use the `StatsMakie` package to help us with this

df = DataFrame(a = 1:10, b = 10*rand(10), c = 10 * rand(10))
scene = Scene()
lines!(scene, Data(df), :a, :b, color=:orange)
lines!(scene, Data(df), :a, :c, color=:blue)
display(scene)

# ### OUTPUT:
# ()[https://github.com/Akshat-mehrotra/codein/blob/master/JULIA/Makie%20Images/plot6.png]
#
# Scatter plotting the data from DataFrame
scene = Scene()
scatter!(scene, Data(df), :a, :b, color=:blue)
display(scene)
# ### OUTPUT:
# ()[https://github.com/Akshat-mehrotra/codein/blob/master/JULIA/Makie%20Images/plot7.png]
#

# Violin and Box Plots of data from DataFrames
using RDatasets
scene = Scene()

d = dataset("Ecdat", "Fatality");

violin!(scene,Data(d), :Year, :Perinc, color=:gray)

display(scene)

# ### OUTPUT:
# ()[https://github.com/Akshat-mehrotra/codein/blob/master/JULIA/Makie%20Images/plot8.png]
#
boxplot!(d.Year, d.Perinc, color = :black)
display(scene)
# ### OUTPUT:
# ()[https://github.com/Akshat-mehrotra/codein/blob/master/JULIA/Makie%20Images/plot9.png]

ext: jl


text:
## Setup
import Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using Measures, LaTeXStrings

using Plots
pgfplotsx()
default(size=(350, 240))

cd(@__DIR__)
include("readdata.jl")
mkpath("./fig")

## Load 2v2 data
#exp_path = joinpath(homedir(), "servicemesh_results", "simpleflipsplit2v2", "HO_220201_194358_new_rand_alg")
exp_path = joinpath(homedir(), "servicemesh_results", "simpleflipsplit2v2uneven", "custom", "2v2uneven_min1scale")
microservices = 4

data_paths = read_all(exp_path)

datas = [loaddata(data_paths[k]) for k in keys(data_paths)]
data_names = ["Kubernetes" "Mesh aware" "RL agent" "Reactive"]

smoothing_window = 1000
downsampling = 1000
days = 100
final_sample = days * 24 * 60 # We sample every minute
sample_range = 1:downsampling:final_sample

t = gettime(datas[1])[sample_range]

## Reward plot
reward_data = getvalues.(datas, "training/reward"; smoothing_window, sample_range)
reward_plot = plot(t, reward_data, title="Reward", xaxis="Time [Days]", ylim=[20, 40], label=data_names, legend=:bottomright, legend_column=2)
savefig(reward_plot, "fig/double_reward.tex")

## Missed jobs
missed_data = getvalues.(datas, "Job1/missed"; smoothing_window, sample_range)
missed_data .+= getvalues.(datas, "Job2/missed"; smoothing_window, sample_range)
missed_plot = plot(t, missed_data, title="Missed deadlines", yaxis="Missed [jobs/s]", xaxis="Time [Days]", ylim=[0, 0.2], label=data_names, legend=:topright)
savefig(missed_plot, "fig/double_missed.tex")

## Arrivals
arrivals_data = [getvalues(datas[1], "Job$i/arrivals"; smoothing_window, sample_range) for i in 1:2]
arrivals_data_full = [getvalues(datas[1], "Job$i/arrivals"; smoothing_window=1, sample_range) for i in 1:2]
arrivals_plot = plot(t, arrivals_data, color=[:black :purple], linestyle=[:solid :dash], title="Job arrivals", yaxis="[jobs/s]", xaxis="Time [Days]", ylim=[0, 4], label=[L"W_1" L"W_2"], legend=:topright, size=(230,170))
plot!(arrivals_plot, t, arrivals_data_full, color=[:black :purple], linestyle=[:solid :dash], alpha=0.3, label="")
savefig(arrivals_plot, "fig/double_arrivals.tex")

## Combined scale and util but separate ms 4
begin
    ## M1
    p = plot(layout=(2, 1), show=true, size=(250, 200), link=:x, bottommargin=-7mm, topmargin=-7mm)

    data = getvalues.(datas, "MS1/scale"; smoothing_window, sample_range)
    plot!(p[1], t, data, yaxis="Scale", label="", xformatter=_->"", title=L"M_1", yrange=(0.5, 7.0))

    data = getvalues.(datas, "MS1/utilization"; smoothing_window, sample_range)
    plot!(p[2], t, [min.(1, d) for d in data], xformatter=_->"", label="", yaxis="Utilization", yrange=(0.3, 1.0))

    savefig(p, "fig/double_scale_util_1.tex")
    display(p)

    ## M2
    p = plot(layout=(2, 1), show=true, size=(250, 200), link=:x, bottommargin=-7mm, topmargin=-7mm)

    data = getvalues.(datas, "MS2/scale"; smoothing_window, sample_range)
    plot!(p[1], t, data, label="", xformatter=_->"", yformatter=_->"", title=L"M_2", yrange=(0.5, 7.0))

    data = getvalues.(datas, "MS2/utilization"; smoothing_window, sample_range)
    plot!(p[2], t, [min.(1, d) for d in data], xformatter=_->"", yformatter=_->"", label="", yrange=(0.3, 1.0))

    savefig(p, "fig/double_scale_util_2.tex")
    display(p)

    ## M3
    p = plot(layout=(2, 1), show=true, size=(250, 200), link=:x, bottommargin=-7mm, topmargin=-7mm)

    data = getvalues.(datas, "MS3/scale"; smoothing_window, sample_range)
    plot!(p[1], t, data, label=data_names, xformatter=_->"", yaxis="Scale", title=L"M_3", yrange=(0.5, 7.0), legend=:outerleft, legend_column=1)

    data = getvalues.(datas, "MS3/utilization"; smoothing_window, sample_range)
    plot!(p[2], t, [min.(1, d) for d in data], xaxis="Time [days]", label="", yaxis="Utilization", yrange=(0.3, 1.0))

    savefig(p, "fig/double_scale_util_3.tex")
    display(p)

    ## M4
    p = plot(layout=(2, 1), show=true, size=(250, 200), link=:x, bottommargin=-7mm, topmargin=-7mm)

    data = getvalues.(datas, "MS4/scale"; smoothing_window, sample_range)
    plot!(p[1], t, data, label="", xformatter=_->"", yformatter=_->"", title=L"M_4", yrange=(0.5, 7.0))

    data = getvalues.(datas, "MS4/utilization"; smoothing_window, sample_range)
    plot!(p[2], t, [min.(1, d) for d in data], yformatter=_->"", xaxis="Time [days]", label="", yrange=(0.3, 1.0))

    savefig(p, "fig/double_scale_util_4.tex")
    display(p)
end

## Load simple data 2
exp_path = joinpath(homedir(), "servicemesh_results", "simpleflip", "HO_220224_154551", "2longsimpleflip")
microservices = 2

data_paths = read_all(exp_path)

datas = [loaddata(data_paths[k]) for k in keys(data_paths)]

## Reward plot
reward_data = getvalues.(datas, "training/reward"; smoothing_window, sample_range)
reward_plot = plot(t, reward_data, title="Reward", xaxis="Time [Days]", ylim=[7, 13], label=data_names, legend=:topleft, legend_column=2)
savefig(reward_plot, "fig/simple2_reward.tex")

## Arrival jobs
arrival_data = getvalues(datas[1], "Job1/arrivals"; smoothing_window, sample_range)
arrival_full_data = getvalues(datas[1], "Job1/arrivals"; smoothing_window=1, sample_range)
arrival_plot = plot(t, arrival_full_data, title="Job arrivals", yaxis="[jobs/s]", xaxis="Time [Days]", linecolor=:black, alpha=0.1, ylim=[0, 4], label="Sampled data", size=(250, 200), legend=:outertop)
plot!(arrival_plot, t, arrival_data, label="Smoothed data", linecolor=:black)
savefig(arrival_plot, "fig/simple2_arrivals.tex")

## Lost jobs
missed_data = getvalues.(datas, "Job1/missed"; smoothing_window, sample_range)
missed_plot = plot(t, missed_data, title="Missed deadlines", yaxis="[jobs/s]", xaxis="Time [Days]", ylim=[0, 0.1], label=data_names, legend=:topright)
savefig(missed_plot, "fig/simple2_missed.tex")
dropped_data = getvalues.(datas, "Job1/dropped"; smoothing_window, sample_range)
lost_plot = plot(t, missed_data .+ dropped_data, title="Lost jobs", yaxis="[jobs/s]", xaxis="Time [Days]", ylim=[0, 0.1], label=data_names, legend=:topright)
savefig(lost_plot, "fig/simple2_lost.tex")

## Lost/arrival jobs
lost_arrival_data = [100 .* lost_data[i] ./ (arrival_data[i] .+ 0.0001) for i in 1:4]
lost_arrival_plot = plot(t, lost_arrival_data, title="Percentage of lost jobs", yaxis="[%]", xaxis="Time [Days]", ylim=[0, 5], label=data_names, legend=:topright)
savefig(lost_arrival_plot, "fig/simple2_lost_arrival.tex")

## Combined scale and util but separate ms 2
begin
    p = plot(layout=(2, 1), yrange=(1.0, 6.0), show=true, size=(250, 200), link=:x, bottommargin=-5mm, topmargin=-5mm)

    data = getvalues.(datas, "MS1/scale"; smoothing_window, sample_range)
    plot!(p[1], t, data, yaxis="Scale", label=data_names, xformatter=_->"", legend=:outertop, title=L"M_1", yrange=(1.0, 6.0), legend_column=1)

    # util
    data = getvalues.(datas, "MS1/utilization"; smoothing_window, sample_range)
    plot!(p[2], t, [min.(1, d) for d in data], xaxis="Time [days]", label="", yaxis="Utilization", yrange=(0.4, 1.0))

    savefig(p, "fig/simple2_scale_util_1.tex")
    display(p)

    p = plot(layout=(2, 1), yrange=(1.0, 6.0), show=true, size=(250, 200), link=:x, bottommargin=-5mm, topmargin=-5mm)

    data = getvalues.(datas, "MS2/scale"; smoothing_window, sample_range)
    plot!(p[1], t, data, label="", xformatter=_->"", yformatter=_->"", legend=:outertop, title=L"M_2", yrange=(1.0, 6.0), legend_column=1)

    # util
    data = getvalues.(datas, "MS2/utilization"; smoothing_window, sample_range)
    plot!(p[2], t, [min.(1, d) for d in data], yformatter=_->"", xaxis="Time [days]", label="", yrange=(0.4, 1.0))

    savefig(p, "fig/simple2_scale_util_2.tex")
    display(p)
end

## Load simple data 3
exp_path = joinpath(homedir(), "servicemesh_results", "simpleflipL3min1", "custom", "len3min1nofracval")
microservices = 3

data_paths = read_all(exp_path)

datas = [loaddata(data_paths[k]) for k in keys(data_paths)]

## Reward plot
reward_data = getvalues.(datas, "training/reward"; smoothing_window, sample_range)
reward_plot = plot(t, reward_data, title="Reward", xaxis="Time [Days]", ylim=[8, 19], label=data_names, legend=:bottomright, legend_column=2)
savefig(reward_plot, "fig/simple_reward.tex")

## Arrival jobs
arrival_data = getvalues(datas[1], "Job1/arrivals"; smoothing_window, sample_range)
arrival_full_data = getvalues(datas[1], "Job1/arrivals"; smoothing_window=1, sample_range)
arrival_plot = plot(t, arrival_full_data, title="Job arrivals", yaxis="[jobs/s]", xaxis="Time [Days]", linecolor=:black, alpha=0.1, ylim=[0, 4], label="Sampled data", size=(250, 200), legend=:top)
plot!(arrival_plot, t, arrival_data, label="Smoothed data", linecolor=:black)
savefig(arrival_plot, "fig/simple_arrivals.tex")

## Missed jobs
missed_data = getvalues.(datas, "Job1/missed"; smoothing_window, sample_range)
missed_plot = plot(t, missed_data, title="Missed deadlines", yaxis="Missed [jobs/s]", xaxis="Time [Days]", ylim=[0, 0.1], label=data_names, legend=:topright)
savefig(missed_plot, "fig/simple_missed.tex")

## Combined scale and util but separate ms 2
begin
    ## M1
    p = plot(layout=(2, 1), yrange=(1.0, 6.0), show=true, size=(250, 200), link=:x, bottommargin=-5mm, topmargin=-5mm)

    data = getvalues.(datas, "MS1/scale"; smoothing_window, sample_range)
    plot!(p[1], t, data, yaxis="Scale", label="", xformatter=_->"", title=L"M_1", yrange=(1.0, 5.0))

    data = getvalues.(datas, "MS1/utilization"; smoothing_window, sample_range)
    plot!(p[2], t, [min.(1, d) for d in data], xaxis="Time [days]", label="", yaxis="Utilization", yrange=(0.3, 0.9))

    savefig(p, "fig/simple_scale_util_1.tex")
    display(p)

    ## M2
    p = plot(layout=(2, 1), yrange=(1.0, 6.0), show=true, size=(250, 200), link=:x, bottommargin=-5mm, topmargin=-5mm)

    data = getvalues.(datas, "MS2/scale"; smoothing_window, sample_range)
    plot!(p[1], t, data, label="", xformatter=_->"", yformatter=_->"", title=L"M_2", yrange=(1.0, 5.0))

    data = getvalues.(datas, "MS2/utilization"; smoothing_window, sample_range)
    plot!(p[2], t, [min.(1, d) for d in data], yformatter=_->"", xaxis="Time [days]", label="", yrange=(0.3, 0.9))

    savefig(p, "fig/simple_scale_util_2.tex")
    display(p)

    ## M3
    p = plot(layout=(2, 1), yrange=(1.0, 6.0), show=true, size=(250, 200), link=:x, bottommargin=-5mm, topmargin=-5mm)

    data = getvalues.(datas, "MS3/scale"; smoothing_window, sample_range)
    plot!(p[1], t, data, label=data_names, xformatter=_->"", yformatter=_->"", title=L"M_3", yrange=(1.0, 5.0), legend=:right, legend_column=1)

    data = getvalues.(datas, "MS3/utilization"; smoothing_window, sample_range)
    plot!(p[2], t, [min.(1, d) for d in data], yformatter=_->"", xaxis="Time [days]", label="", yrange=(0.3, 0.9))

    savefig(p, "fig/simple_scale_util_3.tex")
    display(p)
end

## Read ho data
data, labels = load_ho_values(joinpath(homedir(), "servicemesh_results", "simpleflipL3", "HO_220309_105934_len3", "ho_search", "ho_data.csv"))
idxs = [3, 5]
ho_plot = scatter(data[:, idxs], data[:, end], xlabel=reshape(labels[idxs], 1, :), layout=(1, 2), size=(320, 200), label="", alpha=0.2, show=true, margin=-5mm)
plot!(ho_plot[1], ylabel="Relative reward", title="Hyperparameter tuning")
savefig(ho_plot, "fig/ho.tex")

nothing

max_stars_repo_name: albheim/ServiceMeshControl

max_issues_repo_name: albheim/ServiceMeshControl

max_forks_repo_name: albheim/ServiceMeshControl


text:
function strategy_positive_mom_binary(t::Returns, lookback::Integer)::Weights
    return @pipe moving(mean, t, lookback) |>
        TimeSeries.map((timestamp, value) -> (timestamp, value >= 0 ? 1 : 0), _) |>
        TimeSeries.lag(_, 1)
end

function strategy_prive_above_sma_binary(t::Returns, lookback::Integer)::Weights
    rets = TimeSeries.cumprod(1 .+ t, dims = 1)
    sma = moving(mean, rets, lookback)
    function islarger(lhs, rhs) lhs > rhs ? 1.0 : 0.0 end
    return TimeSeries.lag(islarger.(rets, sma), 1)
end

function strategy_vol_target_mean(t::Returns, lookback::Integer)::Weights
    target = values(std(t))[1] * sqrt(252)
    return @pipe moving(std, t, lookback) |>
        TimeSeries.map((timestamp, value) -> (timestamp, target / (value * sqrt(252))), _) |>
        TimeSeries.lag(_, 1)
end

function strategy_vol_target(target::Float64)
    function a(t::Returns, lookback::Integer)::Weights
        return @pipe moving(std, t, lookback) |>
            TimeSeries.map((timestamp, value) -> (timestamp, target / (value * sqrt(252))), _) |>
            TimeSeries.lag(_, 1)
    end
end

function strategy_vol_threshold(threshold::Float64)
    function a(t::Returns, lookback::Integer)::Weights
        return @pipe TimeSeries.lag(moving(std, t, lookback), 1) |>
            TimeSeries.map((timestamp, value) -> (timestamp, value <= threshold ? 1 : 0), _)
    end
end

function strategy_semideviation_target_mean(t::Returns, lookback::Integer)::Weights
    target = values(semideviation(t)) * sqrt(252)
    return @pipe moving(semideviation, t, lookback) |>
        TimeSeries.map((timestamp, value) -> (timestamp, target / (value * sqrt(252))), _) |>
        TimeSeries.lag(_, 1)
end

function strategy_semideviation_threshold(threshold::Float64)
    function a(t::Returns, lookback::Integer)::Weights
        return @pipe TimeSeries.lag(moving(semideviation, t, lookback), 1) |>
            TimeSeries.map((timestamp, value) -> (timestamp, value <= threshold ? 1 : 0), _)
    end
end


text:
"""
	grdvector(cmd0::String="", arg1=nothing, arg2=nothing, kwargs...)

Takes two 2-D grid files which represents the x- and y-components of a vector field and produces
a vector field plot by drawing vectors with orientation and length according to the information
in the files. Alternatively, polar coordinate r, theta grids may be given instead.

Full option list at [`grdvector`]($(GMTdoc)grdvector.html)

Parameters
----------

- **A** | **polar** :: [Type => Bool]

    The grid contain polar (r, theta) components instead of Cartesian (x, y) [Default is Cartesian components].
    ($(GMTdoc)grdvector.html#a)
- $(GMT.opt_B)
- $(GMT.opt_C)
- **G** | **fill** :: [Type => Str | Number]

    Sets color or shade for vector interiors [Default is no fill].
    ($(GMTdoc)grdvector.html#g)
- **I** | **inc** :: [Type => Sytr | Number]	``Arg=[x]dx[/dy]``

    Only plot vectors at nodes every x_inc, y_inc apart (must be multiples of original grid spacing).
    ($(GMTdoc)grdvector.html#i)
- **N** | **noclip** | **no_clip** :: [Type => Bool]

    Do NOT clip symbols that fall outside map border
    ($(GMTdoc)grdvector.html#n)
- **Q** | **vec** | **vector** | **arrow** :: [Type => Str]

    Modify vector parameters. For vector heads, append vector head size [Default is 0, i.e., stick-plot].
    ($(GMTdoc)grdvector.html#q)
- $(GMT.opt_P)
- $(GMT.opt_R)
- **S** | **vec_scale** :: [Type => Str | Number]		``Arg = [i|l]scale[unit]``

    Sets scale for vector plot length in data units per plot distance measurement unit [1].
    ($(GMTdoc)grdvector.html#s)
- **T** | **sign_scale** :: [Type => Bool]

    Means the azimuths of Cartesian data sets should be adjusted according to the signs of the
    scales in the x- and y-directions [Leave alone].
    ($(GMTdoc)grdvector.html#t)
- $(GMT.opt_U)
- $(GMT.opt_V)
- **W** | **pen** :: [Type => Str | Number]

    Sets the attributes for the particular line.
    ($(GMTdoc)grdvector.html#w)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- **Z** | **azimuth** :: [Type => Bool]

    The theta grid provided contains azimuths rather than directions (implies -A).
    ($(GMTdoc)grdvector.html#z)
- $(GMT.opt_V)
- $(GMT.opt_f)
"""
function grdvector(cmd0::String="", arg1=nothing, arg2=nothing; first=true, kwargs...)

	length(kwargs) == 0 && return monolitic("grdvector", cmd0, arg1, arg2)

	d, K, O = init_module(first, kwargs...)		# Also checks if the user wants ONLY the HELP mode

	cmd::String = parse_BJR(d, "", "", O, " -JX12c/0")[1]
	cmd = parse_common_opts(d, cmd, [:I :UVXY :f :p :t :params], first)[1]
	cmd = parse_these_opts(cmd, d, [[:A :polar], [:N :noclip :no_clip], [:S :vec_scale], [:T :sign_scale], [:Z :azimuth]])

    # Check case in which the two grids were transmitted by name.
    (cmd0 != "" && isa(arg1, String)) && (cmd0 *= " " * arg1; arg1 = nothing)

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, arg1)	# Find how data was transmitted

	N_used = got_fname == 0 ? 1 : 0		# To know whether a cpt will go to arg1 or arg2
	cmd, arg1, arg2, = add_opt_cpt(d, cmd, CPTaliases, 'C', N_used, arg1, arg2)
	opt_Q = parse_Q_grdvec(d, [:Q :vec :vector :arrow])
	!occursin(" -G", opt_Q) && (cmd = add_opt_fill(cmd, d, [:G :fill], 'G'))	# If fill not passed in arrow, try from regular option
	cmd *= add_opt_pen(d, [:W :pen], "W", true)									# TRUE to also seek (lw,lc,ls)
	(!occursin(" -C", cmd) && !occursin(" -W", cmd) && !occursin(" -G", opt_Q)) && (cmd *= " -W0.5")	# If still nothing, set -W.
	(opt_Q != "") && (cmd *= opt_Q)

    return finish_PS_module(d, "grdvector " * cmd, "", K, O, true, arg1, arg2)
end

# ---------------------------------------------------------------------------------------------------
function parse_Q_grdvec(d::Dict, symbs::Array{<:Symbol})::String
	(show_kwargs[1]) && return print_kwarg_opts(symbs, "NamedTuple | String")
	cmd::String = ""
    if ((val = find_in_dict(d, symbs)[1]) !== nothing)
		if (isa(val, String))  cmd *= " -Q" * val		# An hard core GMT string directly with options
		else                   cmd *= " -Q" * vector_attrib(val)
		end
		if ((ind = findfirst("+g", cmd)) !== nothing)   # -Q0.4+e+gred+n0.4+pcyan+h0
			cmd *= " -G" * split(cmd[ind[1]+2:end], "+")[1]	# Add a -G (does the same) to have at least one of the -G, -W, -C
		end
	end
	return cmd
end

# ---------------------------------------------------------------------------------------------------
grdvector(arg1, arg2=nothing, cmd0::String="";  kw...) = grdvector(cmd0, arg1, arg2; first=true, kw...)
grdvector!(cmd0::String="", arg1=nothing, arg2=nothing; kw...) = grdvector(cmd0, arg1, arg2; first=false, kw...)
grdvector!(arg1, arg2=nothing, cmd0::String=""; kw...) = grdvector(cmd0, arg1, arg2; first=false, kw...)


text:
module SymbolServer

using Serialization, Pkg
include("from_static_lint.jl")

const storedir = abspath(joinpath(@__DIR__, "..", "..", "store"))
const c = Pkg.Types.Context()
const depot = Dict("manifest" => c.env.manifest,
                    "installed" => (VERSION < v"1.1.0-DEV.857" ? c.env.project["deps"] : c.env.project.deps),
                    "packages" => Dict{String,Any}())

if Sys.isunix()
    global const nullfile = "/dev/null"
elseif Sys.iswindows()
    global const nullfile = "nul"
else
    error("Platform not supported")
end

while true
    message, payload = deserialize(stdin)

    try
        if message == :debugmessage
            @info(payload)
            serialize(stdout, (:success, nothing))
        elseif message == :close
            break
        elseif message == :get_installed_packages_in_env
            pkgs = (VERSION < v"1.1.0-DEV.857" ? c.env.project["deps"] : c.env.project.deps)
            serialize(stdout, (:success, pkgs))
        elseif message == :get_core_packages
            core_pkgs = load_core()["packages"]
            SymbolServer.save_store_to_disc(core_pkgs["Base"], joinpath(storedir, "Base.jstore"))
            SymbolServer.save_store_to_disc(core_pkgs["Core"], joinpath(storedir, "Core.jstore"))
            serialize(stdout, (:success, nothing))
        elseif message == :get_all_packages_in_env
            pkgs = if VERSION < v"1.1.0-DEV.857"
                 Dict{String,Vector{String}}(n=>(p->get(p, "uuid", "")).(v) for (n,v) in c.env.manifest)
            else
                 Dict{String,Vector{String}}(v.name=>[string(n),] for (n,v) in c.env.manifest)
            end
            serialize(stdout, (:success, pkgs))
        elseif message == :load_package
            redirect_stdout(open(nullfile, "w")) do # seems necessary incase packages print on startup
                SymbolServer.import_package(payload, depot)
            end
            for  (uuid, pkg) in depot["packages"]
                SymbolServer.save_store_to_disc(pkg, joinpath(storedir, "$uuid.jstore"))
            end
            serialize(stdout, (:success, collect(keys(depot["packages"]))))
        elseif message == :load_all
            for pkg in (VERSION < v"1.1.0-DEV.857" ? c.env.project["deps"] : c.env.project.deps)
                SymbolServer.import_package(pkg, depot)
            end
            for  (uuid, pkg) in depot["packages"]
                SymbolServer.save_store_to_disc(pkg, joinpath(storedir, "$uuid.jstore"))
            end
            serialize(stdout, (:success, collect(keys(depot["packages"]))))
        else
            serialize(stdout, (:failure, nothing))
        end
    catch err
        serialize(stdout, (:failure, err))
    end
end

end


text:
"""
    FAM.jl

Description:
    Options, structures, and logic for the Fuzzy ARTMAP (FAM) module.
"""

"""
    opts_FAM()

Implements a Fuzzy ARTMAP learner's options.

# Examples
```julia-repl
julia> my_opts = opts_FAM()
```
"""
@with_kw mutable struct opts_FAM <: ARTOpts @deftype Float
    # Vigilance parameter: [0, 1]
    rho = 0.6; @assert rho >= 0.0 && rho <= 1.0
    # Choice parameter: alpha > 0
    alpha = 1e-7; @assert alpha > 0.0
    # Match tracking parameter
    epsilon = 1e-3; @assert epsilon > 0.0 && epsilon < 1.0
    # Learning parameter: (0, 1]
    beta = 1.0; @assert beta > 0.0 && beta <= 1.0
    # Uncommitted node flag
    uncommitted::Bool = true
    # Display flag
    display::Bool = true
    # Maximum number of epochs during training
    max_epochs::Int = 1
end # opts_FAM()

"""
    FAM <: ARTMAP

Fuzzy ARTMAP struct.
"""
mutable struct FAM <: ARTMAP
    opts::opts_FAM
    config::DataConfig
    W::RealMatrix
    labels::IntegerVector
    n_categories::Int
    epoch::Int
end # FAM <: ARTMAP

"""
    FAM()

Implements a Fuzzy ARTMAP learner.

# Examples
```julia-repl
julia> FAM()
FAM
    opts: opts_FAM
    ...
```
"""
function FAM()
    opts = opts_FAM()
    FAM(opts)
end # FAM()

"""
    FAM(;kwargs...)

Implements a Fuzzy ARTMAP learner with keyword arguments.

# Examples
```julia-repl
julia> FAM(rho=0.7)
FAM
    opts: opts_FAM
    ...
```
"""
function FAM(;kwargs...)
    opts = opts_FAM(;kwargs...)
    FAM(opts)
end # FAM(;kwargs...)

"""
    FAM(opts)

Implements a Fuzzy ARTMAP learner with specified options.

# Examples
```julia-repl
julia> opts = opts_FAM()
julia> FAM(opts)
FAM
    opts: opts_FAM
    ...
```
"""
function FAM(opts::opts_FAM)
    FAM(opts,                       # opts_FAM
        DataConfig(),               # config
        Array{Float}(undef, 0,0),   # W
        Array{Int}(undef, 0),       # labels
        0,                          # n_categories
        0                           # epoch
    )
end # FAM(opts::opts_FAM)


text:
"""
$(SIGNATURES)
"""
average(m::AbstractMatrix) = (sum(m) / length(m))

"""
$(SIGNATURES)
"""
average_absolute(m::AbstractMatrix) = (sum(abs.(m)) / length(m))

"""
$(SIGNATURES)
"""
test_average_absolute(m::AbstractMatrix) = (average_absolute(m)) < 1e6

"""
$(SIGNATURES)
"""
function smallest_positive(x1::Float64, x2::Float64)
    (x1 > 0.0 && x2 < 0.0) && return x1
    (x1 < 0.0 && x2 > 0.0) && return x2
    (x1 > 0.0 && x2 > 0.0) && return min(x1, x2)
    return nothing
end

"""
$(SIGNATURES)

Approximates the gradient of a scalar field.
"""
function ∇(ϕ::AbstractMatrix{Float64})
    height, width = size(ϕ)

    # x, y only represents the first and the second variables. No reference to graphical representations.
    ∇ϕx = zeros(Float64, height, width)   # divergence on the right edge will be filled with zeros
    ∇ϕy = zeros(Float64, height, width)   # divergence on bottom edge will be filled with zeros

    ∇ϕx[begin:end-1, :] = ϕ[begin+1:end, :] - ϕ[begin:end-1, :]
    ∇ϕy[:, begin:end-1] = ϕ[:, begin+1:end] - ϕ[:, begin:end-1]

    fill_borders!(∇ϕx, 0.0)
    fill_borders!(∇ϕy, 0.0)

    return ∇ϕx, ∇ϕy
end

"""
$(SIGNATURES)

Calculate the second-order Laplace operator by convolution of a kernel.
"""
function laplacian(ϕ::AbstractMatrix{Float64})
    height, width = size(ϕ)

    # Embed matrix within a larger matrix for better vectorization and avoid duplicated code
    # Padded matrix adds 1 row/col after the size of ϕ.
    # ϕ is inserted in padded matrix within  1:1+height+1 x 1:1+width+1.
    # The rest of the padded matrix (its borders) are set at 0.0.
    padded_ϕ = zeros(Float64, 1 + height + 1, 1 + width + 1)
    padded_ϕ[begin+1:end-1, begin+1:end-1] .= ϕ[1:end, 1:end]

    # Convolution = up + down + left + right
    ∇²ϕ =
        padded_ϕ[begin+1-1:end-1-1, begin+1:end-1] +
        padded_ϕ[begin+1+1:end-1+1, begin+1:end-1] +
        padded_ϕ[begin+1:end-1, begin+1-1:end-1-1] +
        padded_ϕ[begin+1:end-1, begin+1+1:end-1+1] -
        4.0 * padded_ϕ[begin+1:end-1, begin+1:end-1]

    return ∇²ϕ[1:end-1, 1:end-1]
end


text:
module NodesAndModes

using DelimitedFiles # to read quadrature node data
import LinearAlgebra: diagm, eigen, Symmetric
import SpecialFunctions: gamma
include("meshgrid.jl")

# each type of element shape - used for dispatch only
abstract type AbstractElemShape end
struct Line <: AbstractElemShape end
struct Quad <: AbstractElemShape end
struct Tri <: AbstractElemShape end
struct Hex <: AbstractElemShape end
struct Wedge <: AbstractElemShape end
struct Pyr <: AbstractElemShape end
struct Tet <: AbstractElemShape end

export AbstractElemShape
export Line
export Tri, Quad
export Hex, Wedge, Pyr, Tet

# shared functions
export vandermonde, grad_vandermonde, basis
export nodes, equi_nodes, quad_nodes
include("common_functions.jl") # VDM, grad_VDM, and docs

# interp nodes for Tet(), Tri()
include("warpblend_interp_nodes.jl")

# export some convenient 1D routines by default
include("nodes_and_modes_1D.jl")
export gauss_lobatto_quad, gauss_quad
export jacobiP, grad_jacobiP

##################
# 1D elements
##################
include("line_element.jl")

# ##################
# # 2D elements
# ##################
include("tri_element.jl")
include("quad_element.jl")

##################
# 3D elements
##################
include("hex_element.jl")
include("wedge_element.jl")
include("pyr_element.jl")
include("tet_element.jl")

end # module


text:
# Autogenerated wrapper script for OpenEXR_jll for i686-w64-mingw32-cxx03
export libIex, libIlmThread, libOpenEXR, libOpenEXRCore, libOpenEXRUtil

using Imath_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("OpenEXR")
JLLWrappers.@declare_library_product(libIex, "libIex-3_1.dll")
JLLWrappers.@declare_library_product(libIlmThread, "libIlmThread-3_1.dll")
JLLWrappers.@declare_library_product(libOpenEXR, "libOpenEXR-3_1.dll")
JLLWrappers.@declare_library_product(libOpenEXRCore, "libOpenEXRCore-3_1.dll")
JLLWrappers.@declare_library_product(libOpenEXRUtil, "libOpenEXRUtil-3_1.dll")
function __init__()
    JLLWrappers.@generate_init_header(Imath_jll, Zlib_jll)
    JLLWrappers.@init_library_product(
        libIex,
        "bin\\libIex-3_1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libIlmThread,
        "bin\\libIlmThread-3_1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libOpenEXR,
        "bin\\libOpenEXR-3_1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libOpenEXRCore,
        "bin\\libOpenEXRCore-3_1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libOpenEXRUtil,
        "bin\\libOpenEXRUtil-3_1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#
using DiffMesh, DiffEqFlux, LinearAlgebra, Plots
using GalacticOptim, ParameterSchedulers
#using GLMakie
#----------------------------------------------------------------------------#
using Descartes
ci = translate([0.0,0.0])Circle(1.0)
sq = translate([0.0,0.0])Square([3.0,3.0];center=true)
ss = diff(sq,ci)
sdf = x -> FRep(ss,x)

#sdf = (x,y) -> @. (x^2 + y^2) - 1

#nprt = 2
#shape = x -> @. sin(pi/2*x)
#pou = femPOU2D(nprt,shape,sdf)

nprt = 6
pou  = nnPOU2D(nprt,sdf;depth=4,act=Flux.tanh)

nx = ny = 4
bsp = TensorProdCosines(nprt,nx,ny)
net = MeshNet2D(pou,bsp)

p,re = Flux.destructure(net)
#--------------------------------------#
function utrue(x,y)

    u  = similar(x)
    r2 = @. x*x + y*y

    u .= @. 1 - r2
    u .= @. sin(3*pi*x) + cos(3*pi*y)*exp(x*y)
    u .= @. x*y^3 + sin(r2)  + exp(x^2*y)
    return u
end

# train
ne = 1000
Re = rand(1,ne)
θe = 2*pi*rand(1,ne)

xe = @. Re * cos(θe)
ye = @. Re * sin(θe)

xe = rand(1,1000) .* 2.5
ye = rand(1,1000) .* 2.5

ue = utrue(xe,ye)

# test
np = 50
xp = Array(range(-2f5,2f5,length=np))
yp = Array(range(-2f5,2f5,length=np))

xp = xp * ones(1,np)
yp = ones(np) * yp'
up = utrue(xp,yp)
Ωp = sdf(xp,yp) .< 1e-8

xpp = reshape(xp,np*np)' # flattened
ypp = reshape(yp,np*np)'
upp = reshape(up,np*np)'
Ωpp = reshape(Ωp,np*np)'

#----------------------------------------------------------------------------#
it = 0 # iter
#fig = Figure()
#ax1 = Axis(fig[1,1])
#ax2 = Axis(fig[1,2])
#ax3 = Axis(fig[1,3])
#ax1.title = "Pred"
#ax2.title = "True"
#ax3.title = "Err"
#fig

#--------------------------------------#
function loss(p,xe)
    phi = re(p)

    fe = phi(xe,ye)
    ee = ue .- fe
    loss = Flux.Losses.mse(fe,ue)
    return loss
end
#--------------------------------------#
function makeplot(phi)

    fpp = phi(xpp,ypp)
    fp  = reshape(fpp,np,np)
    ep  = @. abs(up-fp)

#   heatmap!(ax1,Ωp.*fp)
#   heatmap!(ax2,Ωp.*up)
#   heatmap!(ax3,Ωp.*ep)
#   fig

    p1 = heatmap(Ωp.*fp)
    p2 = heatmap(Ωp.*up)
    p3 = heatmap(Ωp.*ep)

    return
end
#--------------------------------------#
function cb1(p,l)
    phi = re(p)

    leastSq!(phi,xe,ye,ue)

    fpp = phi(xpp,ypp)
    epp = (upp .- fpp) .* Ωpp
    er  = norm(epp,Inf)

    global it += 1
    println("it: $it, loss: $l, pointwise er: $er")
    (it % 10) == 0 && makeplot(phi)
    (l < 1e-6) && return true
    return false
end
#--------------------------------------#
function cb2(p,l)
    phi = re(p)

    fpp = phi(xpp,ypp)
    epp = (upp .- fpp) .* Ωpp
    er  = norm(epp,Inf)

    global it += 1
    println("it: $it, loss: $l, pointwise er: $er")
    (it % 10) == 0 && makeplot(phi)
    (l < 1e-6) && return true
    return false
end
#----------------------------------------------------------------------------#
# training

func = OptimizationFunction(loss,GalacticOptim.AutoZygote())
prob = OptimizationProblem(func,p,xe)

sch  = ParameterSchedulers.SinExp(λ0=1e-0,λ1=1e-2,γ=5e-1,period=50)
opt1 = ParameterSchedulers.Scheduler(sch,ADAM())
opt1 = ADAM(1e-2)
opt2 = BFGS()

res  = GalacticOptim.solve(prob,opt1,cb=cb1,maxiters=200)
prob = remake(prob,u0=res.minimizer)
res  = GalacticOptim.solve(prob,opt2,cb=cb1,maxiters=200)
prob = remake(prob,u0=res.minimizer)
res  = GalacticOptim.solve(prob,opt2,cb=cb2,maxiters=100)

mzr = res.minimizer
phi = re(mzr)
#----------------------------------------------------------------------------#
return


text:
export maxcol, myround, clip, clip!, trim, trim!, MINWEIGHT, MINPROB, sortcols

const MINWEIGHT = 1e-10 # trim() deletes prior points with less than this weight.
const MINPROB = 1e-100  # clip() sets to zero entries smaller than this.
const DIGITS = 5 ## digits to round off to when calling unique(), by default

using Statistics

"""
    mean(x::Weighted)
Takes the mean along `dims=ndims(x)` taking `weights(x)` into account. Now returns a Weighted.
"""
function Statistics.mean(x::Weighted)
    arr = sum( x.array .* reshape(weights(x), Iterators.repeated(1,ndims(x)-1)...,:) ; dims=ndims(x))
    Weighted(arr, [1.0], x.opt)
end

"""
    maxcol(x)
    maxcol(x, weights)
Returns the column of `x` with the largest weight, as the same type.
If `x::Weighted{Matrix}` then weights need not be given.
"""
function maxcol(x::Weighted{<:AbsMat}, λ=x.weights)
    (m,i) = findmax(λ)
    Weighted(x.array[:, i:i], [x.weights[i]], unnormalise(x.opt))
end
maxcol(x::AbsMat, λ::AbsVec) = begin (m,i) = findmax(λ); x[:, i:i] end

"""    sort(x::Weighted)
Sorts according to `weights(x)`, with `rev=true` thus biggest first by default. See `sortcols()` to sort by `array(x)` instead. """
function Base.sort(x::Weighted; rev=true)
    ind = sortperm(weights(x), rev=rev)
    ndims(x)==1 && return Weighted(x.array[ind], x.weights[ind], x.opt)
    ndims(x)==2 && return Weighted(x.array[:,ind], x.weights[ind], x.opt)
    ndims(x)==3 && return Weighted(x.array[:,:,ind], x.weights[ind], x.opt)
    error("sort doesn't understand >=4-index Weighted yet, sorry")
end

"""    sortcols(x::WeightedMatrix)
Sorts according to `θ₁`, i.e. the first row of `array(x)`."""
function sortcols(x::WeightedMatrix; kw...)
    perm = sortperm(x.array[1,:]; kw...)
    Weighted( x.array[:,perm], x.weights[perm], x.opt)
end

"""
    round(x::Weighted; digits)
Rounds off θ `x.array`... digits need not be an integer. Always rounds -0.0 to 0.0, too.
    myround(x, digits)
Works without keywords!
"""
Base.round(x::Weighted; digits::Real=DIGITS) =
    Weighted( myround.(x.array, digits), x.weights, x.opt )

myround(x::Weighted, digits::Real=DIGITS) = round(x; digits=digits)

function myround(num::Real, dig::Real=DIGITS) ## this allows me to round(π, 1.5)
    factor = 10^(dig-floor(dig))
    round(num*factor; digits=Int(floor(dig)))/factor |> flipzero
end
myround(num::Real, dig::Int) = round(num; digits=dig) |> flipzero

flipzero(x::Real) = ifelse(x==-0.0, zero(x), x)

using GroupSlices
using EllipsisNotation

"""
    unique(x) = unique!(copy(x))
    unique(f, x)
Removes duplicate points while combining their weights. Decided up to `digits=$DIGITS` digits,
and after applying function `f` if given. Now works for any `ndims(x)`.

Note BTW that both of these return `x` partially sorted, done so that among nearly co-incident points,
the position of the heaviest is kept. But the result is not sure to be completely sorted.
"""
Base.unique(x::Weighted, y::Weighted=x; digits=DIGITS) = Weighted(unique_(x,y;digits=digits)..., x.opt)
Base.unique(x::Weighted, f::Function; kw...) = unique(x, f(x); kw...)
Base.unique(f::Function, x::Weighted; kw...) = unique(x, f(x); kw...)

@inline function unique_(x::Weighted, y::Weighted=x; digits=DIGITS)

    if ndims(x)==2 && size(x.array,1)==0 && size(y.array,1)==0 ## special case: cannot alter order
        return x.array, x.weights
    end

    perm = sortperm(y.weights, rev=true)

    ind = groupslices( myround.(y.array, digits)[..,perm] ,ndims(x))

    k = size(x, ndims(x))
    w = zeros(eltype(x.weights), k)

    if ndims(x)==2 && size(x,1)==0 && size(y,1)>0 ## special case: cannot alter order, but can move weight around
        for i=1:k
            w[perm[ind[i]]] += x.weights[perm[i]]
        end
        return x.array, w
    end

    for i=1:k
        w[ind[i]] += x.weights[perm[i]]
    end

    uind = unique(ind)

    x.array[..,perm[uind]], w[uind] ## .. is EllipsisNotation
end

"""
    trim(x::Weighted)
Removes points with weight less than cutoff `cut=$MINWEIGHT`, either copying or mutating.

    trim(Π::Weighted, L::Weighted)
Removes the same columns from both, using the first one's weights. Returns tuple `(Π,L)`.
"""
trim(x::Weighted; kw...) = Weighted(trim_(x; kw...)..., x.opt)

@inline function trim_(x::Weighted; cut=MINWEIGHT)
    keep = weights(x) .> cut
    if size(x,1)==0 ## special case used for MultiModel... do not move or delete things.
        keep[:] .= true
    end
    x.array[.., keep], x.weights[keep] ## .. is EllipsisNotation
end

function trim(x::Weighted, y::Weighted; cut=MINWEIGHT)
    keep = weights(x) .> cut
    x_array = x.array[..,keep]
    x_weights = x.weights[keep]
    y_array = y.array[..,keep]
    y_weights = y.weights[keep]
    Weighted(x_array, x_weights, x.opt), Weighted(y_array, y_weights, y.opt)
end

"""
    mapslices(f, x::Weighted)
If no dimensions are given, then `f` acts on slices `x.array[:,...:, c]` for `c=1:lastlength(x)`.
`x.weights` are untouched. But you should use `mapcols` really.
"""
function Base.mapslices(f::Function, x::Weighted; dims=collect(1:ndims(x)-1))
    array = mapslices(f, x.array; dims=dims)
    @assert size(array, ndims(array)) == length(x.weights) "mapslices must preserve lastlength(array)"
    Weighted(array, x.weights, addlname(x.opt, "map-") |> unclamp)
end

using SliceMap

slicedoc = """
    MapCols(f, x::Weighted)
    MapCols{d}(f, x::Weighted)
    ThreadMapCols{d}(f, x::Weighted)
Like `mapslices(f,x)` but for SVector column slices, each of length `d`.
`x.weights` are untouched.
"""

@doc slicedoc
SliceMap.MapCols(f::Function, M::WeightedMatrix, args...) =
    MapCols{size(M,1)}(f, M, args...)

SliceMap.MapCols{d}(f::Function, M::WeightedMatrix, args...) where {d} =
    Weighted(MapCols{d}(f, M.array, args...), M.weights, addlname(M.opt, "map-") |> unclamp)

@doc slicedoc
SliceMap.ThreadMapCols{d}(f::Function, M::WeightedMatrix, args...) where {d} =
    Weighted(ThreadMapCols{d}(f, M.array, args...), M.weights, addlname(M.opt, "map-") |> unclamp)

@doc slicedoc
SliceMap.mapcols(f::Function, M::WeightedMatrix, args...) =
    Weighted(mapcols(f, M.array, args...), M.weights, addlname(M.opt, "map-") |> unclamp)

"""    log(Π)
Log all entries, approximately `== log.(Π)` but with nice labels etc.
"""
function Base.log(Π::Weighted)
    if Π.opt.clamp==false || Π.opt.lo <0
        @warn "taking log of Weighted Array which isn't clamped to positive numbers" maxlog=3
    end
    Weighted(log.(Π.array), Π.weights, addlname(Π.opt, "log-") |> unclamp)
end

"""    sqrt(Π)
Sqrt all entries, approximately `== sqrt.(Π)` but with nice labels etc.
"""
function Base.sqrt(Π::Weighted)
    hi = sqrt(Π.opt.hi)
    if Π.opt.clamp==false || Π.opt.lo <0
        @warn "taking sqrt of Weighted Array which isn't clamped to positive numbers"  maxlog=3
        hi = Inf
    end
    Weighted(sqrt.(Π.array), Π.weights, clamp(addlname(Π.opt, "sqrt-"),0,hi))
end

"""    tanh(Π)
Compactifies all entries to [-1,1], approximately `== tanh.(Π)` but with nice labels etc.
"""
Base.tanh(Π::Weighted) = Weighted(tanh.(Π.array), Π.weights, clamp(addlname(Π.opt, "tanh-"),-1,1) )

"""    sigmoid(Π)
Compactifies all entries to [0,1], approximately `== sigmoid.(Π)` but with nice labels etc.
"""
sigmoid(x::Number) = 1 / (1 + exp(-x))
sigmoid(Π::Weighted) = Weighted(sigmoid.(Π.array), Π.weights, clamp(addlname(Π.opt, "sigmoid-"),0,1))


text:
abstract type AbstractRenewableFormulation <: AbstractDeviceFormulation end
abstract type AbstractRenewableDispatchFormulation <: AbstractRenewableFormulation end
struct RenewableFixed <: AbstractRenewableFormulation end
struct RenewableFullDispatch <: AbstractRenewableDispatchFormulation end
struct RenewableConstantPowerFactor <: AbstractRenewableDispatchFormulation end

########################### renewable generation variables #################################
function activepower_variables!(canonical::Canonical,
                               devices::IS.FlattenIteratorWrapper{R}) where R<:PSY.RenewableGen
    add_variable(canonical,
                 devices,
                 Symbol("P_$(R)"),
                 false,
                 :nodal_balance_active;
                 lb_value = x -> 0.0,
                 ub_value = x -> PSY.get_rating(PSY.get_tech(x)))
    return
end

function reactivepower_variables!(canonical::Canonical,
                                 devices::IS.FlattenIteratorWrapper{R}) where R<:PSY.RenewableGen
    add_variable(canonical,
                 devices,
                 Symbol("Q_$(R)"),
                 false,
                 :nodal_balance_reactive)
    return
end

####################################### Reactive Power Constraints #########################
function reactivepower_constraints!(canonical::Canonical,
                                    devices::IS.FlattenIteratorWrapper{R},
                                    device_formulation::Type{RenewableFullDispatch},
                                    system_formulation::Type{<:PM.AbstractPowerModel}) where R<:PSY.RenewableGen
    range_data = Vector{NamedMinMax}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        tech = PSY.get_tech(d)
        name = PSY.get_name(d)
        if isnothing(PSY.get_reactivepowerlimits(tech))
            limits = (min = 0.0, max = 0.0)
            range_data[ix] = (PSY.get_name(d), limits)
            @warn("Reactive Power Limits of $(name) are nothing. Q_$(name) is set to 0.0")
        else
            range_data[ix] = (name, PSY.get_reactivepowerlimits(tech))
        end
    end
    device_range(canonical,
                range_data,
                Symbol("reactiverange_$(R)"),
                Symbol("Q_$(R)"))
    return
end

function reactivepower_constraints!(canonical::Canonical,
                                    devices::IS.FlattenIteratorWrapper{R},
                                    device_formulation::Type{RenewableConstantPowerFactor},
                                    system_formulation::Type{<:PM.AbstractPowerModel}) where R<:PSY.RenewableGen
    names = (PSY.get_name(d) for d in devices)
    time_steps = model_time_steps(canonical)
    p_variable_name = Symbol("P_$(R)")
    q_variable_name = Symbol("Q_$(R)")
    constraint_name = Symbol("reactiverange_$(R)")
    canonical.constraints[constraint_name] = JuMPConstraintArray(undef, names, time_steps)
    for t in time_steps, d in devices
        name = PSY.get_name(d)
        pf = sin(acos(PSY.get_powerfactor(PSY.get_tech(d))))
        canonical.constraints[constraint_name][name, t] = JuMP.@constraint(canonical.JuMPmodel,
                                canonical.variables[q_variable_name][name, t] ==
                                canonical.variables[p_variable_name][name, t] * pf)
    end
    return
end

######################## output constraints without Time Series ############################
function _get_time_series(canonical::Canonical,
                          devices::IS.FlattenIteratorWrapper{<:PSY.RenewableGen})
    initial_time = model_initial_time(canonical)
    use_forecast_data = model_uses_forecasts(canonical)
    parameters = model_has_parameters(canonical)
    time_steps = model_time_steps(canonical)
    device_total = length(devices)
    ts_data_active = Vector{Tuple{String, Int64, Float64, Vector{Float64}}}(undef, device_total)
    ts_data_reactive = Vector{Tuple{String, Int64, Float64, Vector{Float64}}}(undef, device_total)
    for (ix, device) in enumerate(devices)
        bus_number = PSY.get_number(PSY.get_bus(device))
        name = PSY.get_name(device)
        tech = PSY.get_tech(device)
        pf = sin(acos(PSY.get_powerfactor(PSY.get_tech(device))))
        active_power = use_forecast_data ? PSY.get_rating(tech) : PSY.get_activepower(device)
        if use_forecast_data
            forecast = PSY.get_forecast(PSY.Deterministic,
                                        device,
                                        initial_time,
                                        "get_rating",
                                        length(time_steps))
            ts_vector = TS.values(PSY.get_data(forecast))
        else
            ts_vector = ones(time_steps[end])
        end
        ts_data_active[ix] = (name, bus_number, active_power, ts_vector)
        ts_data_reactive[ix] = (name, bus_number, active_power * pf, ts_vector)
    end
    return ts_data_active, ts_data_reactive
end

function activepower_constraints!(canonical::Canonical,
                                devices::IS.FlattenIteratorWrapper{R},
                                device_formulation::Type{<:AbstractRenewableDispatchFormulation},
                                system_formulation::Type{<:PM.AbstractPowerModel}) where R<:PSY.RenewableGen
    parameters = model_has_parameters(canonical)
    use_forecast_data = model_uses_forecasts(canonical)
    if !parameters && !use_forecast_data
        range_data = [(PSY.get_name(d), (min = 0.0, max = PSY.get_rating(PSY.get_tech(d)))) for d in devices]
        device_range(canonical,
                    range_data,
                    Symbol("activerange_$(R)"),
                    Symbol("P_$(R)"))
        return
    end
    ts_data_active, _ = _get_time_series(canonical, devices)
    if parameters
        device_timeseries_param_ub(canonical,
                            ts_data_active,
                            Symbol("activerange_$(R)"),
                            UpdateRef{R}("get_rating"),
                            Symbol("P_$(R)"))
    else
        device_timeseries_ub(canonical,
                            ts_data_active,
                            Symbol("activerange_$(R)"),
                            Symbol("P_$(R)"))
    end
    return
end

########################## Addition of to the nodal balances ###############################
function nodal_expression!(canonical::Canonical,
                           devices::IS.FlattenIteratorWrapper{R},
                           system_formulation::Type{<:PM.AbstractPowerModel}) where R<:PSY.RenewableGen
    parameters = model_has_parameters(canonical)
    ts_data_active, ts_data_reactive = _get_time_series(canonical, devices)
    if parameters
        include_parameters(canonical,
                           ts_data_active,
                           UpdateRef{R}("get_rating"),
                           :nodal_balance_active)
        include_parameters(canonical,
                           ts_data_reactive,
                           UpdateRef{R}("get_rating"),
                           :nodal_balance_reactive)
        return
    end
    for t in model_time_steps(canonical)
        for device_value in ts_data_active
            _add_to_expression!(canonical.expressions[:nodal_balance_active],
                            device_value[2],
                            t,
                            device_value[3]*device_value[4][t])
        end
        for device_value in ts_data_reactive
            _add_to_expression!(canonical.expressions[:nodal_balance_reactive],
                            device_value[2],
                            t,
                            device_value[3]*device_value[4][t])
        end
    end
    return
end

function nodal_expression!(canonical::Canonical,
                           devices::IS.FlattenIteratorWrapper{R},
                           system_formulation::Type{<:PM.AbstractActivePowerModel}) where R<:PSY.RenewableGen
    parameters = model_has_parameters(canonical)
    ts_data_active, ts_data_reactive = _get_time_series(canonical, devices)
    if parameters
        include_parameters(canonical,
                           ts_data_active,
                           UpdateRef{R}("get_rating"),
                           :nodal_balance_active)
        return
    end
    for t in model_time_steps(canonical)
        for device_value in ts_data_active
            _add_to_expression!(canonical.expressions[:nodal_balance_active],
                            device_value[2],
                            t,
                            device_value[3]*device_value[4][t])
        end
    end
    return
end

##################################### renewable generation cost ############################
function cost_function(canonical::Canonical,
                       devices::IS.FlattenIteratorWrapper{PSY.RenewableDispatch},
                       device_formulation::Type{D},
                       system_formulation::Type{<:PM.AbstractPowerModel}) where D<:AbstractRenewableDispatchFormulation
    add_to_cost(canonical,
                devices,
                Symbol("P_RenewableDispatch"),
                :fixed,
                -1.0)
    return
end


text:
__precompile__(true)
module ChaoticMaps

#using Reexport
########################################
# Common structures used by all the maps.
########################################
include("InformationStruct.jl")
export Info

##################################################################
# Individual maps (each map module may contain several submodules)
##################################################################
include("Logistic.jl")
export Logistic

end


text:
"""
    Carries training information for a given tree node
"""
mutable struct TrainNode{T<:AbstractFloat}
    gain::T
    𝑖::Union{Nothing, AbstractVector{UInt32}}
    ∑::Vector{T}
    h::Vector{Vector{T}}
    hL::Vector{Vector{T}}
    hR::Vector{Vector{T}}
    gains::Matrix{T}
end

function TrainNode(nvars, nbins, K, T)
    node = TrainNode{T}(
            zero(T),
            nothing,
            zeros(T, 2*K+1),
            [zeros(T, (2*K+1) * nbins) for j in 1:nvars],
            [zeros(T, (2*K+1) * nbins) for j in 1:nvars],
            [zeros(T, (2*K+1) * nbins) for j in 1:nvars],
            zeros(T, nbins, nvars))
    return node
end

# single tree is made of a vectors of length num nodes
struct Tree{T<:AbstractFloat}
    feat::Vector{Int}
    cond_bin::Vector{UInt8}
    cond_float::Vector{T}
    gain::Vector{T}
    pred::Matrix{T}
    split::Vector{Bool}
end

Tree(x::Vector{T}) where T <: AbstractFloat = Tree(zeros(Int, 1), zeros(UInt8, 1), zeros(T, 1), zeros(T, 1), reshape(x, :, 1), zeros(Bool, 1))
Tree(depth, K, ::T) where {T <: AbstractFloat} = Tree(zeros(Int, 2^depth-1), zeros(UInt8, 2^depth-1), zeros(T, 2^depth-1), zeros(T, 2^depth-1), zeros(T, K, 2^depth-1), zeros(Bool, 2^depth-1))

# eval metric tracking
mutable struct Metric
    iter::Int
    metric::Float32
end
Metric() = Metric(0, Inf)

# gradient-boosted tree is formed by a vector of trees
struct GBTree{T<:AbstractFloat}
    trees::Vector{Tree{T}}
    params::EvoTypes
    metric::Metric
    K::Int
    levels
end


text:
function equals(obj::BigFractionField, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_instance(::Type{BigFractionField})
    return jcall(BigFractionField, "getInstance", BigFractionField, ())
end

function get_one(obj::BigFractionField)
    return jcall(obj, "getOne", BigFraction, ())
end

function get_runtime_class(obj::BigFractionField)
    return jcall(obj, "getRuntimeClass", Class, ())
end

function get_zero(obj::BigFractionField)
    return jcall(obj, "getZero", BigFraction, ())
end

function hash_code(obj::BigFractionField)
    return jcall(obj, "hashCode", jint, ())
end


text:
function [J,dXdx,X]=JacobiDG(DG,F,Grid,Trans)
ksi=DG.xw;
eta=DG.xw;
n=DG.OrdPoly+1;
XE1=Trans(DG,Grid.Nodes(F.N(1)).P,Grid.Nodes(F.N(2)).P);
XE2=Trans(DG,Grid.Nodes(F.N(2)).P,Grid.Nodes(F.N(3)).P);
XE3=Trans(DG,Grid.Nodes(F.N(4)).P,Grid.Nodes(F.N(3)).P);
XE4=Trans(DG,Grid.Nodes(F.N(1)).P,Grid.Nodes(F.N(4)).P);
X=zeros(n,n,2);
for j=1:n
  for i=1:n
    X(i,j,:)=0.5*((1-ksi(i))*XE4(j,:)+(1+ksi(i))*XE2(j,:)+(1-eta(j))*XE1(i,:)+(1+eta(j))*XE3(i,:))-...
      0.25*((1-ksi(i))*((1-eta(j))*Grid.Nodes(F.N(1)).P(1:2)+(1+eta(j))*Grid.Nodes(F.N(4)).P(1:2))+...
      (1+ksi(i))*((1-eta(j))'*Grid.Nodes(F.N(2)).P(1:2)+(1+eta(j))*Grid.Nodes(F.N(3)).P(1:2)));
  end
end
dXdx=zeros(n,n,2,2);
dXdx(:,:,1,1)=DG.DS*X(:,:,1);
dXdx(:,:,2,1)=DG.DS*X(:,:,2);
dXdx(:,:,1,2)=X(:,:,1)*DG.DS';
dXdx(:,:,2,2)=X(:,:,2)*DG.DS';
J=dXdx(:,:,1,1).*dXdx(:,:,2,2)-dXdx(:,:,1,2).*dXdx(:,:,2,1);
end


text:
using DeepDiffs
using TimerOutputs

verify = "verify" in ARGS
num_threads = 4

struct Benchmark
    name::String
    benchmark::Union{Number, String}
    verify::Union{Number, String}
    requires_fasta::Bool
end
Benchmark(name, benchmark, verify) = Benchmark(name, benchmark, verify, false)

const BENCHMARKS = [
    Benchmark("binarytrees", 21, 10),
    Benchmark("fannkuchredux", 12, 7),
    Benchmark("fasta", 25000000, 1000),
    Benchmark("knucleotide", "fasta.txt", "knucleotide/knucleotide-input.txt", true),
    Benchmark("mandelbrot", 16000, 200),
    Benchmark("nbody", 50000000, 1000),
    Benchmark("pidigits", 10000, 27),
    Benchmark("regexredux", "fasta.txt", "regexredux/regexredux-input.txt", true),
    Benchmark("revcomp", "fasta.txt", "revcomp/revcomp-input.txt", true),
    Benchmark("spectralnorm", 5500, 100),
]

function run_benchmarks()
    verify && println("VERIFYING!")
    errored = false
    result_file = "result.bin"
    TimerOutputs.reset_timer!()
    for benchmark in BENCHMARKS
        if !verify && benchmark.requires_fasta
            if !isfile(joinpath(@__DIR__, "fasta.txt"))
                fasta_gen = joinpath(@__DIR__, "fasta", "fasta.jl")
                @info "Generating fasta file"
                run(pipeline(`$(Base.julia_cmd()) $fasta_gen 25000000` ;stdout = "fasta.txt"))
            end
        end
        dir = benchmark.name
        _arg = verify ? benchmark.verify : benchmark.benchmark
        println("Running $dir")
        bdir = joinpath(@__DIR__, dir)
        arg, input = _arg isa String ? ("", "$(_arg)") : (string(_arg), "")
        @timeit dir begin
            for file in readdir(bdir)
                endswith(file, ".jl") || continue
                println("    $file:")
                withenv("JULIA_NUM_THREADS" => num_threads) do
                    if !isempty(input)
                        cmd = pipeline(`$(Base.julia_cmd()) $(joinpath(bdir, file)) `; stdin=input, stdout = result_file)
                    else
                        cmd = pipeline(`$(Base.julia_cmd()) $(joinpath(bdir, file)) $(arg)`; stdout = result_file)
                    end
                    @timeit file run(cmd)
                end
                if verify
                    bench_output = read(result_file, String)
                    correct_output = read(joinpath(bdir, string(dir, "-output.txt")), String)
                    if bench_output != correct_output
                        println(deepdiff(correct_output, bench_output))
                        errored = true
                    end
                end
            end
        end
    end
    TimerOutputs.print_timer(; compact=true, allocations=false)

    if errored
        println()
        error("Some verification failed")
    end
end

run_benchmarks()


text:
module Runtests_withPlot

# Run all tests with the activated plot package

import ModiaResult
using Test

const  test_title = "Tests with plot package " * ModiaResult.currentPlotPackage()

@testset "$test_title" begin
    include("include_all.jl")
end

end

max_stars_repo_name: ModiaSim/ModiaResult

max_issues_repo_name: ModiaSim/ModiaResult

max_forks_repo_name: ModiaSim/ModiaResult


text:
###############################################################################
#
# Abstract original data type
#
###############################################################################
"""
    abstract type AbstractUngriddedData{FT}

Hierachy of AbstractUngriddedData
- [`AbstractMODIS500m`](@ref)
- [`AbstractMODIS1km`](@ref)
"""
abstract type AbstractUngriddedData{FT} end

###############################################################################
#
# MODIS raw data at 500 m resolution
#
###############################################################################
"""
    abstract type AbstractMODIS500m{FT}

Hierachy of AbstractMODIS500m
- [`MOD15A2Hv006LAI`](@ref)
"""
abstract type AbstractMODIS500m{FT} <: AbstractUngriddedData{FT} end

"""
    struct MOD09A1v006LAI

Terra Surface Reflectance
"""
struct MOD09A1v006NIRv{FT} <: AbstractMODIS500m{FT} end

"""
    struct MOD15A2Hv006LAI

Leaf area index
"""
struct MOD15A2Hv006LAI{FT} <: AbstractMODIS500m{FT} end

###############################################################################
#
# MODIS raw data at 1 km resolution
#
###############################################################################
"""
    abstract type AbstractMODIS1km{FT}

Hierachy of AbstractMODIS1km
"""
abstract type AbstractMODIS1km{FT} <: AbstractUngriddedData{FT} end


text:
struct Ray{T<:Real}
    origin::Point{3, T}
    direction::Vec{3, T}
end

origin(ray::Ray) = ray.origin
direction(ray::Ray) = ray.direction
at(ray::Ray{T}, t::T) where {T<:Real} = origin(ray) + t * direction(ray)

import Base.≈
≈(ray1::Ray, ray2::Ray; kwargs...) = ≈(origin(ray1), origin(ray2); kwargs...) && ≈(direction(ray1), direction(ray2); kwargs...)

struct HitRecord{T<:Real}
    point::Point{3, T}
    normal::Vec{3, T}
    t::T
    front_face::Bool
end

point(hitrecord::HitRecord) = hitrecord.point
normal(hitrecord::HitRecord) = hitrecord.normal
t(hitrecord::HitRecord) = hitrecord.t
front_face(hitrecord::HitRecord) = hitrecord.front_face

function face_normal(ray::Ray{T}, outward_normal::Vec{3, T}) where {T<:Real}
    front_face = dot(direction(ray), outward_normal) < 0
    normal = front_face ? outward_normal : -outward_normal
    (front_face, normal)
end

function HitRecord(ray::Ray{T}, outward_normal::Vec{3, T}, t::T) where {T<:Real}
    (front_face, normal) = face_normal(ray, outward_normal)
    HitRecord(at(ray, t), normal, t, front_face)
end


text:
"""
    NeuralNetwork(X, y, w₁, w₂, output)

Creates a container to hold all the necessary components for the network that other package functions refer to.

    Usage:
```julia-repl
NeuralNetwork(X, y, w₁, w₂, output)
```
"""
struct NeuralNetwork
    X::Array
    y::Array
    w₁::Array
    w₂::Array
    output::Array
end


text:
#Find and return the variable names in the expression x. This function is called in the main function mean_formula.
find_variables(x) = find_variables!(Symbol[], x) #this is so that we can call the function without the exclamation
function find_variables!(var_names, x::Number) #if the variable name is a number then just return it without doing anything
    return var_names
end
function find_variables!(var_names, x::Symbol) # if the variable name is a symbol then push it to the list of var_names because its a name of a column
    push!(var_names, x)
end

function find_variables!(var_names, x::Expr) # if the variable is a expression object then we have to crawl through each argument of the expression
    # safety checking
    if x.head == :call  # check for + symbol bc we are summing linear combinations within each expression argument
      # pass the remaining expression
      for argument in x.args[2:end] # since the first argument is the :+ call
        find_variables!(var_names, argument) #recursively find the names given in each argument so check if number, if symbol, if expression etc.again agian
    end
end
return var_names
end

#Search the expression x for the variable name in var, and then create the julia interpretable expression. We extend this function to allow for recursion on the variable vars... This function is called in the main function mean_formula.
function search_variables!(x::Expr, var::Symbol)
    for i in eachindex(x.args)
        if x.args[i] == var # if the argument is one of the variables given then just put it in the right format df[:x1]
            x.args[i] = Meta.parse(string(:input_data_from_user,"[!, ", ":", var, "]"))
        elseif x.args[i] isa Expr # else if the argument is an expression (i.e not a varaible (symbol) or a number) then
            search_variables!(x.args[i], var) #go through this function recursively on each of the arguments of the expression object
        end
    end
    return x
end

function search_variables!(x::Expr, vars...) # this is for when you have more than one variable name found in the string
    for var in vars #goes through each of the variables in the vector vars
        x = search_variables!(x, var) #runs the recursion on each variable in vars
    end
    return x
end

"""
    mean_formula(user_formula_string::String, df::DataFrame)
Construction of the evaluated mean vector, given formula string and named dataframe of covariates as input. This function allows for transformations of any functions supported by julia.
"""
function mean_formula(user_formula_string::String, df::DataFrame)
    global input_data_from_user = df #this is so we can call whatever name the user has for the dataframe
    users_formula_expression = Meta.parse(user_formula_string)
    if(users_formula_expression isa Expr)
        found_markers = TraitSimulation.find_variables(users_formula_expression) #store the vector of symbols of the found variables
        #X = Matrix(df[:, found_markers)
        dotted_args = map(Base.Broadcast.__dot__, users_formula_expression.args) # adds dots to the arguments in the expression
        dotted_expression = Expr(:., dotted_args[1], Expr(:tuple, dotted_args[2:end]...)) #reformats the exprssion arguments by changing the variable names to tuples of the variable names to keep the dot structure of julia
        julia_interpretable_expression = TraitSimulation.search_variables!(dotted_expression, found_markers...) #gives me the julia interpretable exprsesion with the dataframe provided
        mean_vector = eval(Meta.parse(string(julia_interpretable_expression))) #evaluates the julia interpretable expression on the dataframe provided
    else
        mean_vector = [users_formula_expression for i in 1:size(df, 1)]
		found_markers = TraitSimulation.find_variables(users_formula_expression)
    end
    mean_vector, found_markers
end

"""
    FixedEffectTerms(effectsizes::AbstractVecOrMat, snps::AbstractVecOrMat)
Construction of the proper String expression for evaluation in the simulation process,
 using the specified vectors of regression coefficients (Effect Sizes) as a vector of numbers and snp names as a vector of strings.
"""
function FixedEffectTerms(effectsizes::AbstractVecOrMat, snps::AbstractVecOrMat)
 # implementation
    fixed_terms = ""
    for i in 1:length(effectsizes)
        expression = " + " * string(effectsizes[i]) * "(" * snps[i] * ")"
        fixed_terms = fixed_terms * expression
    end
    return String(fixed_terms)
end

"""
append_terms
Allows us to append terms to create a VarianceComponent type
"""
function append_terms!(AB, summand)
	A_esc = esc(summand.args[2])	# elements in args are symbols,
	B_esc = esc(summand.args[3])
	push!(AB.args, :(VarianceComponent($A_esc, $B_esc)))
end

"""
this vc macro allows us to create a vector of VarianceComponent objects for simulation so with_bigfloat_precis, precision::Integer)
so that the user can type out @vc V[1] ⊗ Σ[1] + V[2] ⊗ Σ[2] + .... + V[m] ⊗ Σ[m]
"""
macro vc(expression)
	n = length(expression.args)
	AB = :(VarianceComponent[]) # AB is an empty vector of variance components list of symbols
	if expression.args[1] != :+ #if first argument is not plus (only one vc)
		summand = expression
		append_terms!(AB, summand)
	else #MULTIPLE VARIANCE COMPONENTS if the first argument is a plus (Sigma is a sum multiple variance components)
		for i in 2:n
			summand = expression.args[i]
			append_terms!(AB, summand)
		end
	end
	return(:($AB))
end

"""
vcobjectuple(vcobject)
This function creates a tuple of Variance Components, given a vector of variancecomponent objects to be compatible with VarianceComponentModels.jl
"""
function  vcobjtuple(vcobject::Vector{VarianceComponent})
	m = length(vcobject)
	d = size(vcobject[1].Σ, 1)
	n = size(vcobject[1].V, 1)
	Σ = ntuple(x -> zeros(d, d), m)
	V = ntuple(x -> zeros(n, n), m)
	for i in eachindex(vcobject)
		copyto!(V[i], vcobject[i].V)
		copyto!(Σ[i], vcobject[i].Σ)
	end
	return(Σ, V)
end