text:
function remove_id_from_name(str)

    if startswith(str,"id_")
        return str[4:length(str)]
    end

    if endswith(str,"_id")
        return str[1:length(str)-3]
    end

    return str

end

function build_module_name(str::String
                          ;module_name_for_all_schemas::Union{String,Missing} = missing)
   if ismissing(module_name_for_all_schemas)
      return StringCases.classify(str)
   else
      return module_name_for_all_schemas
   end
end

function build_enum_module_name(str::String)
   StringCases.classify(str)
end

function build_enum_type_name(str::String)
   uppercase(StringCases.snakecase(str))
end

function build_enum_name_w_module(str::String)
   string(build_enum_module_name(str),
         ".",
         build_enum_type_name(str))
end

function build_enum_value(str::String)
   StringCases.snakecase(str)
end

function build_struct_name(str::String)
   StringCases.classify(str)
end

function build_struct_abstract_name(str::String)
   return "I$(build_struct_name(str))"
end

function build_field_name(str::String,
                          lang_code::String
                         ;replace_ids = false,
                          is_onetomany = false)
   return build_field_name(tovector(str),
                           lang_code::String
                          ;replace_ids = replace_ids,
                           is_onetomany = is_onetomany)
end

function build_field_name(str_arr::Vector{String},
                          lang_code::String
                         ;replace_ids = false,
                          is_onetomany = false)

    if (replace_ids)
       str_arr = map(x -> remove_id_from_name(x),
                     str_arr)
    end

    field_name = join(str_arr, '_')
    field_name = StringCases.camelize(field_name)
    if is_onetomany
      field_name = PostgresORMUtil.pluralize(field_name,lang_code)
    end

    return field_name
end

function get_fieldtype_from_coltype(coltype::String,
                                    elttype::String,
                                    customtypes::Dict,
                                    ;tablename::String = "",
                                     colname::String = "")

   attrtype = missing
   customtypes_names = keys(customtypes) |> collect |> n -> string.(n)

   if (coltype == "character"
      || coltype == "character varying"
      || coltype == "text"
      || coltype == "uuid")
      attrtype = "String"
   elseif (coltype == "boolean")
      attrtype = "Bool"
   elseif (coltype == "smallint")
      attrtype = "Int16"
   elseif (coltype == "integer" || coltype == "interval")
      attrtype = "Int32"
   elseif (coltype == "bigint")
      attrtype = "Int64"
   elseif (coltype == "numeric")
      attrtype = "Float64"
   elseif (coltype == "bytea")
      attrtype = "Vector{UInt8}"
   elseif (coltype == "date")
      attrtype = "Date"
   elseif (coltype == "time without time zone")
      attrtype = "Time"
   elseif (coltype == "timestamp without time zone")
      attrtype = "DateTime"
   elseif (coltype == "timestamp with time zone")
      attrtype = "ZonedDateTime"
   elseif (coltype == "ARRAY")
      if (elttype == "_text" || elttype == "_varchar")
        attrtype = "Vector{String}"
      elseif elttype[2:end] in customtypes_names
        elttype = elttype[2:end] # remove the leading underscore
        attrtype = "Vector{$(build_enum_name_w_module(elttype))}"
      else
        error("Unknown array type[$elttype] for table[$tablename] column[$colname]")
      end
   elseif (coltype == "tsvector")
      attrtype = missing
   elseif (coltype == "bytea")
      attrtype = missing
   elseif (coltype == "USER-DEFINED")
      attrtype = build_enum_name_w_module(elttype)
   else
      error("Unknown type[$coltype] for table[$tablename] column[$colname]")
   end

   return attrtype
end

function generate_julia_code(dbconn::LibPQ.Connection,
                             outdir::String
                            ;lang_code = "eng",
                             module_name_for_all_schemas::Union{String,Missing} = "Model")
   object_model =
      generate_object_model(dbconn,
                            lang_code,
                            module_name_for_all_schemas = module_name_for_all_schemas)
   generate_structs_from_object_model(object_model, outdir)
   generate_orms_from_object_model(object_model, outdir)
   generate_enums_from_object_model(object_model, outdir)

end

function generate_object_model(
   dbconn::LibPQ.Connection,
   lang_code::String
  ;ignored_columns::Vector{String} = Vector{String}(),
   camelcase_is_default::Bool = true,
   exceptions_to_default::Vector{String} = Vector{String}(),
   module_name_for_all_schemas::Union{String,Missing} = missing)

   db_analysis = SchemaInfo.analyse_db_schema(dbconn)
   custom_types = SchemaInfo.get_custom_types(dbconn)

   # Initialize result as Dict
   object_model = Dict()

   # Initialize the various components of the result
   modules = []
   structs = []
   fields = []
   enums = []

   # Deal with the custom types
   for (k,v) in custom_types
      enum_name_w_module = build_enum_name_w_module(k)
      enum_values = []
      for val in v[:possible_values]
         enum_value = build_enum_value(val)
         push!(enum_values, enum_value)
      end
      push!(enums, Dict(:module_name => build_enum_module_name(k),
                        :type_name => build_enum_type_name(k),
                        :values => enum_values))
   end # ENDOF for (k,v) in custom_types

   # Deal with the schemas in ordre to fill in modules, structs and fields
   for (schema, schemadef) in db_analysis

      module_name = if ismissing(module_name_for_all_schemas)
                       build_module_name(schema
                                        ;module_name_for_all_schemas = module_name_for_all_schemas)
                    else
                       module_name_for_all_schemas
                    end

      _module = Dict(:name => module_name,
                                  :schema => schema)
      push!(modules, _module)

      for (table,tabledef) in schemadef
         # Ignore table partitions (the partitioned table is the one of interest)
         if tabledef[:is_partition]
            continue
         end
         struct_name = build_struct_name(table)
         struct_abstract_name = build_struct_abstract_name(table)
         _struct = Dict(:name => struct_name,
                        :module => _module,
                        :abstract_name => struct_abstract_name,
                        :table => table,
                        :schema => schema)
         push!(structs, _struct)

         # Initialize some temporary arrays for the different types of fields
         struct_manytoone_fields = []
         struct_id_fields = []
         struct_basic_fields = []

         # Loop over the FKs of the table to add the complex field
         for (fkname,fkdef) in tabledef[:fks]
            manytoone_field = Dict()
            manytoone_field[:struct] = _struct
            manytoone_field[:referenced_table] = fkdef[:referenced_table]
            manytoone_field[:referenced_cols] = fkdef[:referenced_cols]
            manytoone_field[:cols] = fkdef[:referencing_cols]
            manytoone_field[:is_id] = false
            if occursin("onetoone",fkname)
               manytoone_field[:is_onetoone] = true
               manytoone_field[:is_manytoone] = false
            else
               manytoone_field[:is_onetoone] = false
               manytoone_field[:is_manytoone] = true
            end
            manytoone_field[:is_onetomany] = false
            manytoone_field[:is_enum] = false

            # Build a field name by one of the following options:
            # Case 1: The FK is composed of one column only. In this case we use
            #           the name of the column.
            # Case 2.1: The FK is composed of several columns and it is the only
            #             FK in this table pointing to the given target table.
            #             In this case we use the name of the target table.
            # case 2.2 The FK is composed of several columns and it there are
            #             several FKs in this table pointing to the given target
            #             table.  In this case we concatenate the referencing columns
            nb_fks_same_table_same_targeted_table = tabledef[:fks] |>
               n -> filter(x -> x[2][:referenced_table] ==
                                 manytoone_field[:referenced_table],
                           n) |> length
            manytoone_field[:name] = if length(manytoone_field[:cols]) == 1
                  build_field_name(manytoone_field[:cols],
                                   lang_code,
                                  ;replace_ids = true)
               else
                  if nb_fks_same_table_same_targeted_table == 1
                     build_field_name(manytoone_field[:referenced_table][:table],
                                      lang_code)
                  else
                     build_field_name(manytoone_field[:cols],
                                      lang_code
                                     ;replace_ids = true)
                  end
               end

            # Build the referenced module.struct name
            # NOTE: We cannot get the information from 'structs' because the
            #         struct of interest is probably not available yet

            struct_name_w_module = (
               build_module_name(manytoone_field[:referenced_table][:schema]
                 ;module_name_for_all_schemas = module_name_for_all_schemas)
             * "."
             * build_struct_name(manytoone_field[:referenced_table][:table])
            )

            struct_abstract_name_w_module = (
               build_module_name(manytoone_field[:referenced_table][:schema]
                 ;module_name_for_all_schemas = module_name_for_all_schemas)
             * "."
             * build_struct_abstract_name(manytoone_field[:referenced_table][:table])
            )

            manytoone_field[:field_type] = struct_name_w_module
            manytoone_field[:field_abstract_type] = struct_abstract_name_w_module

            push!(struct_manytoone_fields, manytoone_field)

            # Check whether all the referencing columns are in the PK columns, if
            #   yes also add the property to the id property of the struct
            if all(map(x -> x in tabledef[:pk],
                       manytoone_field[:cols]))
                manytoone_field[:is_id] = true
            end

         end # for (fkname,fkdef) in tabledef[:fks]

         # Store the referencing columns for control
         referencing_cols = []
         for v in (collect(values(tabledef[:fks])) |>
                   vect -> map(x -> x[:referencing_cols],vect))
           push!(referencing_cols,v...)
         end

         # Loop over the PKs of the table, if the pk is found in a FKs, skip
         #   because we already added it when looping over the FKs
         for pkcol in tabledef[:pk]

            # Skip if we already handled this column
            if pkcol in referencing_cols continue end

            id_field = Dict()
            id_field[:struct] = _struct
            id_field[:cols] = tovector(pkcol)
            id_field[:is_id] = true
            id_field[:is_manytoone] = false
            id_field[:is_onetoone] = false
            id_field[:is_onetomany] = false
            id_field[:is_enum] = false
            field_name = build_field_name(pkcol,lang_code)
            id_field[:name] = field_name

            field_type =
               get_fieldtype_from_coltype(tabledef[:cols][pkcol][:type],
                                          tabledef[:cols][pkcol][:elttype_if_array],
                                          custom_types)

            # Check if it is an enum
            if tabledef[:cols][pkcol][:type] == "USER-DEFINED"
                id_field[:is_enum] = true
            end

            id_field[:field_type] = field_type
            push!(struct_id_fields, id_field)
         end # for pk in tabledef[:pk]

         # Loop over all the columns of the table and skip the ones that have
         #   already been mapped to complex properties or id properties
         for (colname, coldef) in tabledef[:cols]
            if colname in referencing_cols continue end
            if colname in tabledef[:pk] continue end

            basic_field = Dict()
            basic_field[:struct] = _struct
            basic_field[:cols] = tovector(colname)
            basic_field[:is_id] = false
            basic_field[:is_manytoone] = false
            basic_field[:is_onetoone] = false
            basic_field[:is_onetomany] = false
            basic_field[:is_enum] = false
            field_name = build_field_name(colname,lang_code)
            basic_field[:name] = field_name

            field_type =
               get_fieldtype_from_coltype(coldef[:type],
                                          coldef[:elttype_if_array],
                                          custom_types
                                          ;tablename = table, colname = colname)

            # Check if it is an enum
            if coldef[:type] == "USER-DEFINED"
                basic_field[:is_enum] = true
            end

            basic_field[:field_type] = field_type
            push!(struct_basic_fields, basic_field)

         end

         # Enrich the arrays for all structs with this struct
         push!(fields, struct_manytoone_fields...)
         push!(fields, struct_id_fields...)
         push!(fields, struct_basic_fields...)

      end # ENDOF for (table,tabledef) in schemadef

   end # ENDOF for (schema, schemadef) in db_analysis

   # Now that we have all the struct known we can loop though the  manytoone
   #  fieds and build the corresponding onetomany fields.
   # Eg. For the manytoone (complex prop) 'Rental.staff' we build the onetomany
   #       'Staff.rentals of type Vector{IRental}'
   for manytoone_field in filter(x -> x[:is_manytoone] == true, fields)
      onetomany_field = Dict()
      onetomany_field[:is_id] = false
      onetomany_field[:is_manytoone] = false
      onetomany_field[:is_onetoone] = false
      onetomany_field[:is_onetomany] = true
      onetomany_field[:is_enum] = false
      onetomany_type_name_w_module = manytoone_field[:field_type] # Public.Staff

      # Build a field_name using one of the following options:
      # Case1: The class holding the manytoone field has only one manytoone field
      #          pointing to the target class. In that case the name of the
      #          onetomany field in the target class will simply be the plural
      #          of the name of the manytoone field.
      #          Eg. Rental.staff -> rentals
      # Case2: There are several manytoone fields in the class pointing to the
      #          the same target class. In that case we concatenate the manytoone
      #          field name and the onetoname table name and we make it plural
      #          Eg. Sentence.judge => CivilServant.judgeSentences
      #              Sentence.executionner => CivilServant.executionnerSentences
      #
      nb_manytoone_same_class_and_same_target_class = fields |>
           n -> filter(x -> (x[:struct][:name] == manytoone_field[:struct][:name]),n) |>
           n -> filter(x -> x[:is_manytoone],n) |>
           n -> filter(x -> x[:field_type] == manytoone_field[:field_type],n) |>
           length

      onetomany_field_name = if nb_manytoone_same_class_and_same_target_class == 1
            manytoone_field[:struct][:table]
         else
            string(manytoone_field[:name],"_",manytoone_field[:struct][:table])
         end
      onetomany_field_name = build_field_name(onetomany_field_name,
                                              lang_code
                                             ;is_onetomany = true)
      onetomany_field_type =
        "Vector{$(manytoone_field[:struct][:module][:name]).$(manytoone_field[:struct][:name])}"
      onetomany_field_abstract_type =
         "Vector{$(manytoone_field[:struct][:module][:name]).$(manytoone_field[:struct][:abstract_name])}"

      onetomany_field[:name] = onetomany_field_name
      onetomany_field[:field_type] = onetomany_field_type
      onetomany_field[:field_abstract_type] = onetomany_field_abstract_type
      onetomany_field[:manytoone_field] = manytoone_field

      # Loop over the structs to link it to the
      for _struct in structs
         if "$(_struct[:module][:name]).$(_struct[:name])" == onetomany_type_name_w_module
            onetomany_field[:struct] = _struct
         end
      end
      if !haskey(onetomany_field,:struct)
         error("Unable to find a struct $onetomany_type_name_w_module")
      end

      # Add the fields to the other fields
      # NOTE: There is no risk to creaet an infinite loop because we loop over
      #           the filtered array and not the array itself
      push!(fields,onetomany_field)
   end

   return Dict(:modules => modules,
               :structs => structs,
               :fields => fields,
               :enums => enums)

end

function generate_enums_from_object_model(object_model::Dict, outdir::String)

   outdir = joinpath(outdir,"enum")
   enums = object_model[:enums]

   result = ""
   for e in enums
      result *= "module $(e[:module_name])\n"
      result *= "  export $(e[:type_name])\n"
      result *= "  @enum $(e[:type_name]) begin\n"

      counter = 0
      for v in e[:values]
         counter += 1
         result *= "    $v = $counter \n"
      end

      result *= "  end\n" # close type
      result *= "end\n\n" # close module
   end

   if !isdir(outdir)
      mkpath(outdir)
   end
   # @info file_path
   file_path = joinpath(outdir,"enums.jl")
   write(file_path,result)

end

function generate_structs_from_object_model(object_model::Dict, outdir::String)

   outdir = joinpath(outdir,"structs")

   modules = object_model[:modules]
   structs = object_model[:structs]
   fields = object_model[:fields]

   indent = " "

   # Reset content
   for _struct in structs
      _struct[:struct_content] = ""
   end

   # ################ #
   # Open the struct #
   # ################ #
   for _struct in structs
      str = "mutable struct $(_struct[:name]) <: $(_struct[:abstract_name]) \n\n"
      _struct[:struct_content] *= str
   end

   # ################## #
   # Declare the fields #
   # ################## #
   for f in fields

      if !haskey(f,:struct)
         @error f
         return
      end

      _struct = f[:struct]

      field_name = f[:name]
      field_type = if (f[:is_manytoone] || f[:is_onetoone] || f[:is_onetomany])
         f[:field_abstract_type]
      else
         f[:field_type]
      end
      str = "  $field_name::Union{Missing,$field_type}\n"

      if f[:name] == "policeOfficerRequests"
        @info "$(f[:name]) -> $str"
      end

      _struct[:struct_content] *= str
   end

   # ############################ #
   # Create the first constructor #
   # ############################ #
   for _struct in structs
      str = "\n"
      str *= "  $(_struct[:name])(args::NamedTuple) = $(_struct[:name])(;args...)"
      _struct[:struct_content] *= str
   end

   # ########################### #
   # Open the second constructor #
   # ########################### #
   for _struct in structs
      str = "\n"
      str *= "  $(_struct[:name])(;\n"
      _struct[:struct_content] *= str
   end

   # ####################################################### #
   # Add the arguments declaration of the second constructor #
   # ####################################################### #
   for f in fields

      _struct = f[:struct]

      field_name = f[:name]
      str = "    $field_name = missing,\n"
      _struct[:struct_content] *= str
   end

   # ######################################################### #
   # Close the arguments declaration of the second constructor #
   # ######################################################### #
   for _struct in structs
      str = ""
      str *= "  ) = (\n"
      str *= "    x = new("
      _struct[:struct_content] *= str
   end

   # #################################################################################### #
   # Add the 'new(missing, missing, ...)' arguments assignment of the second constructor #
   # #################################################################################### #
   for f in fields

      _struct = f[:struct]

      str = "missing,"
      _struct[:struct_content] *= str
   end

   # ############################################################ #
   # Close 'new(missing, missing, ...)' of the second constructor #
   # ############################################################ #
   for _struct in structs
      str = ");\n"
      _struct[:struct_content] *= str
   end

   # ###################################################### #
   # Add the arguments assignment of the second constructor #
   # ###################################################### #
   for f in fields

      _struct = f[:struct]

      str = "    x.$(f[:name]) = $(f[:name]);\n"
      _struct[:struct_content] *= str
   end

   # ######################################################## #
   # Close the arguments assignment of the second constructor #
   # ######################################################## #
   for _struct in structs
      str = "    return x\n"
      str *= "  )\n"
      _struct[:struct_content] *= str
   end

   # ################ #
   # Close the struct #
   # ################ #
   for _struct in structs
      str = "\nend "
      _struct[:struct_content] *= str
   end

   # ############## #
   # Write to files #
   # ############## #

   # Empty the modules dirs
   for _module in modules
      module_dir = joinpath(outdir,_module[:name])
      rm(module_dir, recursive=true, force = true)
      mkpath(module_dir)
   end

   # Write abstract types
   for _struct in structs
      str = "abstract type $(_struct[:abstract_name]) <: IEntity end\n"
      module_dir = joinpath(outdir,_struct[:module][:name])
      file_path = joinpath(module_dir,"abstract-types.jl")
      io = open(file_path, "a");
      write(io,str)
      close(io);
   end

   # Write the content of structs to files
   for _struct in structs
      module_dir = joinpath(outdir,_struct[:module][:name])
      if !isdir(module_dir)
         mkpath(module_dir)
      end
      file_path = joinpath(module_dir,"$(_struct[:name]).jl")
      write(file_path,_struct[:struct_content])
   end

end

function generate_orms_from_object_model(object_model::Dict, outdir::String)

   outdir = joinpath(outdir,"orms")

   modules = object_model[:modules]
   structs = object_model[:structs]
   fields = object_model[:fields]

   # Check how may different modules we have, if only one, there is no need
   #   to put the ORMs' modules in separate submodules of the ORM root module
   #   (i.e. we don't need a module ORM.Public, ORM.SchmeName1, ...)
   nb_different_modules = length(unique(x -> x[:name],modules))

   orm_root_module = "ORM"

   indent = " "

   # ####################################### #
   # Data type, ORM, schema name, table name #
   # ####################################### #

   for _struct in object_model[:structs]

      struct_name = "$(_struct[:module][:name]).$(_struct[:name])"
      orm_name = if nb_different_modules > 1
            "$(orm_root_module).$(_struct[:module][:name]).$(_struct[:name])ORM"
         else
            "$(orm_root_module).$(_struct[:name])ORM"
         end

      table = _struct[:table]
      schema = _struct[:schema]

      orm_content = ("
data_type = $struct_name
PostgresORM.get_orm(x::$struct_name) = return($(orm_name))
get_schema_name() = \"$schema\"
get_table_name() = \"$table\"
")

      _struct[:orm_content] = orm_content
      _struct[:orm_name] = orm_name

   # write(filename_for_orm_module,orm_content)

   end #ENDOF `for _struct in structs`

   # ############################# #
   # Columns selection and mapping #
   # ############################# #

   # Declare the function and open the Dict
   for _struct in object_model[:structs]
      str = "\n\n"
      str *= "# Declare the mapping between the properties and the database columns\n"
      str *= "get_columns_selection_and_mapping() = return columns_selection_and_mapping"
      str *= "\nconst columns_selection_and_mapping = Dict(\n"
      _struct[:orm_content] *= str
   end #ENDOF `for _struct in structs`

   mapping_arr = []
   for f in fields

      # Skip the onetomany fields because they do not have any correson column
      if f[:is_onetomany] continue end

      field_name = f[:name]
      # colnames =  "[$(join(f[:cols],", "))]"

      colnames =  if length(f[:cols]) > 1
                     "[" * join( string.("\"", f[:cols], "\""), ", ") * "]"
                  else
                     "\"$(f[:cols][1])\""
                  end
      # colnames = "r4r4W"
      str = "$(repeat(indent,1)) :$(field_name) => $colnames, \n"
      f[:struct][:orm_content] *= str

   end # ENDOF `for id_field in object_model[:id_fields]`

   # Close the Dict for 'columns_selection_and_mapping'
   for _struct in object_model[:structs]
      _struct[:orm_content] *= ")\n\n"
   end

   # ############# #
   # ID properties #
   # ############# #

   # Open the function
   for _struct in object_model[:structs]
      _struct[:orm_content] *= "\n"
      _struct[:orm_content] *= "# Declare which properties are used to uniquely identify an object\n"
      _struct[:orm_content] *= "get_id_props() = return ["
   end

   # Add the fieds
   for f in filter(x -> x[:is_id], fields)
      f[:struct][:orm_content] *= ":$(f[:name]),"
   end

   # Close the function
   for _struct in object_model[:structs]
      _struct[:orm_content] *= "]"
   end

   # ###################### #
   # onetomany_counterparts #
   # ###################### #
   # Declare the function and open the Dict
   for _struct in object_model[:structs]
      _struct[:orm_content] *= "\n\n"
      _struct[:orm_content] *= "# Associate the onetomany properties to the corresponding manytoone peroperties in the other classes \n"
      _struct[:orm_content] *= "get_onetomany_counterparts() = return onetomany_counterparts\n"
      _struct[:orm_content] *= "const onetomany_counterparts = Dict(\n"
   end

   # Add the fieds
   for f in filter(x -> x[:is_onetomany], fields)

      field_name = f[:name]
      field_type = f[:field_type]
      manytoone = f[:manytoone_field]
      manytoone_struct =
         "$(manytoone[:struct][:module][:name]).$(manytoone[:struct][:name])"
      manytoone_field_name = manytoone[:name]

      tip_for_data_type = "# The struct where the associated manytoone property is"
      tip_for_property = "# The name of the associated manytoone property"
      tip_for_action_on_remove =
         "# Change this to 'PostgresORM.CRUDType.delete' if the object doesn't make sense when orphaned"

      str = "
  :$(field_name) => (
    data_type = $manytoone_struct, $tip_for_data_type
    property = :$manytoone_field_name, $tip_for_property
    action_on_remove = PostgresORM.CRUDType.update), $tip_for_action_on_remove \n"

      f[:struct][:orm_content] *= str
   end

   # Close the Dict
   for _struct in object_model[:structs]
      _struct[:orm_content] *= "\n)"
   end

   # ############## #
   # Types override #
   # ############## #
   # Declare the function and open the Dict
   for _struct in object_model[:structs]
      _struct[:orm_content] *= "\n\n"
      _struct[:orm_content] *= "# Override the abstract types \n"
      _struct[:orm_content] *= "get_types_override() = return types_override\n"
      _struct[:orm_content] *= "const types_override = Dict(\n"
   end

   # Add the fieds for the manytoone and onetomany fields
   for f in filter(x -> (x[:is_manytoone] || x[:is_onetoone] || x[:is_onetomany]),
                        fields)
      field_name = f[:name]
      field_type = f[:field_type]
      str = "$(repeat(indent,1)) :$(field_name) => $field_type, \n"
      f[:struct][:orm_content] *= str
   end

   # Close the Dict
   for _struct in object_model[:structs]
      _struct[:orm_content] *= "\n)"
   end

   # ############# #
   # Track changes #
   # ############# #
   # Declare the function and open the Dict
   for _struct in object_model[:structs]
      _struct[:orm_content] *= "\n\n"
      _struct[:orm_content] *= "# Specify whether we want to track the changes to the objects of this class \n"
      _struct[:orm_content] *= "# get_track_changes() = false # Uncomment and modify if needed \n"
      _struct[:orm_content] *= "# get_creator_property() = :a_property_symbol # Uncomment and modify if needed \n"
      _struct[:orm_content] *= "# get_editor_property() = :a_property_symbol # Uncomment and modify if needed \n"
      _struct[:orm_content] *= "# get_creation_time_property() = :a_property_symbol # Uncomment and modify if needed \n"
      _struct[:orm_content] *= "# get_update_time_property() = :a_property_symbol # Uncomment and modify if needed \n"
   end

   # ############## #
   # Write to files #
   # ############## #

   # Empty the modules dirs
   for _module in modules
      module_dir = joinpath(outdir,_module[:name])
      rm(module_dir, recursive=true, force = true)
      mkpath(module_dir)
   end

   # Write abstract types
   for _struct in structs
      module_dir = joinpath(outdir,_struct[:module][:name])
      orm_name_wo_module =
         _struct[:orm_name][
            findlast(".",_struct[:orm_name])[1]+1:length(_struct[:orm_name])
            ]
      file_path = joinpath(module_dir,"$orm_name_wo_module.jl")
      io = open(file_path, "w");
      write(io,_struct[:orm_content])
      close(io);
   end

   return object_model

end

function generate_julia_struct_from_table(
    dbconn::LibPQ.Connection,
    schema_name::String,
    table_name::String,
    struct_name::String,
    filename_for_struct::String,
    filename_for_orm_module::String,
   ;ignored_columns::Vector{String} = Vector{String}(),
    camelcase_is_default::Bool = true,
    exceptions_to_default::Vector{String} = Vector{String}())

    @info "generate_julia_struct_from_table"

    # Retrieve the columns names and types

    query_string = "SELECT column_name,
                           data_type AS column_type,
                           udt_name AS element_type
                    from information_schema.columns
                    WHERE table_schema = \$1 AND table_name = \$2"

    cols = execute_plain_query(query_string,
                               [schema_name,table_name],
                               dbconn)

    # For convenience, we make sure that 'colnames_to_camelcase' and
    #   'colnames_left_as_it_is' are vectors
    if camelcase_is_default
        colnames_to_camelcase = cols[:,:column_name]
        filter!(x -> !(x in exceptions_to_default),colnames_to_camelcase)
    else
        colnames_to_camelcase = exceptions_to_default
    end

    mapping_arr = []
    attrs_declarations = []
    attrsnames = []
    colsnames = []

    # Declare the indentation unit
    indent = "    "

    for c in eachrow(cols)

       colname = c.column_name
       coltype = c.column_type
       elttype = c.element_type
       attrname = colname

       # Skip ignored columns
       if (colname in ignored_columns)
          continue
       end

       if (colname in colnames_to_camelcase)
          attrname = StringCases.camelize(colname)
       end

       attrtype = Union{Missing,String}
       push!(colsnames,colname)
       push!(attrsnames,attrname)
       push!(mapping_arr,":$attrname => \"$colname\"")

       # Deduce the corresponding Julia type
       if (coltype == "character"
          || coltype == "character varying"
          || coltype == "text"
          || coltype == "uuid"
          || coltype == "USER-DEFINED")
          attrtype = "String"
       elseif (coltype == "boolean")
          attrtype = "Bool"
       elseif (coltype == "smallint")
          attrtype = "Int16"
       elseif (coltype == "integer" || coltype == "interval")
          attrtype = "Int32"
       elseif (coltype == "bigint")
          attrtype = "Int64"
       elseif (coltype == "numeric")
          attrtype = "Float64"
       elseif (coltype == "date")
          attrtype = "Date"
       elseif (coltype == "time without time zone")
          attrtype = "Time"
       elseif (coltype == "timestamp without time zone")
          attrtype = "DateTime"
       elseif (coltype == "timestamp with time zone")
          attrtype = "ZonedDateTime"
       elseif (coltype == "ARRAY")
          if (elttype == "_text" || elttype == "_varchar")
             attrtype = "Vector{String}"
          else
             error("Unknown array type[$elttype] for column[$colname]")
          end
       else
          error("Unknown type[$coltype] for column[$colname]")
       end

       # Declare the attribute
       push!(attrs_declarations,
            "$attrname::Union{Missing,$attrtype}")

    end # ENDOF for-loop on columns

    #
    # Prepare the constructors
    #
    constructor1 = "$struct_name(args::NamedTuple) = $struct_name(;args...)"
    constructor2 = "$struct_name(;\n"
    constructor2 *= join(string.(repeat(indent,2),attrsnames," = missing,\n"))
    constructor2 *= " ) = (\n"

    constructor2 *= "$(repeat(indent,3))x = new(" * join(repeat(["missing"],length(attrsnames)),", ") * ");\n"
    constructor2 *=
      join(string.(repeat(indent,3),"x.",attrsnames," = " , attrsnames),"; \n") * ";"

    constructor2 *= "\n$(repeat(indent,3))return x"
    constructor2 *= "\n$(repeat(indent,3)))"

    #
    # Prepare the mapping
    #
    mapping = ""

    #
    # Build the strings
    #
    struct_content = ("
abstract type I$struct_name <: IEntity end \n
mutable struct $struct_name <: I$struct_name

  $(join(attrs_declarations,"\n  "))

  $constructor1
  $constructor2

end
")

   orm_content = ("
data_type = $struct_name
PostgresORM.get_orm(x::$struct_name) = return($(struct_name)ORM)
gettablename() = \"$schema_name.$table_name\"
const columns_selection_and_mapping =
  Dict(
     $(join(mapping_arr,",\n$(repeat(indent,1))"))
  )
")

    write(filename_for_struct,struct_content)
    write(filename_for_orm_module,orm_content)

    cols
end


text:
using QuasiGeostrophy, FFTW, Test, BenchmarkTools
using LinearAlgebra
using Plots, GraphRecipes
using SymbolicUtils
import SymbolicUtils: Chain, Postwalk
include(pwd() * "/test/test_utils.jl")
const tol = 1e1
gr(size = (300,300))
## Define 1D Test
Ω = S¹(0, 2π)
Nx = Ny = Nz = 2^3;
fourier_grid = FourierGrid(Nx, Ω)
fieldnames = ("u", "σ", "v")
u, σ, v = create_fields(names = fieldnames, grid = fourier_grid)
∂x = create_operators(fourier_grid)
# initialize FourierField data
x = fourier_grid.grid[1]
u(sin.(x) .+ 1.0); σ(cos.(x));
v( sin.(x));
starting_energy = (u*u).data[1]
# Wrap Around Impero Objects
@wrapper u=u v=v σ=σ
Δ  = Operator(∂x^2)
Δ²  = Operator(∂x^4)
∂x = Operator(∂x)
κ = 10.0/Nx

# This is just syntax for now
∂t = Operator(nothing, OperatorMetaData(nothing, "∂t"))

@pde_system pde_system = [
    ∂t(u) = -∂x(u^2) + κ*Δ(u) - 1e-6*Δ²(u),
]

pde_plot = plot(pde_system[1]);
abstract type TimeSteppingMethod end
struct  RK1 <: TimeSteppingMethod end #include Δt as a part of RK1
struct  RK2 <: TimeSteppingMethod end

function evolve_pde(pde_system, Δt, ::RK1)
    data = pde_system[1].lhs.operand.data.data
    rhs = evaluate(pde_system[1].rhs)
    data .+= Δt .* rhs
    return nothing
end

function evolve_pde(pde_system, Δt, ::RK2)
    data = pde_system[1].lhs.operand.data.data
    olddata = copy(data)
    stage1 = evaluate(pde_system[1].rhs)
    data .+= Δt .* stage1
    stage2 = evaluate(pde_system[1].rhs)
    data .= olddata .+ Δt/2 .* (stage1 + stage2)
    return nothing
end

const Δt = 0.1 / Nx
T = 3
steps = round(Int, T / Δt)
plotsteps = round(Int, steps / 10)
for i in 1:steps
    evolve_pde(pde_system, Δt, RK2())
    if i%plotsteps==0
        p1 = plot(u.data, ylims = (-2,2))
        display(p1)
        sleep(0.01)
    end
end

ending_energy = evaluate(u * u)[1]
energy_loss = (starting_energy - ending_energy)/starting_energy * 100

ext: jl

max_stars_repo_name: sandreza/QuasiGeostrophy

max_issues_repo_name: sandreza/QuasiGeostrophy

max_forks_repo_name: sandreza/QuasiGeostrophy


text:
__precompile__()
module FinancialToolbox

	using SpecialFunctions: erf
	using Dates,Requires
	include("dates.jl");
	function __init__()
		@require DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74" include("financial_dual.jl")
		@require HyperDualNumbers = "50ceba7f-c3ee-5a84-a6e8-3ad40456ec97" include("financial_hyper.jl")
	end
	include("financial.jl")
	export
	    normcdf,
		normpdf,
		blsprice,
		blkprice,
		blsdelta,
		blsgamma,
		blsvega,
		blsrho,
		blstheta,
		blslambda,
		blsimpv,
		blkimpv,
		##	ADDITIONAL Function
		blspsi,
		blsvanna,
		#Dates
		fromExcelNumberToDate,
		daysact,
		yearfrac

end#End Module


text:
# Random samples from determinantal point processes

"""
Computes a random sample from the determinantal point process defined by the
spectral factorization object `L`.

Inputs:

    `L`: `Eigen` factorization object of an N x N matrix

Output:

    `Y`: A `Vector{Int}` with entries in [1:N].

References:

    Algorithm 18 of \\cite{HKPV05}, as described in Algorithm 1 of \\cite{KT12}.

    @article{HKPV05,
        author = {Hough, J Ben and Krishnapur, Manjunath and Peres, Yuval and Vir\'{a}g, B\'{a}lint},
        doi = {10.1214/154957806000000078},
        journal = {Probability Surveys},
        pages = {206--229},
        title = {Determinantal Processes and Independence},
        volume = {3},
        year = {2005}
        archivePrefix = {arXiv},
        eprint = {0503110},
    }

    @article{KT12,
        author = {Kulesza, Alex and Taskar, Ben},
        doi = {10.1561/2200000044},
        journal = {Foundations and Trends in Machine Learning},
        number = {2-3},
        pages = {123--286},
        title = {Determinantal Point Processes for Machine Learning},
        volume = {5},
        year = {2012},
        archivePrefix = {arXiv},
        eprint = {1207.6083},
    }

    TODO Check loss of orthogonality - a tip from Zelda Mariet
"""
function rand(L::LinearAlgebra.Eigen{S,T}) where {S<:Real,T}
    N = length(L.values)
    J = Int[]
    for n=1:N
        λ = L.values[n]
        rand() < λ/(λ+1) && push!(J, n)
    end

    V = L.vectors[:, J]
    Y = Int[]
    nV = size(V, 2)
    while true
        # Select i from 𝒴=[1:N] (ground set) with probabilities
        # Pr(i) = 1/|V| Σ_{v∈V} (v⋅eᵢ)²

        #Compute selection probabilities
        Pr = zeros(N)
        for i=1:N
            for j=1:nV #TODO this loop is a bottleneck - why?
                Pr[i] += (V[i,j])^2 #ith entry of jth eigenvector
            end
            Pr[i] /= nV
        end
        @assert abs(1-sum(Pr)) < N*eps() #Check normalization

        #Simple discrete sampler
        i, ρ = N, rand()
        for j=1:N
            if ρ < Pr[j]
                i = j
                break
            else
                ρ -= Pr[j]
            end
        end
        push!(Y, i)
        nV == 1 && break #Done
        #V = V⊥ #an orthonormal basis for the subspace of V ⊥ eᵢ
        V[i, :] = 0 #Project out eᵢ
        V = full(qrfact!(V)[:Q])[:, 1:nV-1]
        nV = size(V, 2)
    end
    return Y
end


text:
# Show cancellation errors in stresslet when evaluated directly

push!(LOAD_PATH, string(pwd(),"/src"))
using PyPlot
using ModifiedStokesSolver

for ifprime = [false, true]
    N = 100
    z = linspace(eps(), 3, N)
    E = zeros(N, 6)
    for i=1:N
        Tp = ModifiedStokesSolver.Ti_split_pow(z[i], ifprime)
        Td = ModifiedStokesSolver.Ti_split_direct(z[i], ifprime)
        for j=1:6
            E[i, j] = abs(Tp[j]-Td[j]) / maximum(abs.(Tp))
        end
    end

    figure(ifprime ? 1 : 2)
    clf()
    semilogy(z, E[:,1], label="T1S")
    semilogy(z, E[:,2], label="T1L")
    semilogy(z, E[:,3], label="T2S")
    semilogy(z, E[:,4], label="T2L")
    semilogy(z, E[:,5], label="T3S")
    semilogy(z, E[:,6], label="T3L")
    ylim(1e-16, 1)

    if ifprime
        stem(ModifiedStokesSolver.POWERSERIES_ZMAX_PRIME*[1.0, 1.0], [1e-16, 1], label="cutoff")
    else
        stem(ModifiedStokesSolver.POWERSERIES_ZMAX*[1.0, 1.0], [1e-16, 1], label="cutoff")
    end
    title("ifprime=$ifprime")

    legend()
    grid("on")
end
show()


text:
#We generate random eternity2 feasible puzzles

#First, we create the pieces, randomly simpling the color of the edges. Then, we shuffle the pieces

struct Eternity2Generator <: AbstractModelGenerator
    n::Int
    m::Int
    k::Int
end

"""
    fill_with_generator!(cpmodel::CPModel, gen::GraphColoringGenerator)
Fill a CPModel with the variables and constraints generated. We fill it directly instead of
creating temporary files for efficiency purpose.

This generator create graps for the Eternity2 problem.
"""

"""
function fill_with_generator!(cpmodel::CPModel, gen::Eternity2Generator; seed=123)
    cpmodel.limit.numberOfSolutions = 1
    if !isnothing(seed)
        Random.seed!(seed)
    end
    rng = MersenneTwister(seed)
    n=gen.n
    m=gen.m
    k=gen.k
    colors=1:k

    src_v = Matrix{Int}(undef,n,m+1) #horizontal edges
    src_h = Matrix{Int}(undef,n+1,m) #vertical edges
    fill!(src_h,0)
    fill!(src_v,0)

    shuff = shuffle(rng,1:n*m)

    for i in 2:n
        for j in 1:m
            src_h[i,j]=rand(rng,colors)
        end
    end

    for i in 1:n
        for j in 2:m
            src_v[i,j]=rand(rng,colors)
        end
    end

    pieces = Matrix{Int}(undef, n*m, 4)

    for i = 1:n
        for j = 1:m
            pieces[shuff[(i-1)*m + j],1] = src_h[i,j]
            pieces[shuff[(i-1)*m + j],2] = src_v[i,j+1]
            pieces[shuff[(i-1)*m + j],3] = src_h[i+1,j]
            pieces[shuff[(i-1)*m + j],4] = src_v[i,j]
        end
    end

    table = Matrix{Int}(undef, 6, 4*n*m) #n*m pieces with four different orientations, 6 for orientation  + u + r + d + l + id

    for i = 1:n*m
        for k = 1:4
            table[1,4*(i-1) + k] = i
            table[6,4*(i-1) + k] = 4*(i-1) + k
            for j = 2:5
                table[j,4*(i-1) + k] = pieces[i, (j+k+1)%4+1]
            end
        end
    end

    id = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)
    u = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)#up
    r = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)#right
    d = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)#down
    l = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)#left

    src_v = Matrix{SeaPearl.AbstractIntVar}(undef,n,m+1) #vertical edges
    src_h = Matrix{SeaPearl.AbstractIntVar}(undef,n+1,m) #horizontal edges
    orientation = Matrix{SeaPearl.AbstractIntVar}(undef, n,m)

    for i = 1:n
        src_v[i,m+1] = SeaPearl.IntVar(0, k, "src_v"*string(i)*string(m+1), cpmodel.trailer)
        SeaPearl.addVariable!(cpmodel, src_v[i,m+1]; branchable=false)
    end
    for j =1:m
        src_h[n+1,j] = SeaPearl.IntVar(0, k, "src_h"*string(n+1)*string(j), cpmodel.trailer)
        SeaPearl.addVariable!(cpmodel, src_h[n+1,j]; branchable=false)
    end

    for i = 1:n, j=1:m
        orientation[i,j] = SeaPearl.IntVar(1, 4*n*m, "orientation_"*string(i)*string(j), cpmodel.trailer)
        SeaPearl.addVariable!(cpmodel, orientation[i,j]; branchable=true)
        id[i,j] = SeaPearl.IntVar(1, n*m, "id_"*string(i)*string(j), cpmodel.trailer)
        SeaPearl.addVariable!(cpmodel, id[i,j]; branchable=false)
        src_v[i,j] = SeaPearl.IntVar(0, k, "src_v"*string(i)*string(j), cpmodel.trailer)
        SeaPearl.addVariable!(cpmodel, src_v[i,j]; branchable=false)
        src_h[i,j] = SeaPearl.IntVar(0, k, "src_h"*string(i)*string(j), cpmodel.trailer)
        SeaPearl.addVariable!(cpmodel, src_h[i,j]; branchable=false)
    end

    for i = 1:n, j=1:m
        u[i,j] = src_h[i,j]
        d[i,j] = src_h[i+1,j]
        l[i,j] = src_v[i,j]
        r[i,j] = src_v[i,j+1]

        vars = SeaPearl.AbstractIntVar[id[i,j], u[i,j], r[i,j],d[i,j], l[i,j], orientation[i,j]]
        SeaPearl.addConstraint!(cpmodel, SeaPearl.TableConstraint(vars, table, cpmodel.trailer))

        if (j==m) SeaPearl.addConstraint!(cpmodel, SeaPearl.EqualConstant(r[i,j], 0, cpmodel.trailer)) end
        if (j==1) SeaPearl.addConstraint!(cpmodel, SeaPearl.EqualConstant(l[i,j], 0, cpmodel.trailer)) end
        if (i==1) SeaPearl.addConstraint!(cpmodel, SeaPearl.EqualConstant(u[i,j], 0, cpmodel.trailer)) end
        if (i==n) SeaPearl.addConstraint!(cpmodel, SeaPearl.EqualConstant(d[i,j], 0, cpmodel.trailer)) end
    end

    return nothing
end
"""

function fill_with_generator!(cpmodel::CPModel, gen::Eternity2Generator; seed=nothing)
    cpmodel.limit.numberOfSolutions = 1
    if !isnothing(seed)
        Random.seed!(seed)
    end
    rng = MersenneTwister(seed)
    n=gen.n
    m=gen.m
    k=gen.k
    colors=1:k

    table = Matrix{Int}(undef, 5, 4*n*m)

    src_v = Matrix{Int}(undef,n,m+1) #horizontal edges
    src_h = Matrix{Int}(undef,n+1,m) #vertical edges
    fill!(src_h,0)
    fill!(src_v,0)

    shuff = shuffle(rng,1:n*m)

    for i in 2:n
        for j in 1:m
            src_h[i,j]=rand(rng,colors)
        end
    end

    for i in 1:n
        for j in 2:m
            src_v[i,j]=rand(rng,colors)
        end
    end
    pieces = Matrix{Int}(undef, n*m, 4)

    for i = 1:n
        for j = 1:m
            pieces[shuff[(i-1)*m + j],1] = src_h[i,j]
            pieces[shuff[(i-1)*m + j],2] = src_v[i,j+1]
            pieces[shuff[(i-1)*m + j],3] = src_h[i+1,j]
            pieces[shuff[(i-1)*m + j],4] = src_v[i,j]
        end
    end

    for i = 1:n*m
        for k = 1:4
            table[1,4*(i-1) + k] = i
            for j = 2:5
                table[j,4*(i-1) + k] = pieces[i, (j+k+1)%4+1]
            end
        end
    end

    id = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)
    u = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)#up
    r = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)#right
    d = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)#down
    l = Matrix{SeaPearl.AbstractIntVar}(undef, n, m)#left
    src_v = Matrix{SeaPearl.AbstractIntVar}(undef,n,m+1) #horizontal edges
    src_h = Matrix{SeaPearl.AbstractIntVar}(undef,n+1,m) #vertical edges

    for i = 1:n
        src_v[i,m+1] = SeaPearl.IntVar(0, k, "src_v"*string(i)*string(m+1), cpmodel.trailer)
        SeaPearl.addVariable!(cpmodel, src_v[i,m+1]; branchable=false)
    end
    for j =1:m
        src_h[n+1,j] = SeaPearl.IntVar(0, k, "src_h"*string(n+1)*string(j), cpmodel.trailer)
        SeaPearl.addVariable!(cpmodel, src_h[n+1,j]; branchable=false)
    end

    for i = 1:n, j=1:m
        id[i,j] = SeaPearl.IntVar(1, n*m, "id_"*string(i)*string(j), cpmodel.trailer)
        SeaPearl.addVariable!(cpmodel, id[i,j]; branchable=true)
        src_v[i,j] = SeaPearl.IntVar(0, k, "src_v"*string(i)*string(j), cpmodel.trailer)
        SeaPearl.addVariable!(cpmodel, src_v[i,j]; branchable=false)
        src_h[i,j] = SeaPearl.IntVar(0, k, "src_h"*string(i)*string(j), cpmodel.trailer)
        SeaPearl.addVariable!(cpmodel, src_h[i,j]; branchable=false)
    end

    for i = 1:n, j=1:m
        u[i,j] = src_h[i,j]
        d[i,j] = src_h[i+1,j]
        l[i,j] = src_v[i,j]
        r[i,j] = src_v[i,j+1]

        vars = SeaPearl.AbstractIntVar[id[i,j], u[i,j], r[i,j],d[i,j], l[i,j]]
        SeaPearl.addConstraint!(cpmodel, SeaPearl.TableConstraint(vars, table, cpmodel.trailer))

        if (j==m) SeaPearl.addConstraint!(cpmodel, SeaPearl.EqualConstant(r[i,j], 0, cpmodel.trailer)) end
        if (j==1) SeaPearl.addConstraint!(cpmodel, SeaPearl.EqualConstant(l[i,j], 0, cpmodel.trailer)) end
        if (i==1) SeaPearl.addConstraint!(cpmodel, SeaPearl.EqualConstant(u[i,j], 0, cpmodel.trailer)) end
        if (i==n) SeaPearl.addConstraint!(cpmodel, SeaPearl.EqualConstant(d[i,j], 0, cpmodel.trailer)) end
    end

    #breaking some symmetries

    #if count(==(2),count(==(0),pieces,dims=2))==4
        #index = findfirst(==(2),vec(count(==(0),pieces,dims=2)))
        #SeaPearl.addConstraint!(model,SeaPearl.EqualConstant(id[1,1], index, trailer))
        #SeaPearl.assign!(id[1,1].domain, index)
    #end

    SeaPearl.addConstraint!(cpmodel, SeaPearl.AllDifferent(id, cpmodel.trailer))
    cpmodel.adhocInfo = pieces
    return nothing

end


text:
function render!(fbuffer, hbuffer, cmap, hmap::Matrix{T}, pos::Tuple{T,T,T}, θ::T) where T<:Real
    # draw background (light blue sky)
    fill!(fbuffer, RGB24(0.529f0, 0.808f0, 0.980f0))
    fill!(hbuffer, WIN_HEIGHT)
    # precompute some constants
    w = Base.multiplicativeinverse(MAP_WIDTH%UInt)
    h = Base.multiplicativeinverse(MAP_HEIGHT%UInt)
    # compute orientation angle θ sine and cosine
    sθ, cθ = sincos(θ)
    # draw from front to back (low z coordinate to high z coordinate)
    @inbounds for z in 1 : VIEW_DIST
        # find view-line on map, this calculation corresponds to a field of view of 90°
        px = pos[1] - (cθ+sθ)*z; py = pos[2] - (cθ-sθ)*z
        # segment the line in equal increments
        dx = 2z*cθ/WIN_WIDTH; dy = -2z*sθ/WIN_WIDTH
        # compute local z-scale factor
        scale = ZSCALE/z
        # raster horizontal view-line and draw the vertical for each segment
        for j in 1 : WIN_WIDTH
            # compute map indices from px & py
            pxi = UInt(1) + rem(unsafe_trunc(Int, px)%UInt, w)
            pyi = UInt(1) + rem(unsafe_trunc(Int, py)%UInt, h)
            # compute the on-screen height of the feature
            height = unsafe_trunc(Int, muladd(pos[3]-hmap[pxi,pyi], scale, HORIZON))
            # assign color to column in framebuffer
            color = cmap[pxi,pyi]
            for i in max(height,1) : hbuffer[j]
                fbuffer[j,i] = color
            end
            # keep track of feature height
            hbuffer[j] = min(height, hbuffer[j])
            # move to next view-line feature position
            px += dx; py += dy
        end
    end
end

# "fisheye camera"
#ϕ = range(-π/4,π/4,length=WIN_WIDTH)
# compute position of feature
#if z > VIEW_DIST÷2
#    sθ, cθ = sincos(θ-ϕ[j])
#    px = pos[1] - sθ*z; py = pos[2] - cθ*z
#end

max_stars_repo_name: jebej/VoxelSpace

max_issues_repo_name: jebej/VoxelSpace

max_forks_repo_name: jebej/VoxelSpace


text:
using Merlin
using MLDatasets
using HDF5
using ProgressMeter

const h5file = ".data/wordembeds_nyt100.h5"

function main()
    traindata, testdata, worddict, chardict = setup_data()
    println(length(testdata))
    traindata = traindata[1:10000]
    wordembeds = h5read(h5file, "v")
    charembeds = rand(Float32, 10, 100)
    model = Model2(wordembeds)

    opt = SGD(0.0001, momentum=0.99, nesterov=true)
    #opt = SGD(0.001)
    for epoch = 1:30
        println("epoch: $epoch")
        prog = Progress(length(traindata))
        loss = 0.0
        for (w,cs,y) in shuffle(traindata)
            z = model(w, y)
            loss += sum(z.data)
            vars = gradient!(z)
            for v in vars
                #clipnorm!(v.grad, 0.1)
                opt(v.data, v.grad)
            end
            next!(prog)
        end
        loss /= length(traindata)
        println("loss: $loss")

        ys = cat(1, map(x -> x[3].data[2:end], testdata)...) # remove ROOT head
        zs = cat(1, map(x -> model(x[1])[2:end], testdata)...)
        @assert length(ys) == length(zs)
        acc = mean(i -> ys[i] == zs[i] ? 1.0 : 0.0, 1:length(ys))
        println("test acc.: $acc")
        println("")
    end
end

include("data.jl")
include("model.jl")
main()

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples


text:
using Distributions

# Simulator step! function uses no memory

# optimize! allocates 264.447 MB
# Doubling the time limit roughly doubles memory usage
# Initializing the simulator distributions to actual distributions (as opposed to undef) removes simulation allocation
# After optimizations (primarily moving the internal functions outside) the memory usage for 2s is 256 bytes
# (was about 500 MB)

# Channels allocate the first few times they're used, but then stop allocating
# So let's leave the channels alone for now

# I think the profiler is the main culprit when it comes to memory usage
# Let's check if using locks uses any memory
# Locks use no memory, and seem to be extremely fast
# I could change the profiler-load-balancer system to use locks instead of channels
# The coordinator acquires a lock whenever it has latency data to write
# The latency data is written into a ring-buffer specific for each worker
# The coordinator releases the lock and signals the profiler that there's new data to be processed
# (using another lock, or a semaphore)

# The profiler could be an object instead of a thread
# This object contains a lock and a ring-buffer for each worker
# The coordinator acquires the lock and pushes new samples into the ring-buffers
# The load-balancer acquires the lock and computes the mean and variance for each worker, after which it runs its optimizer
# The load-balancer could also be an object, containing a lock and a partition vector
# The load-balancer writes into this vector whenever it has changes
# At the start of each iteration, the coordinator looks for changes in this vector
# (or it uses some other signaling mechanism to detect changes)

# Channel memory usage seems to be fine
# It allocates whenever it needs to increase its buffer size, but the next time you don't need to re-allocate
# So let's have a look at the profiler memory usage, and let's leave the channels alone for now

# The moving time window allocates
# pushing 100 values allocates 14 KiB

# Let's replace the moving time window with my own
# Yep. Let's replace the moving time window with a ring buffer

function lb_benchmark()
    nworkers = 100
    nwait = round(Int, 100 / 2)
    comp_distributions = [Gamma() for _ in 1:nworkers]
    comm_distributions = [Gamma() for _ in 1:nworkers]
    sim = EventDrivenSimulator(;nwait, nworkers, comm_distributions, comp_distributions)
    @btime step!($sim, 100)
end

function profiler_push(chin)
    for i in 1:100
        timestamp = Time(now())
        v = CodedComputing.ProfilerInput(1, 1/2, 1/10, timestamp, rand(), rand())
        push!(chin, v)
    end
    return
end

function profiler_take(chout)
    while !isready(chout)
        sleep(0.1)
    end
    while isready(chout)
        take!(chout)
    end
    return
end

function profiler_main()
    # setup the profiler
    # put some values into it
    # take the result out
    # meanwhile, measure the memory usage

    chin, chout = CodedComputing.setup_profiler_channels()
    nworkers = 9
    qlower, qupper = 0.0, 1.0
    windowsize = Second(1)
    minsamples = 1
    task = Threads.@spawn CodedComputing.latency_profiler(chin, chout; nworkers, qlower, qupper, windowsize, minsamples)

    println("pushing")
    profiler_push(chin)
    sleep(0.1)

    println("pushing")
    profiler_push(chin)
    sleep(0.1)

    println("pushing")
    profiler_push(chin)
    sleep(0.1)

    println("taking")
    profiler_take(chout)

    println("pushing")
    profiler_push(chin)
    sleep(0.1)

    println("pushing")
    profiler_push(chin)
    sleep(0.1)

    println("pushing")
    profiler_push(chin)
    sleep(0.1)

    println("taking")
    profiler_take(chout)

    println("closing")
    close(chin)
    wait(task)
    close(chout)
    return
end

function cb_foo!(w)
    for i in 1:100
        timestamp = Time(now())
        v = CodedComputing.ProfilerInput(i, 1/2, 1/10, timestamp, 1.0, 1.0)
        pushfirst!(w, v)
    end
    w
end

function cb_clear!(w, windowsize)
    t = Time(now())
    rv = 0
    while length(w) > 0 && (t - w[end].timestamp) > windowsize
        pop!(w)
        rv += 1
    end
    rv
end

function cb_main()
    w = CircularBuffer{CodedComputing.ProfilerInput}(100)
    @time cb_foo!(w)
    println(length(w))

    @time var((x)->getfield(x, :comm_delay), w)
    return

    sleep(2)
    # @time cb_foo!(w)
    # println(length(w))
    windowsize = Second(1)
    @time cb_clear!(w, windowsize)
    println(length(w))
end

function window_foo!(w)
    for i in 1:100
        timestamp = Time(now())
        v = CodedComputing.ProfilerInput(i, 1/2, 1/10, timestamp, 1.0, 1.0)
        fit!(w, (timestamp, v))
    end
    w
end

using Dates, OnlineStats
function window_main()
    windowsize = Second(1)
    w = MovingTimeWindow(windowsize, valtype=CodedComputing.ProfilerInput, timetype=Time)
    sizehint!(w.values, 1000)
    # sizehint!(w, 100)
    @time window_foo!(w)
    sleep(1)
    @time window_foo!(w)
    sleep(1)
    @time window_foo!(w)
    sleep(1)
    @time window_foo!(w)
end

function lock_main()
    l = ReentrantLock()
    @time lock(l)
    @time unlock(l)
    @time lock(l)
    @time unlock(l)
    @time lock(l)
    @time unlock(l)
end

function ch_foo(ch, n)
    for i in 1:n
        v = Foo(i, i+1, float(i))
        push!(ch, v)
    end
    ch
end

function ch_bar(ch)
    while !isempty(ch)
        take!(ch)
    end
    ch
end

struct Foo
    a::Int
    b::Int
    c::Float64
end

function channel_main()
    ch = Channel{Foo}(Inf)
    # ns = [1000, 500, 2000, 1000]
    for _ in 1:10
        n = ceil(Int, rand() * 100)
        @time ch_foo(ch, n)
        @time ch_bar(ch)
        sleep(0.1)
    end
    sleep(10)
    for _ in 1:10
        n = ceil(Int, rand() * 100)
        @time ch_foo(ch, n)
        @time ch_bar(ch)
        sleep(0.1)
    end
    # @time ch_foo(ch, 500)
    # @time ch_bar(ch)
    # @time ch_foo(ch, 2000)
    # @time ch_bar(ch)
    # @time ch_foo(ch, 1000)
    # @time ch_bar(ch)
end

function sim_main()

    # PriorityQueue
    # @btime step!($sim, 100)
    # 302.834 μs (1105 allocations: 54.22 KiB)

    # BinaryMinHeap
    # @btime step!($sim, 100)
    # 208.338 μs (1000 allocations: 31.25 KiB)

    # Integer => Int
    # @btime step!($sim, 100)
    # 110.700 μs (0 allocations: 0 bytes)

    nworkers = 9
    nwait = 6
    nsubpartitions = 160.0
    θs = fill(1/nworkers, nworkers)
    ps = fill(nsubpartitions, nworkers)

    # values recorded on eX3 with low variance and very low comm. latency relative to comp.
    comp_mcs = [255.07814185599995, 235.33940168554452, 253.72451957702958, 116.82766857029705, 116.75048824871284, 116.63288112475249, 126.99890071128709, 126.4382596942574, 117.12455760475244]
    comp_vcs = [0.0027557071631392205, 0.003592251480945663, 0.002281988669487589, 0.0007248118462346298, 0.0005408877863892769, 0.0005133790987450482, 0.0005865578758950107, 0.0006378358097673953, 0.0005586345559702721]
    comm_mcs = [3.373097029702779e-5, 3.189267326732839e-5, 3.427457425742762e-5, 3.50826633663372e-5, 3.664775247524789e-5, 3.760737623762412e-5, 3.5680287128712425e-5, 3.737996039603971e-5, 3.6443900990098125e-5]
    comm_vcs = [1.4755830783262175e-12, 1.5582646160146614e-12, 1.3494372741908877e-12, 2.0233815302431064e-12, 2.7092806020993726e-12, 2.554836195077552e-12, 1.7668490363708934e-12, 1.5867281370452314e-12, 1.8507092575245513e-12]

    comp_distributions = CodedComputing.distribution_from_mean_variance.(Gamma, comp_mcs .* θs ./ ps, comp_vcs .* θs ./ ps)
    comm_distributions = CodedComputing.distribution_from_mean_variance.(Gamma, comm_mcs, comm_vcs)
    sim = EventDrivenSimulator(;nwait=nwait, nworkers, comm_distributions, comp_distributions)

    @time step!(sim, 100)
end

function lb_main()
    nworkers = 2
    nwait = 1
    min_processed_fraction = 0.1
    time_limit = 1.0 # must be floating-point
    θs = [0.3, 0.7]
    qs = 1 ./ [2, 3]
    ps = round.(Int, 1 ./ qs)
    chin, chout = CodedComputing.setup_loadbalancer_channels()

    # put some random values into the load-balancer input
    Random.seed!(123)
    worker = 1
    v1 = CodedComputing.ProfilerOutput(worker, θs[worker], qs[worker], rand(), rand(), rand(), rand())
    push!(chin, v1)
    v1 = CodedComputing.ProfilerOutput(worker, θs[worker], qs[worker], rand(), rand(), rand(), rand())
    push!(chin, v1)
    v1 = CodedComputing.ProfilerOutput(worker, θs[worker], qs[worker], rand(), rand(), rand(), rand())
    push!(chin, v1)

    worker = 2
    v2 = CodedComputing.ProfilerOutput(worker, θs[worker], qs[worker], rand(), rand(), rand(), rand())
    push!(chin, v2)

    # start the load-balancer
    task = Threads.@spawn CodedComputing.load_balancer(chin, chout; min_processed_fraction, nwait, nsubpartitions=ps, nworkers, time_limit, min_improvement=1)

    t0 = time_ns()
    while (time_ns() - t0)/1e9 < 10 && isready(chin)
        sleep(0.1)
    end

    while isready(chout)
        take!(chout)
    end

    worker = 1
    v1 = CodedComputing.ProfilerOutput(worker, θs[worker], qs[worker], rand(), rand(), rand(), rand())
    push!(chin, v1)
    v1 = CodedComputing.ProfilerOutput(worker, θs[worker], qs[worker], rand(), rand(), rand(), rand())
    push!(chin, v1)
    v1 = CodedComputing.ProfilerOutput(worker, θs[worker], qs[worker], rand(), rand(), rand(), rand())
    push!(chin, v1)

    t0 = time_ns()
    while (time_ns() - t0)/1e9 < 10 && isready(chin)
        sleep(0.1)
    end

    while isready(chout)
        take!(chout)
    end

    close(chin)
    wait(task)
    close(chout)
    return
end

function main(ps, ps_prev)
    nworkers = 72
    nwait = 36
    # nworkers = 9
    # nwait = 1
    nsubpartitions = 160.0
    # ps = copy(ps0)
    θs = fill(1/nworkers, nworkers)
    # ps = copy(ps_prev)

    # nslow = 3
    # ps_obs = fill(nsubpartitions, nworkers)
    # comp_ms = fill(1.0, nworkers)
    # comp_ms[1:nslow] .*= 2
    # comp_vs = comp_ms ./ 100
    # comp_mcs = comp_ms ./ (θs ./ ps_obs)
    # comp_vcs = comp_vs ./ (θs ./ ps_obs)

    # comm_mcs = fill(1e-2, nworkers)
    # comm_vcs = comm_mcs ./ 100

    # eX3 values with low variance and very small comm. latency
    # comp_mcs = [255.07814185599995, 235.33940168554452, 253.72451957702958, 116.82766857029705, 116.75048824871284, 116.63288112475249, 126.99890071128709, 126.4382596942574, 117.12455760475244]
    # comp_vcs = [0.0027557071631392205, 0.003592251480945663, 0.002281988669487589, 0.0007248118462346298, 0.0005408877863892769, 0.0005133790987450482, 0.0005865578758950107, 0.0006378358097673953, 0.0005586345559702721]
    # comm_mcs = [3.373097029702779e-5, 3.189267326732839e-5, 3.427457425742762e-5, 3.50826633663372e-5, 3.664775247524789e-5, 3.760737623762412e-5, 3.5680287128712425e-5, 3.737996039603971e-5, 3.6443900990098125e-5]
    # comm_vcs = [1.4755830783262175e-12, 1.5582646160146614e-12, 1.3494372741908877e-12, 2.0233815302431064e-12, 2.7092806020993726e-12, 2.554836195077552e-12, 1.7668490363708934e-12, 1.5867281370452314e-12, 1.8507092575245513e-12]

    # comm. latency equal to comp. latency
    # comm_mcs = copy(comp_mcs .* θs ./ nsubpartitions)
    # comm_vcs = copy(comp_vcs .* θs ./ nsubpartitions)

    # comm. latency equal to half of comp. latency
    # comm_mcs = copy(comp_mcs .* θs ./ nsubpartitions ./ 2)
    # comm_vcs = copy(comp_vcs .* θs ./ nsubpartitions ./ 2)

    comp_mcs = [83.83719400727274, 84.56425873454548, 83.24744768000001, 85.09492904727274, 86.62222161454545, 86.56492212363639, 84.18081268363638, 83.2964639418182, 83.63603345454547, 85.52723904, 85.89885661090909, 78.39971525818183, 93.71899182545458, 85.06894481454546, 85.09361448727273, 86.73630545454546, 85.32156596363637, 84.40723630545455, 84.44267933090909, 82.6533564509091, 84.94617053090911, 84.90497064727273, 84.32782423272727, 84.80634338909091, 84.74789800727272, 90.32809192727272, 84.78158516363638, 79.42350749090909, 86.52342138181817, 84.51160738909093, 84.64746606545455, 84.11714967272728, 84.18197172363638, 84.91479365818182, 83.04630289454546, 86.9430754327273, 84.56391953454548, 82.88443828363637, 86.16268706909094, 84.96690315636366, 85.49810461090908, 84.27044596363636, 84.09113727999998, 84.44407127272729, 85.21930240000002, 84.55733469090909, 85.48306798545454, 85.43646283636365, 84.59088482909092, 89.27814429090911, 85.69318813090909, 82.89300008727272, 84.66070068363636, 84.22216052363638, 84.26186385454547, 84.2088694109091, 87.36550842181819, 84.30621387636366, 83.68611566545457, 84.92868846545456, 83.99821399272729, 84.27077440000001, 91.8733388218182, 84.28995083636363, 85.13140084363637, 86.28002333090912, 85.47884549818181, 85.11117451636365, 83.85590574545455, 84.5384607418182, 83.28374976, 84.23413847272728]
    comp_vcs = [0.0011421443791165577, 0.001965937010903701, 0.001080289134005863, 0.0015147010681869218, 0.0019543853740587075, 0.0018525280119243218, 0.0010717424262526945, 0.0013552023930692842, 0.0013218279306052023, 0.001646853834402909, 0.0011971988786598198, 0.0008171234908152827, 0.0011180002317882714, 0.0014849310589688388, 0.0013503309833066961, 0.002003399670388999, 0.002229045828245602, 0.0016310696849770677, 0.0013699264700789177, 0.001313373340831927, 0.001440510953185439, 0.0015190398374354923, 0.0012164160375029052, 0.0008306335177574833, 0.0011758994786901626, 0.0010003010861745945, 0.0013473881742654802, 0.0007348181667276463, 0.0023673479141235, 0.0012086931203448797, 0.0008219273744405475, 0.0012246601989724356, 0.0011504895403956735, 0.0013126575204895218, 0.0010350639226623624, 0.001785326509537193, 0.0013303620153629623, 0.0009981376354634575, 0.0017964488424443016, 0.001630187729733738, 0.0017519484032818982, 0.001787923696131851, 0.0014579672492872793, 0.001153181321664806, 0.0013350055705060714, 0.0011140041528969352, 0.0011822385440088466, 0.001246437813973071, 0.0009866541269556797, 0.0007728917443144192, 0.0014659113197203564, 0.0013141219658085656, 0.0013415545268761183, 0.001504927399968425, 0.000906065910986086, 0.0010870860165416388, 0.0012550196169302695, 0.0008419391399846764, 0.0010460085177690004, 0.0015293742748779598, 0.0009529932101843303, 0.0014254519016226395, 0.0008692140107794807, 0.0008258379861804196, 0.0017921203526416393, 0.0013239507762105058, 0.001927217880471291, 0.0016084685805719378, 0.0008582012743286067, 0.0014071052024702247, 0.0012597131371129017, 0.0013191295808503373]
    comm_mcs = [0.0003480254545454544, 0.0003319341919191919, 0.00026587455555555564, 0.00034138355555555556, 0.00040316833333333324, 0.00039426245454545454, 0.0003415080000000001, 0.0003219556666666666, 0.0003220089191919191, 0.00038495763636363644, 0.0004099111616161616, 0.0002096090606060607, 0.0002259690404040403, 0.00036263934343434346, 0.00039680437373737374, 0.00037868467676767685, 0.00036656339393939405, 0.0003677447070707071, 0.00035856226262626247, 0.0003042973535353536, 0.0003993650808080807, 0.00041439392929292927, 0.00040552959595959587, 0.00040939215151515157, 0.00041946872727272706, 0.0003035104040404039, 0.00041239234343434355, 0.0002445649797979797, 0.0003739147373737374, 0.00042839323232323243, 0.0004370764949494949, 0.00038970366666666657, 0.0004284080202020202, 0.00043861960606060595, 0.00041204406060606067, 0.00039340965656565655, 0.00044712492929292936, 0.0004082409797979797, 0.00039524609090909083, 0.00043600814141414135, 0.0004214690303030302, 0.00044303437373737376, 0.0004123091818181818, 0.000435883898989899, 0.0004593682323232322, 0.0004801075656565655, 0.00043758155555555554, 0.0004502059090909091, 0.0004367451515151514, 0.00026511119191919186, 0.000407111090909091, 0.00044593664646464643, 0.00044228803030303015, 0.000437117595959596, 0.0004675319797979798, 0.0004789977070707069, 0.00035643654545454537, 0.0004665113030303031, 0.00044163787878787893, 0.00040496073737373753, 0.0004637384646464644, 0.0004160973030303031, 0.00020500094949494935, 0.0004884162121212122, 0.0004058337979797981, 0.0003695217171717172, 0.00039872561616161625, 0.00041293092929292917, 0.0004977244747474749, 0.0004509595252525253, 0.00046671063636363655, 0.00047462048484848474]
    comm_vcs = [5.740965332390451e-8, 4.345452720940765e-8, 2.0897826589721686e-8, 3.9728614195903475e-8, 4.783874768185863e-8, 4.746413291931858e-8, 4.628794348258591e-8, 4.656202443670707e-8, 2.7572550612902546e-8, 4.70896132374031e-8, 4.235644613027693e-8, 1.8075039649862175e-9, 8.830231602806428e-9, 3.322531477545785e-8, 5.105769057712298e-8, 3.608274134141073e-8, 3.91950921956125e-8, 4.6208766182308144e-8, 4.083966421352698e-8, 2.4841435527885144e-8, 4.3822700146215706e-8, 4.882312347214655e-8, 4.320883324624073e-8, 3.7677578237199255e-8, 5.280018543405695e-8, 1.924283053094787e-8, 5.917481436683162e-8, 6.488923716141019e-9, 3.19733472787795e-8, 5.665397354431981e-8, 4.2398035149946876e-8, 3.4169793846282844e-8, 4.8652965888928915e-8, 4.50717038498953e-8, 6.389193203478423e-8, 4.8686624690387e-8, 4.9491465989298095e-8, 5.2208408375938995e-8, 3.4346688861133144e-8, 4.8051166583576e-8, 4.31302555697668e-8, 4.834922571694111e-8, 5.623120069057301e-8, 4.2799384903242254e-8, 5.000090975007734e-8, 5.391810740503352e-8, 4.341524975719643e-8, 5.440975129410279e-8, 4.470642305651242e-8, 1.1782371344862134e-8, 4.785839480238568e-8, 6.184400377152143e-8, 5.0141780705342516e-8, 5.53503587409074e-8, 4.327395670434303e-8, 5.5518941104974865e-8, 3.039298803533887e-8, 4.6083321766938514e-8, 4.0201333080389357e-8, 3.993063539980984e-8, 5.314513615560225e-8, 4.481117759217081e-8, 9.494668103913932e-10, 5.5252223501763073e-8, 4.202337081874707e-8, 3.7601071229576554e-8, 4.208200520001426e-8, 3.801127667739909e-8, 5.5409340272451336e-8, 4.8944886495683677e-8, 5.7064170423605156e-8, 5.956194455121953e-8]

    ms = comp_mcs .* θs ./ ps .+ comp_mcs
    p = sortperm(ms, rev=true)
    comp_mcs .= comp_mcs[p]
    comp_vcs .= comp_vcs[p]
    comm_mcs .= comm_mcs[p]
    comm_vcs .= comm_vcs[p]

    # comm_mcs = copy(comp_mcs)
    # comm_vcs = copy(comp_vcs)

    # return CodedComputing.distribution_from_mean_variance.(Gamma, comm_mcs, comm_vcs)

    # sim_nwait = ceil(Int, nworkers/2)
    # comp_distributions = Vector{Gamma}(undef, nworkers)
    # comm_distributions = Vector{Gamma}(undef, nworkers)
    comp_distributions = [Gamma() for _ in 1:nworkers]
    comm_distributions = [Gamma() for _ in 1:nworkers]
    sim = EventDrivenSimulator(;nwait, nworkers, comm_distributions, comp_distributions)
    min_processed_fraction = nwait / nworkers / nsubpartitions

    # ls = zeros(nworkers)
    # ∇s = zeros(nworkers)
    contribs = zeros(nworkers)

    ps, latency0, contrib0, latency, contrib = CodedComputing.optimize2!(ps, ps_prev, sim; θs, contribs, comp_mcs, comp_vcs, comm_mcs, comm_vcs, min_contribution=min_processed_fraction, time_limit=1.0)

    @info "ps: $ps"
    @info "contribs: $contribs"
    ms = comp_mcs .* θs ./ ps .+ comm_mcs
    @info "ms: $ms"
    @info "latency $(latency / latency0), contrib. $(contrib/ contrib0) ($(contrib / min_processed_fraction))"

    if (!isapprox(latency, latency0, rtol=1e-3) && latency < latency0) && (isapprox(contrib, min_processed_fraction, rtol=1e-2) || min_processed_fraction < contrib)
        # latency improved and contribution is high enough
        @info "improved latency; accepting solution"
        ps_prev .= ps
    elseif (!isapprox(contrib, contrib0, rtol=1e-2) && contrib0 < contrib) && ((isapprox(latency, latency0, rtol=1e-3) || latency < latency0))
        # contribution improved and latency remained unchanged or improved
        @info "improved contrib; accepting solution"
        ps_prev .= ps
    else
        @info "no improvement; rejecting solution"
    end

    # if contrib0 < contrib && latency < latency0
    #     @info "accepting new solution"
    #     ps_prev .= ps
    # else
    #     @info "rejecting new solution"
    # end

    return

    ls = zeros(nworkers)
    simulation_niterations = 100
    simulation_nsamples = 10
    function simulate(ps)

        # setup compute latency distributions
        for i in 1:nworkers
            m = comp_mcs[i] * θs[i] / ps[i]
            v = comp_vcs[i] * θs[i] / ps[i]
            sim.comp_distributions[i] = CodedComputing.distribution_from_mean_variance(Gamma, m, v)
        end

        # run nsamples simulations, each consisting of niterations steps
        ls .= 0
        for _ in 1:simulation_nsamples
            sim = EventDrivenSimulator(sim)
            step!(sim, simulation_niterations)
            ls .+= sim.nfresh ./ simulation_niterations
        end
        ls ./= simulation_nsamples
    end

    # setup communication latency distributions
    for i in 1:nworkers
        m = comm_mcs[i]
        v = comm_vcs[i]
        sim.comm_distributions[i] = CodedComputing.distribution_from_mean_variance(Gamma, m, v)
    end

    function finite_diff(ps, i, δ=ceil(Int, ps[i]/10))
        pi0 = ps[i]
        ps[i] = pi0 + δ
        ls = simulate(ps)

        # while isapprox(ls[i], 0)
        #     δ *= 2
        #     ps[i] = pi0 + δ
        #     ls = simulate(1 ./ ps)
        # end

        forward = ls[i] * θs[i] / ps[i]
        ps[i] = pi0 - δ
        ps[i] = max(1, ps[i])
        ls = simulate(ps)
        backward = ls[i] * θs[i] / ps[i]
        rv = (forward - backward) / (pi0 + δ - ps[i])
        ps[i] = pi0
        rv
    end

    function finite_diff1(ps, i, δ=ceil(Int, ps[i]/10))
        ps[1] += δ
        ls = simulate(ps)
        forward = ls[i] * θs[i] / ps[i]
        ps[1] -= 2δ
        ls = simulate(ps)
        backward = ls[i] * θs[i] / ps[i]
        ps[1] += δ
        println("forward: $forward, backward: $backward")
        (forward - backward) / 2δ
    end

    # println(finite_diff(ps, 1))
    # println(finite_diff(ps, 1))
    # println(finite_diff(ps, 1))

    min_processed_fraction = sim_nwait / nworkers / nsubpartitions

    ls = simulate(ps)
    losses = ls .* θs ./ ps
    f0 = var(losses)
    println("prev. loss: $f0, contribution: $(sum(losses))")
    println(sort(losses))

    Is = collect(1:nworkers)
    for _ in 1:10
        shuffle!(Is)
        # println("start")

        ls = simulate(ps)
        losses = ls .* θs ./ ps
        μ = mean(losses)

        for i in Is

            δ = finite_diff(ps, i)

            if isapprox(δ, 0)
                x = ps[i]*1.1
            else
                # solve for x in losses[i] + (x - ps[i]) * δ - μ = 0
                x = (μ - losses[i]) / δ + ps[i]
                x = min(x, ps[i]*1.1)
                x = max(x, ps[i]*0.9)
                x = max(x, 1.0)
            end
            ps[i] = x
            # println("x: $x")

            # solve for x in losses[i] + δ*x - μ = 0
            # x = (μ - losses[i]) / δ
            # ps[i] +=

        end

        # println("finish")

        ls = simulate(ps)
        losses = ls .* θs ./ ps
        ps ./= min_processed_fraction / sum(losses)
    end

    ps .= round.(ps)

    ls = simulate(ps)
    losses = ls .* θs ./ ps
    f1 = var(losses)
    # println("$(sum(losses)) / $min_processed_fraction")
    println("new loss: $f1, contribution: $(sum(losses))")
    println(sort(losses))

    println("f0: $f0, f1: $f1")
    # println(ps)
    if f1 < f0 * 0.99
        println("found a sufficiently better solution")
        return ps
    else
        println("new solution not sufficiently better; discarding it")
        return ps0
    end
end

"""

Return a lower bound on the log of the probability of `z1 + y1 <= z2 + y2`, where `z1`, `y1`,
`z2`, and `y2` are drawn from the probability distributions `dz1`, `dy1`, `dz2`, and `dy2`,
respectively.
"""
function less_than_lower_bound(dz1, dy1, dz2, dy2)

    # compute the midpoint of the means
    cz = (mean(dz1) + mean(dz2)) / 2
    cy = (mean(dy1) + mean(dy2)) / 2

    # prob. of z1 <= cz and z2 >= cz
    v1 = logcdf(dz1, cz) + logccdf(dz2, cz)

    # prob. of y1 <= cy and y2 >= cy
    v2 = logcdf(dy1, cy) + logccdf(dy2, cy)

    v1 + v2
end

"""

For the `i`-th random variable, return a lower bound on the log of the probability of a sample drawn
from this variable being smaller than all other variables.
"""
function less_than_lower_bound(i::Integer, dzs, dys)

    # compute the midpoint of the means
    cz = mean(mean, dzs)
    cy = mean(mean, dys)

    # prob. of zi <= cz and all others >= cz
    pz = logcdf(dzs[i], cz)
    for j in 1:length(dzs)
        if j != i
            pz += logccdf(dzs[j], cz)
        end
    end

    # prob. of yi <= cy and all others >= cy
    py = logcdf(dys[i], cy)
    for j in 1:length(dys)
        if j != i
            py += logccdf(dys[j], cy)
        end
    end

    # prob. of both being true
    pz + py
end

"""

For each random variable, return a lower bound on the log of the probability of a sample drawn
from this variable being smaller than all other variables.
"""
function less_than_lower_bounds!(rv, dzs, dys)
    length(rv) == length(dzs) || throw(DimensionMismatch("rv has dimension $(length(rv)), but dzs has dimension $(length(dzs))"))
    length(rv) == length(dys) || throw(DimensionMismatch("rv has dimension $(length(rv)), but dys has dimension $(length(dys))"))
    n = length(rv)
    rv .= 0

    # compute the midpoint of the means
    cz = mean(mean, dzs)
    cy = mean(mean, dys)

    # for each i, prob. of zi <= cz and all others >= cz
    pz = 0.0
    for i in 1:n
        pz += logccdf(dzs[i], cz)
    end
    for i in 1:n
        v = logccdf(dzs[i], cz)
        pz -= v
        rv[i] += pz + logcdf(dzs[i], cz)
        pz += v
    end

    # for each i, prob. of yi <= cy and all others >= cy
    py = 0.0
    for i in 1:n
        py += logccdf(dys[i], cy)
    end
    for i in 1:n
        v = logccdf(dys[i], cy)
        py -= v
        rv[i] += py + logcdf(dys[i], cy)
        py += v
    end

    rv
end

less_than_lower_bounds(dzs, dys) = less_than_lower_bounds!(zeros(length(dzs)), dzs, dys)

# The lower bound looks great
# Next, let's change the load-balancer to evaluate contribution in the log domain

function bound_main()

    nworkers = 9
    nsubpartitions = 160.0
    θs = fill(1/nworkers, nworkers)
    ps = fill(nsubpartitions, nworkers)

    # values recorded on eX3 with low variance and very low comm. latency relative to comp.
    comp_mcs = [255.07814185599995, 235.33940168554452, 253.72451957702958, 116.82766857029705, 116.75048824871284, 116.63288112475249, 126.99890071128709, 126.4382596942574, 117.12455760475244]
    comp_vcs = [0.0027557071631392205, 0.003592251480945663, 0.002281988669487589, 0.0007248118462346298, 0.0005408877863892769, 0.0005133790987450482, 0.0005865578758950107, 0.0006378358097673953, 0.0005586345559702721]
    comm_mcs = [3.373097029702779e-5, 3.189267326732839e-5, 3.427457425742762e-5, 3.50826633663372e-5, 3.664775247524789e-5, 3.760737623762412e-5, 3.5680287128712425e-5, 3.737996039603971e-5, 3.6443900990098125e-5]
    comm_vcs = [1.4755830783262175e-12, 1.5582646160146614e-12, 1.3494372741908877e-12, 2.0233815302431064e-12, 2.7092806020993726e-12, 2.554836195077552e-12, 1.7668490363708934e-12, 1.5867281370452314e-12, 1.8507092575245513e-12]

    dzs = CodedComputing.distribution_from_mean_variance.(Gamma, comp_mcs .* θs ./ ps, comp_vcs .* θs ./ ps)
    dys = CodedComputing.distribution_from_mean_variance.(Gamma, comm_mcs, comm_vcs)

    println([less_than_lower_bound(i, dzs, dys) for i in 1:nworkers])
    println(less_than_lower_bounds(dzs, dys))
    return

    dz1 = Gamma(1)
    dz2 = Gamma(1)
    dy1 = Gamma(2)
    dy2 = Gamma(1)

    # bound prob. of 1 being faster than 2

    # take the midpoint of the means
    cz = (mean(dz1) + mean(dz2)) / 2
    cy = (mean(dy1) + mean(dy2)) / 2

    @info "cz: $cz, cy: $cy"

    # prob. of z1 <= cz and z2 >= cz
    v1 = cdf(dz1, cz) * (1 - cdf(dz2, cz))

    # prob. of y1 <= cy and y2 >= cy
    v2 = cdf(dy1, cy) * (1 - cdf(dy2, cy))

    dzs = [dz1, dz2]
    dys = [dy1, dy2]
    v1 * v2, exp(less_than_lower_bound(dz1, dy1, dz2, dy2)), exp(less_than_lower_bound(1, dzs, dys))
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

plain(x) = sprint(plain, x)

function plain(io::IO, content::Vector)
    isempty(content) && return
    for md in content[1:end-1]
        plain(io, md)
        println(io)
    end
    plain(io, content[end])
end

plain(io::IO, md::MD) = plain(io, md.content)

function plain{l}(io::IO, header::Header{l})
    print(io, "#"^l*" ")
    plaininline(io, header.text)
    println(io)
end

function plain(io::IO, code::Code)
    println(io, "```", code.language)
    println(io, code.code)
    println(io, "```")
end

function plain(io::IO, p::Paragraph)
    plaininline(io, p.content)
    println(io)
end

function plain(io::IO, list::List)
    for (i, item) in enumerate(list.items)
        print(io, list.ordered ? "$i. " : "  * ")
        plaininline(io, item)
        println(io)
    end
end

function plain(io::IO, md::HorizontalRule)
    println(io, "–" ^ 3)
end

# Inline elements

plaininline(x) = sprint(plaininline, x)

function plaininline(io::IO, md...)
    for el in md
        plaininline(io, el)
    end
end

plaininline(io::IO, md::Vector) = !isempty(md) && plaininline(io, md...)

plaininline(io::IO, md::Image) = plaininline(io, "![", md.alt, "](", md.url, ")")

plaininline(io::IO, s::String) = print(io, s)

plaininline(io::IO, md::Bold) = plaininline(io, "**", md.text, "**")

plaininline(io::IO, md::Italic) = plaininline(io, "*", md.text, "*")

plaininline(io::IO, md::Code) = print(io, "`", md.code, "`")

plaininline(io::IO, br::LineBreak) = println(io)

plaininline(io::IO, x) = writemime(io, MIME"text/plain"(), x)

# writemime

Base.writemime(io::IO, ::MIME"text/plain", md::MD) = plain(io, md)


text:
using Revise
using SpineOpt
#using Gurobi
using JuMP
using Cbc

optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
m = run_spineopt(ARGS...; mip_solver=optimizer, lp_solver=optimizer)

#optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
#optimizer = optimizer_with_attributes(Gurobi.Optimizer, "OutputFlag" => 0, "MIPGap" => 1e-2)
#m = run_spineopt(ARGS...)


text:
isbanded(A::AbstractTriangular) = isbanded(parent(A))
bandwidths(A::Union{UpperTriangular,UnitUpperTriangular}) =
    (min(0,bandwidth(parent(A),1)), bandwidth(parent(A),2))
bandwidths(A::Union{LowerTriangular,UnitLowerTriangular}) =
    (bandwidth(parent(A),1), min(0,bandwidth(parent(A),2)))

triangularlayout(::Type{Tri}, ::ML) where {Tri,ML<:BandedColumns} = Tri{ML}()
triangularlayout(::Type{Tri}, ::ML) where {Tri,ML<:BandedRows} = Tri{ML}()
triangularlayout(::Type{Tri}, ::ML) where {Tri,ML<:ConjLayout{<:BandedRows}} = Tri{ML}()

function tribandeddata(::TriangularLayout{'U'}, A)
    B = triangulardata(A)
    u = bandwidth(B,2)
    D = bandeddata(B)
    view(D, 1:u+1, :)
end

function tribandeddata(::TriangularLayout{'L'}, A)
    B = triangulardata(A)
    l,u = bandwidths(B)
    D = bandeddata(B)
    view(D, u+1:l+u+1, :)
end

tribandeddata(A) = tribandeddata(MemoryLayout(typeof(A)), A)

Base.replace_in_print_matrix(A::Union{UpperTriangular{<:Any,<:AbstractBandedMatrix},
                                      UnitUpperTriangular{<:Any,<:AbstractBandedMatrix}}, i::Integer, j::Integer, s::AbstractString) =
    -bandwidth(A,1) ≤ j-i ≤ bandwidth(A,2) ? s : Base.replace_with_centered_mark(s)

Base.replace_in_print_matrix(A::Union{LowerTriangular{<:Any,<:AbstractBandedMatrix},
                                      UnitLowerTriangular{<:Any,<:AbstractBandedMatrix}}, i::Integer, j::Integer, s::AbstractString) =
    -bandwidth(A,1) ≤ j-i ≤ bandwidth(A,2) ? s : Base.replace_with_centered_mark(s)

# Mul
@lazylmul UpperTriangular{T, <:AbstractBandedMatrix{T}} where T
@lazylmul UnitUpperTriangular{T, <:AbstractBandedMatrix{T}} where T
@lazylmul LowerTriangular{T, <:AbstractBandedMatrix{T}} where T
@lazylmul UnitLowerTriangular{T, <:AbstractBandedMatrix{T}} where T

@inline function materialize!(M::BlasMatLmulVec{<:TriangularLayout{'U',UNIT,<:BandedColumnMajor},
                                                <:AbstractStridedLayout}) where UNIT
    A,x = M.A,M.B
    tbmv!('U', 'N', UNIT, size(A,1), bandwidth(A,2), tribandeddata(A), x)
end

@inline function materialize!(M::BlasMatLmulVec{<:TriangularLayout{'L',UNIT,<:BandedColumnMajor},
                                                <:AbstractStridedLayout}) where UNIT
    A,x = M.A,M.B
    tbmv!('L', 'N', UNIT, size(A,1), bandwidth(A,1), tribandeddata(A), x)
end

@inline function materialize!(M::BlasMatLmulVec{<:TriangularLayout{UPLO,UNIT,BandedRowMajor},
                                                <:AbstractStridedLayout}) where {UPLO,UNIT}
    A,x = M.A,M.B
    tbmv!(UPLO, 'T', UNIT, transpose(tribandeddata(A)), x)
end

@inline function materialize!(M::BlasMatLmulVec{<:TriangularLayout{UPLO,UNIT,ConjLayout{BandedRowMajor}},
                                                <:AbstractStridedLayout}) where {UPLO,UNIT}
    A,x = M.A,M.B
    tbmv!(UPLO, 'C', UNIT, tribandeddata(A)', dest)
end

# Ldiv
@lazyldiv UpperTriangular{T, <:AbstractBandedMatrix{T}} where T
@lazyldiv UnitUpperTriangular{T, <:AbstractBandedMatrix{T}} where T
@lazyldiv LowerTriangular{T, <:AbstractBandedMatrix{T}} where T
@lazyldiv UnitLowerTriangular{T, <:AbstractBandedMatrix{T}} where T

for UNIT in ('N', 'U')
    @eval begin
        @inline function materialize!(M::BlasMatLdivVec{<:TriangularLayout{'U',$UNIT,<:BandedColumnMajor},
                                        <:AbstractStridedLayout})
            A,x = M.A,M.B
            tbsv!('U', 'N', $UNIT, size(A,1), bandwidth(A,2), tribandeddata(A), x)
        end

        @inline function materialize!(M::BlasMatLdivVec{<:TriangularLayout{'L',$UNIT,<:BandedColumnMajor},
                                                        <:AbstractStridedLayout})
            A,x = M.A,M.B
            tbsv!('L', 'N', $UNIT, size(A,1), bandwidth(A,1), tribandeddata(A), x)
        end
    end
    for UPLO in ('U', 'L')
        @eval begin
            @inline function materialize!(M::BlasMatLdivVec{<:TriangularLayout{$UPLO,$UNIT,BandedRowMajor},
                                                        <:AbstractStridedLayout})
                A,x = M.A,M.B
                tbsv!($UPLO, 'T', $UNIT, transpose(tribandeddata(A)), x)
            end

            @inline function materialize!(M::BlasMatLdivVec{<:TriangularLayout{$UPLO,$UNIT,ConjLayout{BandedRowMajor}},
                                                        <:AbstractStridedLayout})
                A,x = M.A,M.B
                tbsv!($UPLO, 'C', $UNIT, tribandeddata(A)', x)
            end
        end
    end
end


text:
import NbodyGradient: kepler_init
# Runs a simple test of kepler_init.jl

# Not sure why this test can't see these from the module...
const GNEWT = 39.4845 / 365.242^2
const third = 1.0 / 3.0

@testset "kepler_init" begin

ntrial = 10
for itrial = 1:ntrial
    t0 = 2.4
    mass = 1.0
    period = 1.5
    elements = ones(Float64, 6)
    while elements[3]^2 + elements[4]^2 >= 0.2^2
        elements = [period,rand() * period,randn(),randn(),rand() * pi,rand() * pi]
    end
    elements_diff = zeros(Float64, 6)
    ecc = sqrt(elements[3]^2 + elements[4]^2)
    ntime = 1000
    time = linearspace(t0, t0 + period, ntime)
    xvec = zeros(Float64, 3, ntime)
    vvec = zeros(Float64, 3, ntime)
    vfvec = zeros(Float64, 3, ntime)
    timebig = big.(time)
    # dt = 1e-8
    dt = period / ntime
    dtbig = big(dt)
    jac_init = zeros(Float64, 7, 7)
    # Compute Jacobian in BigFloat precision with analytic formulae:
    jac_init_big = zeros(BigFloat, 7, 7)
    # Compute Jacobian in BigFloat precision with finite differences:
    jac_init_num = zeros(BigFloat, 7, 7)
    # Variations in the elements for finite differences:
    delements = big.([1e-15,1e-15,1e-15,1e-15,1e-15,1e-15,1e-15])
    elements_name = ["P","t0","ecosom","esinom","inc","Omega","Mass"]
    cartesian_name = ["x","y","z","vx","vy","vz","mass"]
    # println("Elements: ",elements)
    massbig = big(mass)
    for i = 1:ntime
        x, v = kepler_init(timebig[i], massbig, big.(elements))
        if i == 1
            x_ekep, v_ekep = kepler_init(time[i], mass, elements, jac_init)
            # Redo in BigFloat
            elements_big = big.(elements)
            x_ekep_big, v_ekep_big = kepler_init(timebig[i], massbig, elements_big, jac_init_big)
            # Check that these agree:
            #    println("x-x_ekep: ",maximum(abs.(x-x_ekep))," v-v_ekep: ",maximum(abs.(v-v_ekep)))
            # Now take some derivatives:
            for j = 1:7
                elements_diff = big.(elements)
                if j <= 6
                    elements_diff[j] -= delements[j]
                else
                    massbig -= delements[7]
                end
                x_minus, v_minus = kepler_init(timebig[i], massbig, elements_diff)
                elements_diff .= elements
                if j <= 6
                    elements_diff[j] += delements[j]
                else
                    massbig += 2.0 * delements[7]
                end
                x_plus, v_plus = kepler_init(timebig[i], massbig, elements_diff)
                if j == 7
                    massbig -= delements[7]
                end
                for k = 1:3
                    jac_init_num[  k,j] = .5 * (x_plus[k] - x_minus[k]) / delements[j]
                    jac_init_num[3 + k,j] = .5 * (v_plus[k] - v_minus[k]) / delements[j]
                end
            end
            jac_init_num[7,7] = 1.0
            #= Dont need for CI Testing
            #for j=1:7, k=1:7
            if abs(jac_init[k,j]-jac_init_num[k,j]) > 1e-8
                println(elements_name[j]," ",cartesian_name[k]," ",jac_init[k,j]," ",jac_init_num[k,j]," ",jac_init[k,j]-jac_init_num[k,j])
            end
            if abs(jac_init[k,j]/jac_init_num[k,j]-1) > 1e-12 && jac_init_num[k,j] != 0.0
                println(elements_name[j]," ",cartesian_name[k]," ",jac_init[k,j]," ",jac_init_num[k,j]," ",jac_init[k,j]/jac_init_num[k,j]-1)
            end
        end
        println("Jacobians: difference ",maximum(abs.(jac_init-jac_init_num)))
        println("Jacobians: log diff   ",maximum(abs.(jac_init[jac_init_num .!= 0.0]./jac_init_num[jac_init_num .!= 0.0].-1))) =#
        end
        xvec[:,i] = x
        vvec[:,i] = v
        # Compute finite difference velocity:
        xf, vf = kepler_init(timebig[i] + dtbig, massbig, big.(elements))
        vfvec[:,i] = (xf - x) / dt
    end
    period = elements[1]
    # omega=atan2(elements[4],elements[3])
    omega = atan(elements[4], elements[3])
    semi = (GNEWT * mass * period^2 / 4 / pi^2)^third
    # Compute the focus:
    focus = [semi * ecc,0.,0.]
    xfocus = semi * ecc * (-cos(omega) - sin(omega))
    zfocus = semi * ecc * sin(omega)
    # Check that we have an ellipse (we're assuming that the motion is in the x-z plane; no longer true):
    Atot = 0.0
    dAdt = zeros(Float64, ntime)
    dx = xvec[1,1] - xvec[1,ntime - 1]
    dy = xvec[2,1] - xvec[2,ntime - 1]
    dz = xvec[3,1] - xvec[3,ntime - 1]
    dAvec = 0.5 * [xvec[2,1] * dz - xvec[3,1] * dy;xvec[3,1] * dx - xvec[1,1] * dz;xvec[1,1] * dy - xvec[2,1] * dx]
    # println("dA: ",norm(dAvec)," ",dAvec/norm(dAvec))
    # dA = 0.5*(xvec[3,1]*dx-xvec[1,1]*dz)
    dAdt[1] = norm(dAvec) / (time[2] - time[1])
    Atot += norm(dAvec)
    for i = 2:ntime
        dx = xvec[1,i] - xvec[1,i - 1]
        dy = xvec[2,i] - xvec[2,i - 1]
        dz = xvec[3,i] - xvec[3,i - 1]
        dAvec = 0.5 * [xvec[2,i] * dz - xvec[3,i] * dy;xvec[3,i] * dx - xvec[1,i] * dz;xvec[1,i] * dy - xvec[2,i] * dx]
        #  println("dA: ",norm(dAvec)," ",dAvec/norm(dAvec))
        # dA = 0.5*(xvec[3,i]*dx-xvec[1,i]*dz)
        dAdt[i] = norm(dAvec) / (time[i] - time[i - 1])
        Atot += norm(dAvec)
    end
    # println("Total area: ",Atot," ",pi*sqrt(1-ecc^2)*semi^2," ratio: ",Atot/pi/semi^2/sqrt(1-ecc^2))
    # using PyPlot
    # plot(time,dAdt)
    # axis([minimum(time),maximum(time),0.,1.5*maximum(dAdt)])
    # println("Scatter in dAdt: ",std(dAdt))

    # plot(time,vvec[1,:])
    # plot(time,vfvec[1,:],".")
    # plot(time,vvec[1,:]-vfvec[1,:],".")
    # plot(time,vvec[3,:])
    # plot(time,vfvec[3,:],".")
    # plot(time,vvec[3,:]-vfvec[3,:],".")
    # Check that velocities match finite difference values

    @test isapprox(jac_init, jac_init_num;norm=maxabs)
    @test isapprox(jac_init, convert(Array{Float64,2}, jac_init_big);norm=maxabs)
end
end


text:
# A collection of extensions to the DomainSets package.

using DomainSets: inverse_map, forward_map, MappedDomain

##########################################################################
### Distances and Normals
##########################################################################
export normal, distance
# Function that returns the normal vector when evaluated at the domain boundary
# This vector is of unit length and points to the OUTSIDE.
function normal(d::Domain, x)
    if hascanonicaldomain(d)
        normal(convert(MappedDomain, d), x)
    else
        error("Normal not available for this domain type")
    end
end

# Auxiliary function that returns distance from the boundary in some metric
function distance(x, d::Domain)
    if hascanonicaldomain(d)
        distance(x, convert(MappedDomain, d))
    else
        error("Domain distance not available for this domain type")
    end
end

distance(x, ::UnitSimplex) = min(minimum(x),1-sum(x))

distance(x, d::AbstractInterval) = min(abs(maximum(d)-x),abs(x-minimum(d)))

distance(x, d::UnitBall) = 1-norm(x)

distance(x, d::UnitSphere) = 1-norm(x)

distance(x,d::UnionDomain) = indomain(x,d) ? sum(map(di->max(0,distance(x,di)),components(d))) : maximum(map(di->distance(x,di),components(d)))

normal(d::UnionDomain, x) = normal(components(d)[findmin(map(di->abs(distance(x,di)),components(d)))[2]], x)

distance(x,d::IntersectDomain) = minimum(map(di->distance(x,di),components(d)))

normal(d::IntersectDomain, x) = normal(components(d)[findmin(map(di->distance(x,di),components(d)))[2]], x)

distance(x,d::SetdiffDomain) = indomain(x,d) ? min(abs(distance(x,d.domains[1])),abs(distance(x,d.domains[2]))) : -1*min(abs(distance(x,d.domains[1])),abs(distance(x,d.domains[2])))

normal(d::SetdiffDomain, x) = abs(distance(x,d.domains[1]))<abs(distance(x,d.domains[2])) ? normal(d.domains[1], x) : -1*normal(d.domains[2], x)

distance(x, t::ProductDomain) = minimum(map(distance, x, components(t)))

distance(x, d::MappedDomain) = distance(inverse_map(d)(x),superdomain(d))

function normal(d::MappedDomain, x)
    x = applymap(forward_map(d),normal(superdomain(d),inverse_map(d)(x)))
    x0 = inverse(inverse_map(d))(zeros(size(x)))
   (x-x0)/norm(x-x0)
end
function normal(t::ProductDomain, x)
    index = findmin(map(distance, x, components(t)))[2]
    [(i==index)*normal(component(t,i),x[i]) for i =1:length(components(t))]
end