text:
# `read_file` will read the entire file as a single string. From this string we can
# start parsing the model description
function read_file(filename::String)
    file = open(filename)
    contents = readall(file)
    close(file)
    return contents
end

# Look a pattern consisting of a `#` followed by 0 or more character of any type except newline
# and ending in newline
strip_comments(source::String) = replace(source, r"\#.*\n"," ")

# Keywords that id the type of sentence.
# Note that the "export" version of a start keyword must go first in order for the
# store_logical_lines function to work
const start_keywords = ["qsspecies", "qsstate", "species", "state", "reaction",
                       "variable", "derivative", "observed", "parameter", "constant",
                       "forcing", "component", "function", "compartment"]

# Tokenize the string by whitespace, tab or newline
tokenize(source::String) = split(source, r"\s|\n|\t")

# Separate into logical lines knowing that a valid sentence must begin with a start keyword
# Whitespace is ignored. To make sure there are no problems, the reconstructed sentences are always
# appended surrounded by whitespace
function logical_lines(source::Array{SubString{ASCIIString},1})
 lines = ASCIIString[]
 c = 0
 for i in source
   if i == ""
     continue
   elseif in(i, start_keywords)
     c += 1
     push!(lines, i)
   else
     lines[c] *= " $i "
   end
 end
 return lines
end

# Store the lines in a dictionary where the key indicates the type of sentence
function store_logical_lines(source::Array{ASCIIString, 1})
 lines = Dict{ASCIIString, Array{ASCIIString, 1}}()
 [lines[i] = ASCIIString[] for i = start_keywords]
 for i in source
   for j in start_keywords
     if contains(i, j)
       push!(lines[j], i)
       break
     end
   end
 end
 return lines
end

#### Parsers for each type of sentence ####

## Parse species
# species <name> in <comparment> = <value> in <units>
# species <name> in <comparment> = <value>
# species <name> in <compartment>
# species <name>
function parse_species(s::ASCIIString)
  sentence = tokenize(strip(s))
  c = 1
  name = ""
  value = 0.0
  compartment = "__default_comp__"
  unit = ""
  comp_given = false
  # Process non-blank tokens
  for i in sentence
    i == "" && continue
    # The first token should be the keyword
    if c == 1
      c = 2
    # The second token should be the name of the parameter
    elseif c == 2
      name = i
      c = 3
    # The third token could be the keyword "in" if there is a compartment
    elseif (c == 3) & (i == "in")
      comp_given = true
    # The fourth token could be the compartment (only if preceded by "in")
    elseif (c == 3) & (comp_given == true) & (i != "=")
      compartment = i
    # The 3rd or 5th token should be the = sign
    elseif (c == 3) & (i == "=")
      c = 4
    # The 4th or 6th token should be the value
    elseif (c == 4) & (i != "in")
      try value = parse(i) catch error("Error in species $name: You must supply a value") end
    # If there is an "in" keyword after the value, it indicates the units
    elseif (c == 4) & (i == "in")
      c = 5
    else
      unit *= i
    end
  end
  return name, dt.Species(value, compartment, dt.Unit(unit), {}, false)
end

## Parse qsspecies
# qsspecies <name> in <comparment> = <value> in <units>
# qsspecies <name> in <comparment> = <value>
# qsspecies <name> in <compartment>
# qsspecies <name>
function parse_qsspecies(s::ASCIIString)
  name, species = parse_species(s)
  species.QSS = true
  return name, species
end

## Parse parameter
# parameter <name> = <value> in <unit>
function parse_parameter(s::ASCIIString)
  sentence = tokenize(strip(s))
  c = 1
  name = ""
  value = 0.0
  unit = ""
  # Process non-blank tokens
  for i in sentence
    i == "" && continue
    # The first token should be the keyword
    if c == 1
      c = 2
    # The second token should be the name of the parameter
    elseif c == 2
      name = i
      c = 3
    # The third token should be the equal sign, otherwise we have an error
    elseif c == 3
      i != "=" ? error("Error in parameter $name: I was expecting an equal = sign.") : (c = 4)
    # The fourth (and posterior) token should be the expression
    elseif (c == 4) & (i != "in")
      try value = parse(i) catch error("Error in parameter $name: You must supply a value") end
    elseif (c == 4) & (i == "in")
      c = 5
    elseif (c == 5)
      unit *= i
    end
  end
  !isa(value, Number) && error("Error in parameter $name: You must supply a value")
  return name, dt.Parameter(value, dt.Unit(unit))
end

## Parse constant
# constant <name> = <value> in <unit>
function parse_constant(s::ASCIIString)
  return parse_parameter(s)
end

## Parse state
# state <name> = <value> in <units>
# state <name> = <value>
# state <name>
function parse_state(s::ASCIIString)
  name, state_as_param =  parse_parameter(s)
  state = dt.State(state_as_param.Value, state_as_param.Units, false)
  return name, state
end

## Parse qsstate
# qsstate <name> = <value> in <units>
# qsstate <name> = <value>
# qsstate <name>
function parse_qsstate(s::ASCIIString)
  name, state_as_param = parse_parameter(s)
  state = dt.State(state_as_param.Value, state_as_param.Units, true)
  return name, state
end

## Parse reaction
# reaction <name> in <compartment> = <expression> for <stoichiometry> in <unit>
# reaction <name> in <compartment> = <expression> for <stoichiometry>
# reaction <name> = <expression> for <stoichiometry> in <unit>
# reaction <name> = <expression> for <stoichiometry>
function parse_reaction(s::ASCIIString)
  sentence = tokenize(strip(s))
  c = 1
  name = ""
  compartment = "__default_comp__"
  comp_given = false
  expr = ""
  unit = ""
  substrates = ""
  products = ""
  # Process non-blank tokens
  for i in sentence
    i == "" && continue
    # The first token should be the keyword
    if c == 1
      c = 2
    # The second token should be the name of the variable
    elseif c == 2
      name = i
      c = 3
    # The third token may be the in keyword
    elseif (c == 3) && (i == "in") && (comp_given == false)
      comp_given = true
    # The fourth token may be the compartment
    elseif (c == 3) && (comp_given == true) && (i != "=")
      compartment = i
    # The 3rd/5th token should be the = sign
    elseif c == 3
      i != "=" ?
          error("Error in reaction $name: I was expecting an equal = sign.") :
          (c = 4)
    # The 4th/6th and posterior tokens until "for" should be the expression
    elseif (c == 4) & (i != "for")
      expr *= i
    # After the mathematical expression the "for" is used to indicate the start of chemical stoichiometry
    elseif (c == 4) & (i == "for")
      c = 5
    # All tokens between "for" and "->" indicate substrates
    elseif (c == 5) & (i != "->")
      substrates *= i
    # The token "->" is used to separate substrates from products
    elseif (c == 5) && (i == "->")
      c = 6
    # All tokens after "->" represent products until the "in" keyword
    elseif (c == 6) & (i != "in")
      products *= i
    elseif (c == 6) & (i == "in")
      c = 7
    elseif (c == 7)
      unit *= i
    end
  end
  # Try to catch situations where the stoichiometry was not given or when the was a parsing error...
  (c == 4) &&
          error("Error in reaction $name: I am missing the stoichiometry of the reaction.")
  (c == 5) &&
          error("Error in reaction $name: I could not parse the stoichiometry correctly.")
  # Make sure that we have a valid mathematical expression. If not raise an error.
  # We don't store the parsed expression as it makes manipulation more difficult, but it is good to check
  # now if we can later parse it
  rhs = try parse(expr) catch
            error("Error in reaction $name: I don't understand the equation: $expr") end
  # Create the arrays of substrates and products (first make sure that all coefficients are provided)
  Substrates = components_in_stoich(substrates)
  Products = components_in_stoich(products)
  # Construct the Reaction object and return it
  return name, dt.Reaction(Substrates, Products, expr, compartment, dt.Unit(unit))
end

# Parse a combination of chemical species with or without pattern-matching
# n species[[component[form]...]]
# n species
function components_in_stoich(reactants::ASCIIString)
    output = dt.Reactant[]
    # Strip by + sign
    chunks = split(reactants, "+")
    (length(chunks[1]) == 0) && return output
    # Iterate over chunks
    for i in chunks
      # temporary strings to store stoichiometry and species description
      stoich_string = ""
      species_string = ""
      parse_stoich = true
      # Iterate over characters
        for j in i
          # If we are parsing stoich and we hit a letter, we are finished with stoich
          parse_stoich && ismatch(r"[a-zA-Z]", string(j)) && (parse_stoich = false)
          # Put the character in the right susbtring depending on whether we are parsing stoich or not
          parse_stoich ? stoich_string *= string(j) : species_string *= string(j)
        end
      push!(output, dt.Reactant(species_string, float(stoich_string)))
    end
    return output
end

## Parse variable
# variable <name> = <expression> in <unit>
function parse_variable(s::ASCIIString)
  sentence = tokenize(strip(s))
  c = 1
  name = ""
  expr = ""
  unit = ""
  # Process non-blank tokens
  for i in sentence
    i == "" && continue
    # The first token should be the keyword
    if c == 1
      c = 2
    # The second token should be the name of the variable
    elseif c == 2
      name = i
      c = 3
    # The third token should be the equal sign, otherwise we have an error
    elseif c == 3
      i != "=" ? error("Error in equation $name: I was expecting an equal = sign.") : (c = 4)
    # The fourth until the in keyword should be the mathematical expression
    elseif (c == 4) & (i != "in")
      expr *= i
    elseif (c == 4) & (i == "in")
      c = 5
    else
      unit *= i
    end
  end
  c == 3 && error("Error in variable $name: You have to specify an expression bound to the variable.")
  # Make sure that we have a valid mathematical expression. If not raise an error.
  # We don't store the parsed expression as it makes manipulation more difficult, but it is good to check
  # now if we can later parse it
  rhs = try parse(expr) catch
            error("Error in variable $name: I don't understand the equation...") end

  return name, dt.Variable(expr, dt.Unit(unit))
end

## Parse observed
# observed <name> = <expression> in <unit>
function parse_observed(s::ASCIIString)
  name, variable = parse_variable(s)
  return name, dt.Observed(variable.Expr, variable.Units)
end

## Parse derivative
# observed <name> of <state> = <expression> in <unit>
# observed <name> of <state> = <expression>
function parse_derivative(s::ASCIIString)
  sentence = tokenize(strip(s))
  c = 1
  name = ""
  state = ""
  expr = ""
  unit = ""
  # Process non-blank tokens
  for i in sentence
    i == "" && continue
    # The first token should be the keyword
    if c == 1
      c = 2
    # The second token should be the name of the variable
    elseif c == 2
      name = i
      c = 3
    elseif c == 3
      i != "of" ? error("Error in equation $name: I was expecting the keyword of to indicate the state the derivative is associated to.") : (c = 4)
# The fifth token should be the equal sign, otherwise we have an error
    elseif c == 4
      state = i
      c = 5
    elseif c == 5
      i != "=" ? error("Error in equation $name: I was expecting an equal = sign.") : (c = 6)
# The sixth until the in keyword should be the mathematical expression
    elseif (c == 6) & (i != "in")
      expr *= i
    elseif (c == 6) & (i == "in")
      c = 7
    else
      unit *= i
    end
  end
  c == 5 && error("Error in derivative $name: You need to supply a function.")

  rhs = try parse(expr) catch
            error("Error in derivative $name: I don't understand the equation...") end

  return name, dt.Derivative(expr, state, dt.Unit(unit))
end

## Parse forcing
# forcing <name> = <values> at <times> in <unit>
# forcing <name> = <values> at <times>
# forcing <name>
function parse_forcing(s::ASCIIString)
  sentence = tokenize(strip(s))
  c = 1
  name = ""
  values = ""
  times = ""
  unit = ""
  # Process non-blank tokens
  for i in sentence
    i == "" && continue
    # The first token should be the keyword
    if c == 1
      c = 2
    # The second token should be the name of the parameter
    elseif c == 2
      name = i
      c = 3
    # The third token should be the equal sign, otherwise we have an error
    elseif c == 3
      i != "=" ? error("Error in forcing $name: I was expecting an equal = sign.") : (c = 4)
    # The tokens between the = sign and "at" keyword should be the values
    elseif (c == 4) & (i != "at")
      values *= i
    elseif (c == 4) & (i == "at")
      c = 5
    elseif (c == 5) & (i != "in")
      times *= i
    elseif (c == 5) & (i == "in")
      c = 6
    else
      unit *= i
    end
  end
  c == 4 && error("Error in forcing $name: I could not parse the times correctly.")
  c == 3 && (values = "[0.0, 0.0]"; times = "[0.0,0.1]")
  values_vec = eval(parse(values))
  times_vec = eval(parse(times))
  length(values_vec) != length(times_vec) && error("Error in forcing $name: Default values and times must have the same length")
  return name, dt.Forcing(times_vec, values_vec, dt.Unit(unit))
end

# Parse component
# component <name> <labels> of <species> = <values>
function parse_component(s::ASCIIString)
  sentence = tokenize(strip(s))
  c = 1
  name = ""
  values = ""
  species = ""
  forms = ""
  # Process non-blank tokens
  for i in sentence
    i == "" && continue
    # The first token should be the keyword
    if c == 1
      c = 2
    # The second token should be the name of the parameter
    elseif c == 2
      name = i
      c = 3
    # The tokens between second and the keyword "of" should be an array of forms
    elseif c == 3
      forms *= i
      ismatch(r"\]" ,i) && (c = 4)
    # The token after the array of forms should be the keyword "of"
    elseif c == 4
      i != "of" ?
          error("Error in component $name: I was expecting the \"of\" keyword to specify the species.") :
          (c = 5)
    # After than we get the species
    elseif c == 5
      species = i
      c = 6
    # Then an = sign
    elseif c == 6
      i != "=" ?
          error("Error in component $name: I was expecting an equal = sign.") :
          (c = 7)
    # The remaining tokens are the values
    elseif (c == 7)
      values *= i
    end
  end
  c == 3 && error("Error in component $name: I could not parse it correctly. Did you keep spaces between the brackets and contents?")
  values_vec = eval(parse(values))
  forms_vec = eval(parse(forms))
  forms_vec = ASCIIString[string(i) for i in forms_vec]
  abs(sum(values_vec) - 1.0) > eps(1.0)  && error("Error in component $name: The sum of default initial values for the different forms must be equal to 1")
  return name, dt.Component(name, species, values_vec, forms_vec)
end

## Parse compartment
# compartment <name> is <constant/variable> = <value> in <unit>
function parse_compartment(s::ASCIIString)
  sentence = tokenize(strip(s))
  c = 1
  name = ""
  value = 0.0
  unit = ""
  ctype = ""
  # Process non-blank tokens
  for i in sentence
    i == "" && continue
    # The first token should be the keyword
    if c == 1
      c = 2
    # The second token should be the name of the parameter
    elseif c == 2
      name = i
      c = 3
    elseif c == 3
      i != "is" ? error("Error in compartment $name: I was expecting they keyword 'is' followed by const/var") : (c = 4)
    elseif c == 4
      ctype = i
      c = 5
    # The fifth token should be the equal sign, otherwise we have an error
    elseif c == 5
      i != "=" ? error("Error in compartment $name: I was expecting an equal = sign.") : (c = 6)
    # The fourth (and posterior) token should be the expression
    elseif (c == 6) & (i != "in")
      try value = parse(i) catch error("Error in compartment $name: You must supply a value") end
    elseif (c == 6) & (i == "in")
      c = 7
    elseif (c == 7)
      unit *= i
    end
  end
  if ctype == "const"
    return name, dt.Parameter(value, dt.Unit(unit))
  elseif ctype == "var"
    return name, dt.State(value, dt.Unit(unit), false)
  else
    error("Error in compartment $name: You must specify the type of compartment with the keywords const or var")
  end
end

#### Create a dictionary that stores all the elements of the model description ####
function process_file(file::ASCIIString)
    source_file = read_file(file);
    process_text(source_file)
end

# Parse the model and store all the ouputs in a ReactionSource object
function process_text(source_file::ASCIIString)
    # OrderedDicts where the sentences will be stored according to typed structures
    sentences = {"constant" => OrderedDict{ASCIIString, dt.Parameter}(),
                 "parameter" => OrderedDict{ASCIIString, dt.Parameter}(),
                 "forcing" => OrderedDict{ASCIIString, dt.Forcing}(),
                 "state" => OrderedDict{ASCIIString, dt.State}(),
                 "qsstate" => OrderedDict{ASCIIString, dt.State}(),
                 "species" => OrderedDict{ASCIIString, dt.Species}(),
                 "qsspecies" => OrderedDict{ASCIIString, dt.Species}(),
                 "component" => OrderedDict{ASCIIString, dt.Component}(),
                 "variable" => OrderedDict{ASCIIString, dt.Variable}(),
                 "reaction" => OrderedDict{ASCIIString, dt.Reaction}(),
                 "compartment" => OrderedDict{ASCIIString, Union(dt.Parameter, dt.State)}(),
                 "derivative" => OrderedDict{ASCIIString, dt.Derivative}(),
                 "observed" => OrderedDict{ASCIIString, dt.Observed}()}

    # Parsers for every type of line
    parsers = {"constant" => parse_constant,
               "parameter" => parse_parameter,
               "forcing" => parse_forcing,
               "state" => parse_state,
               "qsstate" => parse_qsstate,
               "species" => parse_species,
               "qsspecies" => parse_qsspecies,
               "component" => parse_component,
               "variable" => parse_variable,
               "reaction" => parse_reaction,
               "compartment" => parse_compartment,
               "derivative" => parse_derivative,
               "observed" => parse_observed}

    # Read the file and parse it into logical lines
    source_file = strip_comments(source_file)
    source_file = tokenize(source_file)
    source_file = logical_lines(source_file)
    source_lines = store_logical_lines(source_file)

    # Process each type of input and store in the dictionaries
    for i in collect(keys(sentences))
      for j in source_lines[i]
        name, object = parsers[i](j)
        sentences[i][name] = object
      end
    end

    # Components are actually part of species.
    # Iterate over components and add them to the corresponding species
    for (key,val) in sentences["component"]
      push!(sentences["species"][val.Species].Components, val)
    end

    # Merge quasi-steady and dynamic variables
    merge!(sentences["species"], sentences["qsspecies"])
    merge!(sentences["state"], sentences["qsstate"])

    # Put compartments in the right compartments
    for (key,val) in sentences["compartment"]
      isa(val, dt.Parameter) ?
        sentences["parameter"][key] = val :
        sentences["state"][key] = val
    end

    # Check if we have any quasi steady state chemical species or state variable
    if (length(sentences["qsspecies"]) > 1) | (length(sentences["qsstate"]) > 1)
      DynamicType = "DAE"
    else
      DynamicType = "ODE"
    end

    # Create the ReactionSource object
    dt.ReactionSource(DynamicType,
                  sentences["constant"],
                  sentences["parameter"],
                  sentences["forcing"],
                  sentences["species"],
                  sentences["reaction"],
                  sentences["state"],
                  sentences["derivative"],
                  sentences["variable"],
                  sentences["observed"])

end


text:
export greedyfit, update_component, update_coeffs

"""
$(TYPEDSIGNATURES)

Performs a greedy feature selection of an `HermiteMapComponent` up to `maxterms` features, based on the ensemble matrix `X`.

The following optionial settings can be tuned:
* `α::Float64 = αreg`: a Tikkonov regularization parameter of the cost function
* `withconstant::Bool = false`: the option to remove the constant feature in the greedy optimization routine, if the zero feature is the constant function for the basis of interest.
* `withqr::Bool = false`: improve the conditioning of the optimization problem with a QR factorization of the feature basis (recommended option)
* `maxpatience::Int64 = 10^5`: for `optimkind = split`, the maximum number of extra terms that can be added without decreasing the validation error before the greedy optimmization get stopped.
* `verbose::Bool = false`: prints details of the optimization procedure, current component optimize, training and validation errors, number of features added.
* `hessprecond::Bool=true`: use a preconditioner based on the Gauss-Newton of the Hessian of the loss function to accelerate the convergence.
* `b::String="CstProHermiteBasis"`: several bases for the feature expansion are available, see src/hermitemap/basis.jl for more details
* `ATMcriterion::String="gradient"`: sensitivty criterion used to select the feature to add to the expansion. The default uses the derivative of the cost function with respect to the coefficient
   of the features in the reduced margin of the current set of features.
"""
function greedyfit(m::Int64, Nx::Int64, X, maxterms::Int64; α::Float64 = αreg, withconstant::Bool = false, withqr::Bool = false,
                   maxpatience::Int64 = 10^5, verbose::Bool = true, hessprecond::Bool=true,
                   b::String="CstProHermiteBasis", ATMcriterion::String="gradient")

    @assert maxterms >=1 "maxterms should be >= 1"
    best_valid_error = Inf
    patience = 0

    train_error = Float64[]

    # Initialize map C to identity
    C = HermiteMapComponent(m, Nx; α = α, b = b);

    # Compute storage # Add later storage for validation S_valid
    S = Storage(C.I.f, X)

    # Remove or not the constant i.e the multi-index [0 0 0]
    if withconstant == false && iszerofeatureactive(C.I.f.B.B) == false
        # Compute the reduced margin
        reduced_margin = getreducedmargin(getidx(C))
        f = ExpandedFunction(C.I.f.B, reduced_margin, zeros(size(reduced_margin,1)))
        C = HermiteMapComponent(f; α = αreg)
        S = Storage(C.I.f, X)
        coeff0 = getcoeff(C)
        dJ = zero(coeff0)

        negative_log_likelihood!(nothing, dJ, coeff0, S, C, X)
        _, opt_dJ_coeff_idx = findmax(abs.(dJ))

        opt_idx = reduced_margin[opt_dJ_coeff_idx:opt_dJ_coeff_idx,:]

        f = ExpandedFunction(C.I.f.B, opt_idx, zeros(size(opt_idx,1)))
        C = HermiteMapComponent(f; α = α)
        S = Storage(C.I.f, X)
    end

    if withqr == true
        F = QRscaling(S)
    end

    # Compute initial loss on training set
    push!(train_error, negative_log_likelihood!(0.0, nothing, getcoeff(C), S, C, X))

    if verbose == true
        println(string(ncoeff(C))*" terms - Training error: "*string(train_error[end]))
    end

    # Compute the reduced margin
    reduced_margin = getreducedmargin(getidx(C))

    while ncoeff(C) < maxterms
        coeff_old = getcoeff(C)
        idx_old = getidx(C)

        idx_new, reduced_margin = update_component!(C, X, reduced_margin, S; ATMcriterion = ATMcriterion)

        # Update storage with the new feature
        S = update_storage(S, X, idx_new[end:end,:])

        # Update C
        C = HermiteMapComponent(IntegratedFunction(S.f); α = C.α)
        # Optimize coefficients
        if withqr == false
            coeff0 = getcoeff(C)
            if hessprecond == true
                precond = zeros(ncoeff(C), ncoeff(C))
                precond!(precond, coeff0, S, C, X)
                precond_chol = cholesky(Symmetric(precond); check = false)

                if issuccess(precond_chol) == true
                    res = Optim.optimize(Optim.only_fg!(negative_log_likelihood(S, C, X)), coeff0,
                          Optim.LBFGS(; m = 10, P = Preconditioner(Symmetric(precond), precond_chol)))
                elseif cond(Diagonal(precond)) < 10^6
                    res = Optim.optimize(Optim.only_fg!(negative_log_likelihood(S, C, X)), coeff0,
                          Optim.LBFGS(; m = 10, P = Diagonal(precond)))
                else
                    res = Optim.optimize(Optim.only_fg!(negative_log_likelihood(S, C, X)), coeff0,
                          Optim.LBFGS(; m = 10))
                end
            else

            res = Optim.optimize(Optim.only_fg!(negative_log_likelihood(S, C, X)), coeff0,
                                 Optim.LBFGS(; m = 10))
            end

            if Optim.converged(res) == true
                setcoeff!(C, Optim.minimizer(res))
                # Compute new loss on training set
                push!(train_error, negative_log_likelihood!(0.0, nothing, getcoeff(C), S, C, X))
            else
                # if optimization wasn't successful, return map
                setcoeff!(C, Optim.minimizer(res))
                append!(train_error, Inf*ones(maxterms - size(getcoeff(C), 1) + 1))
                println("Optimization wasn't successful")
                break
            end
        else

            coeff0 = getcoeff(C)
            # F = QRscaling(S)
            F = updateQRscaling(F, S)

            # Check condition number of basis evaluations: cond(ψoff ⊗ ψd) = cond(Q R) = cond(R) as Q is an orthogonal matrix
            if cond(F.R) > 1e9 || size(F.R,1) < size(F.R,2)
                println("Warning: Condition number reached")
                append!(train_error, Inf*ones(maxterms - size(getcoeff(C), 1) + 1))
                # Revert to the previous map component
                C = HermiteMapComponent(ExpandedFunction(C.I.f.B, idx_old, coeff_old); α = C.α)
                break
            end
            mul!(coeff0, F.U, coeff0)

            if hessprecond == true
                qrprecond = zeros(ncoeff(C), ncoeff(C))
                qrprecond!(qrprecond, coeff0, F, S, C, X)
                qrprecond_chol = cholesky(Symmetric(qrprecond); check = false)

                if issuccess(qrprecond_chol) == true
                    res = Optim.optimize(Optim.only_fg!(qrnegative_log_likelihood(F, S, C, X)), coeff0,
                                         Optim.LBFGS(; m = 10, P = Preconditioner(Symmetric(qrprecond), qrprecond_chol)))
                elseif cond(Diagonal(qrprecond)) < 10^6
                    res = Optim.optimize(Optim.only_fg!(qrnegative_log_likelihood(F, S, C, X)), coeff0,
                                         Optim.LBFGS(; m = 10, P = Diagonal(qrprecond)))
                else
                    res = Optim.optimize(Optim.only_fg!(qrnegative_log_likelihood(F, S, C, X)), coeff0,
                                         Optim.LBFGS(; m = 10))
                end
                # mul!(S.ψoffψd0, S.ψoffψd0, F.Uinv)
                # mul!(S.ψoffdψxd, S.ψoffdψxd, F.Uinv)
            else
                res = Optim.optimize(Optim.only_fg!(qrnegative_log_likelihood(F, S, C, X)), coeff0,
                                     Optim.LBFGS(; m = 10))
            end

            # Reverse to the non-QR space and update in-place the coefficients
            C.I.f.coeff .= F.Uinv*Optim.minimizer(res)

            if Optim.converged(res) == true
                # The computation is the non-QR space is slightly faster,
                # even if we prefer the QR form to optimize the coefficients
                push!(train_error, negative_log_likelihood!(0.0, nothing, getcoeff(C), S, C, X))
            else
                println("Optimization wasn't successful")
                append!(train_error, Inf*ones(maxterms - size(getcoeff(C), 1) + 1))
                break
            end
        end

        if verbose == true
            println(string(ncoeff(C))*" terms - Training error: "*string(train_error[end]))
        end
    end

    return C, train_error
end

"""
$(TYPEDSIGNATURES)

Compute the sensitivity of the loss function with respect to the multi-indices in the reduced-margin,
determines the feature to add to the exsiting set based on the `ATMcriterion`,
and finally returns the new set of multi-indices and the updated redcued margin.
"""
function update_component!(C::HermiteMapComponent, X, reduced_margin::Array{Int64,2}, S::Storage; ATMcriterion::String="gradient")
    m = C.m
    Nψ = C.Nψ
    idx_old = getidx(C)

    idx_new = vcat(idx_old, reduced_margin)

    # Define updated map
    f_new = ExpandedFunction(C.I.f.B, idx_new, vcat(getcoeff(C), zeros(size(reduced_margin,1))))
    C_new = HermiteMapComponent(f_new; α = αreg)

    # Set coefficients based on previous optimal solution
    coeff_new, coeff_idx_added, idx_added = update_coeffs(C, C_new)

    # Compute gradient after adding the new elements
    S = update_storage(S, X, reduced_margin)
    dJ = zero(coeff_new)
    negative_log_likelihood!(nothing, dJ, coeff_new, S, C_new, X)
    # Find function in the reduced margin most correlated with the residual

    if ATMcriterion == "gradient"
        _, opt_dJ_coeff_idx = findmax(abs.(dJ[coeff_idx_added]))
    elseif ATMcriterion == "normalized gradient"
        # The definition of ψnorm has a rescaling by 1/√Ne,
        # but it doesn't matter as all the entries are multiplied by the same constant
        _, opt_dJ_coeff_idx = findmax(abs.(dJ[coeff_idx_added]).^2 ./S.ψnorm[coeff_idx_added].^2)
    else
        error(ATMcriterion*" is not implemented.")
    end

    opt_idx = idx_added[opt_dJ_coeff_idx,:]

    # Update multi-indices and the reduced margins based on opt_idx
    # reducedmargin_opt_coeff_idx = Bool[opt_idx == x for x in eachslice(reduced_margin; dims = 1)]

    idx_new, reduced_margin = updatereducedmargin(idx_old, reduced_margin, opt_dJ_coeff_idx)

    return idx_new, reduced_margin
end

"""
$(TYPEDSIGNATURES)

Returns the value and index of the new coefficients, as well as their associated multi-indices.
"""
function update_coeffs(Cold::HermiteMapComponent, Cnew::HermiteMapComponent)

    Nψ = Cold.Nψ
    idx_old = getidx(Cold)
    idx_new = getidx(Cnew)

    coeff_old = getcoeff(Cold)

    # Declare vectors for new coefficients and to track added terms
    Nψnew = Cnew.Nψ
    coeff_new = zeros(Nψnew)
    coeff_added = ones(Nψnew)

    # Update coefficients
    @inbounds for i = 1:Nψ
        idx_i = Bool[idx_old[i,:] == x for x in eachslice(idx_new; dims = 1)]
        coeff_new[idx_i] .= coeff_old[i]
        coeff_added[idx_i] .= 0.0
    end

    # Find indices of added coefficients and corresponding multi-indices
    coeff_idx_added = findall(coeff_added .> 0.0)
    idx_added = idx_new[coeff_idx_added,:]

    return coeff_new, coeff_idx_added, idx_added
end

"""
$(TYPEDSIGNATURES)

Performs a greedy feature selection of an `HermiteMapComponent` up to `maxterms` features, based on the pair of ensemble matrices `X` (training set) and `Xvalid` (validation set).

The following optionial settings can be tuned:
* `α::Float64 = αreg`: a Tikkonov regularization parameter of the cost function
* `withconstant::Bool = false`: the option to remove the constant feature in the greedy optimization routine, if the zero feature is the constant function for the basis of interest.
* `withqr::Bool = false`: improve the conditioning of the optimization problem with a QR factorization of the feature basis (recommended option)
* `maxpatience::Int64 = 10^5`: for `optimkind = split`, the maximum number of extra terms that can be added without decreasing the validation error before the greedy optimmization get stopped.
* `verbose::Bool = false`: prints details of the optimization procedure, current component optimize, training and validation errors, number of features added.
* `hessprecond::Bool=true`: use a preconditioner based on the Gauss-Newton of the Hessian of the loss function to accelerate the convergence.
* `b::String="CstProHermiteBasis"`: several bases for the feature expansion are available, see src/hermitemap/basis.jl for more details
* `ATMcriterion::String="gradient"`: sensitivty criterion used to select the feature to add to the expansion. The default uses the derivative of the cost function with respect to the coefficient
   of the features in the reduced margin of the current set of features.
"""
function greedyfit(m::Int64, Nx::Int64, X, Xvalid, maxterms::Int64; α::Float64 = αreg, withconstant::Bool = false,
                   withqr::Bool = false, maxpatience::Int64 = 10^5, verbose::Bool = true, hessprecond::Bool=true,
                   b::String="CstProHermiteBasis", ATMcriterion::String="gradient")

    best_valid_error = Inf
    patience = 0

    train_error = Float64[]
    valid_error = Float64[]

    # Initialize map C to identity
    C = HermiteMapComponent(m, Nx; α = α, b = b);

    # Compute storage # Add later storage for validation S_valid
    S = Storage(C.I.f, X)
    Svalid = Storage(C.I.f, Xvalid)

    # Remove or not the constant i.e the multi-index [0 0 0]
    if withconstant == false && iszerofeatureactive(C.I.f.B.B) == false

        # Compute the reduced margin
        reduced_margin = getreducedmargin(getidx(C))
        f = ExpandedFunction(C.I.f.B, reduced_margin, zeros(size(reduced_margin,1)))
        C = HermiteMapComponent(f; α = α)
        S = Storage(C.I.f, X)
        coeff0 = getcoeff(C)
        dJ = zero(coeff0)

        negative_log_likelihood!(nothing, dJ, coeff0, S, C, X)
        _, opt_dJ_coeff_idx = findmax(abs.(dJ))

        opt_idx = reduced_margin[opt_dJ_coeff_idx:opt_dJ_coeff_idx,:]

        f = ExpandedFunction(C.I.f.B, opt_idx, zeros(size(opt_idx,1)))
        C = HermiteMapComponent(f; α = αreg)
        S = Storage(C.I.f, X)
        Svalid = Storage(C.I.f, Xvalid)
    end

    if withqr == true
        F = QRscaling(S)
    end

    # Compute initial loss on training set
    push!(train_error, negative_log_likelihood!(0.0, nothing, getcoeff(C), S, C, X))
    push!(valid_error, negative_log_likelihood!(0.0, nothing, getcoeff(C), Svalid, C, Xvalid))

    if verbose == true
        println(string(ncoeff(C))*" terms - Training error: "*
        string(train_error[end])*", Validation error: "*string(valid_error[end]))
    end

    # Compute the reduced margin
    reduced_margin = getreducedmargin(getidx(C))

    while ncoeff(C) < maxterms
        coeff_old = getcoeff(C)
        idx_old = getidx(C)
        idx_new, reduced_margin = update_component!(C, X, reduced_margin, S; ATMcriterion = ATMcriterion)

        # Update storage with the new feature
        S = update_storage(S, X, idx_new[end:end,:])
        Svalid = update_storage(Svalid, Xvalid, idx_new[end:end,:])

        # Update C
        C = HermiteMapComponent(IntegratedFunction(S.f); α = C.α)

        # Optimize coefficients
        if withqr == false
            coeff0 = getcoeff(C)
            if hessprecond == true
                precond = zeros(ncoeff(C), ncoeff(C))
                precond!(precond, coeff0, S, C, X)
                precond_chol = cholesky(Symmetric(precond); check = false)

                if issuccess(precond_chol) == true
                    res = Optim.optimize(Optim.only_fg!(negative_log_likelihood(S, C, X)), coeff0,
                          Optim.LBFGS(; m = 10, P = Preconditioner(Symmetric(precond), precond_chol)))
                elseif cond(Diagonal(precond)) < 10^6
                    res = Optim.optimize(Optim.only_fg!(negative_log_likelihood(S, C, X)), coeff0,
                          Optim.LBFGS(; m = 10, P = Diagonal(precond)))
                else #no preconditioner
                    res = Optim.optimize(Optim.only_fg!(negative_log_likelihood(S, C, X)), coeff0,
                          Optim.LBFGS(; m = 10))
                end
            else
                res = Optim.optimize(Optim.only_fg!(negative_log_likelihood(S, C, X)), coeff0,
                      Optim.LBFGS(; m = 10))
            end

            setcoeff!(C, Optim.minimizer(res))

            if Optim.converged(res) == true
                # Compute new loss on training and validation sets
                push!(train_error, negative_log_likelihood!(0.0, nothing, getcoeff(C), S, C, X))
                push!(valid_error, negative_log_likelihood!(0.0, nothing, getcoeff(C), Svalid, C, Xvalid))
            else
                println("Optimization wasn't successful")
                # Compute new loss on training and validation sets
                append!(train_error, Inf*ones(maxterms - size(getcoeff(C), 1) + 1))
                append!(valid_error, Inf*ones(maxterms - size(getcoeff(C), 1) + 1))
                break
            end

        else
            coeff0 = getcoeff(C)
            F = updateQRscaling(F, S)
            # F = QRscaling(S)

            # Check condition number of basis evaluations: cond(ψoff ⊗ ψd) = cond(Q R) = cond(R) as Q is an orthogonal matrix
            if cond(F.R) > 1e9 || size(F.R,1) < size(F.R,2)
                println("Warning: Condition number reached")
                append!(train_error, Inf*ones(maxterms - size(getcoeff(C), 1) + 1))
                append!(valid_error, Inf*ones(maxterms - size(getcoeff(C), 1) + 1))
                # Revert to the previous map component
                C = HermiteMapComponent(ExpandedFunction(C.I.f.B, idx_old, coeff_old); α = C.α)
                break
            end

            mul!(coeff0, F.U, coeff0)
            if hessprecond == true
                qrprecond = zeros(ncoeff(C), ncoeff(C))
                qrprecond!(qrprecond, coeff0, F, S, C, X)
                qrprecond_chol = cholesky(Symmetric(qrprecond); check = false)

                if issuccess(qrprecond_chol) == true
                    res = Optim.optimize(Optim.only_fg!(qrnegative_log_likelihood(F, S, C, X)), coeff0,
                                         Optim.LBFGS(; m = 10, P = Preconditioner(Symmetric(qrprecond), qrprecond_chol)))
                elseif cond(Diagonal(qrprecond)) < 10^6
                    res = Optim.optimize(Optim.only_fg!(qrnegative_log_likelihood(F, S, C, X)), coeff0,
                                         Optim.LBFGS(; m = 10, P = Diagonal(qrprecond)))
                else
                    res = Optim.optimize(Optim.only_fg!(qrnegative_log_likelihood(F, S, C, X)), coeff0,
                                         Optim.LBFGS(; m = 10))
                end
            else
                res = Optim.optimize(Optim.only_fg!(qrnegative_log_likelihood(F, S, C, X)), coeff0,
                                     Optim.LBFGS(; m = 10))
            end
            # Reverse to the non-QR space and update in-place the coefficients
            mul!(view(C.I.f.coeff,:), F.Uinv, Optim.minimizer(res))

            if Optim.converged(res) == true
                # The computation is the non-QR space is slightly faster,
                # even if we prefer the QR form to optimize the coefficients
                push!(train_error, negative_log_likelihood!(0.0, nothing, getcoeff(C), S, C, X))
                push!(valid_error, negative_log_likelihood!(0.0, nothing, getcoeff(C), Svalid, C, Xvalid))
            else
                println("Optimization wasn't successful")
                # Compute new loss on training and validation sets
                append!(train_error, Inf*ones(maxterms - size(getcoeff(C), 1) + 1))
                append!(valid_error, Inf*ones(maxterms - size(getcoeff(C), 1) + 1))
                break
            end

        end

        if verbose == true
            println(string(ncoeff(C))*" terms - Training error: "*
            string(train_error[end])*", Validation error: "*string(valid_error[end]))
        end

        # Update patience
        if valid_error[end] >= best_valid_error
            patience += 1
        else
            best_valid_error = deepcopy(valid_error[end])
            patience = 0
        end

        # Check if patience exceeded maximum patience
        if patience >= maxpatience
            println("Patience has been exceeded")
            break
        end
    end

    return C, (train_error, valid_error)
end

ext: jl


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""SessionAffinityConfig represents the configurations of session affinity.

    IoK8sApiCoreV1SessionAffinityConfig(;
        clientIP=nothing,
    )

    - clientIP::IoK8sApiCoreV1ClientIPConfig : clientIP contains the configurations of Client IP based session affinity.
"""
mutable struct IoK8sApiCoreV1SessionAffinityConfig <: SwaggerModel
    clientIP::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ClientIPConfig } # spec name: clientIP

    function IoK8sApiCoreV1SessionAffinityConfig(;clientIP=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1SessionAffinityConfig, Symbol("clientIP"), clientIP)
        setfield!(o, Symbol("clientIP"), clientIP)
        o
    end
end # type IoK8sApiCoreV1SessionAffinityConfig

const _property_map_IoK8sApiCoreV1SessionAffinityConfig = Dict{Symbol,Symbol}(Symbol("clientIP")=>Symbol("clientIP"))
const _property_types_IoK8sApiCoreV1SessionAffinityConfig = Dict{Symbol,String}(Symbol("clientIP")=>"IoK8sApiCoreV1ClientIPConfig")
Base.propertynames(::Type{ IoK8sApiCoreV1SessionAffinityConfig }) = collect(keys(_property_map_IoK8sApiCoreV1SessionAffinityConfig))
Swagger.property_type(::Type{ IoK8sApiCoreV1SessionAffinityConfig }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1SessionAffinityConfig[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1SessionAffinityConfig }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1SessionAffinityConfig[property_name]

function check_required(o::IoK8sApiCoreV1SessionAffinityConfig)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1SessionAffinityConfig }, name::Symbol, val)
end


text:
# points
convert(::Type{Point2}, p::Loc) = Point2(cx(p), cy(p))
convert(::Type{Point3}, p::Loc) = Point3(cx(p), cy(p), cz(p))

convert(::Type{Loc}, ::Origin) = u0()
@cxxdereference convert(::Type{Loc}, p::Point2) = xy( x(p)[], y(p)[])
@cxxdereference convert(::Type{Loc}, p::Point3) = xyz(x(p)[], y(p)[], z(p)[])
for WP ∈ (:WeightedPoint2, :WeightedPoint3)
    # lossy conversion
    @eval @cxxdereference convert(::Type{Loc}, wp::$WP) =
        convert(Loc, point(wp))
end

for P ∈ (:Point2, :Point3
       , :WeightedPoint2, :WeightedPoint3)
    @eval begin
        convert(::Type{$P}, p::Point) = convert($P, p.position)
        @cxxdereference convert(::Type{Point}, p::$P) =
            Khepri.point(convert(Loc, p))
    end
end

# directions/vectors
convert(T::Type{<:Union{Direction2,Vector2}}, u::Vec) =
    let v = in_world(u)
        T(cx(v), cy(v))
    end
convert(T::Type{<:Union{Direction3,Vector3}}, u::Vec) =
    let v = in_world(u)
        T(cx(v), cy(v), cz(v))
    end

convert(::Type{Vec}, ::NullVector) = vxy(0, 0, world_cs)
@cxxdereference convert(::Type{Vec}, u::Vector2) =
    vxy(x(u)[], y(u)[], world_cs)
@cxxdereference convert(::Type{Vec}, u::Vector3) =
    vxyz(x(u)[], y(u)[], z(u)[], world_cs)
for D ∈ (:Direction2, :Direction3)
    @eval @cxxdereference convert(::Type{Vec}, d::$D) =
        convert(Vec, vector(d))
end

# circles
convert(::Type{Circle2}, c::Union{Circle,SurfaceCircle,Khepri.CircularPath}) =
    Circle2(convert(Point2, c.center), c.radius^2)
convert(::Type{Circle3}, c::Union{Circle,SurfaceCircle,Khepri.CircularPath}) =
    Circle3(convert(Point3, c.center), c.radius^2,
            convert(Vector3, vz(1, c.center.cs)))

@cxxdereference convert(::Type{Path}, c::Circle2) =
    circular_path(convert(Loc, center(c)), √squared_radius(c))
@cxxdereference convert(::Type{Path}, c::Circle3) =
    let cₒ = convert(Loc, center(c)),
        n = convert(Vec, orthogonal_vector(supporting_plane(c)))
        circular_path(loc_from_o_vz(cₒ, n), √squared_radius(c))
    end
for C ∈ (:Circle2, :Circle3)
    @eval begin
        @cxxdereference convert(::Type{Circle}, c::$C) =
            let cp = convert(Path, c)
                circle(cp.center, cp.radius)
            end
        @cxxdereference convert(::Type{Surface}, c::$C) =
            surface(convert(Circle, c))
        @cxxdereference convert(::Type{SurfacePath}, c::$C) =
            surface_path(convert(Path, c))
        @cxxdereference convert(::Type{SurfaceCircle}, c::$C) =
            let cp = convert(Path, c)
                surface_circle(cp.center, cp.radius)
            end
    end
end

# spheres
convert(::Type{Sphere3}, s::Sphere) =
    Sphere3(convert(Point3, s.center), s.radius^2)
@cxxdereference convert(::Type{Sphere}, s::Sphere3) =
    sphere(convert(Loc, center(s)), √squared_radius(s))

# lines/rays/segments
for T ∈ (:Line2, :Line3
       , :Ray2, :Ray3
       , :Segment2, :Segment3)
    @eval convert(::Type{$T}, l::Union{Line,Khepri.PolygonalPath}) =
        $T(convert.(_pointfor($T), l.vertices[1:2])...)
end

for L ∈ (:Line2, :Line3)
    @eval @cxxdereference convert(::Type{Path}, l::$L) =
        polygonal_path(convert.(Loc, point.(Ref(l), [0, 1])))
end
for R ∈ (:Ray2, :Ray3)
    @eval @cxxdereference convert(::Type{Path}, r::$R) =
        let s = source(r)
            polygonal_path(convert.(Loc, [s, s + to_vector(r)]))
        end
end
for S ∈ (:Segment2, :Segment3)
    @eval @cxxdereference convert(::Type{Path}, s::$S) =
        polygonal_path(convert.(Loc, [source(s), target(s)]))
end

for T ∈ (:Line2, :Line3
       , :Ray2, :Ray3
       , :Segment2, :Segment3)
    @eval @cxxdereference convert(::Type{Line}, t::$T) =
        let p = convert(Path, t)
            line(p.vertices)
        end
end

# bbox2/rectangle
convert( ::Type{IsoRectangle2},
        r::Union{Rectangle,SurfaceRectangle,Khepri.RectangularPath}) =
    let p = convert(Point2, r.corner),
        q = convert(Point2, p + Vector2(r.dx, r.dy))

        IsoRectangle2(p, q)
    end
convert( ::Type{Bbox2},
        r::Union{Rectangle,SurfaceRectangle,Khepri.RectangularPath}) =
    bbox(convert(IsoRectangle2, r))

@cxxdereference convert(::Type{Path}, ir::IsoRectangle2) =
    let bl = min(ir),
        Δtr = convert(Vec, max(ir) - bl)
        rectangular_path(convert(Loc, bl), Δtr.x, Δtr.y)
    end
@cxxdereference convert(::Type{Rectangle}, ir::IsoRectangle2) =
    let rp = convert(Path, ir)
        rectangle(rp.corner, rp.dx, rp.dy)
    end
@cxxdereference convert(::Type{Surface}, ir::IsoRectangle2) =
    surface(convert(Rectangle, ir))
@cxxdereference convert(::Type{SurfacePath}, ir::IsoRectangle2) =
    surface_path(convert(Path, ir))
@cxxdereference convert(::Type{SurfaceRectangle}, ir::IsoRectangle2) =
    let rp = convert(Path, ir)
        surface_rectangle(rp.corner, rp.dx, rp.dy)
    end
@cxxdereference convert(T::Type{<:Union{Path,Rectangle,SurfacePath,
                                        SurfaceRectangle}}, b::Bbox2) =
    convert(T, IsoRectangle2(b))

# bbox3/box/cuboid
convert(::Type{Bbox3}, c::Union{Box,Cuboid}) = bbox(convert(IsoCuboid3, c))
convert(::Type{IsoCuboid3}, c::Cuboid) =
    bounding_box(convert.(Point3, [c.b0, c.b1, c.b2, c.b3,
                                   c.t0, c.t1, c.t2, c.t3]))
convert(::Type{IsoCuboid3}, b::Box) =
    let p = convert(Point3, b.c),
        q = convert(Point3, p + vxyz(b.dx, b.dy, b.dz))
        IsoCuboid3(p, q)
    end
@cxxdereference convert(T::Type{<:Union{Box,Cuboid}}, b::Bbox3) =
    convert(T, IsoCuboid3(b))
@cxxdereference convert(::Type{Box}, ic::IsoCuboid3) =
    box(convert(Loc, min(ic)), convert(Loc, max(ic)))
@cxxdereference convert(::Type{Cuboid}, ic::IsoCuboid3) =
    cuboid(convert.(Loc, vertex.(ic, 0:7))...)

# triangles
convert(T::Type{<:Union{Triangle2,Triangle3}},
        p::Union{Polygon,SurfacePolygon,SurfaceRegularPolygon,
                 Khepri.ClosedPolygonalPath}) =
    T(convert.(_pointfor(T), p.vertices[1:3])...)
convert(T::Type{<:Union{Triangle2,Triangle3}}, sp::SurfacePath) =
    convert(T, sp.path)
for T ∈ (:Triangle2, :Triangle3)
    @eval begin
        @cxxdereference convert(::Type{Path}, t::$T) =
            closed_polygonal_path(convert.(Loc, vertex.(Ref(t), 0:2)))
        @cxxdereference convert(::Type{Polygon}, t::$T) =
            polygon(path_vertices(convert(Path, t)))
        @cxxdereference convert(::Type{Surface}, t::$T) =
            surface(convert(Polygon, t))
        @cxxdereference convert(::Type{SurfacePath}, t::$T) =
            surface_path(convert(Path, t))
        @cxxdereference convert(::Type{SurfacePolygon}, t::$T) =
            surface_polygon(path_vertices(convert(Path, t)))
    end
end

# tetrahedron
convert( ::Type{Tetrahedron3},
        p::Union{Polygon,SurfacePolygon,Khepri.ClosedPolygonalPath}) =
    Tetrahedron3(convert.(Point3, p.vertices[1:4])...)
convert(::Type{Tetrahedron3}, sp::SurfacePath) =
    convert(Tetrahedron3, sp.path)
@cxxdereference convert(::Type{Path}, t::Tetrahedron3) =
    closed_polygonal_path(convert.(Loc, vertex.(Ref(t), 0:3)))
@cxxdereference convert(::Type{Polygon}, t::Tetrahedron3) =
    polygon(path_vertices(convert(Path, t)))
@cxxdereference convert(::Type{Surface}, t::Tetrahedron3) =
    surface(convert(Polygon, t))
@cxxdereference convert(::Type{SurfacePath}, t::Tetrahedron3) =
    surface_path(convert(Path, t))
@cxxdereference convert(::Type{SurfacePolygon}, t::Tetrahedron3) =
    surface_polygon(path_vertices(convert(Path, t)))


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiAutoscalingV2beta1ResourceMetricSource <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    targetAverageUtilization::Any # spec type: Union{ Nothing, Int32 } # spec name: targetAverageUtilization
    targetAverageValue::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApiResourceQuantity } # spec name: targetAverageValue

    function IoK8sApiAutoscalingV2beta1ResourceMetricSource(;name=nothing, targetAverageUtilization=nothing, targetAverageValue=nothing)
        o = new()
        validate_property(IoK8sApiAutoscalingV2beta1ResourceMetricSource, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiAutoscalingV2beta1ResourceMetricSource, Symbol("targetAverageUtilization"), targetAverageUtilization)
        setfield!(o, Symbol("targetAverageUtilization"), targetAverageUtilization)
        validate_property(IoK8sApiAutoscalingV2beta1ResourceMetricSource, Symbol("targetAverageValue"), targetAverageValue)
        setfield!(o, Symbol("targetAverageValue"), targetAverageValue)
        o
    end
end # type IoK8sApiAutoscalingV2beta1ResourceMetricSource

const _property_map_IoK8sApiAutoscalingV2beta1ResourceMetricSource = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("targetAverageUtilization")=>Symbol("targetAverageUtilization"), Symbol("targetAverageValue")=>Symbol("targetAverageValue"))
const _property_types_IoK8sApiAutoscalingV2beta1ResourceMetricSource = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("targetAverageUtilization")=>"Int32", Symbol("targetAverageValue")=>"IoK8sApimachineryPkgApiResourceQuantity")
Base.propertynames(::Type{ IoK8sApiAutoscalingV2beta1ResourceMetricSource }) = collect(keys(_property_map_IoK8sApiAutoscalingV2beta1ResourceMetricSource))
Swagger.property_type(::Type{ IoK8sApiAutoscalingV2beta1ResourceMetricSource }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiAutoscalingV2beta1ResourceMetricSource[name]))}
Swagger.field_name(::Type{ IoK8sApiAutoscalingV2beta1ResourceMetricSource }, property_name::Symbol) =  _property_map_IoK8sApiAutoscalingV2beta1ResourceMetricSource[property_name]

function check_required(o::IoK8sApiAutoscalingV2beta1ResourceMetricSource)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAutoscalingV2beta1ResourceMetricSource }, name::Symbol, val)
end


text:
list_test_str =
"3,4,3,1,2"

list_str = readline("day6_input.txt")

# Test function

function testit(name, result, answer)
    if result == answer
        @info "$name successful, with result $result"
    else
        @warn "$name failed, with result $result"
    end
end

# Part 1

function part1(list, n)
    school = parse.(Int, split(list, ','))
    while n > 0
        m = min(minimum(school), n - 1)
        num = count(school .== m)
        n -= m + 1
        @inbounds for j in eachindex(school)
            school[j] -= m + 1
            if school[j] < 0
                school[j] = 6
            end
        end
        append!(school, repeat([8], num))
    end
    length(school)
end

testit("Test Part 1", part1(list_test_str, 18), 26)
testit("Test Part 1", part1(list_test_str, 80), 5934)
testit("Challenge Part 1", part1(list_str, 80), 373378)

# Part 2

function part2(list, n)
    school = parse.(Int, split(list, ','))
    new = [count(school .== j) for j in 0:8]
    old = zeros(Int, 7)
    while n > 0
        offsprings = new[1] + old[1]
        new[1:end-1] .= @view new[2:end]
        old[1:end-1] .= @view old[2:end]
        new[end] = offsprings
        old[7] = offsprings
        n -= 1
    end
    return sum(new) + sum(old)
end

testit("Test Part 2", part2(list_test_str, 18), 26)
testit("Test Part 2", part2(list_test_str, 80), 5934)
testit("Test Part 2", part2(list_test_str, 256), 26984457539)
testit("Challenge Part 2", part2(list_str, 80), 373378)
testit("Challenge Part 2", part2(list_str, 256), 1682576647495)

function part2_alt1(list, n)
    school = parse.(Int, split(list, ','))
    bycount = [count(school .== j) for j in 0:8]
    @inbounds while n > 0
        offsprings = bycount[1]
        for j in 1:8
            bycount[j] = bycount[j+1]
        end
        bycount[7] += offsprings
        bycount[9] = offsprings
        n -= 1
    end
    return sum(bycount)
end

testit("Test Part 2 alt", part2_alt1(list_test_str, 18), 26)
testit("Test Part 2 alt", part2_alt1(list_test_str, 80), 5934)
testit("Test Part 2 alt", part2_alt1(list_test_str, 256), 26984457539)
testit("Challenge Part 2 alt", part2_alt1(list_str, 80), 373378)
testit("Challenge Part 2 alt", part2_alt1(list_str, 256), 1682576647495)

function part2_alt2(list, n)
    school = parse.(Int, split(list, ','))
    mat = zeros(Int, 9, n + 1)
    for j in 1:6
        mat[j+1, 1] = count(school .== j)
    end
    for k in 2:n+1
        mat[1:end-1, k] .= @view(mat[2:end, k-1])
        mat[7, k] += mat[1, k-1]
        mat[9, k] += mat[1, k-1]
    end
    return mat, sum(@view(mat[:, end]))
end

testit("Test Part 2 alt 2", part2_alt2(list_test_str, 18)[2], 26)
testit("Test Part 2 alt 2", part2_alt2(list_test_str, 80)[2], 5934)
testit("Test Part 2 alt 2", part2_alt2(list_test_str, 256)[2], 26984457539)
testit("Challenge Part 2 alt 2", part2_alt2(list_str, 80)[2], 373378)
testit("Challenge Part 2 alt 2", part2_alt2(list_str, 256)[2], 1682576647495)

nothing


text:
abstract type AbstractSampleBuf{T, N} <: AbstractArray{T, N} end

"""
Represents a multi-channel regularly-sampled buffer that stores its own sample
rate (in samples/second). The wrapped data is an N-dimensional array. A 1-channel
sample can be represented with a 1D array or an Mx1 matrix, and a C-channel
buffer will be an MxC matrix. So a 1-second stereo audio buffer sampled at
44100Hz with 32-bit floating-point samples in the time domain would have the
type SampleBuf{Float32, 2}.
"""
mutable struct SampleBuf{T, N} <: AbstractSampleBuf{T, N}
    data::Array{T, N}
    samplerate::Float64
end

# define constructor so conversion is applied to `sr`
SampleBuf(arr::Array{T, N}, sr::Real) where {T, N} = SampleBuf{T, N}(arr, sr)

"""
Represents a multi-channel regularly-sampled buffer representing the frequency-
domain spectrum of a `SampleBuf`. The wrapped data is an N-dimensional array. A
1-channel sample can be represented with a 1D array or an Mx1 matrix, and a
C-channel buffer will be an MxC matrix. So a 1-second stereo audio buffer
sampled at 44100Hz with 32-bit floating-point samples in the time domain would
have the type SampleBuf{Float32, 2}.
"""
mutable struct SpectrumBuf{T, N} <: AbstractSampleBuf{T, N}
    data::Array{T, N}
    samplerate::Float64
end

# define constructor so conversion is applied to `sr`
SpectrumBuf(arr::Array{T, N}, sr::Real) where {T, N} = SpectrumBuf{T, N}(arr, sr)

SampleBuf(T::Type, sr, dims...) = SampleBuf(Array{T}(undef, dims...), sr)
SpectrumBuf(T::Type, sr, dims...) = SpectrumBuf(Array{T}(undef, dims...), sr)
SampleBuf(T::Type, sr, len::Quantity) = SampleBuf(T, sr, inframes(Int,len,sr))
SampleBuf(T::Type, sr, len::Quantity, ch) =
    SampleBuf(T, sr, inframes(Int,len,sr), ch)
SpectrumBuf(T::Type, sr, len::Quantity) =
    SpectrumBuf(T, sr, inframes(Int,len, sr))
SpectrumBuf(T::Type, sr, len::Quantity, ch) =
    SpectrumBuf(T, sr, inframes(Int,len, sr))

# terminology:
# sample - a single value representing the amplitude of 1 channel at some point in time (or frequency)
# channel - a set of samples running in parallel
# frame - a collection of samples from each channel that were sampled simultaneously

# audio methods
samplerate(buf::AbstractSampleBuf) = buf.samplerate
nchannels(buf::AbstractSampleBuf{T, 2}) where {T} = size(buf.data, 2)
nchannels(buf::AbstractSampleBuf{T, 1}) where {T} = 1
nframes(buf::AbstractSampleBuf) = size(buf.data, 1)

function samplerate!(buf::AbstractSampleBuf, sr)
    buf.samplerate = sr

    buf
end

# define audio methods on raw buffers as well
nframes(arr::AbstractArray) = size(arr, 1)
nchannels(arr::AbstractArray) = size(arr, 2)

# it's important to define Base.similar so that range-indexing returns the
# right type, instead of just a bare array
Base.similar(buf::SampleBuf, ::Type{T}, dims::Dims) where {T} = SampleBuf(Array{T}(undef, dims), samplerate(buf))
Base.similar(buf::SpectrumBuf, ::Type{T}, dims::Dims) where {T} = SpectrumBuf(Array{T}(undef, dims), samplerate(buf))
domain(buf::AbstractSampleBuf) = range(0.0, stop=(nframes(buf)-1)/samplerate(buf), length=nframes(buf))

# There's got to be a better way to define these functions, but the dispatch
# and broadcast behavior for AbstractArrays is complex and has subtle differences
# between Julia versions, so we basically just override functions here as they
# come up as problems
import Base: +, -, *, /
import Base.broadcast

find_buf(args::Tuple) = find_buf(find_buf(args[1]), Base.tail(args))
find_buf(x) = x

for btype in (:SampleBuf, :SpectrumBuf)

        @eval find_buf(bc::Base.Broadcast.Broadcasted{Broadcast.ArrayStyle{$btype{T,N}}}) where {T,N} = find_buf(bc.args)
        @eval find_buf(::$btype, args::Tuple{$btype}) = args[1]
        @eval find_buf(::Any, args::Tuple{$btype}) = args[1]
        @eval find_buf(a::$btype, rest) = a

        @eval Base.BroadcastStyle(::Type{$btype{T,N}}) where {T,N} = Broadcast.ArrayStyle{$btype{T,N}}()
        @eval function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{$btype{T,N}}}, ::Type{ElType}) where {T,N,ElType}
            A = find_buf(bc)
            $btype(Array{ElType}(undef, length.(axes(bc))), samplerate(A))
        end

end # for btype

typename(::SampleBuf{T, N}) where {T, N} = "SampleBuf{$T, $N}"
unitname(::SampleBuf) = "s"
srname(::SampleBuf) = "Hz"
typename(::SpectrumBuf{T, N}) where {T, N} = "SpectrumBuf{$T, $N}"
unitname(::SpectrumBuf) = "Hz"
srname(::SpectrumBuf) = "s"

# from @mbauman's Sparklines.jl package
const ticks = ['▁','▂','▃','▄','▅','▆','▇','█']
# 3-arg version (with explicit mimetype) is needed because we subtype AbstractArray,
# and there's a 3-arg version defined in show.jl
function show(io::IO, ::MIME"text/plain", buf::AbstractSampleBuf)
    println(io, "$(nframes(buf))-frame, $(nchannels(buf))-channel $(typename(buf))")
    len = nframes(buf) / samplerate(buf)
    ustring = unitname(buf)
    srstring = srname(buf)
    print(io, "$(len)$ustring sampled at $(samplerate(buf))$srstring")
    nframes(buf) > 0 && showchannels(io, buf)
end

function showchannels(io::IO, buf::AbstractSampleBuf, widthchars=80)
    # number of samples per block
    blockwidth = round(Int, nframes(buf)/widthchars, RoundUp)
    nblocks = round(Int, nframes(buf)/blockwidth, RoundUp)
    blocks = Array{Char}(undef, nblocks, nchannels(buf))
    for blk in 1:nblocks
        i = (blk-1)*blockwidth + 1
        n = min(blockwidth, nframes(buf)-i+1)
        peaks = Compat.maximum(abs.(float.(buf[(1:n) .+ i .- 1, :])), dims=1)
        # clamp to -60dB, 0dB
        peaks = clamp.(20log10.(peaks), -60.0, 0.0)
        idxs = trunc.(Int, (peaks.+60)/60 * (length(ticks)-1)) .+ 1
        blocks[blk, :] = ticks[idxs]
    end
    for ch in 1:nchannels(buf)
        println(io)
        print(io, String(blocks[:, ch]))
    end
end

"""Get a pointer to the underlying data for the buffer. Will return a Ptr{T},
where T is the element type of the buffer. This is particularly useful for
passing to C libraries to fill the buffer"""
channelptr(buf::Array, channel, frameoffset=0) =
    pointer(buf) + ((channel-1)*nframes(buf)+frameoffset) * sizeof(eltype(buf))
channelptr(buf::AbstractSampleBuf, channel, frameoffset=0) =
    channelptr(buf.data, channel, frameoffset)

"""Mix the channels of the source array into the channels of the dest array,
using coefficients from the `mix` matrix. To mix an M-channel buffer to a
N-channel buffer, `mix` should be MxN. `src` and `dest` should not share
memory."""
function mix!(dest::AbstractMatrix, src::AbstractMatrix, mix::AbstractArray)
    inchans = nchannels(src)
    outchans = nchannels(dest)
    size(mix) == (inchans, outchans) || error("Mix Matrix should be $(inchans)x$(outchans)")
    mul!(dest, src, mix)
end

function mix!(dest::AbstractVector, src::AbstractVector, mix::AbstractArray)
    mix!(reshape(dest, (length(dest), 1)), reshape(src, (length(src), 1)), mix)
    dest
end

function mix!(dest::AbstractVector, src::AbstractMatrix, mix::AbstractArray)
    mix!(reshape(dest, (length(dest), 1)), src, mix)
    dest
end

function mix!(dest::AbstractMatrix, src::AbstractVector, mix::AbstractArray)
    mix!(dest, reshape(src, (length(src), 1)), mix)
end

"""Mix the channels of the source array into the channels of the dest array,
using coefficients from the `mix` matrix. To mix an M-channel buffer to a
N-channel buffer, `mix` should be MxN. `src` and `dest` should not share
memory."""
function mix(src::AbstractArray, mix::AbstractArray)
    dest = similar(src, (nframes(src), size(mix, 2)))
    mix!(dest, src, mix)
end

"""Mix the channels of the `src` array into the mono `dest` array."""
function mono!(dest::AbstractArray, src::AbstractArray)
    mix!(dest, src, ones(nchannels(src), 1) ./ nchannels(src))
end

"""Mix the channels of the `src` array into a mono array."""
function monoch(src::AbstractArray)
    dest = similar(src, (nframes(src), 1))
    mono!(dest, src)
end

# the index types that Base knows how to handle. Separate out those that index
# multiple results
const BuiltinMultiIdx = Union{Colon,
                        Vector{Int},
                        Vector{Bool},
                        AbstractRange{Int}}
const BuiltinIdx = Union{Int, BuiltinMultiIdx}
# the index types that will need conversion to built-in index types. Each of
# these needs a `toindex` method defined for it
const ConvertIdx{T1 <: Quantity, T2 <: Int} = Union{T1,
                                                # Vector{T1}, # not supporting vectors of Quantities (yet?)
                                                # Range{T1}, # not supporting ranges (yet?)
                                                ClosedInterval{T2},
                                                ClosedInterval{T1}}

"""
    toindex(buf::SampleBuf, I)

Convert the given index value to one that Base knows how to use natively for
indexing
"""
function toindex end

toindex(buf::SampleBuf, t::Number) = t
toindex(buf::SampleBuf, t::FrameQuant) = inframes(Int, t) + 1
toindex(buf::SampleBuf, t::Unitful.Time) = inframes(Int, t, samplerate(buf)) + 1
toindex(buf::SampleBuf, t::Unitful.AbstractQuantity) = throw(Unitful.DimensionError(t, s))
toindex(buf::SpectrumBuf, f::Number) = f
toindex(buf::SpectrumBuf, f::FrameQuant) = inframes(Int, f) + 1
toindex(buf::SpectrumBuf, f::Unitful.Frequency) = inframes(Int, f, samplerate(buf)) + 1
toindex(buf::SpectrumBuf, f::Unitful.AbstractQuantity) = throw(Unitful.DimensionError(f, Hz))

# indexing by vectors of Quantities not yet supported
toindex(buf::AbstractSampleBuf, I::ClosedInterval{Int}) =
    toindex(buf, minimum(I)*frames):toindex(buf, maximum(I)*frames)
toindex(buf::AbstractSampleBuf, I::ClosedInterval{T}) where {T <: Quantity} =
    toindex(buf, minimum(I)):toindex(buf, maximum(I))

# AbstractArray interface methods
Base.size(buf::AbstractSampleBuf) = size(buf.data)
Base.IndexStyle(::Type{T}) where {T <: AbstractSampleBuf} = Base.IndexLinear()
# this is the fundamental indexing operation needed for the AbstractArray interface
Base.getindex(buf::AbstractSampleBuf, i::Int) = buf.data[i];

Base.getindex(buf::AbstractSampleBuf, I::ConvertIdx) = buf[toindex(buf, I)]
Base.getindex(buf::AbstractSampleBuf, I1::ConvertIdx, I2::BuiltinIdx) =
    buf[toindex(buf, I1), I2]
# In Julia 0.5 scalar indices are now dropped, so by default indexing
# buf[5, 1:2] gives you a 2-frame single-channel buffer instead of a 1-frame
# two-channel buffer. The following getindex method defeats the index dropping
Base.getindex(buf::AbstractSampleBuf, I1::Int, I2::BuiltinMultiIdx) = buf[I1:I1, I2]

function Base.setindex!(buf::AbstractSampleBuf, val, i::Int)
    buf.data[i] = val
end

# equality
import Base.==
==(buf1::AbstractSampleBuf, buf2::AbstractSampleBuf) =
    samplerate(buf1) == samplerate(buf2) &&
    buf1.data == buf2.data

FFTW.fft(buf::SampleBuf) = SpectrumBuf(FFTW.fft(buf.data), nframes(buf)/samplerate(buf))
FFTW.ifft(buf::SpectrumBuf) = SampleBuf(FFTW.ifft(buf.data), nframes(buf)/samplerate(buf))

# does a per-channel convolution on SampleBufs
for buftype in (:SampleBuf, :SpectrumBuf)
    @eval function DSP.conv(b1::$buftype{T, 1}, b2::$buftype{T, 1}) where {T}
        if !isapprox(samplerate(b1), samplerate(b2))
            error("Resampling convolution not yet supported")
        end
        $buftype(conv(b1.data, b2.data), samplerate(b1))
    end

    @eval function DSP.conv(b1::$buftype{T, N1}, b2::$buftype{T, N2}) where {T, N1, N2}
        if !isapprox(samplerate(b1), samplerate(b2))
            error("Resampling convolution not yet supported")
        end
        if nchannels(b1) != nchannels(b2)
            error("Broadcasting convolution not yet supported")
        end
        out = $buftype(T, samplerate(b1), nframes(b1)+nframes(b2)-1, nchannels(b1))
        for ch in 1:nchannels(b1)
            out[:, ch] = conv(b1.data[:, ch], b2.data[:, ch])
        end

        out
    end

    @eval function DSP.conv(b1::$buftype{T, 1}, b2::StridedVector{T}) where {T}
        $buftype(conv(b1.data, b2), samplerate(b1))
    end

    @eval DSP.conv(b1::StridedVector{T}, b2::$buftype{T, 1}) where {T} = conv(b2, b1)

    @eval function DSP.conv(b1::$buftype{T, 2}, b2::StridedMatrix{T}) where {T}
        if nchannels(b1) != nchannels(b2)
            error("Broadcasting convolution not yet supported")
        end
        out = $buftype(T, samplerate(b1), nframes(b1)+nframes(b2)-1, nchannels(b1))
        for ch in 1:nchannels(b1)
            out[:, ch] = conv(b1.data[:, ch], b2[:, ch])
        end

        out
    end

    @eval DSP.conv(b1::StridedMatrix{T}, b2::$buftype{T, 2}) where {T} = conv(b2, b1)
end


text:
using Measurements
using Plots
using LsqFit
using ColorSchemes

val(x) = Measurements.value(x)

include("constants.jl")
include("data.jl")

## initialise data

methyls = loadmethyls("data/methyls.txt")

zq, dq, dqq, qq = loadR2("data/R2data.txt")
B0 = [600, 700, 800, 950] * 1e6 * 2π / 2.675e8;
B0² = B0.^2;

S2tc = loadS2tc("data/S2tc.txt")
ΔσH = load1Hcsa("data/etaH.txt", S2tc, 950)
ΔσC = load13Ccsa("data/13Ccsa.txt")

## calculate ξH and ξC contributions from CSA

ξHcsa = @. 2/(3*sqrt(5)) * γH * ΔσH * 1e-6 * sqrt(S2tc * 1e-9)
ξCcsa = @. 2/(3*sqrt(5)) * γC * ΔσC * 1e-6 * sqrt(S2tc * 1e-9)

## define fitting functions
# this needs to be called later as the functions rely on definitions of CSA
include("fitting.jl")
include("plotting.jl")

## select residue and run analysis
resi = 21

# get best-fit slopes (without correcting for CSA)
βZQ = fitβ(B0², zq[resi, :])
βDQ = fitβ(B0², dq[resi, :])
βDQQ = fitβ(B0², dqq[resi, :])
βQQ = fitβ(B0², qq[resi, :])

# compute χ2 surface
hval = LinRange(-.3,.3,100)
cval = LinRange(-.3,.3,120)
χ2 = zeros(length(hval), length(cval))
for j=1:length(hval)
    for k=1:length(cval)
        χ2[j,k] = calcχ2(hval[j], cval[k], resi)
    end
end

# get minimum of χ2 surface (gridsearch)
@show χ2min = minimum(vec(χ2))
@show ξHmin = hval[argmin(χ2)[1]]
@show ξCmin = cval[argmin(χ2)[2]]

# plot results
plt1 = plotlinefits(ξHmin, ξCmin, resi)
plt2 = plotχ2(hval, cval, χ2, βZQ, βDQ, βDQQ, βQQ)

plot(plt1, plt2, layout=(1,2), size=(800,400))

max_stars_repo_name: chriswaudby/qq-hahnecho

max_issues_repo_name: chriswaudby/qq-hahnecho

max_forks_repo_name: chriswaudby/qq-hahnecho


text:
# # A trivial SOS decomposition example

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/sos_decomposition.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/sos_decomposition.ipynb)
# **Contributed by**: votroto
# **Adapted from**: Examples 3.25 of [BPT12]
#
# [BPT12] Blekherman, G.; Parrilo, P. A. & Thomas, R. R.
# *Semidefinite Optimization and Convex Algebraic Geometry*.
# Society for Industrial and Applied Mathematics, **2012**.

using Test #src
using DynamicPolynomials
using SumOfSquares
import CSDP

# The polynomial `p = x^2 - x*y^2 + y^4 + 1` is SOS.
# We can, for example, decompose it as
# `p = 3/4*(x - y^2)^2 + 1/4*(x + y)^2 + 1`,
# which clearly proves that `p` is SOS, and there are infinitely many other ways
# to decompose `p` into sums of squares.

# We can use SumOfSquares.jl to find such decompositions.

# First, setup the polynomial of interest.

@polyvar x y
p = x^2 - x*y^2 + y^4 + 1

# Secondly, constrain the polynomial to be nonnegative.
# SumOfSquares.jl transparently reinterprets polyonmial nonnegativity as the
# appropriate SOS certificate for polynomials nonnegative on semialgebraic sets.

model = SOSModel(CSDP.Optimizer)
@constraint(model, cref, p >= 0)

# Thirdly, optimize the feasibility problem!

optimize!(model)

# Lastly, recover a SOS decomposition.
# In general, SOS decompositions are not unique!

sos_dec = sos_decomposition(cref, 1e-4)

# Converting, rounding, and simplifying - Huzza, Back where we began!

polynomial(sos_dec, Float32)

# ## A deeper explanation and the unexplained `1e-4` parameter

# `p = x^2 - x*y^2 + y^4 + 1` can be represented in terms of its Gram matrix as

gram = gram_matrix(cref)

#-

gram.basis.monomials' * gram.Q * gram.basis.monomials

# where the matrix `gram.Q` is positive semidefinite, because `p` is SOS. If we
# could only get the decomposition `gram.Q = V' * V`, the SOS decomposition would
# simply be `||V * monomials||^2`.

# Unfortunately, we can not use Cholesky decomposition, since gram `Q` is only
# semidefinite, not definite. Hence, SumOfSquares.jl uses SVD decomposition
# instead and discards small singular values (in our case `1e-4`).

max_stars_repo_path: examples/SOS decomposition.jl

max_issues_repo_path: examples/SOS decomposition.jl

max_forks_repo_path: examples/SOS decomposition.jl


text:
@testset "psibound" begin

  @testset "psi_lower" begin
    sp = PrimesSet(2, 100);
    fb = []; for x=sp push!(fb, fmpz(x)); end;
    fb = FactorBase(fb)
    @test length(find(x->issmooth(fb, fmpz(x)), 1:256)) == psi_lower(255, 100)[1][end]
  end
end


text:
using Test, YaoBlocks, YaoArrayRegister

function heisenberg(n::Int; periodic::Bool=true)
    Sx(i) = put(n, i=>X)
    Sy(i) = put(n, i=>Y)
    Sz(i) = put(n, i=>Z)

    return sum(1:(periodic ? n : n-1)) do i
        j = mod1(i, n)
        Sx(i) * Sx(j) + Sy(i) * Sy(j) + Sz(i) * Sz(j)
    end
end

const hm = heisenberg(4)

@testset "constructor:time evolution" begin
    te = TimeEvolution(hm, 0.2)
    # copy
    cte = copy(te)
    @test cte == te
    @test cte !== te

    # dispatch
    dispatch!(cte, [2.0])
    @test cte != te
    @test cte.dt == 2.0
    @test setiparams!(cte, 0.5).dt == 0.5
    @test setiparams!(cte, :random).dt != 0.5
    @test setiparams!(cte, :zero).dt == 0.0
end

@testset "test imaginary time evolution" begin
    tei = TimeEvolution(hm, 0.2im)
    r = rand_state(4)
    r1 = copy(r) |> tei
    @test exp(Matrix(mat(hm)) * 0.2) * r.state ≈ r1.state

    @test applymatrix(tei) ≈ mat(tei)
    @test applymatrix(adjoint(tei)) ≈ applymatrix(tei)'
    @test !isunitary(tei)
    @test !isunitary(tei |> mat)
end

@testset "test time evolution" begin
    te = TimeEvolution(hm, 0.2)

    r = rand_state(4)
    r1 = copy(r) |> te
    @test exp(Matrix(mat(hm)) * -0.2im) * r.state ≈ r1.state
    @test applymatrix(adjoint(te)) ≈ applymatrix(te)'

    @test applymatrix(te) ≈ mat(te)
    @test isunitary(te)
    @test isunitary(te |> mat)
    cte = copy(te)
    @test cte == te
end


text:
# Floating point type for simulations
T = Float64

# Geometrical setup
setup = CylinderSetup{T}(;
    name = "cylinders/somecylinders",
    ncell = 3,
    rmin = 2.0,
    rmax = 6.0,
    dmin = 0.2,
    dmax = 0.3,
    height = 1.0,
    bend = 0.0,
    twist = 0.0,
    include_in = false,
    in_ratio = 0.6,
    ecs_shape = :convex_hull,
    ecs_ratio = 0.5,
    # refinement = 0.5,
)

# Get compartimentalized coefficient vectors
coeffs = coefficients(
    setup;
    D = (; in = 0.002 * I(3), out = 0.002 * I(3), ecs = 0.002 * I(3)),
    T₂ = (; in = Inf, out = Inf, ecs = Inf),
    ρ = (; in = 1.0, out = 1.0, ecs = 1.0),
    κ = (; in_out = 1e-4, out_ecs = 1e-4, in = 0.0, out = 0.0, ecs = 0.0),
    γ = 2.67513e-4,
)


text:
# Genie-aided MMSE filter uses true covariance matrix
mmse_genie(y,h_cov,snr) = begin
    rho = 10^(0.1*snr);
    hest = zeros(y)
    (nAntennas,nCoherence,nBatches) = size(y)
    for b in 1:nBatches
        C = scm.toeplitzHe(h_cov[:,b]) # get full cov matrix
        Cr = C + eye(nAntennas)./rho
        hest[:,:,b] = C*(Cr\y[:,:,b])
    end
    hest
end

# Genie-aided OMP algorithm that stops at optimal order
omp_genie(y,h; evaluate_cost = (y0,yhat) -> sum(abs2,y0-yhat)) = begin
    (nAntennas,nCoherence,nBatches) = size(y)
    nGrid = 4*nAntennas # four times oversampled DFT
    grid = linspace(-1,1, nGrid+1)[1:nGrid]
    A = 1/sqrt(nAntennas) * exp.(1im*pi*(0:nAntennas-1) * grid')
    hest = zeros(h)
    for b in 1:nBatches
        hest[:,:,b] = omp_genie_alg(A,y[:,:,b],h[:,:,b], evaluate_cost = evaluate_cost)[1]
    end
    hest
end

"""
    omp_genie_alg(A::Matrix, y::VecOrMat, x0::VecOrMat[, B=A', find_supp = x -> L, solve_restricted = L -> x, evaluate_cost = (y0,yhat) -> c]) -> x, L
Genie aided Orthogonal Matching Pursuit with optimal support size (based on evaluate_cost criterion)
Denoises y=y0+n by approximating y=Ax with a dictionary A

# Arguments
* `A`: Full sensing matrix
* `y`: Observation vector or matrix
* `x`: True signal vector or matrix
* `B`: Mismatched filter matrix (e.g., `A*`)
* `find_supp`: function that returns a logical vector with the best index of the input
* `solve_restricted`: function that solves `y=Ax` on with support of `x` restricted to current index set `L`
* `evaluate_cost`: function that should be minimized and may also take the true y0 as input

# Output
* `x`: k-sparse solution the linear system `y=Ax`
* `L`: logical vector of nonzero indexes of `x`
"""
function omp_genie_alg{T<:Number}(A::Matrix{T},
                                  y::VecOrMat{T},
                                  y0::VecOrMat{T};
                                  # Optional
                                  k_max = 100,
                                  B=A',
			                      find_supp = x -> large_ind(vec(sum(abs2,x,2))),
	                              solve_restricted = L -> A[:,L] \ y,
                                  evaluate_cost = (y0,yhat) -> sum(abs2,y0-yhat)) # MSE criterion
	L = Int[]
    x    = complex(zeros(size(A,2),size(y,2)))
    yhat = zeros(y)
    ytmp = zeros(y)
    cost = evaluate_cost(y0,yhat)
    for k=1:k_max
        cost_prev = cost

		# Union of support sets
		L = union(L,find_supp(B*(y-A*x)))
		# Restricted LS estimate
        x[L,:] = solve_restricted(L)
        ytmp = A*x

        # Calculate angle between y and ytmp
        cost = evaluate_cost(y0,ytmp)
        if cost > cost_prev # cost increases...
            break
        end

        yhat[:] = ytmp[:]
    end
	return yhat, L
end
function large_ind(x;k=1::Int)
	return sortperm(x,rev=true)[1:k]
end

max_stars_repo_name: tum-msv/learning-mmse-est


text:
cd("INSERT_PATH")

using Distributions
using PolynomialRoots

include("ToolBox.jl")
include("predictionGSBR.jl")

# define map parameters
theta = [0.05, 2.55, 0., -0.99]
x0 = 1.
dataSeed = 2
n = 114

# noise mixture: f_{2,l}
w1 = 0.6
lam1 = 1e05
lam2 = 1e-02 * lam1

x = genData(n,theta,x0,w1,lam1,lam2,dataSeed)
# Simulation parameters, T: prediction horizon [T=0 --> Plain GSBReconstruction]
degree,maxiter,burnin,seed,T = 5,500000,250000,3,0

@time thetas, zp, sx0 = predictionGSBR(x,degree,maxiter,burnin,seed,T);

max_stars_repo_name: KosKal/GSBR-DPR

max_issues_repo_name: KosKal/GSBR-DPR

max_forks_repo_name: KosKal/GSBR-DPR


text:
using MultiStartLH
using Documenter

DocMeta.setdocmeta!(MultiStartLH, :DocTestSetup, :(using MultiStartLH); recursive=true)

makedocs(;
    modules=[MultiStartLH],
    authors="hendri54 <hendricksl@protonmail.com> and contributors",
    repo="https://github.com/hendri54/MultiStartLH.jl/blob/{commit}{path}#{line}",
    sitename="MultiStartLH.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)


text:
module ProtoBufTestCodec
using ProtoBuf
using Test
using Random
import Base: ==

import ProtoBuf.meta

macro _rand_int(T,mx,a)
   esc(quote
       $(T)(round(rand() * $(mx)) + $(a))
   end)
end

print_hdr(tname) = println("testing $tname...")

mutable struct TestType <: ProtoType
    val::Any
end

mutable struct TestStr <: ProtoType
    val::AbstractString
end
==(t1::TestStr, t2::TestStr) = (t1.val == t2.val)

mutable struct TestOptional <: ProtoType
    sVal1::TestStr
    sVal2::TestStr
    iVal2::Array{Int64,1}
end

mutable struct TestNested <: ProtoType
    fld1::TestType
    fld2::TestOptional
    fld3::Array{TestStr}
end

mutable struct TestDefaults <: ProtoType
    iVal1::Int64
    sVal2::AbstractString
    iVal2::Array{Int64,1}

    TestDefaults(f1,f2,f3) = new(f1,f2,f3)
    TestDefaults() = new()
end

mutable struct TestOneofs <: ProtoType
    iVal1::Int64
    iVal2::Int64
    iVal3::Int64

    TestOneofs() = new()
end

mutable struct TestMaps <: ProtoType
    d1::Dict{Int,Int}
    d2::Dict{Int32,String}
    d3::Dict{String,String}
    TestMaps() = new()
end

mutable struct TestFilled <: ProtoType
    fld1::TestType
    fld2::TestType
    TestFilled() = new()
end

mutable struct __enum_TestEnum <: ProtoEnum
    UNIVERSAL::Int32
    WEB::Int32
    IMAGES::Int32
    LOCAL::Int32
    NEWS::Int32
    PRODUCTS::Int32
    VIDEO::Int32
    __enum_TestEnum() = new(0,1,2,3,4,5,6)
end
const TestEnum = __enum_TestEnum()

# disable caching of meta since we manually modify them for the tests
meta(t::Type{TestType})         = meta(t, Symbol[], Int[], Dict{Symbol,Any}(), false)
meta(t::Type{TestOptional})     = meta(t, Symbol[], Int[], Dict{Symbol,Any}(), false)
meta(t::Type{TestNested})       = meta(t, Symbol[], Int[], Dict{Symbol,Any}(), false)
const _t_defaults = Dict{Symbol,Any}(:iVal1 => 10, :iVal2 => [1,2,3])
meta(t::Type{TestDefaults})     = meta(t, Symbol[], Int[], _t_defaults, false)
const _t_oneofs = Int[0,1,1]
const _t_oneof_names = [:optval]
meta(t::Type{TestOneofs})       = meta(t,  Symbol[], Int[], Dict{Symbol,Any}(), false, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, _t_oneofs, _t_oneof_names)
meta(t::Type{TestFilled})       = meta(t, Symbol[:fld1], Int[], Dict{Symbol,Any}())

function mk_test_nested_meta(o1::Bool, o2::Bool, o21::Bool, o22::Bool)
    meta1 = mk_test_meta(1, :int64)
    meta2 = mk_test_optional_meta(o21, o22)

    m = meta(TestNested)
    m.symdict[:fld1].occurrence = o1 ? 0 : 1
    m.symdict[:fld2].occurrence = o2 ? 0 : 1
    m.symdict[:fld1].meta = meta1
    m.symdict[:fld2].meta = meta2
    m
end

function mk_test_optional_meta(opt1::Bool, opt2::Bool)
    m = meta(TestOptional)
    m.symdict[:sVal1].occurrence = opt1 ? 0 : 1
    m.symdict[:sVal2].occurrence = opt2 ? 0 : 1
    m
end

function mk_test_meta(fldnum::Int, ptyp::Symbol)
    m = meta(TestType)
    attrib = m.symdict[:val]
    attrib.fldnum = fldnum
    attrib.ptyp = ptyp
    m.numdict = Dict{Int,ProtoMetaAttribs}()
    m.numdict[fldnum] = attrib
    m
end

assert_equal(::Type{Array{T,1}}, ::Type{Array{U,1}}) where {T,U} = @test (T <: U) || (U <: T)
assert_equal(T::Type, U::Type) = @test (T <: U) || (U <: T)
assert_equal(val1::Bool, val2::Bool) = @test val1 == val2
assert_equal(val1::AbstractString, val2::AbstractString) = @test val1 == val2
assert_equal(val1::Number, val2::Number) = @test val1 == val2
assert_equal(val1::Array{T,1}, val2::Array{T,1}) where {T<:Number} = @test val1 == val2
assert_equal(val1::Array{T,1}, val2::Array{T,1}) where {T<:AbstractString} = @test val1 == val2
assert_equal(val1::Dict{K,V}, val2::Dict{K,V}) where {K,V} = @test val1 == val2
function assert_equal(val1, val2)
    typ1 = typeof(val1)
    typ2 = typeof(val2)
    assert_equal(typ1, typ2)

    n = fieldnames(typ1)
    t = typ1.types
    for fld in n
        fldfill1 = isfilled(val1, fld)
        fldfill2 = isfilled(val2, fld)
        @test fldfill1 == fldfill2
        fldfill1 && assert_equal(getfield(val1, fld), getfield(val2, fld))
    end
end

function test_types()
    pb = PipeBuffer()
    testval = TestType(0)
    readval = TestType(0)

    # test enum
    print_hdr("enum")
    testval.val = @_rand_int(Int32, 10^9, 0)
    fldnum = @_rand_int(Int, 100, 1)
    meta = mk_test_meta(fldnum, :enum)
    writeproto(pb, testval, meta)
    readproto(pb, readval, meta)
    assert_equal(testval, readval)

    let typs = [Int32,Int64,UInt32,UInt64,Int32,Int64,UInt64,Int64,UInt32,Int32], ptyps=[:int32,:int64,:uint32,:uint64,:sint32,:sint64,:fixed64,:sfixed64,:fixed32,:sfixed32]
        for (typ,ptyp) in zip(typs,ptyps)
            print_hdr(ptyp)
            for idx in 1:100
                testval.val = convert(typ, @_rand_int(UInt32, 10^9, 0))
                fldnum = @_rand_int(Int, 100, 1)
                meta = mk_test_meta(fldnum, ptyp)
                writeproto(pb, testval, meta)
                readproto(pb, readval, meta)
                assert_equal(testval, readval)
            end
        end
    end

    let typs = [Int32,Int64,Int32,Int64], ptyps=[:int32,:int64,:sint32,:sint64]
        for (typ,ptyp) in zip(typs,ptyps)
            print_hdr(ptyp)
            for idx in 1:100
                testval.val = convert(typ, -1 * @_rand_int(Int32, 10^9, 0))
                fldnum = @_rand_int(Int, 100, 1)
                meta = mk_test_meta(fldnum, ptyp)
                writeproto(pb, testval, meta)
                readproto(pb, readval, meta)
                assert_equal(testval, readval)
            end
        end
    end

    print_hdr("varint overflow...")
    ProtoBuf._write_uleb(pb, -1)
    @test ProtoBuf._read_uleb(pb, Int8) == 0
    ProtoBuf._write_uleb(pb, 1)
    @test ProtoBuf._read_uleb(pb, Int8) == 1
    write(pb, 0xff)
    ProtoBuf._write_uleb(pb, -1)
    @test ProtoBuf._read_uleb(pb, Int32) == 0

    let typs = [Bool], ptyps=[:bool]
        for (typ,ptyp) in zip(typs,ptyps)
            print_hdr(ptyp)
            for idx in 1:100
                testval.val = convert(typ, @_rand_int(UInt32, 1, 0))
                fldnum = @_rand_int(Int, 100, 1)
                meta = mk_test_meta(fldnum, ptyp)
                writeproto(pb, testval, meta)
                readproto(pb, readval, meta)
                assert_equal(testval, readval)
            end
        end
    end

    let typs = [Float64,Float32], ptyps=[:double,:float]
        for (typ,ptyp) in zip(typs,ptyps)
            print_hdr(ptyp)
            for idx in 1:100
                testval.val = convert(typ, @_rand_int(UInt32, 10^9, 0))
                fldnum = @_rand_int(Int, 100, 1)
                meta = mk_test_meta(fldnum, ptyp)
                writeproto(pb, testval, meta)
                readproto(pb, readval, meta)
                assert_equal(testval, readval)
            end
        end
    end

    print_hdr("string")
    for idx in 1:100
        testval.val = randstring(50)
        fldnum = @_rand_int(Int, 100, 1)
        meta = mk_test_meta(fldnum, :string)
        writeproto(pb, testval, meta)
        readproto(pb, readval, meta)
        assert_equal(testval, readval)
    end
end

function test_repeats()
    pb = PipeBuffer()
    testval = TestType(0)
    readval = TestType(0)

    print_hdr("repeated int64")
    for idx in 1:100
        testval.val = collect(Int64, randstring(50))
        readval.val = Int64[]
        fldnum = @_rand_int(Int, 100, 1)
        meta = mk_test_meta(fldnum, :int64)
        meta.ordered[1].occurrence = 2
        writeproto(pb, testval, meta)
        readproto(pb, readval, meta)
        assert_equal(testval, readval)
    end

    print_hdr("repeated and packed int64")
    for idx in 1:100
        testval.val = collect(Int64, randstring(50))
        readval.val = Int64[]
        fldnum = @_rand_int(Int, 100, 1)
        meta = mk_test_meta(fldnum, :int64)
        meta.ordered[1].occurrence = 2
        meta.ordered[1].packed = true
        writeproto(pb, testval, meta)
        readproto(pb, readval, meta)
        assert_equal(testval, readval)
    end

    print_hdr("repeated string")
    for idx in 1:100
        testval.val = [randstring(5) for i in 1:10]
        readval.val = AbstractString[]
        fldnum = @_rand_int(Int, 100, 1)
        meta = mk_test_meta(fldnum, :string)
        meta.ordered[1].occurrence = 2
        writeproto(pb, testval, meta)
        readproto(pb, readval, meta)
        assert_equal(testval, readval)
    end
end

function test_optional()
    print_hdr("optional fields")
    pb = PipeBuffer()
    testval = TestOptional(TestStr(""), TestStr(""), Int64[1,2,3])
    readval = TestOptional(TestStr(""), TestStr(""), Int64[])

    for idx in 1:100
        testval.sVal1 = TestStr(string(@_rand_int(Int, 100, 0)))
        testval.sVal2 = TestStr(randstring(5))
        testval.iVal2 = Int64[@_rand_int(Int,100,0) for i in 1:10]
        sVal1Opt = rand(Bool)
        sVal2Opt = rand(Bool)
        meta = mk_test_optional_meta(sVal1Opt, sVal2Opt)
        fillunset(testval)
        fillset(testval, :iVal2)
        !sVal1Opt && fillset(testval, :sVal1)
        !sVal2Opt && fillset(testval, :sVal2)

        writeproto(pb, testval, meta)
        readval.iVal2 = Int64[]
        readproto(pb, readval, meta)

        assert_equal(testval, readval)
    end
end

function test_nested()
    print_hdr("nested types")
    pb = PipeBuffer()

    testfld1 = TestType(0)
    readfld1 = TestType(0)
    testfld2 = TestOptional(TestStr("1"), TestStr(""), Int64[1,2,3])
    readfld2 = TestOptional(TestStr("1"), TestStr(""), Int64[])
    testval = TestNested(testfld1, testfld2, [TestStr("hello"), TestStr("world")])
    readval = TestNested(readfld1, readfld2, TestStr[])

    for idx in 1:100
        testfld1.val = @_rand_int(Int64, 10^9, 0)
        testfld2.sVal1 = TestStr(string(@_rand_int(Int, 100, 0)))
        testfld2.sVal2 = TestStr(randstring(5))
        testfld2.iVal2 = Int64[@_rand_int(Int, 100, 0) for i in 1:10]

        o1 = rand(Bool)
        o2 = rand(Bool)
        o21 = rand(Bool)
        o22 = rand(Bool)
        meta = mk_test_nested_meta(o1, o2, o21, o22)

        fillunset(testval)
        fillset(testval, :fld3)
        !o1 && fillset(testval, :fld1)
        !o2 && fillset(testval, :fld2)

        fillunset(testfld2)
        fillset(testfld2, :iVal2)
        !o21 && fillset(testfld2, :sVal1)
        !o22 && fillset(testfld2, :sVal2)

        writeproto(pb, testval, meta)
        readfld2.iVal2 = Int64[]
        readval.fld3 = TestStr[]
        readproto(pb, readval, meta)

        assert_equal(testval, readval)
    end
end

function test_defaults()
    print_hdr("default values")
    pb = PipeBuffer()

    testval = TestDefaults()
    readval = TestDefaults()
    testval.iVal1 = @_rand_int(Int, 100, 0)
    writeproto(pb, testval)
    readproto(pb, readval)

    assert_equal(TestDefaults(testval.iVal1, "", [1,2,3]), readval)
end

function test_oneofs()
    print_hdr("oneofs")
    testval = TestOneofs()
    @test isfilled(testval)
    @test isfilled(testval, :iVal1)
    @test !isfilled(testval, :iVal2)
    @test isfilled(testval, :iVal3)
    @test which_oneof(testval, :optval) === :iVal3

    testval.iVal2 = 10
    @test isfilled(testval, :iVal1)
    @test isfilled(testval, :iVal2)
    @test !isfilled(testval, :iVal3)
    @test which_oneof(testval, :optval) === :iVal2

    testval.iVal1 = 10
    @test isfilled(testval, :iVal1)
    @test isfilled(testval, :iVal2)
    @test !isfilled(testval, :iVal3)
    @test which_oneof(testval, :optval) === :iVal2

    testval.iVal3 = 10
    @test isfilled(testval, :iVal1)
    @test !isfilled(testval, :iVal2)
    @test isfilled(testval, :iVal3)
    @test which_oneof(testval, :optval) === :iVal3
end

function test_maps()
    print_hdr("maps")
    pb = PipeBuffer()

    testval = TestMaps()
    readval = TestMaps()
    writeproto(pb, testval)
    readproto(pb, readval)
    assert_equal(testval, readval)

    testval = TestMaps()
    readval = TestMaps()
    testval.d1 = Dict{Int,Int}()
    writeproto(pb, testval)
    readproto(pb, readval)
    @test !isfilled(readval, :d1)

    testval = TestMaps()
    readval = TestMaps()
    testval.d2 = Dict{Int32,String}()
    writeproto(pb, testval)
    readproto(pb, readval)
    @test !isfilled(readval, :d2)

    testval = TestMaps()
    readval = TestMaps()
    testval.d3 = Dict{String,String}()
    writeproto(pb, testval)
    readproto(pb, readval)
    @test !isfilled(readval, :d3)

    testval = TestMaps()
    readval = TestMaps()
    testval.d1 = Dict{Int,Int}()
    testval.d1[1] = 1
    testval.d1[2] = 2
    writeproto(pb, testval)
    readproto(pb, readval)
    @test isfilled(readval, :d1)
    assert_equal(testval, readval)

    testval = TestMaps()
    readval = TestMaps()
    testval.d2 = Dict{Int32,String}()
    testval.d2[Int32(1)] = convert(String, "One")
    testval.d2[Int32(2)] = convert(String, "Two")
    writeproto(pb, testval)
    readproto(pb, readval)
    @test isfilled(readval, :d2)
    assert_equal(testval, readval)

    testval = TestMaps()
    readval = TestMaps()
    testval.d3 = Dict{String,String}()
    testval.d3["1"] = "One"
    testval.d3["2"] = "Two"
    writeproto(pb, testval)
    readproto(pb, readval)
    @test isfilled(readval, :d3)
    assert_equal(testval, readval)
end

function test_misc()
    print_hdr("misc functionality")
    testfld = TestOptional(TestStr("1"), TestStr(""), Int64[1,2,3])
    readfld = TestOptional(TestStr(""), TestStr("1"), Int64[])
    copy!(readfld, testfld)
    assert_equal(readfld, testfld)

    # test add_field!
    readfld = TestOptional(TestStr("1"), TestStr(""), Int64[])
    for iVal2 in testfld.iVal2
        add_field!(readfld, :iVal2, iVal2)
    end
    assert_equal(readfld, testfld)
    @test ProtoBuf.protoisequal(readfld, testfld)

    tf = TestFilled()
    @test !isfilled(tf)
    tf.fld1 = TestType("")
    fillset(tf, :fld1)
    @test isfilled(tf)

    iob = IOBuffer()
    show(iob, meta(TestOptional))
    @test !isempty(take!(iob))
    nothing
end

function test_enums()
    print_hdr("enums")
    @test getfield(TestEnum, lookup(TestEnum, 0)) == TestEnum.UNIVERSAL
    @test getfield(TestEnum, lookup(TestEnum, 1)) == TestEnum.WEB
    @test getfield(TestEnum, lookup(TestEnum, 2)) == TestEnum.IMAGES
    @test getfield(TestEnum, lookup(TestEnum, 3)) == TestEnum.LOCAL
    @test getfield(TestEnum, lookup(TestEnum, 4)) == TestEnum.NEWS
    @test getfield(TestEnum, lookup(TestEnum, 5)) == TestEnum.PRODUCTS
    @test getfield(TestEnum, lookup(TestEnum, 6)) == TestEnum.VIDEO

    @test enumstr(TestEnum, TestEnum.LOCAL) == "LOCAL"
    @test_throws ErrorException enumstr(TestEnum, Int32(12))
end

end # module ProtoBufTestCodec

ProtoBufTestCodec.test_types()
ProtoBufTestCodec.test_enums()
ProtoBufTestCodec.test_oneofs()
ProtoBufTestCodec.test_maps()
ProtoBufTestCodec.test_repeats()
ProtoBufTestCodec.test_optional()
ProtoBufTestCodec.test_nested()
ProtoBufTestCodec.test_defaults()
ProtoBufTestCodec.test_misc()
GC.gc()
println("_metacache has $(length(ProtoBuf._metacache)) entries")
#println(ProtoBuf._metacache)
println("_fillcache has $(length(ProtoBuf._fillcache)) entries")
#println(ProtoBuf._fillcache)


text:
using RasterDataSources, URIs, Test, Dates
using RasterDataSources: rasterpath, zipurl, zipname

@testset "SRTM" begin
    zip_url = URI(scheme = "https", host = "srtm.csi.cgiar.org", path = "/wp-content/uploads/files/srtm_5x5/TIFF/srtm_02_01.zip")
    tile_index1 = CartesianIndex(1, 2)   # [y, x] order
    tile_index2 = CartesianIndex(2, 2)   # [y, x] order
    @test zipurl(SRTM; tile_index = tile_index1) == zip_url
    @test zipurl(SRTM; tile_index = tile_index1) == zip_url
    @test zipname(SRTM; tile_index = tile_index1) == "srtm_02_01.zip"

    raster_path1 = joinpath(ENV["RASTERDATASOURCES_PATH"], "SRTM", "srtm_02_01.tif")
    raster_path2 = joinpath(ENV["RASTERDATASOURCES_PATH"], "SRTM", "srtm_02_02.tif")
    @test rasterpath(SRTM; tile_index = tile_index1) == raster_path1
    @test getraster(SRTM; tile_index = tile_index1) == raster_path1
    @test isfile(raster_path1)

    lon1, lat1 = -175, 60       # Coordinates of [0, 0] pixel of tile x=2, y=1
    @test getraster(SRTM; bounds=((lon1, lon1), (lat1, lat1))) == reshape([raster_path1], 1, 1)
    lon2, lat2 = -171, 55   # Coordinates of [3000, 3000] pixel of tile x=2, y=2
    @test getraster(SRTM; bounds=((lon1, lon2), (lat2, lat1))) == permutedims([raster_path1 raster_path2])
end

max_stars_repo_name: cesaraustralia/RasterDataSources

max_issues_repo_name: cesaraustralia/RasterDataSources

max_forks_repo_name: cesaraustralia/RasterDataSources