More efficient cubic functions (#233)

* More efficient cubic functions

* fix int64

* add test

Author: joaquimg

src/MOI_wrapper.jl

@@ -600,7 +600,7 @@ function MOI.supports(
     # needing to support names for both quadratic and affine constraints.
     # TODO:
     # switch to only check support name for the case of linear
-    # is a solver does not support quadratic constraints it will fain in add_
+    # is a solver does not support quadratic constraints it will fail in add_
     if MOI.supports_constraint(model.optimizer, F, S)
         return MOI.supports(model.optimizer, attr, MOI.ConstraintIndex{F,S}) &&
                MOI.supports(model.optimizer, attr, MOI.ConstraintIndex{G,S})
@@ -1910,7 +1910,7 @@ end
 #
 
 function MOI.optimize!(model::Optimizer)
-    if !isempty(model.updated_parameters)
+    if any(!isnan, values(model.updated_parameters))
         MOI.Utilities.final_touch(model, nothing)
         update_parameters!(model)
     end

src/cubic_objective.jl

@@ -136,40 +136,62 @@ function _try_incremental_cubic_update!(
     # Get the current objective function type from the inner optimizer
     F = MOI.get(model.optimizer, MOI.ObjectiveFunctionType())
 
-    # Compute full new values (not deltas) for robustness
-    # The delta was used to detect changes; now apply full new coefficients
-    new_quad_terms = _parametric_quadratic_terms(model, pf)
-    new_affine_terms = _parametric_affine_terms(model, pf)
-    new_constant = _parametric_constant(model, pf)
-
-    # Apply quadratic coefficient changes
+    # Apply quadratic coefficient changes.
+    # For each changed (var1, var2) pair, recompute its new coefficient from the
+    # base vv (pf.quadratic_data) data plus current pvv contributions
+    # (avoids full copy + full iteration).
     # MOI convention:
     #   - Off-diagonal (v1 != v2): coefficient C means C*v1*v2 (use as-is)
     #   - Diagonal (v1 == v2): coefficient C means (C/2)*v1^2 (multiply by 2)
-    for ((var1, var2), _) in delta_quadratic
-        new_coef = new_quad_terms[(var1, var2)]
-        # Apply MOI coefficient convention
-        moi_coef = var1 == var2 ? new_coef * 2 : new_coef
+    for (var1, var2) in keys(delta_quadratic)
+        new_coef = get(pf.quadratic_data, (var1, var2), zero(T))
+        for term in pf.pvv
+            p = term.index_1
+            first_is_greater = term.index_2.value > term.index_3.value
+            v1 = ifelse(first_is_greater, term.index_3, term.index_2)
+            v2 = ifelse(first_is_greater, term.index_2, term.index_3)
+            if (v1, v2) == (var1, var2)
+                new_coef +=
+                    term.coefficient * _effective_param_value(model, p_idx(p))
+            end
+        end
+        moi_coef = new_coef * ifelse(var1 == var2, 2, 1)
         MOI.modify(
             model.optimizer,
             MOI.ObjectiveFunction{F}(),
             MOI.ScalarQuadraticCoefficientChange(var1, var2, moi_coef),
         )
     end
 
-    # Apply affine coefficient changes (use full new coefficient)
-    for (var, _) in delta_affine
-        new_coef = new_affine_terms[var]
+    # Apply affine coefficient changes.
+    # For each changed variable, recompute its new coefficient from the base
+    # affine_data plus current pv and ppv contributions.
+    for var in keys(delta_affine)
+        new_coef = get(pf.affine_data, var, zero(T))
+        for term in pf.pv
+            if term.variable_2 == var
+                new_coef +=
+                    term.coefficient *
+                    _effective_param_value(model, p_idx(term.variable_1))
+            end
+        end
+        for term in pf.ppv
+            if term.index_3 == var
+                p1_val = _effective_param_value(model, p_idx(term.index_1))
+                p2_val = _effective_param_value(model, p_idx(term.index_2))
+                new_coef += term.coefficient * p1_val * p2_val
+            end
+        end
         MOI.modify(
             model.optimizer,
             MOI.ObjectiveFunction{F}(),
             MOI.ScalarCoefficientChange(var, new_coef),
         )
     end
 
-    # Apply constant change
+    # Apply constant change using the tracked current_constant (no full recompute).
     if !iszero(delta_constant)
-        pf.current_constant = new_constant
+        pf.current_constant += delta_constant
         MOI.modify(
             model.optimizer,
             MOI.ObjectiveFunction{F}(),

src/cubic_parser.jl

@@ -295,24 +295,92 @@ function _expand_power(args, ::Type{T}) where {T}
     return result
 end
 
+"""
+    _sort3(a, b, c) -> (a, b, c) sorted ascending
+
+Sort three integers without heap allocation (in-place bubble sort).
+"""
+function _sort3(a::Real, b::Real, c::Real)
+    if a > b
+        a, b = b, a
+    end
+    if b > c
+        b, c = c, b
+    end
+    if a > b
+        a, b = b, a
+    end
+    return a, b, c
+end
+
+"""
+    _monomial_key(m::_Monomial)::NTuple{4,Int64}
+
+Compute a canonical hash key for a monomial: (degree, sorted_val1, sorted_val2, sorted_val3).
+Uses integer tuple instead of a sorted Vector for faster hashing.
+"""
+function _monomial_key(m::_Monomial)
+    n = length(m.variables)
+    if n == 0
+        return (Int64(0), Int64(0), Int64(0), Int64(0))
+    elseif n == 1
+        a = m.variables[1].value
+        return (Int64(1), Int64(a), Int64(0), Int64(0))
+    elseif n == 2
+        a, b = m.variables[1].value, m.variables[2].value
+        lo, hi = a <= b ? (a, b) : (b, a)
+        return (Int64(2), Int64(lo), Int64(hi), Int64(0))
+    else  # n >= 3; degree > 3 is rejected at classification stage
+        a, b, c = _sort3(
+            m.variables[1].value,
+            m.variables[2].value,
+            m.variables[3].value,
+        )
+        return (Int64(n), Int64(a), Int64(b), Int64(c))
+    end
+end
+
+"""
+    _monomial_vars(key::NTuple{4,Int64})::Vector{MOI.VariableIndex}
+
+Given a monomial key, reconstruct the list of variables.
+"""
+function _monomial_vars(key::NTuple{4,Int64})
+    degree = key[1]
+    if degree == 0
+        return MOI.VariableIndex[]
+    elseif degree == 1
+        return [MOI.VariableIndex(key[2])]
+    elseif degree == 2
+        return [MOI.VariableIndex(key[2]), MOI.VariableIndex(key[3])]
+    else  # degree == 3
+        return [
+            MOI.VariableIndex(key[2]),
+            MOI.VariableIndex(key[3]),
+            MOI.VariableIndex(key[4]),
+        ]
+    end
+end
+
 """
     _combine_like_monomials(monomials::Vector{_Monomial{T}}) where {T}
 
 Combine like monomials (same variables, regardless of order).
+Assumes all monomials have degree ≤ 3.
 """
 function _combine_like_monomials(monomials::Vector{_Monomial{T}}) where {T}
-    # Use a dict keyed by sorted variable tuple
-    combined = Dict{Vector{MOI.VariableIndex},T}()
+    # Key: NTuple{4,Int64} (degree + up to 3 sorted variable indices).
+    combined = Dict{NTuple{4,Int64},T}()
 
     for m in monomials
-        # Sort variables for canonical key
-        key = sort(m.variables, by = v -> v.value)
+        key = _monomial_key(m)
         combined[key] = get(combined, key, zero(T)) + m.coefficient
     end
 
     result = _Monomial{T}[]
-    for (vars, coef) in combined
+    for (key, coef) in combined
         if !iszero(coef)
+            vars = _monomial_vars(key)
             push!(result, _Monomial{T}(coef, vars))
         end
     end
@@ -341,7 +409,7 @@ function _classify_monomial(m::_Monomial)
         else
             return :pp
         end
-    elseif degree == 3
+    else  # degree == 3 (degree > 3 rejected early in _parse_cubic_expression)
         if num_params == 0
             return :vvv  # Invalid - no parameter
         elseif num_params == 1
@@ -351,8 +419,6 @@ function _classify_monomial(m::_Monomial)
         else
             return :ppp
         end
-    else
-        return :invalid  # Degree > 3
     end
 end
 
@@ -377,12 +443,17 @@ function _parse_cubic_expression(
         return nothing
     end
 
+    # Reject any monomial with degree > 3 before combining
+    for m in monomials
+        if _monomial_degree(m) > 3
+            return nothing
+        end
+    end
+
     # Combine like terms
     monomials = _combine_like_monomials(monomials)
 
     # Classify and collect terms
-    cubic_terms = _ScalarCubicTerm{T}[]
-
     cubic_ppp = _ScalarCubicTerm{T}[]
     cubic_ppv = _ScalarCubicTerm{T}[]
     cubic_pvv = _ScalarCubicTerm{T}[]
@@ -399,8 +470,8 @@ function _parse_cubic_expression(
     for m in monomials
         classification = _classify_monomial(m)
 
-        if classification == :invalid || classification == :vvv
-            return nothing  # Invalid degree or no parameter in cubic
+        if classification == :vvv
+            return nothing  # No parameter in cubic term
         elseif classification == :constant
             constant += m.coefficient
         elseif classification == :v

src/cubic_types.jl

@@ -40,17 +40,27 @@ function _normalize_cubic_indices(
     idx2::MOI.VariableIndex,
     idx3::MOI.VariableIndex,
 )
-    params = MOI.VariableIndex[]
-    vars = MOI.VariableIndex[]
-    for idx in (idx1, idx2, idx3)
-        if _is_parameter(idx)
-            push!(params, idx)
-        else
-            push!(vars, idx)
-        end
+    p1 = _is_parameter(idx1)
+    p2 = _is_parameter(idx2)
+    p3 = _is_parameter(idx3)
+    # Place parameters before variables, preserving relative order within each group.
+    if p1 && p2 && p3
+        return idx1, idx2, idx3  # ppp — already ordered
+    elseif p1 && p2             # p p v
+        return idx1, idx2, idx3
+    elseif p1 && p3             # p v p → p p v
+        return idx1, idx3, idx2
+    elseif p2 && p3             # v p p → p p v
+        return idx2, idx3, idx1
+    elseif p1                   # p v v — already ordered
+        return idx1, idx2, idx3
+    elseif p2                   # v p v → p v v
+        return idx2, idx1, idx3
+    elseif p3                   # v v p → p v v
+        return idx3, idx1, idx2
+    else                        # v v v — no parameter (caller validates)
+        return idx1, idx2, idx3
     end
-    all_indices = vcat(params, vars)
-    return all_indices[1], all_indices[2], all_indices[3]
 end
 
 """

src/parametric_cubic_function.jl

@@ -143,12 +143,9 @@ cubic_parameter_parameter_parameter_terms(f::ParametricCubicFunction) = f.ppp
 
 Get the effective parameter value: updated value if available, otherwise current value.
 """
-function _effective_param_value(model, pi::ParameterIndex)
-    if haskey(model.updated_parameters, pi) &&
-       !isnan(model.updated_parameters[pi])
-        return model.updated_parameters[pi]
-    end
-    return model.parameters[pi]
+function _effective_param_value(model, p::ParameterIndex)
+    val = model.updated_parameters[p]
+    return isnan(val) ? model.parameters[p] : val
 end
 
 """