OptimizationProblems.jl/src/ADNLPProblems/variational.jl at 3eef3eb2cd6259c2985f8a4074b7d771635529e7 · arnavk23/OptimizationProblems.jl · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
export variational

function variational(; use_nls::Bool = false, kwargs...)
  model = use_nls ? :nls : :nlp
  return variational(Val(model); kwargs...)
end

function variational(::Val{:nlp}; n::Int = default_nvar, type::Type{T} = Float64, kwargs...) where {T}
  h = one(T) / T(n + 1)

  function f(x; n = length(x))
    S = eltype(x)
    hS = convert(S, h)

    xext = Vector{S}(undef, n + 2)
    xext[1] = zero(S)
    @inbounds for i = 1:n
      xext[i + 1] = x[i]
    end
    xext[n + 2] = zero(S)

    term1 = zero(S)
    @inbounds for i = 1:n
      xi = x[i]
      xip1 = xext[i + 2]
      term1 += xi * (xi - xip1) / hS
    end

    function diffquot(a, b)
      d = b - a
      if d == zero(d)
        return exp(a)
      else
        return (exp(b) - exp(a)) / d
      end
    end

    term2 = zero(S)
    @inbounds for j = 0:n
      a = xext[j + 1]
      b = xext[j + 2]
      term2 += diffquot(a, b)
    end
    return 2 * (term1 + n * (hS / 2) * term2)
  end

  x0 = Vector{T}(undef, n)
  for i = 1:n
    x0[i] = T(i) * h * (one(T) - T(i) * h)
  end

  return ADNLPModels.ADNLPModel(f, x0, name = "variational", minimize = true; kwargs...)
end

function variational(::Val{:nls}; n::Int = default_nvar, type::Type{T} = Float64, kwargs...) where {T}
  h = one(T) / T(n + 1)

  function diffquot(a, b)
    d = b - a
    if d == zero(d)
      return exp(a)
    else
      return (exp(b) - exp(a)) / d
    end
  end

  function F!(r, x)
    S = eltype(x)
    hS = convert(S, h)

    xext = Vector{S}(undef, n + 2)
    xext[1] = zero(S)
    @inbounds for i = 1:n
      xext[i + 1] = x[i]
    end
    xext[n + 2] = zero(S)

    term1 = zero(S)
    @inbounds for i = 1:n
      xi = x[i]
      xip1 = xext[i + 2]
      term1 += xi * (xi - xip1) / hS
    end

    term2 = zero(S)
    @inbounds for j = 0:n
      a = xext[j + 1]
      b = xext[j + 2]
      term2 += diffquot(a, b)
    end

    Fval = 2 * (term1 + n * (hS / 2) * term2)
    r[1] = sqrt(Fval < zero(Fval) ? zero(Fval) : Fval)
    return r
  end

  x0 = Vector{T}(undef, n)
  for i = 1:n
    x0[i] = T(i) * h * (one(T) - T(i) * h)
  end

  return ADNLPModels.ADNLSModel!(F!, x0, 1, name = "variational-nls"; kwargs...)
end