text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id$
#*--------------------------------------------------------------------
type Heatex_Basic
	Heatex_Basic()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of Components"
			)),
			fill(molweight(Dict{Symbol,Any}(
				:Brief=>"Component Mol Weight",
				:Hidden=>true
			)),(NComp)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet Hot Stream",
				:PosX=>0,
				:PosY=>0.508,
				:Symbol=>"^{inHot}"
			)),
			streamPH(Dict{Symbol,Any}(
				:Brief=>"Outlet Hot Stream",
				:PosX=>1,
				:PosY=>0.508,
				:Symbol=>"^{outHot}"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet Cold Stream",
				:PosX=>0.50,
				:PosY=>1,
				:Symbol=>"^{inCold}"
			)),
			streamPH(Dict{Symbol,Any}(
				:Brief=>"Outlet Cold Stream",
				:PosX=>0.50,
				:PosY=>0,
				:Symbol=>"^{outCold}"
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Exchange Surface Area"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Duty",
				:Default=>7000,
				:Lower=>1e-6,
				:Upper=>1e10
			)),
			heat_trans_coeff(Dict{Symbol,Any}(
				:Brief=>"Overall Heat Transfer Coefficient",
				:Default=>1,
				:Lower=>1e-6,
				:Upper=>1e10
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Pressure Drop Hot Side",
				:Default=>0.01,
				:Lower=>0,
				:DisplayUnit=>"kPa" ,
				:Symbol=>"\\Delta P_{hot}"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Pressure Drop Cold Side",
				:Default=>0.01,
				:Lower=>0,
				:DisplayUnit=>"kPa" ,
				:Symbol=>"\\Delta P_{cold}"
			)),
			[
				:(Q = InletHot.F*(InletHot.h-OutletHot.h)),
				:(Q =-InletCold.F*(InletCold.h-OutletCold.h)),
				:(InletHot.F = OutletHot.F),
				:(InletCold.F = OutletCold.F),
				:(OutletHot.z = InletHot.z),
				:(OutletCold.z = InletCold.z),
				:(OutletHot.P = InletHot.P - PdropHotSide),
				:(OutletCold.P = InletCold.P - PdropColdSide),
			],
			[
				"Energy Balance Hot Stream","Energy Balance Cold Stream","Molar Balance Hot Stream","Molar Balance Cold Stream","Hot Stream Molar Fraction Constraint","Cold Stream Molar Fraction Constraint","Pressure Drop Hot Stream","Pressure Drop Cold Stream",
			],
			[:PP,:NComp,:M,],
			[:InletHot,:OutletHot,:InletCold,:OutletCold,:A,:Q,:U,:PdropHotSide,:PdropColdSide,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	M::Array{molweight}
	InletHot::stream
	OutletHot::streamPH
	InletCold::stream
	OutletCold::streamPH
	A::area
	Q::power
	U::heat_trans_coeff
	PdropHotSide::press_delta
	PdropColdSide::press_delta
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Heatex_Basic
function set(in::Heatex_Basic)
	#"Component Molecular Weight"
	M = PP.MolecularWeight()

end
function setEquationFlow(in::Heatex_Basic)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
end
function atributes(in::Heatex_Basic,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="Basic Model for Simplified Heat Exchangers"
	fields[:Info]="Model of a simplified heat exchanger.
This model perform only material and heat balance.

== Assumptions ==
* Steady-State operation;
* No heat loss to the surroundings.

== Specify ==
* The Inlet streams: Hot and Cold;
"
	drive!(fields,_)
	return fields
end
Heatex_Basic(_::Dict{Symbol,Any})=begin
	newModel=Heatex_Basic()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using Test, StaticArrays, DelayEmbeddings

println("\nTesting Dataset...")

@testset "Dataset" begin
  data = Dataset(rand(1001,3))
  @testset "Methods & Indexing" begin
    a = data[:, 1]
    b = data[:, 2]
    c = data[:, 3]

    @test Dataset(a, b, c) == data
    @test size(Dataset(a, b)) == (1001, 2)

    @test data[:, 2:3][:, 1] == data[:, 2]

    @test size(data[1:10,1:2]) == (10,2)
    @test data[1:10,1:2] == Dataset(a[1:10], b[1:10])
    @test data[SVector{10}(1:10), SVector(1, 2)] == data[1:10, 1:2]
  end

  @testset "minmax" begin
    mi = minima(data)
    ma = maxima(data)
    mimi, mama = minmaxima(data)
    @test mimi == mi
    @test mama == ma
    for i in 1:3
      @test mi[i] < ma[i]
    end
  end

  @testset "Conversions" begin
    m = Matrix(data)
    @test Dataset(m) == data

    m = rand(1000, 4)
    @test Matrix(Dataset(m)) == m
  end

end

ext: jl


text:
####
#### Sampling
####
mutable struct Sampling{T<:Function} <: Layer
    options::NamedTuple
end

Sampling{T}(;window=2, padding=0, stride=window, mode=0, maxpoolingNanOpt=0, alpha=1) where T <: Function =
    Sampling{T}((window=window, padding=padding, stride=stride, mode=mode, maxpoolingNanOpt=maxpoolingNanOpt, alpha=alpha))

@inline (m::Sampling{typeof(pool)})(x)   =  pool(x;m.options...)
@inline (m::Sampling{typeof(unpool)})(x) =  unpool(x;m.options...)

Base.show(io::IO,m::Sampling{typeof(pool)})   = print(io,"Pool",m.options)
Base.show(io::IO,m::Sampling{typeof(unpool)}) = print(io,"UnPool",m.options)
"""
    Pool(keywords...)
Creates a pooling layer, `Sampling{typeof(pool)}`, according to given keyword arguments.

    (::Sampling{typeof(pool)})(x)
Compute pooling of input values (i.e., the maximum or average of several adjacent values)
to produce an output with smaller height and/or width.

Currently 4 or 5 dimensional KnetArrays with Float32 or Float64 entries are
supported. If x has dimensions (X1,X2,...,I,N), the result y will have dimensions
(Y1,Y2,...,I,N) where

Yi=1+floor((Xi+2*padding[i]-window[i])/stride[i])

Here I is the number of input channels, N is the number of instances, and Xi,Yi
are spatial dimensions. window, padding and stride are keyword arguments that can
be specified as a single number (in which case they apply to all dimensions), or
an array/tuple with entries for each spatial dimension.

Keywords:

* window=2: the pooling window size for each dimension.
* padding=0: the number of extra zeros implicitly concatenated at the
start and at the end of each dimension.
* stride=window: the number of elements to slide to reach the next pooling
window.
* mode=0: 0 for max, 1 for average including padded values, 2 for average
excluding padded values.
* maxpoolingNanOpt=0: Nan numbers are not propagated if 0, they are
propagated if 1.
* alpha=1: can be used to scale the result.

"""
const Pool = Sampling{typeof(pool)}

"""
    UnPool(kwargs...)
    (::Sampling{typeof(unpool)})(x)

    Reverse of pooling. It has same kwargs with Pool. see `Pool`

    x == pool(unpool(x;o...); o...)
"""
const UnPool = Sampling{typeof(unpool)}
####
#### Filtering
####
mutable struct Filtering{T<:Function,P,A<:ActOrNothing,V<:Bias} <: Layer
    weight::P
    bias::V
    activation::A
    options::NamedTuple
end

function Filtering{T}(;height::Integer, width::Integer, inout::Pair=1=>1,
                       activation::ActOrNothing=NonAct(),
                       winit=xavier, binit=zeros,
                       atype=arrtype,
                       opts...) where T <: Function

    wsize = T===typeof(conv4) ? inout : reverse(inout)
    w = param(height,width,wsize...; init=winit, atype=atype)
    b = binit !== nothing ? Bias(1,1,inout[2],1; init=binit, atype=atype) : Bias(nothing)
    Filtering{T}(w, b, activation; opts...)

end

Filtering{T}(w, b, activation; stride=1, padding=0, mode=0, dilation=1, alpha=1) where T =
    Filtering{T,typeof(w),typeof(activation),typeof(b)}(w, b, activation, (stride=stride, dilation=dilation, mode=mode, alpha=alpha, padding=padding))

@inline (m::Filtering{typeof(conv4)})(x) =
     unmake4D(postConv(m, conv4(m.weight, make4D(x); m.options...)), ndims(x))

Base.show(io::IO,m::Filtering{typeof(conv4),P,A,V}) where {P,A,V} =
    print(io,"Conv{$P,$A,$V}",m.options)

@inline (m::Filtering{typeof(deconv4)})(x) =
    unmake4D(postConv(m, deconv4(m.weight, make4D(x); m.options...)), ndims(x))

Base.show(io::IO,m::Filtering{typeof(deconv4),P,A,V}) where {P,A,V} =
    print(io,"DeConv{$P,$A,$V}",m.options)

"""
    Conv(;height=filterHeight, width=filterWidth, inout = 1 => 1, kwargs...)

Creates and convolutional layer `Filtering{typeof(conv4)}` according to given filter dimensions.

    (m::Filtering{typeof(conv4)})(x) #forward run

If `m.w` has dimensions `(W1,W2,...,I,O)` and
`x` has dimensions `(X1,X2,...,I,N)`, the result `y` will have
dimensions `(Y1,Y2,...,O,N)` where

    Yi=1+floor((Xi+2*padding[i]-Wi)/stride[i])

Here `I` is the number of input channels, `O` is the number of output channels, `N` is the number of instances,
and `Wi,Xi,Yi` are spatial dimensions. `padding` and `stride` are
keyword arguments that can be specified as a single number (in which case they apply to all dimensions),
or an tuple with entries for each spatial dimension.

# Keywords
* `inout=input_channels => output_channels`
* `activation=NonAct()`: nonlinear function applied after convolution, default is identity.
* `pool=nothing`: Pooling layer or window size of pooling
* `winit=xavier`: weight initialization distribution
* `bias=zeros`: bias initialization distribution
* `padding=0`: the number of extra zeros implicitly concatenated at the start and at the end of each dimension.
* `stride=1`: the number of elements to slide to reach the next filtering window.
* `dilation=1`: dilation factor for each dimension.
* `mode=0`: 0 for convolution and 1 for cross-correlation.
* `alpha=1`: can be used to scale the result.
* `handle`: handle to a previously created cuDNN context. Defaults to a Knet allocated handle.

"""
const Conv = Filtering{typeof(conv4)}

"""
    DeConv(;height=filterHeight, width=filterWidth, inout=1=>1, kwargs...)

Creates and deconvolutional layer `Filtering{typeof(deconv4)}`  according to given filter dimensions.

    (m::Filtering{typeof(deconv4)})(x) #forward run

If `m.w` has dimensions `(W1,W2,...,I,O)` and
`x` has dimensions `(X1,X2,...,I,N)`, the result `y` will have
dimensions `(Y1,Y2,...,O,N)` where

    Yi = Wi+stride[i](Xi-1)-2padding[i]

Here `I` is the number of input channels, `O` is the number of output channels, `N` is the number of instances,
and `Wi,Xi,Yi` are spatial dimensions. `padding` and `stride` are
keyword arguments that can be specified as a single number (in which case they apply to all dimensions),
or an tuple with entries for each spatial dimension.

# Keywords
* `inout=input_channels => output_channels`
* `activation=identity`: nonlinear function applied after convolution
* `unpool=nothing`: Unpooling layer or window size of unpooling
* `winit=xavier`: weight initialization distribution
* `bias=zeros`: bias initialization distribution
* `padding=0`: the number of extra zeros implicitly concatenated at the start and at the end of each dimension.
* `stride=1`: the number of elements to slide to reach the next filtering window.
* `dilation=1`: dilation factor for each dimension.
* `mode=0`: 0 for convolution and 1 for cross-correlation.
* `alpha=1`: can be used to scale the result.
* `handle`: handle to a previously created cuDNN context. Defaults to a Knet allocated handle.
"""
const DeConv = Filtering{typeof(deconv4)}

###
### Utils
###
@inline function make4D(x)
    n = ndims(x)
    @assert n < 5 "filtering layers currently supports 4 dimensional arrays"
    n == 4 ? x : reshape(x,size(x)...,ntuple(x->1, 4-n)...)
end

@inline unmake4D(y,dims::Int) = dims>3 ? y : reshape(y,size(y)[1:dims])

@inline postConv(m::Filtering{<:Any,<:Any,<:Activation,<:Bias}, y) = m.activation(m.bias(y))

@inline postConv(m::Filtering{<:Any,<:Any,<:Activation,<:Nothing}, y) = m.activation(y)

@inline postConv(m::Filtering{<:Any,<:Any,<:Nothing,<:Bias}, y) = m.bias(y)


text:
using XDiff

ex = quote
    y[i] = W[i,k] * x[k]
    v[i] = y[i] - x[i]
    z = v[i] * I[i]
end
inputs = [:W => rand(2, 2), :x => rand(2)]
ctx = Dict()
@time dexs = rdiff(ex; ctx=ctx, inputs...)


text:
# tests to check the parts of the whole optimization that includes
# optTopLevel (search in the same level of hybrids) and
# addHybridizationUpdate(net) that will add a new hybrid if needed, and deleteHybridizationUpdate()


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct TapTransformer <: ACBranch
        name::String
        available::Bool
        activepower_flow::Float64
        reactivepower_flow::Float64
        arc::Arc
        r::Float64
        x::Float64
        primaryshunt::Float64
        tap::Float64
        rate::Union{Nothing, Float64}
        services::Vector{Service}
        ext::Dict{String, Any}
        forecasts::InfrastructureSystems.Forecasts
        internal::InfrastructureSystemsInternal
    end

# Arguments
- `name::String`
- `available::Bool`
- `activepower_flow::Float64`
- `reactivepower_flow::Float64`
- `arc::Arc`
- `r::Float64`: System per-unit value, validation range: (-2, 2), action if invalid: error
- `x::Float64`: System per-unit value, validation range: (-2, 4), action if invalid: error
- `primaryshunt::Float64`: System per-unit value, validation range: (0, 2), action if invalid: error
- `tap::Float64`, validation range: (0, 2), action if invalid: error
- `rate::Union{Nothing, Float64}`, validation range: (0, nothing), action if invalid: error
- `services::Vector{Service}`: Services that this device contributes to
- `ext::Dict{String, Any}`
- `forecasts::InfrastructureSystems.Forecasts`: internal forecast storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct TapTransformer <: ACBranch
    name::String
    available::Bool
    activepower_flow::Float64
    reactivepower_flow::Float64
    arc::Arc
    "System per-unit value"
    r::Float64
    "System per-unit value"
    x::Float64
    "System per-unit value"
    primaryshunt::Float64
    tap::Float64
    rate::Union{Nothing, Float64}
    "Services that this device contributes to"
    services::Vector{Service}
    ext::Dict{String, Any}
    "internal forecast storage"
    forecasts::InfrastructureSystems.Forecasts
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function TapTransformer(name, available, activepower_flow, reactivepower_flow, arc, r, x, primaryshunt, tap, rate, services=Device[], ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), )
    TapTransformer(name, available, activepower_flow, reactivepower_flow, arc, r, x, primaryshunt, tap, rate, services, ext, forecasts, InfrastructureSystemsInternal(), )
end

function TapTransformer(; name, available, activepower_flow, reactivepower_flow, arc, r, x, primaryshunt, tap, rate, services=Device[], ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), )
    TapTransformer(name, available, activepower_flow, reactivepower_flow, arc, r, x, primaryshunt, tap, rate, services, ext, forecasts, )
end

# Constructor for demo purposes; non-functional.
function TapTransformer(::Nothing)
    TapTransformer(;
        name="init",
        available=false,
        activepower_flow=0.0,
        reactivepower_flow=0.0,
        arc=Arc(Bus(nothing), Bus(nothing)),
        r=0.0,
        x=0.0,
        primaryshunt=0.0,
        tap=1.0,
        rate=0.0,
        services=Device[],
        ext=Dict{String, Any}(),
        forecasts=InfrastructureSystems.Forecasts(),
    )
end

InfrastructureSystems.get_name(value::TapTransformer) = value.name
"""Get TapTransformer available."""
get_available(value::TapTransformer) = value.available
"""Get TapTransformer activepower_flow."""
get_activepower_flow(value::TapTransformer) = value.activepower_flow
"""Get TapTransformer reactivepower_flow."""
get_reactivepower_flow(value::TapTransformer) = value.reactivepower_flow
"""Get TapTransformer arc."""
get_arc(value::TapTransformer) = value.arc
"""Get TapTransformer r."""
get_r(value::TapTransformer) = value.r
"""Get TapTransformer x."""
get_x(value::TapTransformer) = value.x
"""Get TapTransformer primaryshunt."""
get_primaryshunt(value::TapTransformer) = value.primaryshunt
"""Get TapTransformer tap."""
get_tap(value::TapTransformer) = value.tap
"""Get TapTransformer rate."""
get_rate(value::TapTransformer) = value.rate
"""Get TapTransformer services."""
get_services(value::TapTransformer) = value.services
"""Get TapTransformer ext."""
get_ext(value::TapTransformer) = value.ext

InfrastructureSystems.get_forecasts(value::TapTransformer) = value.forecasts
"""Get TapTransformer internal."""
get_internal(value::TapTransformer) = value.internal


text:
using Revise
using JuMP
using Gurobi
using Potato

path = joinpath(@__DIR__,"abs1n20.xlsx")
extract!(path;f="euclidean")

path = joinpath(DEPOT_PATH[1],"dev","Potato\\test\\B.xlsx")
extract!(path;f="haversine")

GUROBI_ENV = Gurobi.Env()
set_optimizer!(optimizer_with_attributes(() -> Gurobi.Optimizer(GUROBI_ENV)))
set_slack_coeff!(1000.0)
set_surp_coeff!(-1000.0)

prioritize_vertex!()
prioritize_vehicle!()

Threads.nthreads()
test_root = root()
@time buildSub!(test_root)
test_master = master(test_root)
test_duals = getDuals(test_master)
@time test_sub = sub(test_root,test_duals)
@time test_cols = getCols(test_sub)

push!(test_root.columns,test_cols)

@time colGen(test_root,maxCG=Inf,track=true)

ori = origin(test_root)

println(ori.v[:,82,1])

R = Dict(1:length(test_root.columns) .=> test_root.columns)
θ = value.(master(test_root).obj_dict[:θ])

rounded = round.(θ)

o = JuMP.Containers.DenseAxisArray{Float64}(undef,V(),K(),T())
u = JuMP.Containers.DenseAxisArray{Float64}(undef,V(),K(),T())
v = JuMP.Containers.DenseAxisArray{Float64}(undef,V(),K(),T())
x = JuMP.Containers.DenseAxisArray{Float64}(undef,V(),V(),K(),T())
l = JuMP.Containers.DenseAxisArray{Float64}(undef,V(),V(),K(),T())

o .= 0
u .= 0
v .= 0
l .= 0
x .= 0

for k in K(), t in T(), i in K(k).cover
    o[i,k,t] = sum(R[r][(k,t)].o[i] * rounded[r,k,t] for r in keys(R))
    u[i,k,t] = sum(R[r][(k,t)].u[i] * rounded[r,k,t] for r in keys(R))
    v[i,k,t] = sum(R[r][(k,t)].v[i] * rounded[r,k,t] for r in keys(R))
end

for k in K(), t in T(), i in K(k).cover, j in K(k).cover
    x[i,j,k,t] = sum(R[r][(k,t)].x[i,j] * rounded[r,k,t] for r in keys(R))
    l[i,j,k,t] = sum(R[r][(k,t)].l[i,j] * rounded[r,k,t] for r in keys(R))
end

tes = col(u,v,l,o,x)

new = find_branch(R,θ)
println(new)
s(new,R,θ)
to_push = bound(new,"<=",0)
push!(test_root.bounds,to_push)

test_root = root()
@time res = leaf(test_root,Inf,1000.0)

res.upperBound
tes = filter(p -> p.status[end] == "NEW_BOUND",res.visited)

test_master = master(tes[end])
origin(tes[end])

traverse(1000.0)

tes_I = value.(test_master.obj_dict[:I]).data
h_i = [V(i).h for i in V()]

sum(tes_I .* h_i)
dist(1,5) + dist(5,4) + dist(4,1)

[V(i).START for i in V()]

[d(i) for i in V()]

coluna = optimizer_with_attributes(
    Coluna.Optimizer,
    "params" => Coluna.Params(solver = Coluna.Algorithm.TreeSearchAlgorithm()),
    "default_optimizer" => Gurobi.Optimizer
)

colindex = Vector{NamedTuple}()
for k in K(), t in T(), m in 1:sum(K(k).BP[i] for i in K(k).cover)
    push!(colindex,(k=k,t=t,m=m))
end

@axis(𝕂𝕋𝕄, colindex)

𝕂𝕋𝕄[2].indice.k

𝕄 = Dict()
for k in K()
    𝕄[k] = collect(1:sum(K(k).BP[i] for i in K(k).cover))
end

@axis(𝕋, T())

ori = BlockModel(coluna)
@variable(ori, I[V(), vcat(first(T()) - 1, T())])

@variable(ori, u[k = 𝕂, t = 𝕋, m = 𝕄[k], i = V()] >= 0, Int)
@variable(ori, v[k = 𝕂, t = 𝕋, m = 𝕄[k], i = V()] >= 0, Int)
@variable(ori, l[k = 𝕂, t = 𝕋, m = 𝕄[k], i = V(), j = V()] >= 0, Int)
@variable(ori, o[k = 𝕂, t = 𝕋, m = 𝕄[k], i = V()], Bin)
@variable(ori, x[k = 𝕂, t = 𝕋, m = 𝕄[k], i = V(), j = V()], Bin)

#BASIC CONSTRAINTS
@constraint(ori, [k = 𝕂, t = 𝕋, m = 𝕄[k]],
    sum(u[k,t,m,i] for i in K(k).cover) - sum(v[k,t,m,i] for i in K(k).cover) == 0
)

@constraint(ori, [k = 𝕂, t = 𝕋, m = 𝕄[k], i = K(k).cover],
    sum(l[k,t,m,j,i] for j in K(k).cover) - sum(l[k,t,m,i,j] for j in K(k).cover) ==
    u[k,t,m,i] - v[k,t,m,i]
)

@constraint(ori, [k = 𝕂, t = 𝕋, m = 𝕄[k], i = K(k).cover],
    sum(x[k,t,m,j,i] for j in K(k).cover) - sum(x[k,t,m,i,j] for j in K(k).cover) == 0
)

@constraint(ori, [k = 𝕂, t = 𝕋, m = 𝕄[k], i = K(k).cover],
    v[k,t,m,i] <= K(k).Q * o[k,t,m,i]
) #VZ

@constraint(ori, [k = 𝕂, t = 𝕋, m = 𝕄[k], i = K(k).cover, j = K(k).cover],
    l[k,t,m,i,j] <= K(k).Q * x[k,t,m,i,j]
) #XL

@constraint(ori, [k = 𝕂,t = 𝕋, m = 𝕄[k]],
    sum(o[k,t,m,i] for i in K(k).cover) <= 1
) #one start

@constraint(ori, [k = 𝕂, t = 𝕋, i = K(k).cover],
    sum(o[k,t,m,i] for m in 𝕄[k]) <= K(k).BP[i]
)

@objective(ori, Min,
    sum(
        sum(
            dist(i,j) * (
                K(k).vx * x[k,t,m,i,j] +
                K(k).vl * l[k,t,m,i,j]
            )
            for i in K(k).cover, j in K(k).cover
        ) +
        sum(
            K(k).fd * u[k,t,m,i]
            for i in K(k).cover
        ) +
        sum(
            K(k).fp * o[k,t,m,i]
            for i in K(k).cover
        )
        for k in 𝕂, t in 𝕋, m in 𝕄[k]
    ) + #column costs
    sum(
        V(i).h * I[i,t]
        for i in V(), t in T()
    )
)

@constraint(ori, [i = V(), t = T()],
    I[i,t - 1] +
    sum(u[k,t,m,i] for k in 𝕂, m in 𝕄[k]) ==
    sum(v[k,t,m,i] for k in 𝕂, m in 𝕄[k]) +
    d(i,t) + I[i,t]
)

@constraint(ori, [i = V(), t = T()],
    V(i).MIN <= I[i,t] <= V(i).MAX #inventory capacity interval
)

@constraint(ori, [i = V()],
    I[i,first(T())-1] == V(i).START #starting inventory level
)

@dantzig_wolfe_decomposition(ori, decomposition, 𝕂)

master = getmaster(decomposition)
subproblems = getsubproblems(decomposition)

for k in K()
    specify!(
        subproblems[k], lower_multiplicity = 0,
        upper_multiplicity = sum(M[k]) * length(T())
    )
end

optimize!(ori)

value.(u)

println(value.(u))
println("+++++++++++++++")


text:
export marginalrule

@marginalrule Wishart(:out_ν_S) (m_out::Wishart, m_ν::PointMass, m_S::PointMass) = begin
    return (out = prod(ProdAnalytical(), Wishart(mean(m_ν), mean(m_S)), m_out), ν = m_ν, S = m_S)
end


text:
export edges

function edges(p::Polygon{N,T}) where {N,T}
    return Vec{N, LineSegment{T}}(
                ntuple(i->LineSegment(p[ (i-1)%N + 1],
                                      p[  i   %N + 1]
                                     ), Val(N))
               )
end

function edges(p::QuadraticPolygon{N,T}) where {N,T}
    M = N ÷ 2
    return Vec{M, QuadraticSegment{T}}(
                ntuple(i->QuadraticSegment(p[ (i-1)%M + 1],
                                           p[  i   %M + 1],
                                           p[  i      + M]
                                          ), Val(M))
               )
end

function edges(p::Tetrahedron)
    v = vertices(p)
    return Vec(LineSegment(v[1], v[2]),
               LineSegment(v[2], v[3]),
               LineSegment(v[3], v[1]),
               LineSegment(v[1], v[4]),
               LineSegment(v[2], v[4]),
               LineSegment(v[3], v[4]))
end

function edges(p::Hexahedron)
    v = vertices(p)
    return Vec(LineSegment(v[1], v[2]), # lower z
               LineSegment(v[2], v[3]),
               LineSegment(v[3], v[4]),
               LineSegment(v[4], v[1]),
               LineSegment(v[5], v[6]), # upper z
               LineSegment(v[6], v[7]),
               LineSegment(v[7], v[8]),
               LineSegment(v[8], v[5]),
               LineSegment(v[1], v[5]), # lower, upper connections
               LineSegment(v[2], v[6]),
               LineSegment(v[3], v[7]),
               LineSegment(v[4], v[8]))
end

function edges(p::QuadraticTetrahedron)
    v = vertices(p)
    return Vec(QuadraticSegment(v[1], v[2], v[5]),
               QuadraticSegment(v[2], v[3], v[6]),
               QuadraticSegment(v[3], v[1], v[7]),
               QuadraticSegment(v[1], v[4], v[8]),
               QuadraticSegment(v[2], v[4], v[9]),
               QuadraticSegment(v[3], v[4], v[10]))
end

function edges(p::QuadraticHexahedron)
    v = vertices(p)
    return Vec(QuadraticSegment(v[1], v[2], v[ 9]), # lower z
               QuadraticSegment(v[2], v[3], v[10]),
               QuadraticSegment(v[3], v[4], v[11]),
               QuadraticSegment(v[4], v[1], v[12]),
               QuadraticSegment(v[5], v[6], v[13]), # upper z
               QuadraticSegment(v[6], v[7], v[14]),
               QuadraticSegment(v[7], v[8], v[15]),
               QuadraticSegment(v[8], v[5], v[16]),
               QuadraticSegment(v[1], v[5], v[17]), # lower, upper connections
               QuadraticSegment(v[2], v[6], v[18]),
               QuadraticSegment(v[3], v[7], v[19]),
               QuadraticSegment(v[4], v[8], v[20]))
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
# Gives generic rules for generating the strictly new vertices when splitting
#  an E-dimensional simplex with a size reducing factor of k. Subsimplices containing
#  one of the original vertices is not considered.
#
# splits the decomposition part of a SimplexSplitting output. This is a
# (k^E * (E+1) x k) array of integers. (so, ignore the orientations for now)
#
# @return
#  NewVert An array of the vertices of the splitted simplex.
function simplicial_subdivision_single(k, E)
    ##
    # OUTCOME:
    # NewVert is an array of dimension R x k. Each row contains the indices of
    # the original vertices defining the new vertices of the simplicial subdivision.
    # R = (k+E)!/(k! E!)-E-1
    # Subtriangulation is an array of dimension (E+1)k^E x 1. It consists of a
    # stack of k^E columns of dimension E+1, each one of them containing
    # the indices of the vertices generating the corresponding simplex of the subdivision in the order
    # given by [rows(OriginalVertices);rows(NewVert)]
    ##

    # Get simplex splitting rules, ignore orientations for now.
    splits = simplex_split(k, E, orientations = false)

    # Retrieve integer labels from the tensor decomposition. A way of uniquely identifying
    # the rows of splits. Each integer corresponds to a row in
    # splits (but are not necessarily in order).
    integer_labels = (splits .- 1) * ((E+1) .^ collect(0:(k-1)))

    uniquelabels = unique(integer_labels)
    n_rows = size(uniquelabels, 1)

    J = indexin(integer_labels, uniquelabels)

    # IndicesOfUnqiueIntegersinOriginal
    I = Vector{Int}(length(uniquelabels))
    for i = 1:length(uniquelabels)
        I[i] = findfirst(integer_labels, uniquelabels[i])
    end

    Aux = splits[I, :]
    V = repeat(uniquelabels, 1, E + 1)
    aux = repeat((((E+1)^k - 1) ./ E) * transpose(collect(0:E)), size(uniquelabels, 1), 1)

    aux = heaviside0(-abs.(V - aux)) .* repeat(collect(1:n_rows), 1, E + 1)
    aux = aux[find(aux)]

    aux = round.(Int, aux)
    Caux = round.(Int, complementary(aux, n_rows))
    Ipp = [aux; Caux]

    Ip = zeros(Int, size(I))

    # for any i\in Ipp, Ip(i) is the number j\in {1,...,size(Ipp,1)} such that Ipp(j)=i
    Ip[Ipp] = transpose(collect(1:size(I, 1)))

    # Indices of the non-repeated vertices that are not part of the original simplex.
    # These are the truly new vertices generated by the splitting.
    NewVert = Aux[Caux, :]

    allverts = vcat(repeat(collect(1:(E+1)), 1, k), NewVert)

    # Subarrays of dimension E+1 x 1 vertically concatenated. Each subarray contains the
    # indices in the rows of NewVert corresponding to the vertices furnishing the corresponding
    # simplex (analogous to the output of SimplexSplitting, check notes).
    SubTriangulation = round.(Int, Ip[J])
    SubTriangulation = copy(transpose(reshape(SubTriangulation, E + 1, k^E)))

    return allverts, SubTriangulation
end


text:
"""
functions that will enable creation of long String that will be the code for custom fragment shader that will be suited for defined textures
"""
module CustomFragShad

using ModernGL, GeometryTypes, GLFW,  ..ForDisplayStructs, ColorTypes

export createCustomFramgentShader,divideTexteuresToMainAndRest,addSamplerStr

"""
We divide textures into main image texture and the rest
listOfTexturesToCreate - list of textures on the basis of which we will create custom fragment shader code
returns tuple where fist entr is the main image texture specification and second is rest of textures
  """
function divideTexteuresToMainAndRest(listOfTexturesToCreate::Vector{TextureSpec})::Tuple{TextureSpec, Vector{TextureSpec}}
    mainTexture =   filter(it->it.isMainImage,listOfTexturesToCreate )[1] # main image texture
    notMainTextures = filter(it->!it.isMainImage,listOfTexturesToCreate) # textures associated with not main
return (mainTexture, notMainTextures)
end #divideTexteuresToMainAndRest

"""
We will in couple steps create code for fragment shader that will be based on the texture definitions we gave
listOfTexturesToCreate - list of textures on the basis of which we will create custom fragment shader code
maskToSubtrastFrom,maskWeAreSubtracting - texture specifications used in order to generate code needed to diplay diffrence between those two masks - every time we want diffrent diffrence we will need to recreate shader invoking this function
"""
function createCustomFramgentShader(listOfTexturesToCreate::Vector{TextureSpec}
                                    ,maskToSubtrastFrom::TextureSpec
                                    ,maskWeAreSubtracting::TextureSpec )::String
    mainTexture, notMainTextures=divideTexteuresToMainAndRest(listOfTexturesToCreate)
masksConstsants= map( x-> addMasksStrings(x), notMainTextures) |>
(strings)-> join(strings, "")

continuusColorMasks =   filter(it->it.isContinuusMask,listOfTexturesToCreate ) # main image texture

 res =  """
$(initialStrings())

$(addUniformsForNuclearAndSubtr())

$(addMainTextureStrings(mainTexture))
$masksConstsants
$(getNuclearMaskFunctions(continuusColorMasks))
$( getMasksSubtractionFunction(maskToSubtrastFrom,maskWeAreSubtracting))

$(mainFuncString(mainTexture,notMainTextures,maskToSubtrastFrom,maskWeAreSubtracting))
 """
# uncomment for debugging
#   for st in split(res, "\n")
#     @info st
#     end
return res

end #createCustomFramgentShader

"""
some initial constants that are the same irrespective of textures
"""
function initialStrings()::String
return """
#extension GL_EXT_gpu_shader4 : enable    //Include support for this extension, which defines usampler2D
out vec4 FragColor;
in vec3 ourColor;
smooth in vec2 TexCoord0;
uniform int  min_shown_white = 400;//400 ;// value of cut off  - all values above will be shown as white
uniform int  max_shown_black = -200;//value cut off - all values below will be shown as black
uniform float displayRange = 600.0;
uniform float mainImageContribution = 1.0;//controll main image contribution to color
"""
end #initialStrings

"""
managing main texture and on the basis of the type  we will initialize diffrent variables
mainTexture - specification of a texture related to main image
  """
function addMainTextureStrings(mainTexture::TextureSpec)::String
    mainImageName= mainTexture.name
    return """
    uniform $(addSamplerStr(mainTexture, mainImageName)); // main image sampler
    uniform int $(mainImageName)isVisible = 1; // controllin main texture visibility

    $(addWindowingFunc(mainTexture))

    """
end #addMainTextureStrings

"""
setting string representing sampler depending on type
"""
function addSamplerStr(textur::TextureSpec, samplerName::String)::String

    if supertype(Int)== supertype(parameter_type(textur))
        return "isampler2D $(samplerName)"
    elseif supertype(UInt)== supertype(parameter_type(textur))
        return "usampler2D $(samplerName)"
    else # so float ...
        return "sampler2D $(samplerName)"
    end
end #addSamplerStr

"""
giving variable name associated with given type
"""
function addTypeStr(textur::TextureSpec)::String

    if supertype(Int)== supertype(parameter_type(textur))
        return "int"
    elseif supertype(UInt)== supertype(parameter_type(textur))
        return "uint"
    else # so float ...
        return "float"
    end
end #addTypeStr

"""
giving string representing main function defining how windowing of main image should be performed
"""
function addWindowingFunc(textur::TextureSpec)::String
    typeStr = addTypeStr(textur)
    mainImageName= textur.name

   return """

   //in case of $typeStr texture  controlling color display of main image we keep all above some value as white and all below some value as black
   vec4 mainColor(in $typeStr texel)
   {
       if($(mainImageName)isVisible==0){
           return vec4(0.0, 0.0, 0.0, 1.0);
       }
       else if(texel >min_shown_white){
           return vec4(1.0, 1.0, 1.0, 1.0);
           }
       else if (texel< max_shown_black){
           return vec4(0.0, 0.0, 0.0, 1.0);
      }
       else{
         float fl = float(texel-max_shown_black)/displayRange ;
        return vec4(fl, fl, fl, 1.0);
       }
   }

   """

end #addWindowingFunc

"""
Adding string necessary for managing ..Uniforms of masks textures
"""
function addMasksStrings(textur::TextureSpec)
    textName = textur.name
return """

uniform $(addSamplerStr(textur,textName )); // mask image sampler
$(addColorUniform(textur))
uniform int $(textName)isVisible= 0; // controlling visibility

uniform $(addTypeStr(textur))  $(textName)minValue= $(textur.minAndMaxValue[1]); // minimum possible value set in configuration
uniform $(addTypeStr(textur))  $(textName)maxValue= $(textur.minAndMaxValue[2]); // maximum possible value set in configuration
uniform $(addTypeStr(textur))  $(textName)ValueRange= $(textur.minAndMaxValue[2] -textur.minAndMaxValue[1] ); // range of possible values calculated from above
uniform float  $(textName)maskContribution=1.0; //controls contribution of mask to output color

"""
end#addMasksStrings

"""
Adding ..Uniforms resopnsible for colors associated with given mask
"""
function addColorUniform(textur::TextureSpec)
    #in case of multiple colors used by single mask like in case of nuclearm medicine mask
    if(textur.isContinuusMask)
        # colors = textur.colorSet
        # colLen= length(colors)+1
        # colorStrings = join(map(it ->"vec4($(it.r),$(it.g),$(it.b) ,1.0)" , colors) ,",")
        # return "uniform vec4[$(colLen)] $(textur.name)ColorMask  = vec4[$(colLen)](vec4(0.0,0.0,0.0,1.0),$(colorStrings));// we add one so later function operating on this will make easier"
    end
    #in case mask uses only single color
    return "uniform vec4 $(textur.name)ColorMask= vec4(0.0,0.0,0.0,0.0); //controlling colors"

end#addColorUniform

"""
controlling main function - basically we need to return proper FragColor which represents pixel color in given spot
we generete separately r,g and b values by adding contributions from all textures
"""
function mainFuncString( mainTexture::TextureSpec
                        ,notMainTextures::Vector{TextureSpec}
                        ,maskToSubtrastFrom::TextureSpec
                        ,maskWeAreSubtracting::TextureSpec)::String

    mainImageName= mainTexture.name

    notMainTexturesNotCont = filter(it->!it.isContinuusMask ,notMainTextures)#only single color associated
    notMainTexturesCont = filter(it->it.isContinuusMask ,notMainTextures)#multiple colors associated

    masksInfluences= map( x-> setMaskInfluence(x), notMainTextures) |>
                    (strings)-> join(strings, "")

   lll = length(notMainTextures)+1
   #The step function returns 0.0 if x is smaller than edge and otherwise 1.0.

   sumColors = map(letter->
                    map( x-> "  $(x.name)ColorMask.$(letter) * ( $(x.name)Res  * step($(x.name)minValue,$(x.name)Res) * step($(x.name)Res ,$(x.name)maxValue ))/ $(x.name)ValueRange "
                        , notMainTexturesNotCont) |>
   (strings)-> join(strings, " + ")              ,["r", "g", "b"])

   sumColorR= sumColors[1]
   sumColorG= sumColors[2]
   sumColorB= sumColors[3]

   sumColorRCont= map( x->"r$(x.name)getColorForMultiColor($(x.name)Res)", notMainTexturesCont) |>
                    (strings)-> join(strings, " + ")
   sumColorGCont= map( x->"g$(x.name)getColorForMultiColor($(x.name)Res)", notMainTexturesCont) |>
                    (strings)-> join(strings, " + ")
   sumColorBCont= map( x->"b$(x.name)getColorForMultiColor($(x.name)Res)", notMainTexturesCont) |>
                    (strings)-> join(strings, " + ")

   isVisibleList=  map( x-> "$(x.name)isVisible *$(x.name)maskContribution", notMainTextures) |>
   (strings)-> join(strings, " + ")
    return """
    void main()
    {
$(masksInfluences)

float todiv = $(isVisibleList) + $(mainImageName)isVisible*mainImageContribution+ isMaskDiffrenceVis;
 vec4 $(mainImageName)Res = mainColor(texture2D($(mainImageName), TexCoord0).r);
   FragColor = vec4(($(sumColorR)+$(sumColorRCont)
   +$(mainImageName)Res.r*mainImageContribution+rdiffrenceColor(texture2D($(maskToSubtrastFrom.name), TexCoord0).r ,texture2D($(maskWeAreSubtracting.name), TexCoord0).r )
    )/todiv
   ,($(sumColorG)+$(sumColorGCont)
   +$(mainImageName)Res.g*mainImageContribution+gdiffrenceColor(texture2D($(maskToSubtrastFrom.name), TexCoord0).r ,texture2D($(maskWeAreSubtracting.name), TexCoord0).r ))
   /todiv,
   ($(sumColorB)+$(sumColorBCont)
   +$(mainImageName)Res.b*mainImageContribution+bdiffrenceColor(texture2D($(maskToSubtrastFrom.name), TexCoord0).r ,texture2D($(maskWeAreSubtracting.name), TexCoord0).r ) )
   /todiv,
   1.0  ); //long  product, if mask is invisible it just has full transparency

    }

    """
end#mainFuncString

"""
Giving value from texture f the texture is set to be visible otherwise 0
"""
function setMaskInfluence(textur::TextureSpec)
    textName = textur.name

return """

float $(textName)Res = texture2D($(textName), TexCoord0).r * $(textName)isVisible*$(textName)maskContribution ;

"""
end#setMaskInfluence

"""
on the basis of texture type gives proper function controlling color of the mask
"""
function chooseColorFonuction(textur::TextureSpec)::String

    if supertype(Int)== supertype(parameter_type(textur))
        return "imaskColor"
    elseif supertype(UInt)== supertype(parameter_type(textur))
        return "umaskColor"
    else # so float ...
        return "fmaskColor"
    end
end#chooseColorFonuction

# """
# used to display and debug  output - output can be also additionally checked using this tool http://evanw.github.io/glslx/
# """
# function debuggingOutput(listOfTexturesToCreate)
#     strr=  ..CustomFragShad.createCustomFramgentShader(listOfTexturesToCreate)
#     for st in split(strr, "\n")
#     @info st
#     end
# end#

 """
used in order to enable subtracting one mask from the other - hence displaying
pixels where value of mask a is present but mask b not (order is important)
automatically both masks will be set to be invisible and only the diffrence displayed

In order to provide maximum performance and avoid branching inside shader multiple shader programs will be attached and one choosed  that will use diffrence needed
maskToSubtrastFrom,maskWeAreSubtracting - specifications o textures we are operating on
"""
function getMasksSubtractionFunction(maskToSubtrastFrom::TextureSpec
                                    ,maskWeAreSubtracting::TextureSpec)::String

  #letter - r g or b - will create specialized functions for r g and b colors

return map(letter-> """

float $(letter)diffrenceColor($(addTypeStr(maskToSubtrastFrom)) maskkA,$(addTypeStr(maskWeAreSubtracting)) maskkB)
{
  return   max(float(maskkA)-float(maskkB),0.0)*(($(maskToSubtrastFrom.name)ColorMask.$(letter) + $(maskWeAreSubtracting.name)ColorMask.$(letter))/2)*isMaskDiffrenceVis;
}

""",["r","g","b"]  )|>  x->join(x)

end #getMasksSubtractionFunction
"""
Enable displaying for example nuclear medicine data by applying smoothly changing colors
to floating point data
1)in confguration phase we need to have minimum, maximum and range of possible values associated with nuclear mask
2)as uniform we would need to have set of vec4's - those will be used  to display colors
3)colors will be set with algorithm presented below
    a)range will be divided into sections (value ranges) wich numbers will equal  length of colors vector - 1
    b)in each section  the output color will be mix of 2 colors one associated with this section and with next one
        - contribution  of the color associated with given section will  vary from 100% at the begining of the section to 0%
        in the end where 100% of color will be associated with color of next section

"""
function getNuclearMaskFunctions(continuusColorTextSpecs::Vector{TextureSpec} )::String
   #Check in which range we are without if
   #Important first color in color list needs to be doubled in order to make algorithm cleaner - so we will start from index 1 and always there would be some previous index

   tuples= map(x->[ (x.name,[ [a.r,a.g,a.b] for a in x.colorSet  ],"r" ,1),(x.name,[[a.r,a.g,a.b] for a in x.colorSet  ],"g",2),(x.name,[[a.r,a.g,a.b] for a in x.colorSet  ],"b",3)],continuusColorTextSpecs)

   if(!isempty(tuples)) tuples=  reduce(vcat,tuples) end

   return join(map(x->"""

float $(x[3])$(x[1])getColorForMultiColor(float innertexelRes) {

 float texelRes=  clamp(innertexelRes, $(x[1])minValue,$(x[1])maxValue );
        float normalized = (texelRes/float($(x[1])ValueRange))*$(length(x[2])+1);
        uint indexx = uint(floor(normalized)) ;// so we normalize floor  in order to get index of color from color list
        float[$(length(x[2])+1)] colorFloats = float[$(length(x[2])+1)](0.0,$( map(it->it[x[4]],x[2])|> (fls)-> join(fls,",")    )   )  ;
        float normalizedColorPercent= normalized-float(indexx) ;
        return  clamp(colorFloats[indexx +1]*normalizedColorPercent+colorFloats[indexx]*(1- normalizedColorPercent) ,0.0,1.0 );//so we get color from current section and closer we are to the end of this section the bigger influence of this color, the closer we are to the begining the bigger the influence of previous section color
}

        """ ,  tuples  )," ")

end#getNuclearMaskFunction

# if( innertexelRes < $(x[1])minValue  ){
#     texelRes = $(x[1])minValue;
#       }
#       else if(innertexelRes > $(x[1])maxValue ){
#     texelRes = $(x[1])maxValue ;
#           }
#       else{
#     texelRes = innertexelRes ;
#    }

"""
Add ..Uniforms necessary for operation of mask subtraction and  for proper display of nuclear medicine masks
"""
function addUniformsForNuclearAndSubtr()::String
return """
// for mask difference display
uniform int isMaskDiffrenceVis=0 ;//1 if we want to display mask difference
//uniform int maskAIndex=0 ;//marks index of first mask we want to get diffrence visualized
//uniform int maskBIndex=0 ;//marks index of second mask we want to get diffrence visualized
"""
end#add..UniformsForNuclearAndSubtr

end #..CustomFragShad module


text:
# Definition of an orthonormal basis
"""
    OrthonormalBasis{T} <: Basis{T}

A list of vector like objects of type `T` that are mutually orthogonal and normalized to one,
representing an orthonormal basis for some subspace (typically a Krylov subspace). See also
[`Basis`](@ref)

Orthonormality of the vectors contained in an instance `b` of `OrthonormalBasis`
(i.e. `all(dot(b[i],b[j]) == I[i,j] for i=1:length(b), j=1:length(b))`) is not checked when elements are added; it is up to the algorithm that constructs `b` to
guarantee orthonormality.

One can easily orthogonalize or orthonormalize a given vector `v` with respect to a `b::OrthonormalBasis`
using the functions [`w, = orthogonalize(v,b,...)`](@ref orthogonalize) or [`w, = orthonormalize(v,b,...)`](@ref orthonormalize).
The resulting vector `w` of the latter can then be added to `b` using `push!(b, w)`. Note that
in place versions [`orthogonalize!(v, b, ...)`](@ref orthogonalize) or [`orthonormalize!(v, b, ...)`](@ref orthonormalize)
are also available.

Finally, a linear combination of the vectors in `b::OrthonormalBasis` can be obtained by multiplying
`b` with a `Vector{<:Number}` using `*` or `mul!` (if the output vector is already allocated).
"""
struct OrthonormalBasis{T} <: Basis{T}
    basis::Vector{T}
end
OrthonormalBasis{T}() where {T} = OrthonormalBasis{T}(Vector{T}(undef, 0))

# Iterator methods for OrthonormalBasis
Base.IteratorSize(::Type{<:OrthonormalBasis}) = Base.HasLength()
Base.IteratorEltype(::Type{<:OrthonormalBasis}) = Base.HasEltype()

Base.length(b::OrthonormalBasis) = length(b.basis)
Base.eltype(b::OrthonormalBasis{T}) where {T} = T

Base.iterate(b::OrthonormalBasis) = Base.iterate(b.basis)
Base.iterate(b::OrthonormalBasis, state) = Base.iterate(b.basis, state)

Base.getindex(b::OrthonormalBasis, i) = getindex(b.basis, i)
Base.setindex!(b::OrthonormalBasis, i, q) = setindex!(b.basis, i, q)
Base.firstindex(b::OrthonormalBasis) = firstindex(b.basis)
Base.lastindex(b::OrthonormalBasis) = lastindex(b.basis)

Base.first(b::OrthonormalBasis) = first(b.basis)
Base.last(b::OrthonormalBasis) = last(b.basis)

Base.popfirst!(b::OrthonormalBasis) = popfirst!(b.basis)
Base.pop!(b::OrthonormalBasis) = pop!(b.basis)
Base.push!(b::OrthonormalBasis{T}, q::T) where {T} = (push!(b.basis, q); return b)
Base.empty!(b::OrthonormalBasis) = (empty!(b.basis); return b)
Base.sizehint!(b::OrthonormalBasis, k::Int) = (sizehint!(b.basis, k); return b)
Base.resize!(b::OrthonormalBasis, k::Int) = (resize!(b.basis, k); return b)

# Multiplication methods with OrthonormalBasis
function Base.:*(b::OrthonormalBasis, x::AbstractVector)
    y = zero(eltype(x))*first(b)
    mul!(y, b, x)
end
LinearAlgebra.mul!(y, b::OrthonormalBasis, x::AbstractVector) = unproject!(y, b, x, 1, 0)

const BLOCKSIZE = 4096

function project!(y::AbstractVector, b::OrthonormalBasis, x, α::Number = true, β::Number = false, r = Base.OneTo(length(b)))
    # no specialized routine for IndexLinear x because reduction dimension is large dimension
    length(y) == length(r) || throw(DimensionMismatch())
    Threads.@threads for j = 1:length(r)
        @inbounds begin
            if β == 0
                y[j] = α * dot(b[r[j]], x)
            else
                y[j] = β*y[j] + α * dot(b[r[j]], x)
            end
        end
    end
    return y
end

"""
    unproject!(y, b::OrthonormalBasis, x::AbstractVector, α::Number = 1, β::Number = 0, r = Base.OneTo(length(b)))

For a given orthonormal basis `b`, reconstruct the vector-like object `y` that is defined by
expansion coefficients with respect to the basis vectors in `b` in `x`; more specifically
this computes
```
    y = β*y + α * sum(b[r[i]]*x[i] for i = 1:length(r))
```
"""
function unproject!(y, b::OrthonormalBasis, x::AbstractVector, α::Number = true, β::Number = false, r = Base.OneTo(length(b)))
    if y isa AbstractArray && IndexStyle(y) isa IndexLinear && Threads.nthreads() > 1
        return unproject_linear_multithreaded!(y, b, x, α, β, r)
    end
    # general case: using only vector operations, i.e. axpy! (similar to BLAS level 1)
    length(x) == length(r) || throw(DimensionMismatch())
    if β == 0
        rmul!(y, false) # should be hard zero
    elseif β != 1
        rmul!(y, β)
    end
    @inbounds for (i, ri) = enumerate(r)
        y = axpy!(α*x[i], b[ri], y)
    end
    return y
end
function unproject_linear_multithreaded!(y::AbstractArray, b::OrthonormalBasis{<:AbstractArray}, x::AbstractVector, α::Number = true, β::Number = false, r = Base.OneTo(length(b)))
    # multi-threaded implementation, similar to BLAS level 2 matrix vector multiplication
    m = length(y)
    n = length(r)
    length(x) == n || throw(DimensionMismatch())
    for rj in r
        length(b[rj]) == m || throw(DimensionMismatch())
    end
    if n == 0
        return β == 1 ? y : β == 0 ? fill!(y, 0) : rmul!(y, β)
    end
    let m = m, n = n, y = y, x = x, b = b, blocksize = prevpow(2, div(BLOCKSIZE, n))
        Threads.@threads for I = 1:blocksize:m
            unproject_linear_kernel!(y, b, x, I:min(I+blocksize-1, m), α, β, r)
        end
    end
    return y
end
function unproject_linear_kernel!(y::AbstractArray, b::OrthonormalBasis{<:AbstractArray}, x::AbstractVector, I, α::Number, β::Number, r)
    @inbounds begin
        if β == 0
            @simd for i in I
                y[i] = zero(y[i])
            end
        elseif β != 1
            @simd for i in I
                y[i] *= β
            end
        end
        for (j,rj) in enumerate(r)
            xj = x[j]*α
            Vj = b[rj]
            @simd for i in I
                y[i] += Vj[i]*xj
            end
        end
    end
end

"""
    rank1update!(b::OrthonormalBasis, y, x::AbstractVector, α::Number = 1, β::Number = 1, r = Base.OneTo(length(b)))

Perform a rank 1 update of a basis `b`, i.e. update the basis vectors as
```
    b[r[i]] = β*b[r[i]] + α * y * conj(x[i])
```
It is the user's responsibility to make sure that the result is still an orthonormal basis.
"""
@fastmath function rank1update!(b::OrthonormalBasis, y, x::AbstractVector, α::Number = true, β::Number = true, r = Base.OneTo(length(b)))
    if y isa AbstractArray && IndexStyle(y) isa IndexLinear && Threads.nthreads() > 1
        return rank1update_linear_multithreaded!(b, y, x, α, β, r)
    end
    # general case: using only vector operations, i.e. axpy! (similar to BLAS level 1)
    length(x) == length(r) || throw(DimensionMismatch())
    @inbounds for (i, ri) = enumerate(r)
        if β == 1
            b[ri] = axpy!(α*conj(x[i]), y, b[ri])
        elseif β == 0
            b[ri] = mul!(b[ri], α*x[i], y)
        else
            b[ri] = axpby!(α*x[i], y, β, b[ri])
        end
    end
    return b
end
@fastmath function rank1update_linear_multithreaded!(b::OrthonormalBasis{<:AbstractArray}, y::AbstractArray, x::AbstractVector, α::Number, β::Number, r)
    # multi-threaded implementation, similar to BLAS level 2 matrix vector multiplication
    m = length(y)
    n = length(r)
    length(x) == n || throw(DimensionMismatch())
    for rj in r
        length(b[rj]) == m || throw(DimensionMismatch())
    end
    if n == 0
        return b
    end
    blocksize = prevpow(2, div(BLOCKSIZE, n))
    let m = m, n = n, y = y, x = x, b = b, blocksize = prevpow(2, div(BLOCKSIZE, n))
        Threads.@threads for I = 1:blocksize:m
            @inbounds begin
                for (j,rj) in enumerate(r)
                    xj = α*conj(x[j])
                    Vj = b[rj]
                    if β == 0
                        @simd for i = I:min(I+blocksize-1, m)
                            Vj[i] = zero(Vj[i])
                        end
                    elseif β != 1
                        @simd for i = I:min(I+blocksize-1, m)
                            Vj[i] *= β
                        end
                    end
                    if I + blocksize-1 <= m
                        @simd for i = Base.OneTo(blocksize)
                            Vj[I-1+i] += y[I-1+i]*xj
                        end
                    else
                        @simd for i = I:m
                            Vj[i] += y[i]*xj
                        end
                    end
                end
            end
        end
    end
    return b
end

function basistransform!(b::OrthonormalBasis{T}, U::AbstractMatrix) where {T} # U should be unitary or isometric
    if T<:AbstractArray && IndexStyle(T) isa IndexLinear && Threads.nthreads() > 1
        return basistransform_linear_multithreaded!(b, U)
    end
    m, n = size(U)
    m == length(b) || throw(DimensionMismatch())

    b2 = [similar(b[1]) for j = 1:n]
    Threads.@threads for j = 1:n
        mul!(b2[j], b[1], U[1,j])
        for i = 2:m
            axpy!(U[i,j], b[i], b2[j])
        end
    end
    for j = 1:n
        b[j] = b2[j]
    end
    return b
end

function basistransform_linear_multithreaded!(b::OrthonormalBasis{<:AbstractArray}, U::AbstractMatrix) # U should be unitary or isometric
    m, n = size(U)
    m == length(b) || throw(DimensionMismatch())
    K = length(b[1])

    blocksize = prevpow(2, div(BLOCKSIZE, m))
    let b2 = [similar(b[1]) for j = 1:n], K = K, m = m, n = n
        Threads.@threads for I = 1:blocksize:K
            @inbounds for j = 1:n
                b2j = b2[j]
                @simd for i = I:min(I+blocksize-1, K)
                    b2j[i] = zero(b2j[i])
                end
                for k = 1:m
                    bk = b[k]
                    Ukj = U[k,j]
                    @simd for i = I:min(I+blocksize-1, K)
                        b2j[i] += bk[i] * Ukj
                    end
                end
            end
        end
        for j = 1:n
            b[j] = b2[j]
        end
    end
    return b
end

# function basistransform2!(b::OrthonormalBasis, U::AbstractMatrix) # U should be unitary or isometric
#     m, n = size(U)
#     m == length(b) || throw(DimensionMismatch())
#
#     # apply basis transform via householder reflections
#     for j = 1:size(U,2)
#         h, ν = householder(U, j:m, j)
#         lmul!(h, view(U, :, j+1:n))
#         rmul!(b, h')
#     end
#     return b
# end

# Orthogonalization of a vector against a given OrthonormalBasis
orthogonalize(v, args...) = orthogonalize!(copy(v), args...)

function orthogonalize!(v::T, b::OrthonormalBasis{T}, alg::Orthogonalizer) where {T}
    S = promote_type(eltype(v), eltype(T))
    c = Vector{S}(undef, length(b))
    orthogonalize!(v, b, c, alg)
end

function orthogonalize!(v::T, b::OrthonormalBasis{T}, x::AbstractVector, ::ClassicalGramSchmidt) where {T}
    x = project!(x, b, v)
    v = unproject!(v, b, x, -1, 1)
    return (v, x)
end
function reorthogonalize!(v::T, b::OrthonormalBasis{T}, x::AbstractVector, ::ClassicalGramSchmidt) where {T}
    s = similar(x) ## EXTRA ALLOCATION
    s = project!(s, b, v)
    v = unproject!(v, b, s, -1, 1)
    x .+= s
    return (v, x)
end
function orthogonalize!(v::T, b::OrthonormalBasis{T}, x::AbstractVector, ::ClassicalGramSchmidt2) where {T}
    (v, x) = orthogonalize!(v, b, x, ClassicalGramSchmidt())
    return reorthogonalize!(v, b, x, ClassicalGramSchmidt())
end
function orthogonalize!(v::T, b::OrthonormalBasis{T}, x::AbstractVector, alg::ClassicalGramSchmidtIR) where {T}
    nold = norm(v)
    orthogonalize!(v, b, x, ClassicalGramSchmidt())
    nnew = norm(v)
    while eps(one(nnew)) < nnew < alg.η * nold
        nold = nnew
        (v,x) = reorthogonalize!(v, b, x, ClassicalGramSchmidt())
        nnew = norm(v)
    end
    return (v, x)
end

function orthogonalize!(v::T, b::OrthonormalBasis{T}, x::AbstractVector, ::ModifiedGramSchmidt) where {T}
    for (i, q) = enumerate(b)
        s = dot(q, v)
        v = axpy!(-s, q, v)
        x[i] = s
    end
    return (v, x)
end
function reorthogonalize!(v::T, b::OrthonormalBasis{T}, x::AbstractVector, ::ModifiedGramSchmidt) where {T}
    for (i, q) = enumerate(b)
        s = dot(q, v)
        v = axpy!(-s, q, v)
        x[i] += s
    end
    return (v, x)
end
function orthogonalize!(v::T, b::OrthonormalBasis{T}, x::AbstractVector, ::ModifiedGramSchmidt2) where {T}
    (v, x) = orthogonalize!(v, b, x, ModifiedGramSchmidt())
    return reorthogonalize!(v, b, x, ModifiedGramSchmidt())
end
function orthogonalize!(v::T, b::OrthonormalBasis{T}, x::AbstractVector, alg::ModifiedGramSchmidtIR) where {T}
    nold = norm(v)
    (v,x) = orthogonalize!(v, b, x, ModifiedGramSchmidt())
    nnew = norm(v)
    while eps(one(nnew)) < nnew < alg.η * nold
        nold = nnew
        (v,x) = reorthogonalize!(v, b, x, ModifiedGramSchmidt())
        nnew = norm(v)
    end
    return (v, x)
end

# Orthogonalization of a vector against a given normalized vector
function orthogonalize!(v::T, q::T, alg::Union{ClassicalGramSchmidt,ModifiedGramSchmidt}) where {T}
    s = dot(q,v)
    v = axpy!(-s, q, v)
    return (v, s)
end
function orthogonalize!(v::T, q::T, alg::Union{ClassicalGramSchmidt2,ModifiedGramSchmidt2}) where {T}
    s = dot(q,v)
    v = axpy!(-s, q, v)
    ds = dot(q,v)
    v = axpy!(-ds, q, v)
    return (v, s+ds)
end
function orthogonalize!(v::T, q::T, alg::Union{ClassicalGramSchmidtIR,ModifiedGramSchmidtIR}) where {T}
    nold = norm(v)
    s = dot(q,v)
    v = axpy!(-s, q, v)
    nnew = norm(v)
    while eps(one(nnew)) < nnew < alg.η * nold
        nold = nnew
        ds = dot(q,v)
        v = axpy!(-ds, q, v)
        s += ds
        nnew = norm(v)
    end
    return (v, s)
end

"""
    orthogonalize(v, b::OrthonormalBasis, [x::AbstractVector,] algorithm::Orthogonalizer]) -> w, x
    orthogonalize!(v, b::OrthonormalBasis, [x::AbstractVector,] algorithm::Orthogonalizer]) -> w, x

    orthogonalize(v, q, algorithm::Orthogonalizer]) -> w, s
    orthogonalize!(v, q, algorithm::Orthogonalizer]) -> w, s

Orthogonalize vector `v` against all the vectors in the orthonormal basis `b` using the
orthogonalization algorithm `algorithm`, and return the resulting vector `w` and the overlap
coefficients `x` of `v` with the basis vectors in `b`.

In case of `orthogonalize!`, the vector `v` is mutated in place. In both functions, storage
for the overlap coefficients `x` can be provided as optional argument `x::AbstractVector` with
`length(x) >= length(b)`.

One can also orthogonalize `v` against a given vector `q` (assumed to be normalized), in which
case the orthogonal vector `w` and the inner product `s` between `v` and `q` are returned.

Note that `w` is not normalized, see also [`orthonormalize`](@ref).

For algorithms, see [`ClassicalGramSchmidt`](@ref), [`ModifiedGramSchmidt`](@ref), [`ClassicalGramSchmidt2`](@ref),
[`ModifiedGramSchmidt2`](@ref), [`ClassicalGramSchmidtIR`](@ref) and [`ModifiedGramSchmidtIR`](@ref).
"""
orthogonalize, orthogonalize!

# Orthonormalization: orthogonalization and normalization
orthonormalize(v, args...) = orthonormalize!(copy(v), args...)

function orthonormalize!(v, args...)
    out = orthogonalize!(v, args...) # out[1] === v
    β = norm(v)
    v = rmul!(v, inv(β))
    return tuple(v, β, Base.tail(out)...)
end

"""
    orthonormalize(v, b::OrthonormalBasis, [x::AbstractVector,] algorithm::Orthogonalizer]) -> w, β, x
    orthonormalize!(v, b::OrthonormalBasis, [x::AbstractVector,] algorithm::Orthogonalizer]) -> w, β, x

    orthonormalize(v, q, algorithm::Orthogonalizer]) -> w, β, s
    orthonormalize!(v, q, algorithm::Orthogonalizer]) -> w, β, s

Orthonormalize vector `v` against all the vectors in the orthonormal basis `b` using the
orthogonalization algorithm `algorithm`, and return the resulting vector `w` (of norm 1), its
norm `β` after orthogonalizing and the overlap coefficients `x` of `v` with the basis vectors in
`b`, such that `v = β * w + b * x`.

In case of `orthogonalize!`, the vector `v` is mutated in place. In both functions, storage
for the overlap coefficients `x` can be provided as optional argument `x::AbstractVector` with
`length(x) >= length(b)`.

One can also orthonormalize `v` against a given vector `q` (assumed to be normalized), in which
case the orthonormal vector `w`, its norm `β` before normalizing and the inner product `s` between
`v` and `q` are returned.

See [`orthogonalize`](@ref) if `w` does not need to be normalized.

For algorithms, see [`ClassicalGramSchmidt`](@ref), [`ModifiedGramSchmidt`](@ref), [`ClassicalGramSchmidt2`](@ref),
[`ModifiedGramSchmidt2`](@ref), [`ClassicalGramSchmidtIR`](@ref) and [`ModifiedGramSchmidtIR`](@ref).
"""
orthonormalize, orthonormalize!


text: default_block_size() = Bool(VectorizationBase.has_feature(Val(:x86_64_avx512f))) ? 96 : 64


text:
#
# Check whether we are in the "hard case".
#
# Args:
#  H_eigv: The eigenvalues of H, low to high
#  qg: The inner product of the eigenvalues and the gradient in the same order
#
# Returns:
#  hard_case: Whether it is a candidate for the hard case
#  lambda_index: The index of the first lambda not equal to the smallest
#                eigenvalue, which is only correct if hard_case is true.
function check_hard_case_candidate(H_eigv, qg)
    @assert length(H_eigv) == length(qg)
    if H_eigv[1] >= 0
        # The hard case is only when the smallest eigenvalue is negative.
        return false, 1
    end
    hard_case = true
    lambda_index = 1
    hard_case_check_done = false
    while !hard_case_check_done
        if lambda_index > length(H_eigv)
            hard_case_check_done = true
        elseif abs(H_eigv[1] - H_eigv[lambda_index]) > 1e-10
            # The eigenvalues are reported in order.
            hard_case_check_done = true
        else
            if abs(qg[lambda_index]) > 1e-10
                hard_case_check_done = true
                hard_case = false
            end
            lambda_index += 1
        end
    end

    hard_case, lambda_index
end

# Equation 4.38 in N&W (2006)
function calc_p!(lambda::T, min_i, n, qg, H_eig, p) where T
    fill!( p, zero(T) )
    for i = min_i:n
        p[:] -= qg[i] / (H_eig.values[i] + lambda) * H_eig.vectors[:, i]
    end
    return nothing
end

# Choose a point in the trust region for the next step using
# the interative (nearly exact) method of section 4.3 of N&W (2006).
# This is appropriate for Hessians that you factorize quickly.
#
# Args:
#  gr: The gradient
#  H:  The Hessian
#  delta:  The trust region size, ||s|| <= delta
#  s: Memory allocated for the step size, updated in place
#  tolerance: The convergence tolerance for root finding
#  max_iters: The maximum number of root finding iterations
#
# Returns:
#  m - The numeric value of the quadratic minimization.
#  interior - A boolean indicating whether the solution was interior
#  lambda - The chosen regularizing quantity
#  hard_case - Whether or not it was a "hard case" as described by N&W (2006)
#  reached_solution - Whether or not a solution was reached (as opposed to
#      terminating early due to max_iters)
function solve_tr_subproblem!(gr,
                              H,
                              delta,
                              s;
                              tolerance=1e-10,
                              max_iters=5)
    T = eltype(gr)
    n = length(gr)
    delta_sq = delta^2

    @assert n == length(s)
    @assert (n, n) == size(H)
    @assert max_iters >= 1

    # Note that currently the eigenvalues are only sorted if H is perfectly
    # symmetric.  (Julia issue #17093)
    H_eig = eigen(Symmetric(H))
    min_H_ev, max_H_ev = H_eig.values[1], H_eig.values[n]
    H_ridged = copy(H)

    # Cache the inner products between the eigenvectors and the gradient.
    qg = H_eig.vectors' * gr

    # These values describe the outcome of the subproblem.  They will be
    # set below and returned at the end.
    interior = true
    hard_case = false
    reached_solution = true

    # Unconstrained solution
    if min_H_ev >= 1e-8
        calc_p!(zero(T), 1, n, qg, H_eig, s)
    end

    if min_H_ev >= 1e-8 && sum(abs2, s) <= delta_sq
        # No shrinkage is necessary: -(H \ gr) is the minimizer
        interior = true
        reached_solution = true
        lambda = zero(T)
    else
        interior = false

        # The hard case is when the gradient is orthogonal to all
        # eigenvectors associated with the lowest eigenvalue.
        hard_case_candidate, min_i =
            check_hard_case_candidate(H_eig.values, qg)

        # Solutions smaller than this lower bound on lambda are not allowed:
        # they don't ridge H enough to make H_ridge PSD.
        lambda_lb = -min_H_ev + max(1e-8, 1e-8 * (max_H_ev - min_H_ev))
        lambda = lambda_lb

        hard_case = false
        if hard_case_candidate
            # The "hard case". lambda is taken to be -min_H_ev and we only need
            # to find a multiple of an orthogonal eigenvector that lands the
            # iterate on the boundary.

            # Formula 4.45 in N&W (2006)
            calc_p!(lambda, min_i, n, qg, H_eig, s)
            p_lambda2 = sum(abs2, s)
            if p_lambda2 > delta_sq
                # Then we can simply solve using root finding.
            else
                hard_case = true
                reached_solution = true

                tau = sqrt(delta_sq - p_lambda2)

                # I don't think it matters which eigenvector we pick so take
                # the first.
                calc_p!(lambda, min_i, n, qg, H_eig, s)
                s[:] = -s + tau * H_eig.vectors[:, 1]
            end
        end

        if !hard_case
            # Algorithim 4.3 of N&W (2006), with s insted of p_l for consistency
            # with Optim.jl

            reached_solution = false
            for iter in 1:max_iters
                lambda_previous = lambda

                for i=1:n
                    H_ridged[i, i] = H[i, i] + lambda
                end

                R = cholesky(Hermitian(H_ridged)).U
                s[:] = -R \ (R' \ gr)
                q_l = R' \ s
                norm2_s = dot(s, s)
                lambda_update = norm2_s * (sqrt(norm2_s) - delta) / (delta * dot(q_l, q_l))
                lambda += lambda_update

                # Check that lambda is not less than lambda_lb, and if so, go
                # half the way to lambda_lb.
                if lambda < lambda_lb
                    lambda = 0.5 * (lambda_previous - lambda_lb) + lambda_lb
                end

                if abs(lambda - lambda_previous) < tolerance
                    reached_solution = true
                    break
                end
            end
        end
    end

    m = dot(gr, s) + 0.5 * dot(s, H * s)

    return m, interior, lambda, hard_case, reached_solution
end

struct NewtonTrustRegion{T <: Real} <: SecondOrderOptimizer
    initial_delta::T
    delta_hat::T
    eta::T
    rho_lower::T
    rho_upper::T
end

"""
# NewtonTrustRegion
## Constructor
```julia
NewtonTrustRegion(; initial_delta = 1.0,
                    delta_hat = 100.0,
                    eta = 0.1,
                    rho_lower = 0.25,
                    rho_upper = 0.75)
```

The constructor has 5 keywords:
* `initial_delta`, the starting trust region radius
* `delta_hat`, the largest allowable trust region radius
* `eta`, when `rho` is at least `eta`, accept the step
* `rho_lower`, when `rho` is less than `rho_lower`, shrink the trust region
* `rho_upper`, when `rho` is greater than `rho_upper`, grow the trust region

## Description
The `NewtonTrustRegion` method implements Newton's method with a trust region
for optimizing a function. The method is designed to take advantage of the
second-order information in a function's Hessian, but with more stability that
Newton's method when functions are not globally well-approximated by a quadratic.
This is achieved by repeatedly minimizing quadratic approximations within a
dynamically-sized trust region in which the function is assumed to be locally
quadratic. See Wright and Nocedal and Wright (ch. 4, 2006) for a discussion of
trust-region methods in practice.

## References
 - Nocedal, J., & Wright, S. (2006). Numerical optimization. Springer Science & Business Media.
"""
NewtonTrustRegion(; initial_delta::Real = 1.0,
                    delta_hat::Real = 100.0,
                    eta::Real = 0.1,
                    rho_lower::Real = 0.25,
                    rho_upper::Real = 0.75) =
                    NewtonTrustRegion(initial_delta, delta_hat, eta, rho_lower, rho_upper)

Base.summary(::NewtonTrustRegion) = "Newton's Method (Trust Region)"

mutable struct NewtonTrustRegionState{Tx, T, G} <: AbstractOptimizerState
    x::Tx
    x_previous::Tx
    g_previous::G
    f_x_previous::T
    s::Tx
    hard_case::Bool
    reached_subproblem_solution::Bool
    interior::Bool
    delta::T
    lambda::T
    eta::T
    rho::T
end

function initial_state(method::NewtonTrustRegion, options, d, initial_x)
    T = eltype(initial_x)
    n = length(initial_x)
    # Maintain current gradient in gr
    @assert(method.delta_hat > 0, "delta_hat must be strictly positive")
    @assert(0 < method.initial_delta < method.delta_hat, "delta must be in (0, delta_hat)")
    @assert(0 <= method.eta < method.rho_lower, "eta must be in [0, rho_lower)")
    @assert(method.rho_lower < method.rho_upper, "must have rho_lower < rho_upper")
    @assert(method.rho_lower >= 0.)
    # Keep track of trust region sizes
    delta = copy(method.initial_delta)

    # Record attributes of the subproblem in the trace.
    hard_case = false
    reached_subproblem_solution = true
    interior = true
    lambda = NaN

    value_gradient!!(d, initial_x)
    hessian!!(d, initial_x)

    NewtonTrustRegionState(copy(initial_x), # Maintain current state in state.x
                         copy(initial_x), # Maintain previous state in state.x_previous
                         copy(gradient(d)), # Store previous gradient in state.g_previous
                         T(NaN), # Store previous f in state.f_x_previous
                         similar(initial_x), # Maintain current search direction in state.s
                         hard_case,
                         reached_subproblem_solution,
                         interior,
                         T(delta),
                         T(lambda),
                         T(method.eta), # eta
                         zero(T)) # rho
end

function update_state!(d, state::NewtonTrustRegionState, method::NewtonTrustRegion)
    T = eltype(state.x)
    # Find the next step direction.
    m, state.interior, state.lambda, state.hard_case, state.reached_subproblem_solution =
        solve_tr_subproblem!(gradient(d), NLSolversBase.hessian(d), state.delta, state.s)

    # Maintain a record of previous position
    copyto!(state.x_previous, state.x)
    state.f_x_previous  = value(d)
    copyto!(state.g_previous, gradient(d))

    # Update current position
    state.x .+= state.s

    # Update the function value and gradient
    value_gradient!(d, state.x)

    # Update the trust region size based on the discrepancy between
    # the predicted and actual function values.  (Algorithm 4.1 in N&W (2006))
    f_x_diff = state.f_x_previous - value(d)
    if abs(m) <= eps(T)
        # This should only happen when the step is very small, in which case
        # we should accept the step and assess_convergence().
        state.rho = 1.0
    elseif m > 0
        # This can happen if the trust region radius is too large and the
        # Hessian is not positive definite.  We should shrink the trust
        # region.
        state.rho = -1.0
    else
        state.rho = f_x_diff / (0 - m)
    end

    if state.rho < method.rho_lower
        state.delta *= 0.25
    elseif (state.rho > method.rho_upper) && (!state.interior)
        state.delta = min(2 * state.delta, method.delta_hat)
    else
        # else leave delta unchanged.
    end

    if state.rho <= state.eta
        # The improvement is too small and we won't take it.

        # If you reject an interior solution, make sure that the next
        # delta is smaller than the current step. Otherwise you waste
        # steps reducing delta by constant factors while each solution
        # will be the same.
        x_diff = state.x - state.x_previous
        delta = 0.25 * norm(x_diff)

        d.F = state.f_x_previous
        copyto!(state.x, state.x_previous)
        copyto!(gradient(d), state.g_previous)
    end

    false
end

function assess_convergence(state::NewtonTrustRegionState, d, options::Options)
    x_converged, f_converged, g_converged, converged, f_increased = false, false, false, false, false
    if state.rho > state.eta
        # Accept the point and check convergence
        x_converged,
        f_converged,
        g_converged,
        converged,
        f_increased = assess_convergence(state.x,
                                       state.x_previous,
                                       value(d),
                                       state.f_x_previous,
                                       gradient(d),
                                       options.x_abstol,
                                       options.f_reltol,
                                       options.g_abstol)
    end
    x_converged, f_converged, g_converged, converged, f_increased
end

function trace!(tr, d, state, iteration, method::NewtonTrustRegion, options, curr_time=time())
    dt = Dict()
    dt["time"] = curr_time
    if options.extended_trace
        dt["x"] = copy(state.x)
        dt["g(x)"] = copy(gradient(d))
        dt["h(x)"] = copy(NLSolversBase.hessian(d))
        dt["delta"] = copy(state.delta)
        dt["interior"] = state.interior
        dt["hard case"] = state.hard_case
        dt["reached_subproblem_solution"] = state.reached_subproblem_solution
        dt["lambda"] = state.lambda
    end
    g_norm = norm(gradient(d), Inf)
    update!(tr,
            iteration,
            value(d),
            g_norm,
            dt,
            options.store_trace,
            options.show_trace,
            options.show_every,
            options.callback)
end