text:
# See CoreORM.utils.part2.jl for implementation
function util_get_column_type end
function util_is_column_numeric end

# TODO: There's probably a more efficient way than this to retrieve the db name
#         from the LibPQ.Connection
function util_getdbname(dbconn::LibPQ.Connection)
    filter(x -> x.keyword == "dbname",LibPQ.conninfo(dbconn))[1].val
end

function util_getdbhost(dbconn::LibPQ.Connection)
    host = filter(x -> x.keyword == "host",LibPQ.conninfo(dbconn))[1].val
    if ismissing(host)
        host = "localhost"
    end
    return host
end

function util_overwrite_props!(src::T,
                               target::T,
                               ignore_missing_props::Bool) where T <: IEntity
    for psymbol in propertynames(src)
        srcvalue = getproperty(src,psymbol)
        # If the src value is missing and we ignore the missing props we continue
        if (ismissing(srcvalue) && ignore_missing_props)
            continue
        end
        setproperty!(target,psymbol,srcvalue)
    end
end

"""
    util_get_entity_props_for_db_actions(entity::IEntity,
                                         dbconn::LibPQ.Connection,
                                         include_missing::Bool
                                        ;exclude_props::Union{Missing,Vector{Symbol}} = missing)::Dict

Returns a 'flat' dictionary of property names (and column names for the properties
that are IEntity) and values
"""
function util_get_entity_props_for_db_actions(entity::IEntity,
                                              dbconn::LibPQ.Connection,
                                              include_missing::Bool
                                             ;exclude_props::Union{Missing,Vector{Symbol}} = missing)::Dict

    entity_type = typeof(entity)
    orm_module = PostgresORM.get_orm(entity)
    columns_selection_and_mapping = util_get_columns_selection_and_mapping(orm_module)

    props = PostgresORMUtil.getproperties_asdict(entity, include_missing)

    # Exclude the props that are specified in optional argument
    if !ismissing(exclude_props)
        filter!(x -> !(x.first in exclude_props),props)
    end

    # Exclude the non map properties
    props = util_remove_unmapped_properties_from_dict(props,
                                                      columns_selection_and_mapping)

    # Tweak the dictionary of property values by:
    # - replacing the complex types by the IDS (may also replace the property name
    #     the column IDs if the complex type has a composite PK, hence composite FK)
    # - replacing the enums by their IDs
    # - replacing the Dict by a JSON string
    #
    # CAUTION: This functions does not replace all the property names by the
    #            columns names, only for the complex props
    props = util_replace_complex_types_by_id(props,
                                             columns_selection_and_mapping,
                                             entity_type)
    props = util_replace_enums_by_id_if_needed(props, orm_module, dbconn)
    props = util_replace_empty_vector_of_enums_by_missing(props)
    props = util_replace_dict_types(props, orm_module, dbconn)

end

"""
    util_get_entity_props_for_comparison(entity::IEntity,
                                         dbconn::LibPQ.Connection,
                                         include_missing)::Dict

Returns a 'flat' dictionary of property names (and column names for the properties
that are IEntity) and values so that we can compare the values with another
instance of the same type
"""
function util_get_entity_props_for_comparison(entity::IEntity,
                                              dbconn::LibPQ.Connection,
                                              ;drop_trackchanges_properties = true,
                                               drop_unmapped_properties = true)::Dict

    entity_type = typeof(entity)
    orm_module = PostgresORM.get_orm(entity)
    columns_selection_and_mapping = util_get_columns_selection_and_mapping(orm_module)

    props = PostgresORMUtil.getproperties_asdict(entity,
                                              true # include_missing IMPORTANT!!
                                              )

    if drop_unmapped_properties
        props =
           util_remove_unmapped_properties_from_dict(props,
                                                     columns_selection_and_mapping)
    end

    if drop_trackchanges_properties
        props =
            util_remove_trackchanges_properties_from_dict(props, orm_module)
    end

    # Tweak the dictionary of property values
    props = util_replace_complex_types_by_id(props,
                                             columns_selection_and_mapping,
                                             entity_type)
    props = util_replace_enums_by_id_if_needed(props,orm_module,dbconn)
    props = util_replace_empty_vector_of_enums_by_missing(props)
    props = util_replace_dict_types(props,orm_module,dbconn)

    return(props)

end

"""
    util_get_cols_names_and_values(entity::IEntity,
                                   dbconn::LibPQ.Connection)

Transform an entity in a dict with the key being the column name and the value
being the column value
"""
function util_get_cols_names_and_values(entity::IEntity,
                                        dbconn::LibPQ.Connection)

   props_db_actions = util_get_entity_props_for_db_actions(entity,
                                                           dbconn,
                                                           true, # include_missing
                                                           )

   orm_module = get_orm(entity)
   columns_selection_and_mapping = util_get_columns_selection_and_mapping(orm_module)

   props_names = collect(keys(props_db_actions))
   col_names = util_getcolumns(props_names, columns_selection_and_mapping)
   col_values = collect(values(props_db_actions))

   result = Dict{String,Any}()
   for (k,v) in zip(col_names,col_values)
       result[k] = v
   end

   return result

end

"""
    util_get_ids_cols_names_and_values(entity::IEntity)

Transform an entity in a dict with the keys being the ids names and the values
being the ids column values
"""
function util_get_ids_cols_names_and_values(entity::IEntity,
                                            dbconn::LibPQ.Connection)

   cols_names_and_values = util_get_cols_names_and_values(entity,
                                                          dbconn)
   ids_cols_names = util_get_ids_cols_names(entity)

   return filter(x -> x.first in ids_cols_names, cols_names_and_values)
end

function util_diff_entities(previous_state::IEntity, new_state::IEntity,
                            dbconn::LibPQ.Connection,
                            ;drop_trackchanges_properties = true,
                             drop_unmapped_properties = true)::Dict

    previous_state_props =
        util_get_entity_props_for_comparison(previous_state, dbconn
                                             ;drop_trackchanges_properties = drop_trackchanges_properties,
                                             drop_unmapped_properties = drop_unmapped_properties)

    new_state_props =
        util_get_entity_props_for_comparison(new_state,dbconn
                                            ;drop_trackchanges_properties = drop_trackchanges_properties,
                                             drop_unmapped_properties = drop_unmapped_properties)

    props_modified = diff_dict(previous_state_props, new_state_props)

    props_modified
end

function util_clean_and_rename_query_result(data::DataFrame,
                                            mapping::Dict)::DataFrame

    # Clean the mapping from the properties for which we don't have the columns
    #    in the dataframe
    # WARNING: DO NOT USE filter! because that would overwrite 'columns_selection_and_mapping'
    #          in the orm module
    mapping = filter(x -> begin
            cols = PostgresORMUtil.tovector(x.second)
            if all(y -> y in string.(names(data)),cols)
                return true
            else
                return false
            end
            # x.second in string.(names(data))
        end
     ,mapping)

    # Prepare the 2 vectors of properties and columns names
    datamodel_props = Symbol.(collect(keys(mapping)))
    database_cols = Vector{Symbol}()

    for v in collect(values(mapping))
        push!(database_cols,Symbol.(tovector(v))...)
    end

    # Remove duplicates (same columns can be referenced by several properties))
    database_cols = unique(database_cols)

     # Symbol.(collect(values(mapping)))

    # Clean the dataframe from the columns for which we don't have a property
    #    in the mapping
    data = data[:,database_cols]

    # Rename the selected columns
    # rename!(data, f => t for (f, t) = enumerate(mapping)) # does not work
    # rename!(data, f => t for (f, t) = zip(database_cols,
    #                                       datamodel_props))

    # for (f, t) in zip(database_cols,datamodel_props)
    #     rename!(data,f => t)
    # end

    return(data)

end

function util_getcolumns(propertiesnames::Vector,
                         columns_selection_and_mapping::Dict{Symbol,<:Any})

    column_names = String[]

    # NOTE: The mapping may not have the properties because some of the
    #          properties are the columns names of the complex objects
    #          which are added at runtime
    for p in propertiesnames
        if haskey(columns_selection_and_mapping,p)
            push!(column_names,tovector(columns_selection_and_mapping[p])...)
        else
            push!(column_names,string(p))
        end
    end
    column_names

end

function util_id2entity(ptype::DataType,
                        id_values::Vector,
                        load_props_from_db::Bool,
                        load_complex_properties::Bool,
                        dbconn::Union{LibPQ.Connection,Missing})

    # If one of the id values is missing, return missing
    # NOTE: We need to be strict on this and not allow some missing values
    #         because this would cause side effects (eg. when loading the
    #         complex props after a retrieve_entity)
    if any(map(x -> ismissing(x),id_values))
        return missing
    end

    complex_prop = ptype()
    id_properties = util_get_ids_props_names(complex_prop)

    # Check that the ID values that are expected to be IEntity are actually
    #   IEntity, convert them if not
    # @info "util_id2entity for DataType[$ptype] id_properties[$id_properties], id_values[$id_values] "
    for (idx, pair_prop_name_value) in enumerate(zip(id_properties, id_values))
        id_prop_name = pair_prop_name_value[1]
        id_prop_value = pair_prop_name_value[2]
        id_prop_type = util_get_property_real_type(ptype, Symbol(id_prop_name))
        if (id_prop_type <: IEntity && !isa(id_prop_value, IEntity))
            id_values[idx] = util_id2entity(id_prop_type,
                                            tovector(id_prop_value),
                                            false, # load_props_from_db::Bool,
                                            false, #load_complex_properties::Bool,
                                            dbconn)
        end
    end

    # If we retrieve the details of the object we set the
    #   ID property of a filter object and we retrieve the object
    if load_props_from_db
        # Prepare the filter object with the ID property set
        filter_object = ptype()
        # Set the ID property
        PostgresORMUtil.setpropertiesvalues!(filter_object,
                                               id_properties,
                                               id_values)

        queryresult = retrieve_entity(
            filter_object,
            PostgresORM.get_orm(filter_object),
            load_complex_properties,
            dbconn)
        if length(queryresult) == 0
            @warn "Unable to find an entry in database for DataType[$ptype] with ID[$id_values] => return missing."
            return missing
        end
        complex_prop = queryresult[1]

    # If we don't want to retrieve the details of the object we just
    #   create an object and we set it's ID property
    else
        PostgresORMUtil.setpropertiesvalues!(complex_prop,
                                               id_properties,
                                               id_values)
    end # if load_properties_from_db

    return complex_prop

end

function util_dict2entity(props_dict::Dict{Symbol,T},
                          object_type::DataType,
                          building_from_database_result::Bool,
                          retrieve_complex_props::Bool,
                          dbconn::Union{LibPQ.Connection,Missing}) where T <: Any

  # Clone the Dict to make sure that it is of type Dict{Symbol, Any} so that we
  #   can for example replace a value by another one of different type
  props_dict = Dict{Symbol,Any}(props_dict)

  # Drop the entries that do not correspond to a known property
  keys_to_drop = filter(x -> !(x in fieldnames(object_type)),
                        keys(props_dict))
  for k in keys_to_drop
    delete!(props_dict,k)
  end

  # @info props_dict

  # To make it simple, this function only go get some additional information in
  #   the DB if it's handling the result of a query. In other words, when we
  #   create an entity from a dictionnary not originating from a database result
  #   we won't use more information than the one contained in the dictionary.
  if !building_from_database_result
      retrieve_complex_props = false
      dbconn = missing
  end

  # Loop over the properties of the struct and get the values from the dictionary
  for fsymbol in fieldnames(object_type)

      # Get the 'non-missing' type of the property
      ftype = fieldtype(object_type,fsymbol)
      if (typeof(ftype) == Union)
          ftype = get_nonmissing_typeof_uniontype(ftype);
      end

      # Some attributes of the type may not be present
      if !haskey(props_dict,fsymbol)
          continue
      end

      # If attribute is already of the right type do nothing
      if isa(props_dict[fsymbol],ftype)
          continue
      # Override some properties
      elseif ftype <: UUID
          if ismissing(props_dict[fsymbol])
              continue
          end
          props_dict[fsymbol] = UUID(props_dict[fsymbol])

      elseif ftype <: DateTime
          if (ismissing(props_dict[fsymbol]))
              continue
          end
          # "2019-07-24T00:41:49.732Z" becomes "2019-07-24T00:41:49.732"
          date_match =
            match(r"^([0-9]{4}-[0-9]{2}-[0-9]{2}T[0-9]{2}:[0-9]{2}:[0-9]{2}.[0-9]{3})",
                  props_dict[fsymbol])
          props_dict[fsymbol] = DateTime(date_match.match)

      elseif ftype <: ZonedDateTime
          if (ismissing(props_dict[fsymbol]))
              continue
          end
          zdt = string2zoneddatetime(props_dict[fsymbol])
          if isnothing(zdt)
              throw(
                DomainError("Cannot convert string[$(props_dict[fsymbol])] to $ftype for object "
                          * " of type[$object_type]. Please check 'string2zoneddatetime' "))
          else
              props_dict[fsymbol] = zdt
          end

      elseif ftype <: Date
          if (ismissing(props_dict[fsymbol]))
              continue
          end
          # "2019-07-24T00:41:49.732Z" becomes "2019-07-24"
          date_match =
            match(r"^([0-9]{4}-[0-9]{2}-[0-9]{2})",
                  props_dict[fsymbol])
          props_dict[fsymbol] = Date(date_match.match)
      elseif ftype <: Time
          if (ismissing(props_dict[fsymbol]))
              continue
          end
          # "2019-07-24T00:41:49.732Z" becomes "00:41:49"
          time_match =
            match(r"([0-9]{2}:[0-9]{2}:[0-9]{2})",
                  props_dict[fsymbol])
          props_dict[fsymbol] = Time(time_match.match)
      elseif ftype <: Vector{String}
          if (ismissing(props_dict[fsymbol]))
              continue
          end
          # @info props_dict[fsymbol]
          props_dict[fsymbol] =
            PostgresORMUtil.postgresql_string_array_2_string_vector(props_dict[fsymbol])
      elseif ftype <: Vector{T} where T <: Base.Enums.Enum
          props_dict[fsymbol] = string2vector_of_enums(ftype,props_dict[fsymbol])
      elseif ftype <: Enum

          if ismissing(props_dict[fsymbol])
              continue
          end

          # Convert the int or the string to the correct enum instance
          if isa(props_dict[fsymbol],Number)
              props_dict[fsymbol] = int2enum(ftype,
                                             convert(Int,props_dict[fsymbol]))
          else
              props_dict[fsymbol] = string2enum(ftype,
                                                props_dict[fsymbol])
          end
      # If the property is a Dict, we assume that it is serialized as a JSON
      #   string in the database
      elseif ftype <: Dict
          if (ismissing(props_dict[fsymbol]))
              continue
          end
          props_dict[fsymbol] = JSON.parse(props_dict[fsymbol])
      elseif ftype <: IEntity

          # Create an empty instance so that we can access the ORM of the type
          dummy_object = object_type()

          # If the property value is set, either create the
          #   instance with the ID property only or with all the properties
          #   that we rerieve in the database
          if !ismissing(props_dict[fsymbol])

              ptype = util_get_property_real_type(object_type, fsymbol)

              # Check that the ptype is compatible
              if !(ptype <: ftype)
                  throw(DomainError("Cannot convert $ptype to $ftype for object "
                                  * " of type[$object_type]. Please check 'types_override' "
                                  * " in the ORM[$(PostgresORM.get_orm(dummy_object))]"))
              end

              # Retrieve the ID property names of the complex property
              complex_prop_idprops = util_get_ids_props_names(ptype())

              # If building from the result of a query we convert the id to an entity
              if building_from_database_result
                  load_props_from_db = retrieve_complex_props

                  complex_prop = util_id2entity(ptype,
                                                getdictvalues(props_dict[fsymbol],complex_prop_idprops),
                                                # [props_dict[fsymbol]], # id_value::Any,
                                                load_props_from_db, # load_props_from_db
                                                false, # load_complex_properties
                                                dbconn)
              else
                  # If we are building from a dictionnary
                  if isa(props_dict[fsymbol],Dict)
                      complex_prop =
                        util_dict2entity(props_dict[fsymbol], # this is a dictionary
                                         ptype,
                                         building_from_database_result,
                                         retrieve_complex_props,
                                         dbconn)
                  else
                      load_props_from_db = retrieve_complex_props
                      # TODO: properly handle the case of composed foreign key
                      complex_prop =
                          util_id2entity(ptype,
                                         getdictvalues(props_dict[fsymbol],complex_prop_idprops),
                                         # [props_dict[fsymbol]], # id_value::Any,
                                         load_props_from_db, # load_props_from_db
                                         false, # load_complex_properties
                                         dbconn)
                  end
              end

              # Update the properties dictionary with the complex object
              props_dict[fsymbol] = complex_prop

          end # if !ismissing(props_dict[fsymbol])

      # This will handle the case of Vectors and Arrays of IEntity.
      elseif eltype(ftype) <: IEntity

          # Create an empty instance so that we can access the ORM of the type
          dummy_object = object_type()

          # This should only happen when dealing with results that are not
          #   coming from a query result. In which case the input should be a
          #   vector of Dict.
          if !ismissing(props_dict[fsymbol])

              # Check if there is a type override, if not we use the type found
              #  in the struct
              ptype = util_get_property_real_type(object_type, fsymbol)

              # Check that the ptype is compatible
              if !(eltype(ptype) <: eltype(ftype))
                  throw(DomainError("Cannot convert $ptype to $ftype for object "
                                  * " of type[$object_type]. Please check 'types_override' "
                                  * " in the ORM[$(PostgresORM.get_orm(dummy_object))]"))
              end

              # Check that the elements type is not abstract
              if isabstracttype(eltype(ptype))
                  throw(DomainError("Cannot instantiate type[$ptype] of property[:$fsymbol] because it's a vector of abstract elements. "
                                  * " Please set a non-abstract vector in 'types_override'"
                                  * " in the ORM[$(PostgresORM.get_orm(dummy_object))]"))
              end

              # Instantiate a new array/vector of the given type
              complex_vector_prop = ptype()

              # For every dictionary in the vector we create an entity and add
              #   it to the vector
              for dct in props_dict[fsymbol]
                  push!(complex_vector_prop,
                        util_dict2entity(dct, # this is a dictionary
                                         eltype(ptype),
                                         building_from_database_result,
                                         retrieve_complex_props,
                                         dbconn)
                  )
              end

              # Update the properties dictionary with the complex object
              props_dict[fsymbol] = complex_vector_prop

          end # if !ismissing(props_dict[fsymbol])
      end # eltype(ftype) <: IEntity

  end # for fsymbol in fieldnames(object_type)

  obj = object_type(dict2namedtuple(props_dict))
  return(obj)
end

function util_get_property_real_type(datatype::DataType, propname::Symbol)

    dummy_object = datatype()
    orm_module = get_orm(dummy_object)

    # Initialize the type with what we find in the
    ptype = get_nonmissing_typeof_uniontype(fieldtype(datatype,propname))

    # Check if there is a type override for this property, if not we
    #   use the property type

    if (isdefined(orm_module,:types_override)
        && propname in collect(keys(orm_module.types_override)))
        ptype = PostgresORM.get_orm(dummy_object).types_override[propname]
    end
    return ptype
end

function util_convert_namedtuple_to_object(props::NamedTuple, object_type::DataType,
                                           retrieve_complex_props::Bool,
                                           dbconn::Union{LibPQ.Connection,Missing})

    # Initialize the properties of the object with what was
    #   retrieved fom database
    # NOTE: we need to force the dictionary to accept any type of value
    props_dict::Dict{Symbol,Any} = namedtuple2dict(props)

    # At this point we have a dict of values of the dataframe as returned by
    #   the query: no renaming no encapuslation for complex properties
    #   (eg.: film_id, film_release_year, actor_id).
    # We transform this dict to another dict that suits the hierarchy of the julia
    #   data model
    #  eg. {film_id, film_release_year, actor_id} becomes {film{id,release_year}, actor{id}}
    # @info keys(props_dict)
    # @info length(props_dict[:possible_values_as_str_arr])
    props_dict =
        util_convert_flatdictfromdb_to_structuredrenameddict(props_dict,
                                                             object_type)

    # @info keys(props_dict)
    # @info length(props_dict[:possibleValuesAsStrArr])
    # Now that the dict follows the structure of the datamodel we transform it
    util_dict2entity(props_dict,
                     object_type,
                     true, # building from database result
                     retrieve_complex_props,
                     dbconn)

end

function util_convert_flatdictfromdb_to_structuredrenameddict(flatdict::Dict,
                                                              object_type::DataType)
    sample_object = object_type()
    orm_module = get_orm(sample_object)

    result = Dict()

    # Loop over the properties of the struct and try to get the corresponding
    #   values in the flat dict
    for (propname,prop_colnames) in util_get_columns_selection_and_mapping(orm_module)

        # Get the proptype in order to know if we are dealing with a complex
        #   property or not
        proptype = util_get_property_real_type(object_type, propname)

        # if the property is a complex prop we build a dict with the IDs of
        #    the complex object as described in the ORM of the complex property
        if proptype <: IEntity

            result[propname] = Dict()

            # Retrieve the type of the complex object
            complexprop_instance = proptype()
            complexprop_idprops = util_get_ids_props_names(complexprop_instance)
            complexprop_orm = get_orm(complexprop_instance)

            # Check that the number of columns for the FK is consistent with
            #   the number of properties used as PK in the targeted object
            if length(complexprop_idprops) != length(tovector(prop_colnames))
                error("There is an inconsistency in the ORMs definition,
                type[$object_type] has property[$propname] of type[$proptype]
                which is mapped with $(length(tovector(prop_colnames))) columns
                but type[$proptype] has only $(length(complexprop_idprops)) IDs
                properties")
            end

            # If there is a composed foreign key we need to correctly map the
            #   column values to the column ids in the target table
            if isa(prop_colnames, Array) && length(prop_colnames) > 1

                # Loop over the FK columns and use the PKs in the same order
                # TODO: Handle the case where the order of the FKs and PKs columns
                #         is not the same
                counter = 0
                for colname in prop_colnames

                    # Get the value from the flatdict (if exists)
                    if haskey(flatdict, Symbol(colname))
                        colvalue = flatdict[Symbol(colname)]
                    else
                        break
                    end

                    complexprop_idprop = complexprop_idprops[counter+=1]
                    result[propname][complexprop_idprop] = colvalue

                end

            # If it is a non componsed foreign key then we can simply use the
            #   id property of the targeted object
            else
                complexprop_idprop = complexprop_idprops[1]
                complexprop_colname = Symbol(tovector(prop_colnames)[1])
                if haskey(flatdict, complexprop_colname)
                    result[propname][complexprop_idprop] = flatdict[complexprop_colname]
                end

            end # ENFOF `if isa(prop_colnames, Array) && length(prop_colnames) > 1`

        # If the property is simple then we just look for the value in the dict
        else
            colname = Symbol(tovector(prop_colnames)[1])
            if haskey(flatdict,colname)
                result[propname] = flatdict[colname]
            end
        end # ENDOF `if proptype <: IEntity`

    end # ENDOF `for (propname,prop_colnames) in util_get_columns_selection_and_mapping(orm_module)`

    return result

end

# Create the selection as expected by the Query.jl macro
# Eg: `{id = i.document_id, lastname = i.lastname}` )
function util_get_curlybrackets_expression_for_linq_select(
        columns_selection_and_mapping::Dict)

    target_properties = string.(keys(columns_selection_and_mapping))
    source_cols =  values(columns_selection_and_mapping)
    columns_selection_and_mapping_tmp = "{" * join(string.(target_properties) .* " = i." .* source_cols,", ") * "}"

    Meta.parse(columns_selection_and_mapping_tmp)
end

function util_remove_unmapped_properties_from_dict(
        props::Dict,
        columns_selection_and_mapping::Dict)
    @from i in props begin
    @where i.first in collect(keys(columns_selection_and_mapping))
    @select i
    @collect Dict
    end
end

function util_remove_trackchanges_properties_from_dict(
        props::Dict,
        orm_module::Module)

    if (isdefined(orm_module,:track_changes)
        && getproperty(orm_module,:track_changes))

        trackchanges_properties = Symbol[]
        if isdefined(orm_module,:creator_property)
            push!(trackchanges_properties,orm_module.creator_property)
        end
        if isdefined(orm_module,:editor_property)
            push!(trackchanges_properties,orm_module.creator_property)
        end
        if isdefined(orm_module,:creation_time_property)
            push!(trackchanges_properties,orm_module.creation_time_property)
        end
        if isdefined(orm_module,:update_time_property)
            push!(trackchanges_properties,orm_module.update_time_property)
        end

        if length(trackchanges_properties) == 0
            return props
        end

        @from i in props begin
        @where !(i.first in trackchanges_properties)
        @select i
        @collect Dict
        end

    else
        return props
    end
end

function util_replace_complex_types_by_id(props::Dict, mapping::Dict, data_type::DataType)

    mapping = deepcopy(mapping)
    orm_module = get_orm(data_type())

    # NOTE: We loop on the result of a filter that does nothing because if not
    #        the pairs that are added dynamically in the loop will be looped
    #        over as well
    for (prop_symbol, prop_val) in filter(x -> true, props)

      _fieldtype = util_get_property_real_type(data_type, prop_symbol)

      # NOTE: We test on the fieldtype because the property value can be missing
      if _fieldtype <: IEntity

          # @info "Replace[$prop_symbol] with ids"

          dummy_propval = _fieldtype() # Used in case the prop val is missing

          prop_orm_module = get_orm(_fieldtype())
          prop_columns_selection_and_mapping =
            util_get_columns_selection_and_mapping(prop_orm_module)

          # Retrieve the names of the ID properties in the target struct
          id_props_in_target_struct = util_get_ids_props_names(_fieldtype())

          # @info "id_props_in_target_struct[$id_props_in_target_struct] ($_fieldtype)"

          # Whatever we get from the mapping we make it a vector
          idcols = tovector(mapping[prop_symbol])

          prop_val_as_dict =
              PostgresORMUtil.getproperties_asdict(
                if ismissing(prop_val) dummy_propval else prop_val end,
                 true # include missing
                 )

          prop_val_ids_values =
              filter!(x -> first(x) in id_props_in_target_struct,
                      prop_val_as_dict)
          prop_val_ids_values = util_replace_complex_types_by_id(
              prop_val_ids_values,
              util_get_columns_selection_and_mapping(
                if ismissing(prop_val) dummy_propval else prop_val end
               ),
              _fieldtype)

          if length(idcols) != length(prop_val_ids_values)
              error_msg = ""
              error_msg *= "Problem while flattening property[$prop_symbol] of"
              error_msg *= " an $data_type. The number of columns[$(length(idcols))]"
              error_msg *= " found in ORM[$orm_module] does not match the number"
              error_msg *= " of columns[$(length(prop_val_ids_values))] found"
              error_msg *= " for the IDs in ORM [$prop_orm_module]"
              error(error_msg)
          end

          ids_props_names = collect(keys(prop_val_ids_values))
          ids_cols_names = util_getcolumns(ids_props_names,
                                           prop_columns_selection_and_mapping)
          ids_props_values = collect(values(prop_val_ids_values))

          # Make sure the values are ordered in the same order as the one given
          #  by the mapping of the target entity's ORM
          idvalues = []

          prop_colnames_in_target_type = util_get_ids_cols_names(prop_orm_module)
          for pcol in prop_colnames_in_target_type
              for (n,c,v) in zip(ids_props_names,ids_cols_names,ids_props_values)
                  if (pcol == c)
                      push!(idvalues,v)
                  end
              end
          end

          for (k,v) in zip(idcols, idvalues)
             # This prevents a property using the same primary columns to
             #   overwrite a non missing value
             if ismissing(v) && haskey(props,Symbol(k)) && !ismissing(props[Symbol(k)])
               continue
             else
               props[Symbol(k)] = v
             end
          end

          # for (n,c,v,cbis) in zip(ids_props_names,ids_cols_names,ids_props_values,idcols)
          #     @info "$_fieldtype -> $n, $c, $v"
          #     props[Symbol(cbis)] = v
          # end

          # counter = 0
          # for (k,v) in prop_val_ids_values
          #     id_column_name = idcols[counter+=1] # the name as used by the FK
          #     props[Symbol(id_column_name)] = v
          # end

          # Remove the previous reference that has been replaced
          delete!(props,prop_symbol)

      # Convert vector of IEntity
      elseif eltype(prop_val) <: IEntity
          props[prop_symbol] = Vector{Any}()
          for e in prop_val

              push!(props[prop_symbol],
              getproperty(e,
                          PostgresORM.get_orm(e).id_property))
          end
      end # ENDOF if isa(prop_val,IEntity)

  end # ENDOF for (prop_symbol, prop_val) in props

  return props
end

function util_replace_enums_by_id_if_needed(props::Dict,
                                            orm_module::Module,
                                            dbconn::LibPQ.Connection)

    table_name = util_get_table_name(orm_module)
    schema = util_get_schema_name(orm_module)
    columns_selection_and_mapping = util_get_columns_selection_and_mapping(orm_module)

    for (prop_symbol, prop_val) in props
      if typeof(prop_val) <: Enum
          colname = columns_selection_and_mapping[prop_symbol]
          # Chech in the database if the enum os stored as a numeric
          if util_is_column_numeric(dbconn,
                                    table_name,
                                    schema,
                                    colname)
              props[prop_symbol] = enum2int(prop_val)
          end

      end
    end
    return props
end

"""
    util_replace_empty_vector_of_enums_by_missing(props::Dict)

Replace empty vector of an enum by missing. This allows developers to avoid testing both for
nullity and emptiness of the property
"""
function util_replace_empty_vector_of_enums_by_missing(props::Dict)
    for (prop_symbol, prop_val) in props
        if typeof(prop_val) <: Vector{T} where T <: Base.Enums.Enum
            if isempty(prop_val)
                props[prop_symbol] = missing
            end
        end
    end
    return props
end

function util_replace_dict_types(props::Dict,
                                 orm_module::Module,
                                 dbconn::LibPQ.Connection)
    for (prop_symbol, prop_val) in props
       # @info "$prop_symbol is a dict: "
       if isa(prop_val,Dict)
           props[prop_symbol] = JSON.json(prop_val)
       end
    end
    return props
end

function util_get_ids_props_names(orm_module::Module)
    if isdefined(orm_module,:get_id_props)
        return orm_module.get_id_props()
    # This is support for the legacy way of declaring the id properties
    elseif isdefined(orm_module,:id_property)
        return PostgresORMUtil.tovector(orm_module.id_property)
    else
        error("orm_module[$orm_module] is missing a definition for the id property")
    end
end

function util_get_ids_props_names(o::IEntity)
    return util_get_ids_props_names(PostgresORM.get_orm(o))
end

function util_get_table_name(orm_module::Module)

    tablename = if isdefined(orm_module,:get_table_name)
            orm_module.get_table_name()
        # Support of legacy naming
        elseif isdefined(orm_module,:gettablename)
            orm_module.gettablename()
        else
            error("orm_module[$orm_module] is missing 'get_table_name()''")
        end

    # Handle the case where the schema has been put in the table name
    if occursin(".",tablename)
        schema_and_table_names = string.(split(tablename,'.'))
        schema = schema_and_table_names[1]
        tablename = schema_and_table_names[2]
        return tablename
    end
    return tablename
end

function util_get_table_name(o::IEntity)
    return util_get_table_name(PostgresORM.get_orm(o))
end

function util_get_table_with_schema_name(orm_module::Module)
    return "$(util_get_schema_name(orm_module)).$(util_get_table_name(orm_module))"
end

function util_get_table_with_schema_name(o::IEntity)
    return util_get_table_with_schema_name(PostgresORM.get_orm(o))
end

function util_get_schema_name(orm_module::Module)

    if isdefined(orm_module,:get_schema_name)
        return orm_module.get_schema_name()
    # This is support for legacy ORMs where the schema was declared together with
    #   the table name
    elseif isdefined(orm_module,:gettablename)
        tablename = orm_module.gettablename()
        if occursin(".",tablename)
            schema_and_table_names = string.(split(tablename,'.'))
            schema = schema_and_table_names[1]
            return schema
        else
            return "public"
        end
    elseif isdefined(orm_module,:get_table_name)
        tablename = orm_module.get_table_name()
        if occursin(".",tablename)
            schema_and_table_names = string.(split(tablename,'.'))
            schema = schema_and_table_names[1]
            return schema
        else
            error("orm_module[$orm_module] is missing 'get_schema_name()''")
        end
    else
        error("orm_module[$orm_module] is missing 'get_schema_name()''")
    end

end

function util_get_schema_name(o::IEntity)
    return util_get_schema_name(PostgresORM.get_orm(o))
end

function util_get_columns_selection_and_mapping(orm_module::Module)
    if isdefined(orm_module,:get_columns_selection_and_mapping)
        return orm_module.get_columns_selection_and_mapping()
    # This is support for the legacy way of declaring the id properties
    elseif isdefined(orm_module,:columns_selection_and_mapping)
        return orm_module.columns_selection_and_mapping
    else
        error("orm_module[$orm_module] is missing 'get_columns_selection_and_mapping'")
    end
end

function util_get_columns_selection_and_mapping(o::IEntity)
    return util_get_columns_selection_and_mapping(PostgresORM.get_orm(o))
end

function util_get_onetomany_counterparts(orm_module::Module)
    if isdefined(orm_module,:get_onetomany_counterparts)
        return orm_module.get_onetomany_counterparts()
    # This is support for the legacy way of declaring the id properties
    elseif isdefined(orm_module,:onetomany_counterparts)
        return orm_module.onetomany_counterparts
    else
        error("orm_module[$orm_module] is missing 'get_onetomany_counterparts'")
    end
end

function util_get_onetomany_counterparts(o::IEntity)
    return util_get_onetomany_counterparts(PostgresORM.get_orm(o))
end

function util_get_types_override(orm_module::Module)
    if isdefined(orm_module,:get_types_override)
        return orm_module.get_types_override()
    # This is support for the legacy way of declaring the id properties
    elseif isdefined(orm_module,:types_override)
        return orm_module.types_override
    else
        error("orm_module[$orm_module] is missing 'get_types_override'")
    end
end

function util_get_types_override(o::IEntity)
    return util_get_types_override(PostgresORM.get_orm(o))
end

function util_get_track_changes(orm_module::Module)
    if isdefined(orm_module,:get_track_changes)
        return orm_module.get_track_changes()
    # This is support for the legacy way of declaring the id properties
    elseif isdefined(orm_module,:track_changes)
        return orm_module.track_changes
    else
        return false # default
    end
end

function util_get_track_changes(o::IEntity)
    return util_get_track_changes(PostgresORM.get_orm(o))
end

function util_get_creator_property(orm_module::Module)
    if isdefined(orm_module,:get_creator_property)
        return orm_module.get_creator_property()
    # This is support for the legacy way of declaring the id properties
    elseif isdefined(orm_module,:creator_property)
        return orm_module.creator_property
    else
        return missing
    end
end

function util_get_creator_property(o::IEntity)
    return util_get_creator_property(PostgresORM.get_orm(o))
end

function util_get_editor_property(orm_module::Module)
    if isdefined(orm_module,:get_editor_property)
        return orm_module.get_editor_property()
    # This is support for the legacy way of declaring the id properties
    elseif isdefined(orm_module,:editor_property)
        return orm_module.editor_property
    else
        return missing
    end
end

function util_get_editor_property(o::IEntity)
    return util_get_editor_property(PostgresORM.get_orm(o))
end

function util_get_creation_time_property(orm_module::Module)
    if isdefined(orm_module,:get_creation_time_property)
        return orm_module.get_creation_time_property()
    # This is support for the legacy way of declaring the id properties
    elseif isdefined(orm_module,:creation_time_property)
        return orm_module.creation_time_property
    else
        return missing
    end
end

function util_get_creation_time_property(o::IEntity)
    return util_get_creation_time_property(PostgresORM.get_orm(o))
end

function util_get_update_time_property(orm_module::Module)
    if isdefined(orm_module,:get_update_time_property)
        return orm_module.get_update_time_property()
    # This is support for the legacy way of declaring the id properties
    elseif isdefined(orm_module,:update_time_property)
        return orm_module.update_time_property
    else
        return missing
    end
end

function util_get_update_time_property(o::IEntity)
    return util_get_update_time_property(PostgresORM.get_orm(o))
end

function util_get_ids_cols_names(orm_module::Module)
    props = util_get_ids_props_names(orm_module)
    columns_selection_and_mapping = util_get_columns_selection_and_mapping(orm_module)
    result = []
    for p in props
        propcols = tovector(columns_selection_and_mapping[p])
        push!(result,propcols...)
    end
    return result
end

function util_get_ids_cols_names(o::IEntity)
    orm_module = PostgresORM.get_orm(o)
    util_get_ids_cols_names(orm_module)
end

function util_prettify_id_values_as_str(object::T) where T <: IEntity

    idprops_with_values = []
    idprops = util_get_ids_props_names(object)

    for idprop in idprops
        idpropvalue = getproperty(object, idprop)
        if idpropvalue isa IEntity
            push!(idprops_with_values,
                  (name = idprop, value = util_prettify_id_values_as_str(idpropvalue)))
        else
            push!(idprops_with_values,
                  (name = idprop, value = string(idpropvalue)))
        end
    end

    # If there is only one value no need to specify the name of the property
    if length(idprops_with_values) == 1
        return string(first(idprops_with_values).value)
    end

    result = ""
    strings = []
    for idprop_with_value in idprops_with_values
        push!(strings,"$(idprop_with_value.name)[$(idprop_with_value.value)]")
    end
    return join(strings,", ")

    # return result

    # return idprops_with_values
    # return string(idprops_with_values)

    # str = string(values)
    # my_match = match(r"^Any\[(.*)\]$",str)
    # if isnothing(my_match)
    #     my_match = match(r"^\[(.*)\]$",str)
    # end
    # str = string(collect(my_match.captures)[1])
    # str = replace(str,", missing"=>"")
    # str = replace(str,"missing, "=>"")
    # return str
end


text:
abstract type AbstractNLSolverAlgorithm end
abstract type AbstractNLSolverCache end

@enum NLStatus::Int8 begin
  FastConvergence     = 2
  Convergence         = 1
  SlowConvergence     = 0
  VerySlowConvergence = -1
  Divergence          = -2
end

# solver

mutable struct NLSolver{iip,uType,rateType,uTolType,kType,gType,cType,du1Type,ufType,jcType,lsType,C1,C<:AbstractNLSolverCache}
  z::uType
  dz::uType
  tmp::uType
  ztmp::uType
  k::rateType
  ηold::uTolType
  κ::kType
  γ::gType
  c::cType
  max_iter::Int
  nl_iters::Int
  status::NLStatus
  fast_convergence_cutoff::C1
  du1::du1Type
  uf::ufType
  jac_config::jcType
  linsolve::lsType
  weight::uType
  cache::C
end

# algorithms

struct NLFunctional{K,C} <: AbstractNLSolverAlgorithm
  κ::K
  fast_convergence_cutoff::C
  max_iter::Int
end

NLFunctional(; κ=1//100, max_iter=10, fast_convergence_cutoff=1//5) = NLFunctional(κ, fast_convergence_cutoff, max_iter)

struct NLAnderson{K,D,C} <: AbstractNLSolverAlgorithm
  κ::K
  fast_convergence_cutoff::C
  max_iter::Int
  max_history::Int
  aa_start::Int
  droptol::D
end

NLAnderson(; κ=1//100, max_iter=10, max_history::Int=5, aa_start::Int=1, droptol=nothing, fast_convergence_cutoff=1//5) =
  NLAnderson(κ, fast_convergence_cutoff, max_iter, max_history, aa_start, droptol)

struct NLNewton{K,C1,C2} <: AbstractNLSolverAlgorithm
  κ::K
  max_iter::Int
  fast_convergence_cutoff::C1
  new_W_dt_cutoff::C2
end

NLNewton(; κ=1//100, max_iter=10, fast_convergence_cutoff=1//5, new_W_dt_cutoff=1//5) = NLNewton(κ, max_iter, fast_convergence_cutoff, new_W_dt_cutoff)

# caches

mutable struct NLNewtonCache{W,J,T,C} <: AbstractNLSolverCache
  new_W::Bool
  W::W
  J::J
  W_dt::T
  new_W_dt_cutoff::C
end

mutable struct NLNewtonConstantCache{W,J,C} <: AbstractNLSolverCache
  W::W
  J::J
  new_W_dt_cutoff::C
end

struct NLFunctionalCache{uType} <: AbstractNLSolverCache
  z₊::uType
end

struct NLFunctionalConstantCache <: AbstractNLSolverCache end

mutable struct NLAndersonCache{uType,gsType,QType,RType,gType,D} <: AbstractNLSolverCache
  z₊::uType
  dzold::uType
  z₊old::uType
  Δz₊s::gsType
  Q::QType
  R::RType
  γs::gType
  aa_start::Int
  droptol::D
end

mutable struct NLAndersonConstantCache{gsType,QType,RType,gType,D} <: AbstractNLSolverCache
  Δz₊s::gsType
  Q::QType
  R::RType
  γs::gType
  aa_start::Int
  droptol::D
end

ext: jl


text:
function eg1()
    x = 2
    y = eg1a(x)
    x = eg1b()
    return x
end
eg1a(x) = x
eg1b() = eg1c()
eg1c() = 7

function eg2()
    x = 2
    y = eg2()
    return x
end
eg3() = 5

function badcake1()
    y=0
    f = ()-> y+=1
    f()
    return y
end

function badcake2()
    y=sum([1 for ii in 1:3])
    return y
end

GLOBAL_VAR=0
function badcake3()
    global GLOBAL_VAR
    GLOBAL_VAR+=1
end

badbread1() = ()->99

string_demo1() = "outer" * string_demo1a()
string_demo1a() = "inner"

fma_demo1(a) = 8*a+2
fma_demo2(a) = (8*a)*a+(sqrt.(a))


text:
# File for plotting Figure 2

using JLD2
using Statistics, Distributions
using StatsBase
using DataFrames
using Gadfly
using CSV

import Cairo, Fontconfig

### Load data
load_dict = load("data/deep_sigma1.jld2")
plot_data = load_dict["plot_data"]

tmp = plot_data[plot_data.depth .< 5, :]

function EmpiricalDistribution(data::Vector{T} where T <: Real)
    sort!(data) #sort the observations
    empirical_cdf = ecdf(data) #create empirical cdf
    data_clean = unique(data) #remove duplicates to avoid allunique error
    cdf_data = empirical_cdf.(data_clean) #apply ecdf to data
    pmf_data = vcat(cdf_data[1],diff(cdf_data)) #create pmf from the cdf
    DiscreteNonParametric(data_clean,pmf_data) #define distribution
end

# Get CDFs for different depths
a1 = EmpiricalDistribution(tmp[tmp.depth.==1,:deviation_from_median])
a2 = EmpiricalDistribution(tmp[tmp.depth.==2,:deviation_from_median])
a3 = EmpiricalDistribution(tmp[tmp.depth.==3,:deviation_from_median])
a4 = EmpiricalDistribution(tmp[tmp.depth.==4,:deviation_from_median])

# Add colors
colors = ["#264653", "#E7C15F", "#2A9D8F", "#C33149"]
depths = ["1", "2", "3", "4"]
dist = [a1, a2, a3, a4]

# Plot
p = plot(Guide.xticks(ticks=collect(0:0.2:1)), Guide.yticks(ticks=collect(0:0.2:1)));

for i in 1:length(dist)
    push!(p, layer(x=dist[i].support, y=cumsum(dist[i].p), Geom.line, Theme(default_color=color(colors[i]), line_width=0.8mm)));
end

push!(p, Guide.manual_color_key("Depth", depths, colors))
push!(p, Guide.title("Error term variance = 1"))
push!(p, Guide.YLabel("CFD"))
push!(p, Guide.XLabel("λ"))
# push!(p, Theme(line_width=1cm))
# push!(p, Theme(key_position=:none))

draw(PNG("figures/graphs/deep_sigma1.png", 20cm, 12cm, dpi=300), p)

max_stars_repo_name: christianhilscher/emma

max_issues_repo_name: christianhilscher/emma

max_forks_repo_name: christianhilscher/emma


text:
#
# Iterative Halo orbit and manifold solvers
#

"""
Returns the Monodromy Matrix for a Halo orbit.
"""
function monodromy(orbit::CR3BPOrbit, T; verify=true, reltol = 1e-14, abstol = 1e-14)
    initial = Orbit(CartesianStateWithSTM(state(orbit)), system(orbit), epoch(orbit); frame=frame(orbit))
    final   = propagate(orbit, T; reltol=reltol, abstol=abstol, save_everystep=false)[end]

    if verify && !(state(initial)[1:6] ≈ final[1:6])
        throw(ErrorException("The provided `orbit` is not periodic!"))
    end

    SMatrix{6,6}(final[7:end]...) |> Matrix
end

"""
Returns true if a `RestrictedThreeBodySystem` is numerically periodic.
"""
function isperiodic(orbit::CR3BPOrbit, T; reltol = 1e-14, abstol = 1e-14)
    final = propagate(orbit, T; reltol=reltol, abstol=abstol, save_everystep=false)[end]
    return state(orbit)[1:6] ≈ final[1:6]
end

"""
Returns a numerical solution for a Halo orbit about `L`.

__Arguments:__
- `μ`: Non-dimensional mass parameter for the CR3BP system.
- `Az`: Desired non-dimensional Z-amplitude for Halo orbit.
- `ϕ`: Desired Halo orbit phase.
- `L`: Lagrange point to orbit (L1 or L2).
- `hemisphere`: Specifies northern or southern Halo orbit.

__Outputs:__
- Tuple of initial states: `(r, v)` where `r::Vector{<:AbstractFloat}`, `v::Vector{<:Abstractfloat}`.
- Throws `ArgumentError` if L is not `:L1` or `:L2`

__References:__

The iterative scheme was pulled from directly from literature
and sample code, including Rund 2018,
and Dr. Mireles' lecture notes and EarthSunHaloOrbit_NewtonMewhod.m
file available on their website.
Specifically, the half-period iterative scheme (the `F` matrix
in the source code, and corresponding "next guess" computation)
was ported __directly__ from Dr. Mireles' public code and notes, which are
available online.
- [Dr. Mireles Notes](http://cosweb1.fau.edu/~jmirelesjames/hw5Notes.pdf)
- [Dr. Mireles Code](http://cosweb1.fau.edu/~jmirelesjames/notes.html)
- [Rund, 2018](https://digitalcommons.calpoly.edu/theses/1853/).
"""
function halo(μ; Az=0.0, L=1, hemisphere=:northern,
              tolerance=1e-8, max_iter=100,
              reltol=1e-14, abstol=1e-14,
              nan_on_fail = true, disable_warnings = false)

    r₀, v₀, Τ = analyticalhalo(μ; Az=Az, ϕ=0.0, L=L, hemisphere=hemisphere)
    r₀ = r₀[1,:]
    v₀ = v₀[1,:]
    τ  = Τ/2
    δu = Vector{typeof(τ)}(undef, 6)
    Id = Matrix(I(6))

    for i ∈ 1:max_iter

        problem = ODEProblem(
            AstrodynamicalModels.CR3BPWithSTMVectorField,
            vcat(r₀, v₀, Id...),
            (0.0, τ),
            (μ,)
        )

        retcode, rₛ, vₛ, Φ = let
            sols  = solve(problem, Vern9(); reltol=reltol, abstol=abstol)
            final = sols.u[end]
            sols.retcode, final[1:3], final[4:6], SMatrix{6,6}(final[7:end]...)
        end

        AstrodynamicalModels.CR3BPVectorField(δu, vcat(rₛ, vₛ), (μ,), NaN)
        ∂vₛ = δu[4:6]

        # All code in this `if, else, end` block is ported from
        # Dr. Mireles' MATLAB code, which is available on his
        # website: http://cosweb1.fau.edu/~jmirelesjames/notes.html.
        # Lecture notes, which describe this algorithm further,
        # are available for reference as well:
        # http://cosweb1.fau.edu/~jmirelesjames/hw5Notes.pdf
        if Az ≉ 0
            F = @SMatrix [
                Φ[4,1] Φ[4,5] ∂vₛ[1];
                Φ[6,1] Φ[6,5] ∂vₛ[3];
                Φ[2,1] Φ[2,5]  vₛ[2]
            ]

            xᵪ = SVector(r₀[1], v₀[2], τ) - inv(F) * SVector(vₛ[1], vₛ[3], rₛ[2])

            r₀[1] = xᵪ[1]
            v₀[2] = xᵪ[2]
            τ     = xᵪ[3]
        else
            F = @SMatrix [
                Φ[4,3] Φ[4,5] ∂vₛ[1];
                Φ[6,3] Φ[6,5] ∂vₛ[3];
                Φ[2,3] Φ[2,5]  vₛ[2]
            ]

            xᵪ = SVector(r₀[3], v₀[2], τ) - inv(F) * SVector(vₛ[1], vₛ[3], rₛ[2])

            r₀[3] = xᵪ[1]
            v₀[2] = xᵪ[2]
            τ     = xᵪ[3]
        end

        if abs(vₛ[1]) ≤ tolerance && abs(vₛ[3]) ≤ tolerance
            break;
        elseif τ > 5 * one(τ)
            disable_warnings || @warn "Unreasonably large halo period, $τ, ending iterations."
            !nan_on_fail || return [NaN, NaN, NaN], [NaN, NaN, NaN], NaN
            break
        elseif i == max_iter
            disable_warnings || @warn "Desired tolerance was not reached, and iterations have hit the maximum number of iterations: $max_iter."
            !nan_on_fail || return [NaN, NaN, NaN], [NaN, NaN, NaN], NaN
            break
        elseif retcode != :Success
            !disable_warnings || @warn "Integrator returned $(string(retcode))."
            !nan_on_fail || return [NaN, NaN, NaN], [NaN, NaN, NaN], NaN
            break
        end
    end

    return r₀, v₀, 2τ

end

"""
A `halo` wrapper! Returns a `CircularRestrictedThreeBodyOrbit`.
Returns a tuple: `halo_orbit, halo_period`.
"""
function halo(sys::CR3BPParameters, epoch=TAIEpoch(now()); Az=0.0, kwargs...)

    if Az isa Unitful.Length
        Az = Az / lengthunit(sys)
    end

    r,v,T = halo(massparameter(sys); Az = Az, kwargs...)
    cart  = CartesianState(vcat(r,v); lengthunit=lengthunit(sys), timeunit=timeunit(sys), angularunit=angularunit(sys))
    orbit = Orbit(cart, sys)
    return orbit, T
end

"""
Given a __periodic__ `CircularRestrictedThreeBodyOrbit`,
returns a nearby `CircularRestrictedThreeBodyOrbit` on the
__stable manifold__ of the initial state.
"""
function manifold_perturbation(orbit::CR3BPOrbit, V::AbstractVector; eps = 1e-8)
    V = normalize(V)
    r = state(orbit).r + eps * V[1:3]
    v = state(orbit).v + eps * V[4:6]

    cart = CartesianState(vcat(r,v); lengthunit=lengthunit(orbit), timeunit=timeunit(orbit), angularunit=angularunit(orbit))

    return Orbit(cart, system(orbit), epoch(orbit))
end

#=
"""
Given a __periodic__ `CircularRestrictedThreeBodyOrbit`,
returns a collection of `Trajectory` instances representing
the stable manifold near orbit.
"""
function stable_manifold(orbit::NormalizedSynodicCR3BPOrbit, T::Real;
                         duration = T, reltol = 1e-14,
                         abstol = 1e-14, eps = 1e-8,
                         saveat = 1e-2, num_trajectories = 10, parallel = false)

    @assert duration > 0 "The provided duration cannot be zero or negative. Duration provided was $duration"
    @assert isperiodic(orbit, T) "The provided orbit is not periodic!"

    M = monodromy(orbit, T)
    V = stable_eigenvector(M)

    orbit = NormalizedSynodicSTMCR3BPOrbit(orbit)
    traj = propagate(orbit, T; reltol = reltol, abstol = abstol, saveat = saveat)
    if length(traj) ÷ num_trajectories < 1
        ind = 1:length(traj)
    else
        ind  = 1 : length(traj) ÷ num_trajectories : length(traj) - (length(traj) % num_trajectories)
    end

    f = step -> reverse(propagate(manifold(step, step.state.stm * V; eps = eps), -duration;
                                  reltol = reltol, abstol = abstol))
    return parallel ? pmap(f, traj[ind]) : map(f, traj[ind])
end

"""
Given a __periodic__ `CircularRestrictedThreeBodyOrbit`,
returns a collection of `Trajectory` instances representing
the unstable manifold near orbit.
"""
function unstable_manifold(orbit::NormalizedSynodicCR3BPOrbit, T::Real;
                           duration = T, reltol = 1e-14,
                           abstol = 1e-14, eps = 1e-8,
                           num_trajectories = 10,
                           saveat = 1e-2, parallel = false)

    @assert duration > zero(duration) "The provided duration cannot be zero or negative."
    @assert isperiodic(orbit, T) "The provided orbit is not periodic!"

    orbit = NormalizedSynodicSTMCR3BPOrbit(orbit)
    M = monodromy(orbit, T)
    V = unstable_eigenvector(M)

    traj = propagate(orbit, T; reltol = reltol, abstol = abstol, saveat = saveat)
    if length(traj) ÷ num_trajectories < 1
        ind = 1:length(traj)
    else
        ind  = 1 : length(traj) ÷ num_trajectories : length(traj) - (length(traj) % num_trajectories)
    end

    f = step -> propagate(manifold(step, step.state.stm * V; eps = eps), duration;
                          reltol = reltol, abstol = abstol, saveat = saveat)
    return parallel ? pmap(f, traj[ind]) : map(f, traj[ind])
end
=#


text:
export iauEe00b
"""
**  Equation of the equinoxes, compatible with IAU 2000 resolutions but
**  using the truncated nutation model IAU 2000B.
**
**  This function is part of the International Astronomical Union's
**  SOFA (Standards Of Fundamental Astronomy) software collection.
**
**  Status:  support function.
**
**  Given:
**     date1,date2  double    TT as a 2-part Julian Date (Note 1)
**
**  Returned (function value):
**                  double    equation of the equinoxes (Note 2)
**
**  Notes:
**
**  1) The TT date date1+date2 is a Julian Date, apportioned in any
**     convenient way between the two arguments.  For example,
**     JD(TT)=2450123.7 could be expressed in any of these ways,
**     among others:
**
**            date1          date2
**
**         2450123.7           0.0       (JD method)
**         2451545.0       -1421.3       (J2000 method)
**         2400000.5       50123.2       (MJD method)
**         2450123.5           0.2       (date & time method)
**
**     The JD method is the most natural and convenient to use in
**     cases where the loss of several decimal digits of resolution
**     is acceptable.  The J2000 method is best matched to the way
**     the argument is handled internally and will deliver the
**     optimum resolution.  The MJD method and the date & time methods
**     are both good compromises between resolution and convenience.
**
**  2) The result, which is in radians, operates in the following sense:
**
**        Greenwich apparent ST = GMST + equation of the equinoxes
**
**  3) The result is compatible with the IAU 2000 resolutions except
**     that accuracy has been compromised for the sake of speed.  For
**     further details, see McCarthy & Luzum (2001), IERS Conventions
**     2003 and Capitaine et al. (2003).
**
**  Called:
**     iauPr00      IAU 2000 precession adjustments
**     iauObl80     mean obliquity, IAU 1980
**     iauNut00b    nutation, IAU 2000B
**     iauEe00      equation of the equinoxes, IAU 2000
**
**  References:
**
**     Capitaine, N., Wallace, P.T. and McCarthy, D.D., "Expressions to
**     implement the IAU 2000 definition of UT1", Astronomy &
**     Astrophysics, 406, 1135-1149 (2003)
**
**     McCarthy, D.D. & Luzum, B.J., "An abridged model of the
**     precession-nutation of the celestial pole", Celestial Mechanics &
**     Dynamical Astronomy, 85, 37-49 (2003)
**
**     McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
**     IERS Technical Note No. 32, BKG (2004)
**
**  This revision:  2008 May 18
**
**  SOFA release 2018-01-30
**
**  Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# double iauEe00b(double date1, double date2)

function iauEe00b(date1::Real, date2::Real)
   return ccall((:iauEe00b, libsofa_c), Cdouble,
               (Cdouble, Cdouble),
               convert(Float64, date1), convert(Float64, date2))
end


text:
using Gadgetron
using Test

@testset  "Meta" begin

	meta = MRD.MetaDict()
	meta["Pie"] = 5

	@test meta["Pie"] == 5

	meta["Penguin"] = [1,2,"Pie"]

	@test meta["Penguin"] == [1,2,"Pie"]

	MRD.push!(meta, "Pie", 6, 7)
	@test meta["Pie"] == [5,6,7]
end

@testset "MetaSerialization" begin

	xmlstring = b"<?xml version=\"1.0\"?>
<ismrmrdMeta>
  <meta>
    <name>pi</name>
    <value>3.14159265</value>
  </meta>
  <meta>
    <name>extra</name>
    <value>Hello, World!</value>
  </meta>
  <meta>
    <name>extra</name>
    <value>654321</value>
  </meta>
  <meta>
    <name>extra</name>
    <value>1.234</value>
    <value>67890</value>
    <value>foobar</value>
  </meta>
  <meta>
    <name>when</name>
    <value>2015-03-21, Sat March 2015</value>
    <value>1426933800</value>
  </meta>
</ismrmrdMeta>"
	buffer = IOBuffer()

  MRD.write_string(buffer,String(xmlstring),UInt64)
  Base.seekstart(buffer)

	meta = MRD.read(buffer,MRD.MetaDict)
	@test meta["pi"] == 3.14159265

  @test meta["extra"][2] == 67890

  @test isempty(meta) == false

  buffer = IOBuffer()
  MRD.write(buffer,meta)
  seekstart(buffer)
  meta2 = MRD.read(buffer,MRD.MetaDict)
  @test meta == meta2

end


text:
using MLDatasets: MNIST
using Plots
using Flux

xtrain, ytrain = MNIST.traindata();
xtest, ytest = MNIST.testdata();

ytrain_OneHot = Flux.onehotbatch(ytrain,0:9);
ytest_OneHot = Flux.onehotbatch(ytest,0:9);

xtrain_vec = zeros(28*28,60000)
xtest_vec = zeros(28*28,10000)

for i in 1:size(xtrain)[3]
    xtrain_vec[:,i] = vec(xtrain[:,:,i])
end

for i in 1:size(xtest)[3]
    xtest_vec[:,i] = vec(xtest[:,:,i])
end

xtrain_vec

flux_model = Chain(Dense(28*28, 30, sigmoid), Dense(30,10, sigmoid))

flux_loss(x,y) = Flux.mse(flux_model(x),y)

ps = params(flux_model)
dataset = [(xtrain_vec, ytrain_OneHot)];
opt = Descent(3.0)

cb = () -> println(flux_loss(xtrain_vec,ytrain_OneHot))

Flux.@epochs 30 Flux.train!(flux_loss, ps, dataset, opt, cb = cb)

acc = 0
for i in 1:10000
    x = Flux.onecold(flux_model(xtest_vec[:,1])).-1
    y = ytest[i]
    if x[1] ≈ y
        acc += 1
    end
end
acc /= 10000
flux_model(xtest_vec[:,2])


text:
using Spectra
using Test

# the x axis
x = collect(0:0.1:150)

# a scale factor that you can change
scale = 0.1

# a fake signal: perfect y
y_tot, y_peaks = gaussiennes([10;20.],[40.;60],[5.;15],x)
bkg, ~ = gaussiennes([10.],[50.],[50.],x)
y_perfect = scale.*(y_tot.+bkg)

# we add noise: observed y
y = scale.*(y_tot .+bkg .+ randn(size(y_tot,1)))

roi = [0 20.; 100. 150]

y1, bas1 = baseline(x,y,roi,"gcvspline",s=0.6)

rmse = sqrt.(sum((bas1 .- bkg).^2))./sum(bkg)

@test rmse <0.05 # error inferior to 5%, should actually be arround 0.026.


text:
"""
    Wrapper to `PetscInitializeNoPointers`. Initialize PETCs with arguments.

# Implementation
I don't know if I am supposed to use PetscInt or not...
"""
function PetscInitialize(args::Vector{String}, filename::String, help::String)
    args2 = ["julia"; args]
    nargs = Cint(length(args2))
    error = ccall( (:PetscInitializeNoPointers, libpetsc),
            PetscErrorCode,
            (Cint,
            Ptr{Ptr{UInt8}},
            Cstring,
            Cstring),
            nargs, args2, filename, help
    )
    @assert iszero(error)
end

"""
    Initialize PETSc with arguments stored in a vector of string
"""
PetscInitialize(args::Vector{String}) = PetscInitialize(args, "", "")

"""
    Initialize PETSc with arguments concatenated in a unique string.
"""
PetscInitialize(args::String) = PetscInitialize(convert(Vector{String}, split(args)), "", "")

"""
    PetscInitialize(cmd_line_args::Bool = true)

Initialize PETSc.

If `cmd_line_args == true`, then command line arguments passed to Julia are used as
arguments for PETSc (leading to a call to `PetscInitializeNoPointers`).

Otherwise, if `cmd_line_args == false`, initialize PETSc without arguments (leading
to a call to `PetscInitializeNoArguments`).
"""
function PetscInitialize(cmd_line_args::Bool = true)
    if (cmd_line_args)
        PetscInitialize(ARGS)
    else
        error = ccall((:PetscInitializeNoArguments, libpetsc), PetscErrorCode, ())
        @assert iszero(error)
    end
end

"""
    Wrapper to PetscFinalize
"""
function PetscFinalize()
    error = ccall( (:PetscFinalize, libpetsc), PetscErrorCode, ())
    @assert iszero(error)
end


text:
@testset "Parsimonious flux balance analysis with StandardModel" begin
    model = load_model(StandardModel, model_paths["e_coli_core.json"])
    d = parsimonious_flux_balance_analysis_dict(
        model,
        Tulip.Optimizer;
        modifications = [
            change_constraint("EX_glc__D_e", -12, -12),
            change_optimizer_attribute("IPM_IterationsLimit", 500),
        ],
        qp_modifications = [
            change_optimizer(OSQP.Optimizer),
            change_optimizer_attribute("verbose", false),
        ],
    )
    v = parsimonious_flux_balance_analysis_vec(
        model,
        Tulip.Optimizer;
        modifications = [
            change_constraint("EX_glc__D_e", -12, -12),
            change_optimizer_attribute("IPM_IterationsLimit", 500),
        ],
        qp_modifications = [
            change_optimizer(OSQP.Optimizer),
            change_optimizer_attribute("verbose", false),
        ],
    )

    # The used optimizer doesn't really converge to the same answer everytime
    # here, we therefore tolerate a wide range of results.
    @test isapprox(d["PGM"], -17.568590034769613, atol = 0.5)
    @test isapprox(v[8], -17.568590034769613, atol = 0.5)
end


text:
open("/dev/tape", "w") do f
    write(f, "Hello tape!")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
typealias NumericalOrCategoricalAesthetic
    Union(Nothing, Vector, DataArray, PooledDataArray)

typealias CategoricalAesthetic
    Union(Nothing, PooledDataArray)

typealias NumericalAesthetic
    Union(Nothing, Matrix, Vector, DataArray)

@varset Aesthetics begin
    x,            Union(NumericalOrCategoricalAesthetic, Distribution)
    y,            Union(NumericalOrCategoricalAesthetic, Distribution)
    z,            Union(Nothing, Function, Matrix)
    size,         Maybe(Vector{Measure})
    color,        Maybe(AbstractDataVector{RGB{Float32}})
    label,        CategoricalAesthetic
    group,        CategoricalAesthetic

    xmin,         NumericalAesthetic
    xmax,         NumericalAesthetic
    ymin,         NumericalAesthetic
    ymax,         NumericalAesthetic

    # hexagon sizes used for hexbin
    xsize,        NumericalAesthetic
    ysize,        NumericalAesthetic

    # fixed lines
    xintercept,   NumericalAesthetic
    yintercept,   NumericalAesthetic

    # boxplots
    middle,       NumericalAesthetic
    lower_hinge,  NumericalAesthetic
    upper_hinge,  NumericalAesthetic
    lower_fence,  NumericalAesthetic
    upper_fence,  NumericalAesthetic
    outliers,     NumericalAesthetic

    # subplots
    xgroup,       CategoricalAesthetic
    ygroup,       CategoricalAesthetic

    # guides
    xtick,        NumericalAesthetic
    ytick,        NumericalAesthetic
    xgrid,        NumericalAesthetic
    ygrid,        NumericalAesthetic
    color_key_colors,     Maybe(Associative)
    color_key_title,      Maybe(String)
    color_key_continuous, Maybe(Bool)
    color_function,       Maybe(Function)
    titles,               Maybe(Dict{Symbol, String})

    # mark some ticks as initially invisible
    xtickvisible,         Maybe(Vector{Bool})
    ytickvisible,         Maybe(Vector{Bool})

    # scale at which ticks should become visible
    xtickscale,           Maybe(Vector{Float64})
    ytickscale,           Maybe(Vector{Float64})

    # pesudo-aesthetics used to indicate that drawing might
    # occur beyond any x/y value.
    xdrawmin,     Any
    xdrawmax,     Any
    ydrawmin,     Any
    ydrawmax,     Any

    # plot viewport extents
    xviewmin,     Any
    xviewmax,     Any
    yviewmin,     Any
    yviewmax,     Any

    # labeling functions
    x_label,      Function, showoff
    y_label,      Function, showoff
    xtick_label,  Function, showoff
    ytick_label,  Function, showoff
    color_label,  Function, showoff
    xgroup_label, Function, showoff
    ygroup_label, Function, showoff
end

# Alternate aesthetic names
const aesthetic_aliases =
    [:colour        => :color,
     :x_min         => :xmin,
     :x_max         => :xmax,
     :y_min         => :ymin,
     :y_max         => :ymax,
     :x_group       => :xgroup,
     :y_group       => :ygroup,
     :x_viewmin     => :xviewmin,
     :x_viewmax     => :xviewmax,
     :y_viewmin     => :yviewmin,
     :y_viewmax     => :yviewmax,
     :x_group_label => :xgroup_label,
     :y_group_label => :ygroup_label,
     :x_tick        => :xtick,
     :y_tick        => :ytick,
     :x_grid        => :xgrid,
     :y_grid        => :ygrid]

# Index as if this were a data frame
function getindex(aes::Aesthetics, i::Integer, j::String)
    getfield(aes, symbol(j))[i]
end

# Return the set of variables that are non-nothing.
function defined_aesthetics(aes::Aesthetics)
    vars = Set{Symbol}()
    for name in Aesthetics.names
        if !is(getfield(aes, name), nothing)
            push!(vars, name)
        end
    end
    vars
end

# Checking aesthetics and giving reasonable error messages.

# Raise an error if any of the given aesthetics are not defined.
#
# Args:
#   who: A string naming the caller which is printed in the error message.
#   aes: An Aesthetics object.
#   vars: Symbol that must be defined in the aesthetics.
#
function undefined_aesthetics(aes::Aesthetics, vars::Symbol...)
    setdiff(Set(vars), defined_aesthetics(aes))
end

function assert_aesthetics_defined(who::String, aes::Aesthetics, vars::Symbol...)
    undefined_vars = undefined_aesthetics(aes, vars...)
    if !isempty(undefined_vars)
        error(@sprintf("The following aesthetics are required by %s but are not defined: %s\n",
                       who, join(undefined_vars, ", ")))
    end
end

function assert_aesthetics_equal_length(who::String, aes::Aesthetics, vars::Symbol...)
    defined_vars = Symbol[]
    for var in filter(var -> !(getfield(aes, var) === nothing), vars)
        push!(defined_vars, var)
    end

    n = length(getfield(aes, vars[1]))
    for i in 2:length(vars)
        if length(getfield(aes, vars[1])) != n
            error(@sprintf("The following aesthetics are required by %s to be of equal length: %s\n",
                           who, join(vars, ", ")))
        end
    end
end

# Create a shallow copy of an Aesthetics instance.
#
# Args:
#   a: aesthetics to copy
#
# Returns:
#   Copied aesthetics.
#
copy(a::Aesthetics) = Aesthetics(a)

# Replace values in a with non-nothing values in b.
#
# Args:
#   a: Destination.
#   b: Source.
#
# Returns: nothing
#
# Modifies: a
#
function update!(a::Aesthetics, b::Aesthetics)
    for name in Aesthetics.names
        if issomething(getfield(b, name))
            setfield(a, name, getfield(b, name))
        end
    end

    nothing
end

# Serialize aesthetics to JSON.

# Args:
#  a: aesthetics to serialize.
#
# Returns:
#   JSON data as a string.
#
function json(a::Aesthetics)
    join([string(a, ":", json(getfield(a, var))) for var in aes_vars], ",\n")
end

# Concatenate aesthetics.
#
# A new Aesthetics instance is produced with data vectors in each of the given
# Aesthetics concatenated, nothing being treated as an empty vector.
#
# Args:
#   aess: One or more aesthetics.
#
# Returns:
#   A new Aesthetics instance with vectors concatenated.
#
function concat(aess::Aesthetics...)
    cataes = Aesthetics()
    for aes in aess
        for var in Aesthetics.names
            setfield!(cataes, var,
                      cat_aes_var!(getfield(cataes, var), getfield(aes, var)))
        end
    end
    cataes
end

cat_aes_var!(a::Nothing, b::Nothing) = a
cat_aes_var!(a::Nothing, b) = copy(b)
cat_aes_var!(a, b::Nothing) = a
cat_aes_var!(a::Function, b::Function) = a === string || a == showoff ? b : a

function cat_aes_var!(a::Dict, b::Dict)
    merge!(a, b)
    a
end

function cat_aes_var!{T <: Base.Callable}(a::AbstractArray{T}, b::AbstractArray{T})
    return append!(a, b)
end

function cat_aes_var!{T <: Base.Callable, U <: Base.Callable}(a::AbstractArray{T}, b::AbstractArray{U})
    return append!(a, b)
end

# Let arrays of numbers clobber arrays of functions. This is slightly odd
# behavior, comes up with with function statistics applied on a layer-wise
# basis.
function cat_aes_var!{T <: Base.Callable, U}(a::AbstractArray{T}, b::AbstractArray{U})
    return b
end

function cat_aes_var!{T, U <: Base.Callable}(a::AbstractArray{T}, b::AbstractArray{U})
    return a
end

function cat_aes_var!{T}(a::AbstractArray{T}, b::AbstractArray{T})
    return append!(a, b)
end

function cat_aes_var!{T, U}(a::AbstractArray{T}, b::AbstractArray{U})
    V = promote_type(T, U)
    if isa(a, DataArray) || isa(b, DataArray)
        ab = DataArray(V, length(a) + length(b))
    else
        ab = Array(V, length(a), length(b))
    end
    i = 1
    for x in a
        ab[i] = x
        i += 1
    end
    for x in b
        ab[i] = x
        i += 1
    end

    return ab
end

function cat_aes_var!(a, b)
    a
end

function cat_aes_var!{T}(xs::PooledDataVector{T}, ys::PooledDataVector{T})
    newpool = T[x for x in union(Set(xs.pool), Set(ys.pool))]
    newdata = vcat(T[x for x in xs], T[y for y in ys])
    PooledDataArray(newdata, newpool, [false for _ in newdata])
end

# Summarizing aesthetics

# Produce a matrix of Aesthetic or Data objects partitioning the original
# Aesthetics or Data object by the cartesian product of xgroup and ygroup.
#
# This is useful primarily for drawing facets and subplots.
#
# Args:
#   aes: Aesthetics or Data objects to partition.
#
# Returns:
#   A Array{Aesthetics} of size max(1, length(xgroup)) by
#   max(1, length(ygroup))
#
function by_xy_group{T <: Union(Data, Aesthetics)}(aes::T, xgroup, ygroup)
    @assert !is(xgroup, nothing) || !is(ygroup, nothing)
    @assert xgroup === nothing || ygroup === nothing ||
            length(xgroup) == length(ygroup)

    n = ygroup === nothing ? 1 : maximum(ygroup)
    m = xgroup === nothing ? 1 : maximum(xgroup)

    xrefs = xgroup === nothing ? [1] : xgroup
    yrefs = ygroup === nothing ? [1] : ygroup

    aes_grid = Array(T, n, m)
    staging = Array(DataArray, n, m)
    for i in 1:n, j in 1:m
        aes_grid[i, j] = T()
    end

    function make_pooled_data_array{T, U, V}(::Type{PooledDataArray{T,U,V}},
                                          arr::AbstractArray)
        PooledDataArray(convert(Array{T}, arr))
    end

    for var in T.names
        # Skipped aesthetics. Don't try to partition aesthetics for which it
        # makes no sense to pass on to subplots.
        if var == :xgroup || var == :ygroup||
           var == :xtick || var == :ytick ||
           var == :xgrid || var == :ygrid ||
           var == :x_drawmin || var == :y_drawmin ||
           var == :x_drawmax || var == :y_drawmax ||
           var == :x_viewmin || var == :y_viewmin ||
           var == :x_viewmax || var == :y_viewmax ||
           var == :color_key_colors
            continue
        end

        vals = getfield(aes, var)
        if typeof(vals) <: AbstractArray
            if !is(xgroup, nothing) && length(vals) != length(xgroup) ||
               !is(ygroup, nothing) && length(vals) != length(ygroup)
                error("Aesthetic $(var) must be the same length as xgroup or ygroup")
            end

            for i in 1:n, j in 1:m
                staging[i, j] = DataArray(eltype(vals), 0)
            end

            for (i, j, v) in zip(Iterators.cycle(yrefs), Iterators.cycle(xrefs), vals)
                push!(staging[i, j], v)
            end

            for i in 1:n, j in 1:m
                if typeof(vals) <: PooledDataArray
                    setfield!(aes_grid[i, j], var,
                              make_pooled_data_array(typeof(vals), staging[i, j]))
                else
                    if !applicable(convert, typeof(vals), staging[i, j])
                        T = eltype(vals)
                        if T <: ColorValue T = ColorValue end
                        da = DataArray(T, length(staging[i, j]))
                        copy!(da, staging[i, j])
                        setfield!(aes_grid[i, j], var, da)
                    else
                        setfield!(aes_grid[i, j], var,
                                  convert(typeof(vals), copy(staging[i, j])))
                    end
                end
            end
        else
            for i in 1:n, j in 1:m
                setfield!(aes_grid[i, j], var, vals)
            end
        end
    end

    aes_grid
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using LinearAlgebra, Random
using Distributions, PDMats, StatsBase, ValueShapes, ArraysOfArrays

@testset "initvals" begin
    @testset "bat_initval" begin
        du = MvNormal(Diagonal(fill(1.0, 2)))
        @test @inferred(bat_default(bat_initval, Val(:algorithm), du)) == InitFromTarget()
        @test @inferred(bat_initval(du)).result isa Vector{<:AbstractFloat}
        @test bat_initval(du).optargs.algorithm == InitFromTarget()
        @test pdf(du, bat_initval(du).result) > 0
        @test @inferred(bat_initval(du, 10)).result isa VectorOfSimilarVectors{<:AbstractFloat}
        @test all(x -> x > 0, pdf.(Ref(du), bat_initval(du, 10).result))

        ds = NamedTupleDist(a = Normal(), b = Exponential(), c = Uniform(-1, 2))
        @test @inferred(bat_default(bat_initval, Val(:algorithm), du)) == InitFromTarget()
        @test @inferred(bat_initval(ds)).result isa ShapedAsNT
        @test bat_initval(ds).optargs.algorithm == InitFromTarget()
        @test pdf(ds, bat_initval(ds).result) > 0
        @test @inferred(bat_initval(ds, 10)).result isa ShapedAsNTArray
        @test all(x -> x > 0, pdf.(Ref(ds), bat_initval(ds, 10).result))
    end
end


text:
module ReactiveMPModelsLGSSMTest

using Test, InteractiveUtils
using Rocket, ReactiveMP, GraphPPL, Distributions
using BenchmarkTools, Random, Plots, Dates, LinearAlgebra, StableRNGs

## Model definition
## -------------------------------------------- ##
@model function univariate_lgssm_model(n, x0, c_, P_)

    x_prior ~ NormalMeanVariance(mean(x0), cov(x0))

    x = randomvar(n)
    c = constvar(c_)
    P = constvar(P_)
    y = datavar(Float64, n)

    x_prev = x_prior

    for i in 1:n
        x[i] ~ x_prev + c
        y[i] ~ NormalMeanVariance(x[i], P)
        x_prev = x[i]
    end

    return x, y
end
## -------------------------------------------- ##
@model function multivariate_lgssm_model(n, x0, A, B, Q, P)

    # We create constvar references for better efficiency
    cA = constvar(A)
    cB = constvar(B)
    cQ = constvar(Q)
    cP = constvar(P)

    # `x` is a sequence of hidden states
    x = randomvar(n)
    # `y` is a sequence of "clamped" observations
    y = datavar(Vector{Float64}, n)

    x_prior ~ MvNormalMeanCovariance(mean(x0), cov(x0))
    x_prev = x_prior

    for i in 1:n
        x[i] ~ MvNormalMeanCovariance(cA * x_prev, cQ)
        y[i] ~ MvNormalMeanCovariance(cB * x[i], cP)
        x_prev = x[i]
    end

    return x, y
end
## -------------------------------------------- ##
## Inference definition
## -------------------------------------------- ##
function univariate_lgssm_inference(data, x0, c, P)
    n = length(data)

    model, (x, y) = univariate_lgssm_model(n, x0, c, P);

    x_buffer = buffer(Marginal, n)
    fe       = keep(Float64)

    x_sub = subscribe!(getmarginals(x), x_buffer)
    f_sub = subscribe!(score(Float64, BetheFreeEnergy(), model), fe)

    update!(y, data)

    unsubscribe!((x_sub, f_sub))

    return x_buffer, fe
end
## -------------------------------------------- ##
function multivariate_lgssm_inference(data, x0, A, B, Q, P)
    n = length(data)

    model, (x, y) = multivariate_lgssm_model(n, x0, A, B, Q, P, options = (limit_stack_depth = 100, ));

    xbuffer = buffer(Marginal, n)
    fe      = keep(Float64)

    x_sub = subscribe!(getmarginals(x), xbuffer)
    f_sub  = subscribe!(score(Float64, BetheFreeEnergy(), model), fe)

    update!(y, data)

    unsubscribe!((x_sub, f_sub))

    return xbuffer, fe
end
## -------------------------------------------- ##

@testset "Linear Gaussian State Space Model" begin

    @testset "Univariate" begin
        ## -------------------------------------------- ##
        ## Data creation
        ## -------------------------------------------- ##
        rng = StableRNG(123)
        P   = 100.0
        n   = 500
        hidden   = collect(1:n)
        data     = hidden + rand(rng, Normal(0.0, sqrt(P)), n);
        x0_prior = NormalMeanVariance(0.0, 10000.0)
        ## -------------------------------------------- ##
        ## Inference execution
        x_estimated, fe = univariate_lgssm_inference(data, x0_prior, 1.0, P);
        ## -------------------------------------------- ##
        ## Test inference results
        @test length(x_estimated) === n
        @test all((mean.(x_estimated) .- 3 .* std.(x_estimated)) .< hidden .< (mean.(x_estimated) .+ 3 .* std.(x_estimated)))
        @test all(var.(x_estimated) .> 0.0)
        @test length(fe) === 1
        @test abs(last(fe) - 1854.297647) < 0.01
        ## -------------------------------------------- ##
        ## Form debug output
        base_output = joinpath(pwd(), "_output", "models")
        mkpath(base_output)
        timestamp        = Dates.format(now(), "dd-mm-yyyy-HH-MM")
        plot_output      = joinpath(base_output, "lgssm_univariate_plot_$(timestamp)_v$(VERSION).png")
        benchmark_output = joinpath(base_output, "lgssm_univariate_benchmark_$(timestamp)_v$(VERSION).txt")
        ## -------------------------------------------- ##
        ## Create output plots
        subrange = 200:215
        m = mean.(x_estimated)[subrange]
        s = std.(x_estimated)[subrange]
        p = plot()
        p = plot!(subrange, m, ribbon = s, label = "Estimated signal")
        p = plot!(subrange, hidden[subrange], label = "Hidden signal")
        p = scatter!(subrange, data[subrange], label = "Observations")
        savefig(p, plot_output)
        ## -------------------------------------------- ##
        ## Create output benchmarks
        benchmark = @benchmark univariate_lgssm_inference($data, $x0_prior, 1.0, $P) seconds=15
        open(benchmark_output, "w") do io
            show(io, MIME("text/plain"), benchmark)
            versioninfo(io)
        end
        ## -------------------------------------------- ##
    end

    @testset "Multivariate" begin
        ## -------------------------------------------- ##
        ## Data creation
        ## -------------------------------------------- ##
        function generate_data(rng, A, B, Q, P)
            x_prev = [ 10.0, -10.0 ]

            x = Vector{Vector{Float64}}(undef, n)
            y = Vector{Vector{Float64}}(undef, n)

            for i in 1:n
                x[i] = rand(rng, MvNormal(A * x_prev, Q))
                y[i] = rand(rng, MvNormal(B * x[i], P))
                x_prev = x[i]
            end

            return x, y
        end
        ## -------------------------------------------- ##
        # Seed for reproducibility
        rng = StableRNG(1234)
        # We will model 2-dimensional observations with rotation matrix `A`
        # To avoid clutter we also assume that matrices `A`, `B`, `P` and `Q`
        # are known and fixed for all time-steps
        θ = π / 35
        A = [ cos(θ) -sin(θ); sin(θ) cos(θ) ]
        B = diageye(2)
        Q = diageye(2)
        P = 25.0 .* diageye(2)
        n = 1_000 # Number of observations
        x, y = generate_data(rng, A, B, Q, P)
        x0 = MvNormalMeanCovariance(zeros(2), 100.0 * diageye(2));
        ## -------------------------------------------- ##
        ## Inference execution
        xmarginals, fe = multivariate_lgssm_inference(y, x0, A, B, Q, P);
        ## Test inference results
        @test length(xmarginals) === n
        # We use 3.0var instead of 3.0std here for easier dot broadcasting with mean
        @test all((mean.(xmarginals) .- 3.0 .* var.(xmarginals)) .< x .< (mean.(xmarginals) .+ 3.0 .* var.(xmarginals)))
        @test all(isposdef.(cov.(xmarginals)))
        @test length(fe) === 1
        @test abs(last(fe) - 6275.9015944677) < 0.01
        ## -------------------------------------------- ##
        ## Form debug output
        base_output = joinpath(pwd(), "_output", "models")
        mkpath(base_output)
        timestamp        = Dates.format(now(), "dd-mm-yyyy-HH-MM")
        plot_output      = joinpath(base_output, "lgssm_multivariate_plot_$(timestamp)_v$(VERSION).png")
        benchmark_output = joinpath(base_output, "lgssm_multivariate_benchmark_$(timestamp)_v$(VERSION).txt")
        ## -------------------------------------------- ##
        ## Create output plots
        slicedim(dim) = (a) -> map(e -> e[dim], a)

        subrange = 100:500
        px = plot()

        px = plot!(px, x[subrange] |> slicedim(1), label = "Hidden Signal (dim-1)", color = :orange)
        px = plot!(px, x[subrange] |> slicedim(2), label = "Hidden Signal (dim-2)", color = :green)

        px = plot!(px, mean.(xmarginals)[subrange] |> slicedim(1), ribbon = var.(xmarginals)[subrange] |> slicedim(1) .|> sqrt, fillalpha = 0.5, label = "Estimated Signal (dim-1)", color = :teal)
        px = plot!(px, mean.(xmarginals)[subrange] |> slicedim(2), ribbon = var.(xmarginals)[subrange] |> slicedim(2) .|> sqrt, fillalpha = 0.5, label = "Estimated Signal (dim-1)", color = :violet)

        savefig(px, plot_output)
        ## -------------------------------------------- ##
        ## Create output benchmarks
        benchmark = @benchmark multivariate_lgssm_inference($y, $x0, $A, $B, $Q, $P) seconds=15
        open(benchmark_output, "w") do io
            show(io, MIME("text/plain"), benchmark)
            versioninfo(io)
        end
        ## -------------------------------------------- ##
    end

end

end


text:
using IntrinsicPersistence
using Test

using StaticArrays
using LinearAlgebra

# sanitycheck
@testset "Sanity checks" begin
    @testset "regular circle" begin
        n = 12
        pts = [SVector(cos(t), sin(t)) for t in range(0, stop = 2π, length = n+1)[1:end-1]]
        landmarks = 1:2:n

        gc = GeodesicComplex(pts, 0.6, landmarks = landmarks)
        res = persistence(gc)
        @test length(res) == 1
        α = first(res).indices
        @test sort(α) == landmarks
    end

    @testset "circle with contraction" begin
        n = 24
        pts = vcat(
            [SVector(1.0cos(t), 1.0sin(t)) for t in range(0, stop = 2π, length = n+1)[1:end-1]],
            [SVector(0.9cos(t), 0.9sin(t)) for t in range(0, stop = 2π, length = n+1)[1:end-1]],
            [SVector(0.8cos(t), 0.8sin(t)) for t in range(0, stop = 2π, length = n+1)[1:end-1]],
            [SVector(0.7cos(t), 0.7sin(t)) for t in range(0, stop = 2π, length = n+1)[1:end-1]],
            [SVector(0.6cos(t), 0.6sin(t)) for t in range(0, stop = 2π, length = n+1)[1:end-1]],
            [SVector(0.5cos(t), 0.5sin(t)) for t in range(0, stop = 2π, length = n+1)[1:end-1]],
        )
        landmarks = 1:2:n
        gc = GeodesicComplex(pts, 0.3, landmarks = landmarks)
        res = persistence(gc)
        @test length(res) == 1
        @test diameter(first(res)) ≈ 1
    end

    @testset "landmarks form a circle, but the underlying space is contractible" begin
        n = 12
        pts = vcat([SVector(0.5cos(t), 0.5sin(t))
                    for t in range(0, stop = 2π, length = n÷2+1)[1:end-1]],
                   [SVector(cos(t), sin(t))
                    for t in range(0, stop = 2π, length = n+1)[1:end-1]],
                   [SVector(0.0, 0.0)])
        landmarks = n÷2+1:2:n+n÷2
        gc = GeodesicComplex(pts, 0.6, landmarks = landmarks)
        res = persistence(gc)
        @test isempty(res)

        res = persistence(gc, shrinked = false)
        @test !isempty(res)
    end

    @testset "cylinder" begin
        n = 12
        pts = [SVector(cos(t), sin(t), h)
               for t in range(0, stop = 2π, length = n+1)[1:end-1]
               for h in 0:0.5:10]
        landmarks = 1:2:length(pts)
        gc = GeodesicComplex(pts, 0.6, landmarks = landmarks)
        res = persistence(gc)
        @test length(res) == 1
    end

    @testset "random full graph" begin
        n = 100
        pts = [SVector(rand(), rand(), rand()) for _ in 1:n]
        landmarks = 1:n
        gc = GeodesicComplex(pts, 2, landmarks = landmarks)
        res = persistence(gc)
        @test isempty(res)
    end

    @testset "regular cone" begin
        n = 12
        pts = [SVector(h/10*cos(t), h/10*sin(t), h)
               for t in range(0, stop = 2π, length = n+1)[1:end-1]
               for h in 0:0.5:10]
        landmarks = 1:2:length(pts)
        gc = GeodesicComplex(pts, 0.6, landmarks = landmarks)
        res = persistence(gc)
        @test isempty(res)
    end

    @testset "random cone" begin
        n = 100_000
        pts = map(rand() for _ in 1:n) do h
            t = rand()
            SVector(h*sinpi(2t), h*cospi(2t), h)
        end
        gc = GeodesicComplex(pts, 0.2)
        res = persistence(gc)
        @test isempty(res)
    end

    @testset "random truncated cone" begin
        n = 100_000
        pts = map(rand() for _ in 1:n) do h
            t = rand()
            SVector((h+1)*sinpi(2t), (h+1)*cospi(2t), h)
        end
        gc = GeodesicComplex(pts, 0.2)
        res = persistence(gc)
        @test length(res) == 1
        @test diameter(first(res)) ≈ 2 atol = 0.2
        @test perimeter(first(res)) ≈ 2π atol = 0.2
    end
end

max_stars_repo_name: mtsch/IntrinsicPersistence

max_issues_repo_name: mtsch/IntrinsicPersistence

max_forks_repo_name: mtsch/IntrinsicPersistence