text:
export
dense,
dropout,
flatten,
Dense,
ELU,
Flatten,
LeakyReLU,
MaxPooling1D,
MaxPooling2D,
MaxPooling3D,
Reshape,
UpSampling1D,
UpSampling2D,
UpSampling3D,
ZeroPadding1D,
ZeroPadding2D,
ZeroPadding3D,
Conv1D,
Conv2D,
Conv3D,
Conv2DTranspose,
Conv3DTranspose,
BatchNormalization,
Dropout,
ae,
ae_num,
ae_init,
ae_to_code,
fc,
sparse_softmax_cross_entropy_with_logits

# for a keras layer, `training` is a keyword
# dropout = tf.keras.layers.Dropout(0.2, noise_shape=None, seed=None)(dense, training=is_training)

Dense(args...;kwargs...) = tf.keras.layers.Dense(args...;kwargs...)
ELU(args...;kwargs...) = tf.keras.layers.ELU(args...;kwargs...)
Flatten(args...;kwargs...) = tf.keras.layers.Flatten(args...;kwargs...)
LeakyReLU(args...;kwargs...) = tf.keras.layers.LeakyReLU(args...;kwargs...)
MaxPooling1D(args...;kwargs...) = tf.keras.layers.MaxPooling1D(args...;kwargs...)
MaxPooling2D(args...;kwargs...) = tf.keras.layers.MaxPooling2D(args...;kwargs...)
MaxPooling3D(args...;kwargs...) = tf.keras.layers.MaxPooling3D(args...;kwargs...)
Reshape(args...;kwargs...) = tf.keras.layers.Reshape(args...;kwargs...)
UpSampling1D(args...;kwargs...) = tf.keras.layers.UpSampling1D(args...;kwargs...)
UpSampling2D(args...;kwargs...) = tf.keras.layers.UpSampling2D(args...;kwargs...)
UpSampling3D(args...;kwargs...) = tf.keras.layers.UpSampling3D(args...;kwargs...)
ZeroPadding1D(args...;kwargs...) = tf.keras.layers.ZeroPadding1D(args...;kwargs...)
ZeroPadding2D(args...;kwargs...) = tf.keras.layers.ZeroPadding2D(args...;kwargs...)
ZeroPadding3D(args...;kwargs...) = tf.keras.layers.ZeroPadding3D(args...;kwargs...)
Conv1D(args...;kwargs...) = tf.keras.layers.Conv1D(args...;kwargs...)
Conv2D(args...;kwargs...) = tf.keras.layers.Conv2D(args...;kwargs...)
Conv3D(args...;kwargs...) = tf.keras.layers.Conv3D(args...;kwargs...)
Conv2DTranspose(args...;kwargs...) = tf.keras.layers.Conv2DTranspose(args...;kwargs...)
Conv3DTranspose(args...;kwargs...) = tf.keras.layers.Conv3DTranspose(args...;kwargs...)
BatchNormalization(args...;kwargs...) = tf.keras.layers.BatchNormalization(args...;kwargs...)
Dropout(args...;kwargs...) = tf.keras.layers.Dropout(args...;kwargs...)

@doc raw"""
    fc(x::Union{Array{Float64,2},PyObject}, output_dims::Array{Int64,1},
    θ::Union{Array{Float64,1}, PyObject};
    activation::String = "tanh")

Creates a fully connected neural network with output dimension `o` and inputs $x\in \mathbb{R}^{m\times n}$.

$$x \rightarrow o_1 \rightarrow o_2 \rightarrow \ldots \rightarrow o_k$$

`θ` is the weights and biases of the neural network, e.g., `θ = ae_init(output_dims)`.

`fc` outputs two tensors:

- the output of the neural network: $u\in \mathbb{R}^{m\times o_k}$.

- the sensitivity of the neural network per sample: $\frac{\partial u}{\partial x}\in \mathbb{R}^{m \times o_k \times n}$

"""
function fc(x::Union{Array{Float64,2},PyObject}, output_dims::Array{Int64,1},
    θ::Union{Array{Float64,1}, PyObject};
    activation::String = "tanh")
    if !haskey(COLIB, "extended_nn")
        install_adept()
        install("ExtendedNN", force=true)
    end
    extended_nn_ = load_system_op(COLIB["extended_nn"]...; multiple=true)
    config = [size(x,2);output_dims]
    x_,config_,θ_ = convert_to_tensor([x,config,θ], [Float64,Int64,Float64])
    x_ = reshape(x_, (-1,))
    u, du = extended_nn_(x_,config_,θ_,activation)
    reshape(u, (size(x,1), config[end])), reshape(du, (size(x,1), config[end], size(x,2)))
end

"""
    ae(x::PyObject, output_dims::Array{Int64}, scope::String = "default";
        activation::Union{Function,String} = "tanh")

Creates a neural network with intermediate numbers of neurons `output_dims`.
"""
function ae(x::Union{Array{Float64},PyObject}, output_dims::Array{Int64}, scope::String = "default";
        activation::Union{Function,String} = "tanh")
    if isa(x, Array)
        x = constant(x)
    end
    flag = false
    if length(size(x))==1
        x = reshape(x, length(x), 1)
        flag = true
    end
    net = x
    variable_scope(scope, reuse=AUTO_REUSE) do
        for i = 1:length(output_dims)-1
            net = dense(net, output_dims[i], activation=activation)
        end
        net = dense(net, output_dims[end])
    end
    if flag && (size(net,2)==1)
        net = squeeze(net)
    end
    return net
end

"""
    ae(x::Union{Array{Float64}, PyObject}, output_dims::Array{Int64}, θ::Union{Array{Float64}, PyObject};
    activation::Union{Function,String} = "tanh")

Creates a neural network with intermediate numbers of neurons `output_dims`. The weights are given by `θ`

# Example 1: Explicitly construct weights and biases
```julia
x = constant(rand(10,2))
n = ae_num([2,20,20,20,2])
θ = Variable(randn(n)*0.001)
y = ae(x, [20,20,20,2], θ)
```

# Example 2: Implicitly construct weights and biases
```julia
θ = ae_init([10,20,20,20,2])
x = constant(rand(10,10))
y = ae(x, [20,20,20,2], θ)
```

See also [`ae_num`](@ref), [`ae_init`](@ref).
"""
function ae(x::Union{Array{Float64}, PyObject}, output_dims::Array{Int64}, θ::Union{Array{Float64}, PyObject};
    activation::Union{Function,String} = "tanh")
    activation=="tanh" && (activation = tanh)
    activation=="sigmoid" && (activation = sigmoid)
    activation=="selu" && (activation = selu)
    activation=="elu" && (activation = elu)
    activation=="relu" && (activation = relu)
    activation=="leaky_relu" && (activation = leaky_relu)
    x = convert_to_tensor(x)
    θ = convert_to_tensor(θ)
    flag = false
    if length(size(x))==1
        x = reshape(x, length(x), 1)
        flag = true
    end
    offset = 0
    net = x

    output_dims = [size(x,2); output_dims]
    for i = 1:length(output_dims)-2
        m = output_dims[i]
        n = output_dims[i+1]
        net = net * reshape(θ[offset+1:offset+m*n], (m, n))  + θ[offset+m*n+1:offset+m*n+n]
        net = activation(net)
        offset += m*n+n
    end
    m = output_dims[end-1]
    n = output_dims[end]
    net = net * reshape(θ[offset+1:offset+m*n], (m, n))  + θ[offset+m*n+1:offset+m*n+n]

    offset += m*n+n
    if offset!=length(θ)
        error("ADCME: the weights and configuration does not match. Required $offset but given $(length(θ)).")
    end
    if flag && (size(net,2)==1)
        net = squeeze(net)
    end
    return net
end

function ae(x::Union{Array{Float64}, PyObject},
    output_dims::Array{Int64},
    θ::Union{Array{Array{Float64}}, Array{PyObject}};
    activation::Union{Function,String} = "tanh")
    if isa(θ, Array{Array{Float64}})
        val = []
        for t in θ
            push!(val, θ'[:])
        end
        val = vcat(val...)
    else
        val = []
        for t in θ
            push!(val, reshape(θ, (-1,)))
        end
        vcat(val...)
    end
    ae(x, output_dims, θ, activation=activation)
end

@doc raw"""
    ae_init(output_dims::Array{Int64}; T::Type=Float64, method::String="xavier")

Return the initial weights and bias values by TensorFlow as a vector. The neural network architecture is

$$o_1 (Input layer) \rightarrow o_2 \rightarrow \cdots \rightarrow o_n (Output layer)$$

Three types of
random initializers are provided

- `xavier` (default). It is useful for `tanh` fully connected neural network.
```math
W^l_i \sim \sqrt{\frac{1}{n_{l-1}}}
```
- `xavier_avg`. A variant of `xavier`
```math
W^l_i \sim \sqrt{\frac{2}{n_l + n_{l-1}}}
```
- `he`. This is the activation aware initialization of weights and helps mitigate the problem
of vanishing/exploding gradients.
```math
W^l_i \sim \sqrt{\frac{2}{n_{l-1}}}
```
"""
function ae_init(output_dims::Array{Int64}; T::Type=Float64, method::String="xavier")
    N = ae_num(output_dims)
    W = zeros(T, N)
    offset = 0
    for i = 1:length(output_dims)-1
        m = output_dims[i]
        n = output_dims[i+1]
        if method=="xavier"
            W[offset+1:offset+m*n] = randn(T, m*n) * T(sqrt(1/m))
        elseif method=="xavier_normal"
            W[offset+1:offset+m*n] = randn(T, m*n) * T(sqrt(2/(n+m)))
        elseif method=="xavier_uniform"
            W[offset+1:offset+m*n] = rand(T, m*n) * T(sqrt(6/(n+m)))
        elseif method=="he"
            W[offset+1:offset+m*n] = randn(T, m*n) * T(sqrt(2/(m)))
        else
            error("Method $method not understood")
        end
        offset += m*n+n
    end
    W
end

"""
    ae_num(output_dims::Array{Int64})

Estimates the number of weights and biases for the neural network. Note the first dimension
should be the feature dimension (this is different from [`ae`](@ref) since in `ae` the feature
dimension can be inferred), and the last dimension should be the output dimension.
"""
function ae_num(output_dims::Array{Int64})
    offset = 0
    for i = 1:length(output_dims)-2
        m = output_dims[i]
        n = output_dims[i+1]
        offset += m*n+n
    end
    m = output_dims[end-1]
    n = output_dims[end]
    offset += m*n+n
    return offset
end

function _ae_to_code(d::Dict, scope::String; activation::String)
    i = 0
    nn_code = ""
    while true
        si = i==0 ? "" : "_$i"
        Wkey = "$(scope)backslashfully_connected$(si)backslashweightscolon0"
        bkey = "$(scope)backslashfully_connected$(si)backslashbiasescolon0"
        if haskey(d, Wkey)
            if i!=0
                nn_code *= "    isa(net, Array) ? (net = $activation.(net)) : (net=$activation(net))\n"
                nn_code *= "    #-------------------------------------------------------------------\n"
            end
            nn_code *= "    W$i = aedict$scope[\"$Wkey\"]\n    b$i = aedict$scope[\"$bkey\"];\n"
            nn_code *= "    isa(net, Array) ? (net = net * W$i .+ b$i') : (net = net *W$i + b$i)\n"
            i += 1
        else
            break
        end
    end
    nn_code = """  global nn$scope\n  function nn$scope(net)
$(nn_code)    return net\n  end """
    nn_code
end

"""
    ae_to_code(file::String, scope::String; activation::String = "tanh")

Return the code string from the feed-forward neural network data in `file`. Usually we can immediately evaluate
the code string into Julia session by
```julia
eval(Meta.parse(s))
```
If `activation` is not specified, `tanh` is the default.
"""
function ae_to_code(file::String, scope::String; activation::String = "tanh")
    d = matread(file)
    s = "let aedict$scope = matread(\"$file\")\n"*_ae_to_code(d, scope; activation = activation)*"\nend\n"
    return s
end

# tensorflow layers from contrib
for (op, tfop) in [(:avg_pool2d, :avg_pool2d), (:avg_pool3d, :avg_pool3d),
        (:dropout, :dropout), (:flatten, :flatten), (:max_pool2d, :max_pool2d), (:max_pool3d, :max_pool3d)]
    @eval begin
        export $op
        $op(args...; kwargs...) = tf.contrib.layers.$tfop(args...; kwargs...)
    end
end

for (op, tfop) in [(:conv1d, :conv1d), (:conv2d, :conv2d), (:conv2d_in_plane, :conv2d_in_plane),
    (:conv2d_transpose, :conv2d_transpose), (:conv3d, :conv3d), (:conv3d_transpose, :conv3d_transpose)]
    @eval begin
        export $op
        function $op(args...;activation = nothing, bias=false,  kwargs...)
            if isa(activation, String)
                if lowercase(activation)=="relu"
                    activation = relu
                elseif lowercase(activation)=="sigmoid"
                    activation = sigmoid
                elseif lowercase(activation) in ["leakyrelu" , "leaky_relu"]
                    activation = leaky_relu
                else
                    error("Activation function $activation not understood")
                end
            end
            if bias
                biases_initializer = tf.zeros_initializer()
            else
                biases_initializer = nothing
            end
            tf.contrib.layers.$tfop(args...; activation_fn = activation, biases_initializer=biases_initializer, kwargs...)
        end
    end
end

export dense, bn
function dense(inputs, units, args...; activation = nothing, kwargs...)
    string2fn = Dict(
        "relu" => relu,
        "tanh" => tanh,
        "sigmoid" => sigmoid,
        "leakyrelu" => leaky_relu,
        "leaky_relu" => leaky_relu,
        "relu6" => relu6,
        "softmax" => softmax,
        "softplus" => softplus,
        "selu" => selu,
        "elu" => elu
    )
    if isa(activation, String)
        if haskey(string2fn, lowercase(activation))
            activation = string2fn[lowercase(activation)]
        else
            error("Activation function $activation not understood")
        end
    end
    tf.contrib.layers.fully_connected(inputs, units,  args...; activation_fn=activation, kwargs...)
end

"""
    bn(args...;center = true, scale=true, kwargs...)

`bn` accepts a keyword parameter `is_training`.
# Example
```julia
bn(inputs, name="batch_norm", is_training=true)
```

!!! note
    `bn` should be used with `control_dependency`
    ```julia
    update_ops = get_collection(UPDATE_OPS)
    control_dependencies(update_ops) do
        global train_step = AdamOptimizer().minimize(loss)
    end
    ```
"""
function bn(args...;center = true, scale=true, kwargs...)
    @warn  """
`bn` should be used with `control_dependency`
```julia
update_ops = get_collection(UPDATE_OPS)
control_dependencies(update_ops) do
    global train_step = AdamOptimizer().minimize(loss)
end
```
""" maxlog=1
    kwargs = Dict{Any, Any}(kwargs)
    if :is_training in keys(kwargs)
        kwargs[:training] = kwargs[:is_training]
        delete!(kwargs, :is_training)
    end
    if :scope in keys(kwargs)
        kwargs[:name] = kwargs[:scope]
        delete!(kwargs, :scope)
    end
    tf.layers.batch_normalization(args...;center = center, scale=scale, kwargs...)
end
sparse_softmax_cross_entropy_with_logits(args...;kwargs...) = tf.nn.sparse_softmax_cross_entropy_with_logits(args...;kwargs...)

export group_conv2d
function group_conv2d(inputs::PyObject, filters::Int64, args...;  groups = 1, scope=scope, kwargs...)
    if groups==1
        return conv2d(inputs, filters, args...;scope=scope, kwargs...)
    else
        dims = size(inputs)
        if mod(dims[end], groups)!=0 || mod(filters, groups)!=0
            error("channels and outputs must be the multiples of `groups`")
        end
        n = div(filters, groups)
        in_ = Array{PyObject}(undef, groups)
        out_ = Array{PyObject}(undef, groups)
        for i = 1:groups
py"""
temp = $inputs[:,:,:,$((i-1)*n):$(i*n)]
"""
            in_[i] = py"temp"
            out_[i] = conv2d(in_[i], n, args...;scope=scope*"_group$i", kwargs...)
        end
        out = concat(out_, dims=4)
        return out
    end
end

export separable_conv2d
function separable_conv2d(inputs, num_outputs, args...; activation=false, bias=false, kwargs...)
    if isa(activation, String)
        if lowercase(activation)=="relu"
            activation = relu
        elseif lowercase(activation)=="tanh"
            activation = tanh
        elseif lowercase(activation)=="sigmoid"
            activation = sigmoid
        elseif lowercase(activation) in ["leakyrelu" , "leaky_relu"]
            activation = leaky_relu
        else
            error("Activation function $activation not understood")
        end
    end
    if bias
        biases_initializer = tf.zeros_initializer()
    else
        biases_initializer = nothing
    end
    tf.contrib.layers.separable_conv2d(inputs, num_outputs, args...;activation_fn = activation,
                        biases_initializer=biases_initializer, kwargs...)
end

export depthwise_conv2d
function depthwise_conv2d(input, num_outputs, args...;
    kernel_size = 3, channel_multiplier = 1, stride = 1, padding="SAME", bias=false, scope="default", reuse=AUTO_REUSE, kwargs...)
    local res, strides
    if isa(kernel_size, Int64)
        kernel_size = (3,3)
    end
    if isa(stride, Int64)
        strides = (1, stride, stride, 1)
    end
    variable_scope(scope, reuse=reuse) do
        filter = get_variable("dconv2d",
                shape=[kernel_size[1], kernel_size[2], size(input,4), channel_multiplier],
                initializer=tf.contrib.layers.xavier_initializer())
        res = tf.nn.depthwise_conv2d(input, filter, strides, padding, args...; kwargs...)
        if bias
            res = tf.contrib.layers.bias_add(res, scope=scope, reuse=AUTO_REUSE)
        end
    end
    return res
end


text:
function create_profile_sets(number_of_hours, data, df0, df1,ic_mva,owpp_mva)
    pu=data["baseMVA"]
    e2me=1000000/pu#into ME/PU
    #e2me=1
    da=0.835;id=0.165
    #e2me=1#into ME/PU
    extradata = Dict{String,Any}()
    extradata["dim"] = Dict{String,Any}()
    extradata["dim"] = number_of_hours
    extradata["gen"] = Dict{String,Any}()
    for (g, gen) in data["gen"]
        extradata["gen"][g] = Dict{String,Any}()
        extradata["gen"][g]["pmax"] = Array{Float64,2}(undef, 1, number_of_hours)
        extradata["gen"][g]["pmin"] = Array{Float64,2}(undef, 1, number_of_hours)
        extradata["gen"][g]["cost"] = [Vector{Float64}() for i=1:number_of_hours]
    end

    for d in 1:number_of_hours
        #Day ahead BE
        #source generator
            extradata["gen"]["1"]["pmax"][1, d] = ic_mva/pu+owpp_mva/pu*df0.wind[d]
            extradata["gen"]["1"]["pmin"][1, d] = 0
            extradata["gen"]["1"]["cost"][d] = [da*df0.daprice[d]/e2me+id*df0.idprice[d]/e2me,0]
        #load generator
            extradata["gen"]["2"]["pmax"][1, d] = 0
            extradata["gen"]["2"]["pmin"][1, d] = -ic_mva/pu-owpp_mva/pu*df0.wind[d]
            extradata["gen"]["2"]["cost"][d] = [da*df0.daprice[d]/e2me+id*df0.idprice[d]/e2me,0]

            #Day ahead UK
            #source generator
                extradata["gen"]["3"]["pmax"][1, d] = ic_mva/pu+owpp_mva/pu*df0.wind[d]
                extradata["gen"]["3"]["pmin"][1, d] = 0
                extradata["gen"]["3"]["cost"][d] = [da*df1.daprice[d]/e2me+id*df1.idprice[d]/e2me,0]
            #load generator
                extradata["gen"]["4"]["pmax"][1, d] = 0
                extradata["gen"]["4"]["pmin"][1, d] = -ic_mva/pu-owpp_mva/pu*df0.wind[d]
                extradata["gen"]["4"]["cost"][d] = [da*df1.daprice[d]/e2me+id*df1.idprice[d]/e2me,0]
            #Wind generator
                extradata["gen"]["5"]["pmax"][1, d] = owpp_mva/pu*df0.wind[d]
                extradata["gen"]["5"]["pmin"][1, d] = 0
                extradata["gen"]["5"]["cost"][d] = [0.0,0.0]
    end
    return extradata
end

function create_profile_sets_wstrg(number_of_hours, data, df0, df1,ic_mva,owpp_mva)
    pu=data["baseMVA"]
    e2me=1000000/pu#into ME/PU
    #e2me=1
    da=0.835;id=0.165
    #e2me=1#into ME/PU
    extradata = Dict{String,Any}()
    extradata["dim"] = Dict{String,Any}()
    extradata["dim"] = number_of_hours
    extradata["gen"] = Dict{String,Any}()
    for (g, gen) in data["gen"]
        extradata["gen"][g] = Dict{String,Any}()
        extradata["gen"][g]["pmax"] = Array{Float64,2}(undef, 1, number_of_hours)
        extradata["gen"][g]["pmin"] = Array{Float64,2}(undef, 1, number_of_hours)
        extradata["gen"][g]["cost"] = [Vector{Float64}() for i=1:number_of_hours]
    end

    storage=[(i,b) for (i,b) in data["ne_storage"]]
    sort!(storage, by=x->x[2]["energy_rating"])
    extradata["ne_storage"] = Dict{String,Any}()
    for (b, bat) in data["ne_storage"]
        extradata["ne_storage"][b] = Dict{String,Any}()
        extradata["ne_storage"][b]["cost_abs"] = Array{Float64,2}(undef, 1, number_of_hours)
        extradata["ne_storage"][b]["cost_inj"] = Array{Float64,2}(undef, 1, number_of_hours)
        extradata["ne_storage"][b]["charge_rating"] = Array{Float64,2}(undef, 1, number_of_hours)
        extradata["ne_storage"][b]["discharge_rating"] = Array{Float64,2}(undef, 1, number_of_hours)
    end

    for d in 1:number_of_hours
        #Day ahead BE
        #source generator
            extradata["gen"]["1"]["pmax"][1, d] = ic_mva/pu+owpp_mva/pu*df0.wind[d]
            extradata["gen"]["1"]["pmin"][1, d] = 0
            extradata["gen"]["1"]["cost"][d] = [da*df0.daprice[d]/e2me+id*df0.idprice[d]/e2me,0]
        #load generator
            extradata["gen"]["2"]["pmax"][1, d] = 0
            extradata["gen"]["2"]["pmin"][1, d] = -ic_mva/pu-owpp_mva/pu*df0.wind[d]
            extradata["gen"]["2"]["cost"][d] = [da*df0.daprice[d]/e2me+id*df0.idprice[d]/e2me,0]

            #Day ahead UK
            #source generator
                extradata["gen"]["3"]["pmax"][1, d] = ic_mva/pu+owpp_mva/pu*df0.wind[d]
                extradata["gen"]["3"]["pmin"][1, d] = 0
                extradata["gen"]["3"]["cost"][d] = [da*df1.daprice[d]/e2me+id*df1.idprice[d]/e2me,0]
            #load generator
                extradata["gen"]["4"]["pmax"][1, d] = 0
                extradata["gen"]["4"]["pmin"][1, d] = -ic_mva/pu-owpp_mva/pu*df0.wind[d]
                extradata["gen"]["4"]["cost"][d] = [da*df1.daprice[d]/e2me+id*df1.idprice[d]/e2me,0]
            #Wind generator
                extradata["gen"]["5"]["pmax"][1, d] = owpp_mva/pu*df0.wind[d]
                extradata["gen"]["5"]["pmin"][1, d] = 0
                extradata["gen"]["5"]["cost"][d] = [0.0,0.0]

                #up_reg
                up_prices=[df0.regup_price[d],df1.regup_price[d]];
                best_up_price=deepcopy(findmax(up_prices));up_prices[best_up_price[2]]=Inf
                worst_up_price=findmin(up_prices)
                best_up_mwh=[df0.regup_mwh[d],df1.regup_mwh[d]][best_up_price[2]]
                worst_up_mwh=[df0.regup_mwh[d],df1.regup_mwh[d]][worst_up_price[2]]
                #down_reg
                dwn_prices=[df0.regdwn_price[d],df1.regdwn_price[d]];
                best_dwn_price=deepcopy(findmax(dwn_prices));dwn_prices[best_dwn_price[2]]=Inf
                worst_dwn_price=findmin(dwn_prices)
                best_dwn_mwh=[df0.regdwn_mwh[d],df1.regdwn_mwh[d]][best_dwn_price[2]]
                worst_dwn_mwh=[df0.regdwn_mwh[d],df1.regdwn_mwh[d]][worst_dwn_price[2]]

                best_up_mwh = best_up_mwh>0 ? best_up_mwh/pu : 0
                best_dwn_mwh = best_dwn_mwh>0 ? best_dwn_mwh/pu : 0
                worst_up_mwh = worst_up_mwh>0 ? worst_up_mwh/pu : 0
                worst_dwn_mwh = worst_dwn_mwh>0 ? worst_dwn_mwh/pu : 0

                best_up_price = best_up_price[1]>0 ? best_up_price[1]/e2me : 0
                best_dwn_price = best_dwn_price[1]>0 ? best_dwn_price[1]/e2me : 0
                worst_up_price = worst_up_price[1]>0 ? worst_up_price[1]/e2me : 0
                worst_dwn_price = worst_dwn_price[1]>0 ? worst_dwn_price[1]/e2me : 0

                for (b,bat) in storage
                    dcr_temp=deepcopy(data["ne_storage"][b]["discharge_rating"])
                    bup_mwh_temp=deepcopy(best_up_mwh)
                    wdwn_mwh_temp=deepcopy(worst_up_mwh)

                    best_up_mwh = dcr_temp>best_up_mwh ?  best_up_mwh : dcr_temp#can it be from the battery?
                    best_up_mwh = best_up_mwh>0 ?  best_up_mwh : 0
                    worst_up_mwh = dcr_temp>(best_up_mwh+worst_up_mwh) ?  worst_up_mwh : dcr_temp-best_up_mwh#can additional be from the battery?
                    worst_up_mwh = worst_up_mwh>0 ?  worst_up_mwh : 0
                    if ((best_up_mwh+worst_up_mwh)>0.0001)
                        extradata["ne_storage"][b]["discharge_rating"][1, d]=best_up_mwh+worst_up_mwh
                        extradata["ne_storage"][b]["cost_inj"][1, d] = (best_up_mwh/(best_up_mwh+worst_up_mwh))*best_up_price+((worst_up_mwh)/(best_up_mwh+worst_up_mwh))*worst_up_price
                        best_up_mwh=bup_mwh_temp-best_up_mwh
                        worst_up_mwh=wdwn_mwh_temp-worst_up_mwh
                    else
                        extradata["ne_storage"][b]["discharge_rating"][1, d]=0
                        extradata["ne_storage"][b]["cost_inj"][1, d]=0
                    end
                    if (isnan(extradata["ne_storage"][b]["discharge_rating"][1, d]))
                        println("Bad data - detected at "*string(d)*" - "*string(b))
                        extradata["ne_storage"][b]["discharge_rating"][1, d]=0
                        extradata["ne_storage"][b]["cost_inj"][1, d]=0
                    end
                    cr_temp=deepcopy(data["ne_storage"][b]["charge_rating"])
                    bdwn_mwh_temp=deepcopy(best_dwn_mwh)
                    wdwn_mwh_temp=deepcopy(worst_dwn_mwh)

                    best_dwn_mwh = cr_temp>best_dwn_mwh ?  best_dwn_mwh : cr_temp#can it be from the battery?
                    best_dwn_mwh = best_dwn_mwh>0 ?  best_dwn_mwh : 0
                    worst_dwn_mwh = cr_temp>(best_dwn_mwh+worst_dwn_mwh) ?  worst_dwn_mwh : cr_temp-best_dwn_mwh#can additional be from the battery?
                    worst_dwn_mwh = worst_dwn_mwh>0 ?  worst_dwn_mwh : 0
                    if ((best_dwn_mwh+worst_dwn_mwh)>0.0001)
                        extradata["ne_storage"][b]["charge_rating"][1, d]=cr_temp
                        extradata["ne_storage"][b]["cost_abs"][1, d] = (best_dwn_mwh/cr_temp)*best_dwn_price+(worst_dwn_mwh/cr_temp)*worst_dwn_price
                        best_dwn_mwh=bdwn_mwh_temp-best_dwn_mwh
                        worst_dwn_mwh=wdwn_mwh_temp-worst_dwn_mwh
                    else
                        extradata["ne_storage"][b]["charge_rating"][1, d]=cr_temp
                        extradata["ne_storage"][b]["cost_abs"][1, d]=0
                    end
                end
    end
    return extradata
end

function scale_bat_data_cordoba!(data, scenario)
    rescale_hourly = x -> (scenario["hours"] / (8760*scenario["planning_horizon"])) * x # yearly limit on energy absoption
    for (s, strg) in get(data, "ne_storage", Dict{String,Any}())
        _PM._apply_func!(strg, "max_energy_absorption", rescale_hourly)
    end
end

function scale_cost_data_cordoba!(data, scenario)
    rescale_hourly = x -> (8760*scenario["planning_horizon"] / scenario["hours"]) * x # scale hourly costs to the planning horizon
    rescale_total  = x -> (                                1 / scenario["hours"]) * x # scale total costs to the planning horizon
    for (g, gen) in data["gen"]
        _PM._apply_func!(gen, "cost", rescale_hourly)
    end
    for (b, branch) in get(data, "ne_branch", Dict{String,Any}())
        _PM._apply_func!(branch, "construction_cost", rescale_total)
        _PM._apply_func!(branch, "co2_cost", rescale_total)
    end
    for (b, branch) in get(data, "branchdc_ne", Dict{String,Any}())
        _PM._apply_func!(branch, "cost", rescale_total)
        _PM._apply_func!(branch, "co2_cost", rescale_total)
    end
    for (c, conv) in get(data, "convdc_ne", Dict{String,Any}())
        _PM._apply_func!(conv, "cost", rescale_total)
        _PM._apply_func!(conv, "co2_cost", rescale_total)
    end
    for (s, strg) in get(data, "ne_storage", Dict{String,Any}())
        _PM._apply_func!(strg, "eq_cost", rescale_total)
        _PM._apply_func!(strg, "inst_cost", rescale_total)
        _PM._apply_func!(strg, "co2_cost", rescale_total)
        _PM._apply_func!(strg, "cost_abs", rescale_hourly)
        _PM._apply_func!(strg, "cost_inj", rescale_hourly)
    end
end

function get_profile_data_sets(d1,d2,data, n, scenario = Dict{String, Any}())
    data["scenario"] = Dict{String, Any}()
    data["scenario_prob"] = Dict{String, Any}()
    windgenprofile_beuk=[];gencostid_beuk=[];gencost_beuk=[];ic_nflow_losses=[];

    df0=DataFrame("time_stamp"=>[],"daprice"=>[],"idprice"=>[],"wind"=>[],"regup_price"=>[],"regup_mwh"=>[],"regdwn_price"=>[],"regdwn_mwh"=>[])
    df1=DataFrame("time_stamp"=>[],"daprice"=>[],"idprice"=>[],"wind"=>[],"regup_price"=>[],"regup_mwh"=>[],"regdwn_price"=>[],"regdwn_mwh"=>[])#MWh_up,EUR_up,MWh_dwn,EUR_dwn
    z0_data=CSV.read("./test/data/input/"*d1*".csv", DataFrames.DataFrame)
    z1_data=CSV.read("./test/data/input/"*d2*".csv", DataFrames.DataFrame)
    z01_data=innerjoin(z0_data,z1_data, makeunique=true,on=:time_stamp)
    for (s, scnr) in scenario["sc_years"]
        start_idx = (parse(Int, s) - 1) * scenario["hours"]

        #tss=[DataFrame("time_stamp"=>ukbe_data.time_stamp, "price"=>0.835.*ukbe_data.be_eumwh.+0.165.*ukbe_data.be_costid),DataFrame("time_stamp"=>ukbe_data.time_stamp, "price"=>0.835.*ukbe_data.uk_eumwh.+0.165.*ukbe_data.uk_costid)]
        tss=[DataFrame("time_stamp"=>z01_data.time_stamp, "price"=>abs.((0.835.*z01_data.EUR_da.+0.165.*z01_data.EUR_id).-(0.835.*z01_data.EUR_da_1.+0.165.*z01_data.EUR_id_1))),DataFrame("time_stamp"=>z01_data.time_stamp, "price"=>z01_data.Wnd_MWh)]
        tss_bins=cluster_ts(tss,n)
        sc=sample_cluster(tss_bins,tss,n)
        sort!(sc)
        for t in sc
            #Zone 0
            push!(df0,[t,z0_data[!,:EUR_da][findfirst(isequal(t),z0_data[!,:time_stamp])],
            z0_data[!,:EUR_id][findfirst(isequal(t),z0_data[!,:time_stamp])],
            z0_data[!,:Wnd_MWh][findfirst(isequal(t),z0_data[!,:time_stamp])],
            z0_data[!,:MWh_up][findfirst(isequal(t),z0_data[!,:time_stamp])],
            z0_data[!,:EUR_up][findfirst(isequal(t),z0_data[!,:time_stamp])],
            z0_data[!,:MWh_dwn][findfirst(isequal(t),z0_data[!,:time_stamp])],
            z0_data[!,:EUR_dwn][findfirst(isequal(t),z0_data[!,:time_stamp])]])
            #Zone 1
            push!(df1,[t,z1_data[!,:EUR_da][findfirst(isequal(t),z1_data[!,:time_stamp])],
            z1_data[!,:EUR_id][findfirst(isequal(t),z1_data[!,:time_stamp])],
            z1_data[!,:Wnd_MWh][findfirst(isequal(t),z1_data[!,:time_stamp])],
            z1_data[!,:MWh_up][findfirst(isequal(t),z1_data[!,:time_stamp])],
            z1_data[!,:EUR_up][findfirst(isequal(t),z1_data[!,:time_stamp])],
            z1_data[!,:MWh_dwn][findfirst(isequal(t),z1_data[!,:time_stamp])],
            z1_data[!,:EUR_dwn][findfirst(isequal(t),z1_data[!,:time_stamp])]])
        end

        scenario["hours"]=length(df0[!,:time_stamp])
        data["scenario"][s] = Dict()
        data["scenario_prob"][s] = scnr["probability"]
        for h in 1 : scenario["hours"]
            network = start_idx + h
            data["scenario"][s]["$h"] = network
        end

    end
    # Return info
    return data, df0,df1
end

function get_n_profile_data(data, n, scenario = Dict{String, Any}())
    data["scenario"] = Dict{String, Any}()
    data["scenario_prob"] = Dict{String, Any}()
    windgenprofile_beuk=[];gencostid_beuk=[];gencost_beuk=[];ic_nflow_losses=[];
    if haskey(scenario, "mc")
        monte_carlo = scenario["mc"]
    else
        monte_carlo = false
    end
    df=DataFrame("time_stamp"=>[],"z0_daprice"=>[],"z0_idprice"=>[],"z1_daprice"=>[],"z1_idprice"=>[],"z0_wind"=>[],"z1_wind"=>[])
    ukbe_data=CSV.read("./test/data/cordoba/input/ukbe_ts_2.csv", DataFrames.DataFrame)
    for (s, scnr) in scenario["sc_years"]
        start_idx = (parse(Int, s) - 1) * scenario["hours"]

        #tss=[DataFrame("time_stamp"=>ukbe_data.time_stamp, "price"=>0.835.*ukbe_data.be_eumwh.+0.165.*ukbe_data.be_costid),DataFrame("time_stamp"=>ukbe_data.time_stamp, "price"=>0.835.*ukbe_data.uk_eumwh.+0.165.*ukbe_data.uk_costid)]
        tss=[DataFrame("time_stamp"=>ukbe_data.time_stamp, "price"=>abs.((0.835.*ukbe_data.be_eumwh.+0.165.*ukbe_data.be_costid).-(0.835.*ukbe_data.uk_eumwh.+0.165.*ukbe_data.uk_costid))),DataFrame("time_stamp"=>ukbe_data.time_stamp, "price"=>ukbe_data.be_wind)]
        tss_bins=cluster_ts(tss,n)
        sc=sample_cluster(tss_bins,tss,n)
        sort!(sc)
        for t in sc
            push!(df,[t,ukbe_data[!,:be_eumwh][findfirst(isequal(t),ukbe_data[!,:time_stamp])], ukbe_data[!,:be_costid][findfirst(isequal(t),ukbe_data[!,:time_stamp])],ukbe_data[!,:uk_eumwh][findfirst(isequal(t),ukbe_data[!,:time_stamp])], ukbe_data[!,:uk_costid][findfirst(isequal(t),ukbe_data[!,:time_stamp])],ukbe_data[!,:be_wind][findfirst(isequal(t),ukbe_data[!,:time_stamp])], ukbe_data[!,:uk_wind][findfirst(isequal(t),ukbe_data[!,:time_stamp])]])
        end

        scenario["hours"]=length(df[!,:time_stamp])
        data["scenario"][s] = Dict()
        data["scenario_prob"][s] = scnr["probability"]
        for h in 1 : scenario["hours"]
            network = start_idx + h
            data["scenario"][s]["$h"] = network
        end

    end
    # Return info
    return data, df
end

ext: jl


text:
using MCMCChains, Test

ProjDir = joinpath(dirname(@__FILE__), "..", "examples", "Bernoulli")
cd(ProjDir) do

  @testset "Bernoulli" begin

    isdir("tmp") &&
      rm("tmp", recursive=true);

    include(joinpath(ProjDir, "bernoulli.jl"))

    if rc == 0
      s = summarize(chns)
      @test s[:theta, :mean][1] ≈ 0.34 atol=0.1
    end

    isdir("tmp") &&
      rm("tmp", recursive=true);

  end # testset

end # cd


text:
export prepareSyntheticVisDataFiles

function prepareSyntheticVisDataFiles(m,Minv::RegularMesh,MScreen::RegularMesh,pad::Int64,theta_phi_rad::Array{Float64,2},
					XlocScreen::Float64,LocCamera::Array{Float64},filename::String,nShots::Int64,
					noiseCamLoc::Float64, noiseScreenLoc::Float64, noiseAngles::Float64,noiseTrans::Float64,method::String)
n = Minv.n;
m = reshape(m,tuple(n...));
pFor = getSfSParam(Minv,theta_phi_rad + noiseAngles*randn(size(theta_phi_rad)),noiseTrans*randn(size(theta_phi_rad,1),3) ,MScreen,XlocScreen + noiseScreenLoc*randn(),
						LocCamera+noiseCamLoc*randn(size(LocCamera)),method,1);
Dremote,pFor = getData(m[:],pFor);
Data = arrangeRemoteCallDataIntoLocalData(Dremote);

# Data = Data + noiseSample*mean(Dlocal[:])*randn(size(Dlocal));

file = matopen(string(filename,"_data.mat"), "w");
write(file,"domain",Minv.domain);
write(file,"n",Minv.n);
write(file,"noiseScreenLoc",noiseScreenLoc);
write(file,"noiseAngles",noiseAngles);
write(file,"noiseCamLoc",noiseCamLoc);
write(file,"ScreenDomain",MScreen.domain);
write(file,"n_screen",MScreen.n);
write(file,"CameraLoc",LocCamera);
write(file,"XScreenLoc",XlocScreen);
write(file,"Data",Data);
write(file,"theta_phi",theta_phi_rad);
# write(file,"m_true",m);
write(file,"pad",pad);
close(file);
return;
end

export readVisDataFile
function readVisDataFile(filename::String)
file = matopen(string(filename,"_data.mat"), "r");
domain = read(file,"domain");
n = read(file,"n");
Data = read(file,"Data");
# m_true = read(file,"m_true")
theta_phi = read(file,"theta_phi");
XScreenLoc = read(file,"XScreenLoc");
CameraLoc = read(file,"CameraLoc");
ScreenDomain = read(file,"ScreenDomain");
n_screen = read(file,"n_screen");
pad = read(file,"pad");
close(file);
return n,domain,XScreenLoc,CameraLoc,ScreenDomain,n_screen,Data,theta_phi,pad;
end

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
function compute_ERI!(I::IntegralHelper{T, RIFIT, O}, aoints::IntegralHelper{T, RIFIT, AtomicOrbitals}) where {T<:AbstractFloat, O<:AbstractRestrictedOrbitals}
    Bμν = aoints["ERI"]

    C = I.orbitals.C
    if eltype(I.orbitals.C) !== T
        C = T.(C)
    end
    @tensoropt (Q=>50x, p=>x, q=>x, μ=>x, ν=>x) begin
        Bpq[Q,p,q] :=  Bμν[Q, μ, ν]*C[μ, p]*C[ν, q]
    end
    I["ERI"] = Bpq
end

function compute_BOO!(I::IntegralHelper{T, E, O}) where {T<:AbstractFloat, E<:AbstractDFERI, O<:AbstractRestrictedOrbitals}
    # Create AO integral object
    basis = I.orbitals.basis
    aoorbs = AtomicOrbitals(I.molecule, basis)
    aoints = IntegralHelper{T}(molecule=I.molecule, orbitals=aoorbs, basis=basis, eri_type=I.eri_type)
    compute_BOO!(I,aoints)
end

function compute_BVV!(I::IntegralHelper{T, E, O}) where {T<:AbstractFloat, E<:AbstractDFERI, O<:AbstractRestrictedOrbitals}
    # Create AO integral object
    basis = I.orbitals.basis
    aoorbs = AtomicOrbitals(I.molecule, basis)
    aoints = IntegralHelper{T}(molecule=I.molecule, orbitals=aoorbs, basis=basis, eri_type=I.eri_type)
    compute_BVV!(I,aoints)
end

function compute_BOV!(I::IntegralHelper{T, E, O}) where {T<:AbstractFloat, E<:AbstractDFERI, O<:AbstractRestrictedOrbitals}
    # Create AO integral object
    basis = I.orbitals.basis
    aoorbs = AtomicOrbitals(I.molecule, basis)
    aoints = IntegralHelper{T}(molecule=I.molecule, orbitals=aoorbs, basis=basis, eri_type=I.eri_type)
    compute_BOV!(I,aoints)
end

function compute_BOO!(I::IntegralHelper{T, E, O}, aoints::IntegralHelper{T, E, AtomicOrbitals}) where {T<:AbstractFloat, E<:AbstractDFERI, O<:AbstractRestrictedOrbitals}

    Bμν = aoints["ERI"]

    core = Options.get("drop_occ")
    o = (1+core):I.molecule.Nα
    Co = I.orbitals.C[:,o]
    if eltype(I.orbitals.C) !== T
        Co = T.(Co)
    end
    @tensoropt (Q=>50x, i=>x, j=>x, μ=>10x, ν=>10x) begin
        Bij[Q,i,j] :=  Bμν[Q, μ, ν]*Co[μ, i]*Co[ν, j]
    end
    I["BOO"] = Bij
end

function compute_BOV!(I::IntegralHelper{T, E, O}, aoints::IntegralHelper{T, E, AtomicOrbitals}) where {T<:AbstractFloat, E<:AbstractDFERI, O<:AbstractRestrictedOrbitals}

    Bμν = aoints["ERI"]

    inac = Options.get("drop_vir")
    core = Options.get("drop_occ")
    ndocc = I.molecule.Nα
    nbf = size(I.orbitals.C,1)
    o = (1+core):ndocc
    v = (ndocc+1):(nbf - inac)
    Cv = I.orbitals.C[:,v]
    Co = I.orbitals.C[:,o]
    if eltype(I.orbitals.C) !== T
        Cv = T.(Cv)
        Co = T.(Co)
    end
    @tensoropt (Q=>50x, a=>8x, i=>x, μ=>10x, ν=>10x) begin
        Bia[Q,i,a] :=  Bμν[Q, μ, ν]*Co[μ, i]*Cv[ν, a]
    end
    I["BOV"] = Bia
end

function compute_BVV!(I::IntegralHelper{T, E, O}, aoints::IntegralHelper{T, E, AtomicOrbitals}) where {T<:AbstractFloat, E<:AbstractDFERI, O<:AbstractRestrictedOrbitals}

    Bμν = aoints["ERI"]

    inac = Options.get("drop_vir")
    ndocc = I.molecule.Nα
    nbf = size(I.orbitals.C,1)
    v = (ndocc+1):(nbf - inac)
    Cv = I.orbitals.C[:,v]
    if eltype(I.orbitals.C) !== T
        Cv = T.(Cv)
    end
    @tensoropt (Q=>50x, a=>8x, b=>8x, μ=>10x, ν=>10x) begin
        Bab[Q,a,b] :=  Bμν[Q, μ, ν]*Cv[μ, a]*Cv[ν, b]
    end
    I["BVV"] = Bab
end

function compute_OOOO!(I::IntegralHelper{T,E,O}) where {T<:AbstractFloat, E<:AbstractDFERI, O<:AbstractRestrictedOrbitals}
    Boo = I["BOO"]
    @tensoropt (i=>x, j=>x, k=>x, l=>x, a=>10x, b=>10x, c=>10x, d=>10x, Q=>20x) begin
        OOOO[i,j,k,l] :=  Boo[Q, i, j]*Boo[Q, k, l]
    end
    I["OOOO"] = OOOO
end

function compute_OOOV!(I::IntegralHelper{T,E,O}) where {T<:AbstractFloat, E<:AbstractDFERI, O<:AbstractRestrictedOrbitals}
    Boo = I["BOO"]
    Bov = I["BOV"]
    @tensoropt (i=>x, j=>x, k=>x, l=>x, a=>10x, b=>10x, c=>10x, d=>10x, Q=>20x) begin
        OOOV[i,j,k,a] :=  Boo[Q, i, j]*Bov[Q, k, a]
    end
    I["OOOV"] = OOOV
end

function compute_OOVV!(I::IntegralHelper{T,E,O}) where {T<:AbstractFloat, E<:AbstractDFERI, O<:AbstractRestrictedOrbitals}
    Boo = I["BOO"]
    Bvv = I["BVV"]
    @tensoropt (i=>x, j=>x, k=>x, l=>x, a=>10x, b=>10x, c=>10x, d=>10x, Q=>20x) begin
        OOVV[i,j,a,b] :=  Boo[Q, i, j]*Bvv[Q, a, b]
    end
    I["OOVV"] = OOVV
end

function compute_OVOV!(I::IntegralHelper{T,E,O}) where {T<:AbstractFloat, E<:AbstractDFERI, O<:AbstractRestrictedOrbitals}
    Bov = I["BOV"]
    @tensoropt (i=>x, j=>x, k=>x, l=>x, a=>10x, b=>10x, c=>10x, d=>10x, Q=>20x) begin
        OVOV[i,a,j,b] :=  Bov[Q, i, a]*Bov[Q, j, b]
    end
    I["OVOV"] = OVOV
end

function compute_OVVV!(I::IntegralHelper{T,E,O}) where {T<:AbstractFloat, E<:AbstractDFERI, O<:AbstractRestrictedOrbitals}
    Bov = I["BOV"]
    Bvv = I["BVV"]
    @tensoropt (i=>x, j=>x, k=>x, l=>x, a=>10x, b=>10x, c=>10x, d=>10x, Q=>20x) begin
        OVVV[i,a,b,c] :=  Bov[Q, i, a]*Bvv[Q, b, c]
    end
    I["OVVV"] = OVVV
end

function compute_VVVV!(I::IntegralHelper{T,E,O}) where {T<:AbstractFloat, E<:AbstractDFERI, O<:AbstractRestrictedOrbitals}
    Bvv = I["BVV"]
    @tensoropt (i=>x, j=>x, k=>x, l=>x, a=>10x, b=>10x, c=>10x, d=>10x, Q=>20x) begin
        VVVV[a,b,c,d] :=  Bvv[Q, a, b]*Bvv[Q, c, d]
    end
    I["VVVV"] = VVVV
end


text:
const FrontalPole = BilateralStructure("frontal pole")

const FrontopolarCortex = BilateralStructure("frontopolar cortex")

const SuperiorFrontalGyrus = BilateralStructure("superior frontal gyrus")

const MiddleFrontalGyrus = BilateralStructure("middle frontal gyrus")

const RostralMiddleFrontal = BilateralStructure("rostral middle frontal gyrus")

const CaudalMiddleFrontalGyrus = BilateralStructure("caudal middle frontal gyrus")

const ParsOrbitalis = BilateralStructure("pars orbitalis")

const ParsTriangularis = BilateralStructure("pars triangularis")

const ParsOpercularis = BilateralStructure("pars opercularis")

const InferiorFrontalGyrus = BilateralStructure("inferior frontal gyrus")

const PrecentralGyrus = BilateralStructure("precentral gyrus")

const MedialFrontalGyrus = BilateralStructure("medial frontal gyrus")

# update SupplementaryMotorArea name

const ParacentralLobule = BilateralStructure("paracentral lobule")

const OrbitalGyrus = BilateralStructure("orbital gyrus")

const LateralOrbitofrontalGyrus = BilateralStructure("lateral_orbitofrontal gyrus")

const MedialOrbitofrontalGyrus = BilateralStructure("lateral orbitofrontal gyrus")

const FrontoMarginalGyrus = BilateralStructure("fronto-marginal gyrus")

const GyrusRectus = BilateralStructure("gyrus rectus")

const RostralGyrus = BilateralStructure("rostral gyrus")

const FrontalLobe = BilateralStructure("frontal lobe")
#const FrontoMarginalSulcus = Sulcus("fronto-marginal_sulcus")


text:
struct InvertibleScalarFunction{S,T}
    f::S
    f⁻¹::T
end

function Base.:∘(isf1::InvertibleScalarFunction, isf2::InvertibleScalarFunction)
    f = isf1.f ∘ isf2.f
    f⁻¹ = isf2.f⁻¹ ∘ isf1.f⁻¹
    return InvertibleScalarFunction(f, f⁻¹)
end

(isf::InvertibleScalarFunction)(x) = isf.f(x)

function Base.inv(isf::InvertibleScalarFunction)
    f, f⁻¹ = isf.f, isf.f⁻¹
    return InvertibleScalarFunction(f⁻¹, f)
end

const RealVector = AbstractVector{<:Real}

is_constant(column::RealVector) = maximum(column) == minimum(column)

function binner(column::RealVector, bins::RealVector)
    rtol = 1.0e-5
    atol = 1.0e-8
    eps = @. atol + rtol * abs(column)
    return searchsortedfirst.(Ref(bins), column + eps) .- 1
end

function uniform_binner(column::RealVector, num_bins::Integer)
    if is_constant(column)
        return ones(Int46, length(column))
    end
    col_min = minimum(column)
    col_max = maximum(column)
    bins = range(col_min, stop = col_max, length = num_bins)
    return binner(column, bins)
end

function quantile_binner(column::RealVector, num_bins::Integer)
    bins = range(0, stop=1, length=num_bins)
    return binner(column, quantile(column, bins))
end

function log_scaler(_::RealVector)
    return InvertibleScalarFunction(log, exp)
end

function logistic_scaler(_::RealVector)
    sigmoid(x::T) where {T<:Real} = one(T) / (one(T) + exp(-x))
    logit(x::T) where {T<:Real} = log(x / (one(T) - x))
    return InvertibleScalarFunction(sigmoid, logit)
end

function min_max_scaler(column::RealVector, range::RealVector = [0, 1])
    max_val = maximum(column)
    min_val = minimum(column)
    min_range, max_range = range[1], range[2]
    f(x) = (x - min_val) / (max_val - min_val) * (max_range - min_range) + min_range
    f⁻¹(y) = (y - min_range) * (max_val - min_val) / (max_range - min_range) + min_val
    return InvertibleScalarFunction(f, f⁻¹)
end

function max_abs_scaler(column::RealVector)
    max_val = maximum([maximum(column), abs(minimum(column))])
    f(x) = x / max_val
    f⁻¹(y) = y * max_val
    return InvertibleScalarFunction(f, f⁻¹)
end

function standardizer(column::RealVector)
    m = mean(column)
    s = std(column)
    f(x) = (x - m) / s
    f⁻¹(y) = y * s + m
    return InvertibleScalarFunction(f, f⁻¹)
end


text:
# cd("/Volumes/SSD Hans/Github/MacroFall2020/AdvMacro")
# mkpath("Assignment4")
cd("/Volumes/SSD Hans/Github/MacroFall2020/AdvMacro/Assignment4")
# mkpath("graphs")
using Plots
# using LateXStrings # Pkg.add("LaTeXStrings") # https://github.com/stevengj/LaTeXStrings.jl
using Dierckx # Pkg.add("Dierckx") # https://github.com/kbarbary/Dierckx.jl
using Interpolations # Pkg.add("Interpolations") # https://github.com/JuliaMath/Interpolations.jl
using ForwardDiff # Pkg.add("ForwardDiff") # https://github.com/JuliaDiff/ForwardDiff.jl
#Pkg.add("Optim")
using Optim #  # https://julianlsolvers.github.io/Optim.jl/stable/
    using Optim: converged, maximum, maximizer, minimizer, iterations
using Roots # Pkg.add("Roots") # https://github.com/JuliaMath/Roots.jl
using Parameters # Pkg.add("Parameters") # https://github.com/mauro3/Parameters.jl
# Call Scaled Interpolation Functions
    include("./Scaled_Interpolation_Functions.jl")
## Minimum Bracketing Function
    # Minimum Bracketing Function
        # This function comes from Numerical Recipes, Section 10.1
        # Input: numbers a,b that initiliaze bracket and objective function F
            # Optional: Bounds for the brackeet so that a,b,c ∈ [x_min,x_max]
        # Output: numbers (a,b,c) and images (fa,fb,fc) such that a<b<c and fb<fa,fc
        function mnbrak(a,b,F::Function,x_min=-Inf,x_max=Inf)
            # Auxiliariy variables
            Tiny    = 1E-20
            G_limit = 100
            # Define golden ratio
            γ = 1.618034
            # Evaluate function at end points
            fa, fb = F(a), F(b)
            # Swap a and b so that we can go downhill in the direction from a to b.
                # This way we know for certain that fb<fa, we only need to find c such that fb<fc
            if fb>fa
            a , b = b , a
            fa,fb = fb, fa
            end
            # Guess for value c expanding bracket away from b -> (a,b,c) or (c,b,a)
                # Check bounds
            c  = max( min( b + γ*(b-a) , x_max ) , x_min)
            fc = F(c)
            # While fb>fc we need to keep on bracketing
            while fb>fc
                # Compute u (new candidate) by parabolic extrapolation from a, b, c.
                    # TINY is used to prevent any possible division by zero.
                r = (b-a)*(fb-fc)
                q = (b-c)*(fb-fa)
                u = min(max(b-((b-c)*q-(b-a)*r)/(2*sign(q-r)*max(abs(q-r),Tiny)),x_min),x_max)
                u_lim = min(max(b+G_limit*(c-b),x_min),x_max)
                # Test various cases for new candidate
                if ((b-u)*(u-c) > 0) # Parabolic u is between b and c
                    fu=F(u)
                    if (fu < fc) # Got a minimum between b and c so bounds are (b,u,c)
                        a, fa, b, fb = b, fb, u, fu
                        break
                    elseif (fu > fb) # Got a minimum between a and u so bounds are (a,b,u)
                        c, fc = u, fu
                        break
                    end
                    # If you got here is because candidate u failed
                    # Get new candidate by expanding interval with golden ratio
                    # Check bounds
                    u  = max(min( c+γ*(c-b) , x_max ),x_min)
                    fu = F(u)
                elseif ((c-u)*(u-u_lim) > 0.0) # Parabolic u is between c and its limit (ulim)
                    fu=F(u)
                    if (fu < fc) # Drop c and replace it with u, get new u with golden expansion
                        b, c, fb, fc = c, u, fc, fu
                        u  = max(min( c+γ*(c-b) , x_max ),x_min)
                        fu = F(u)
                    end
                elseif ((u-u_lim)*(u_lim-c) >= 0.0) # U is above its limit, reign it in!
                    u  = u_lim
                    fu = F(u)
                else # Nothing worked, use golden expansion
                    u  = max(min( c+γ*(c-b) , x_max ),x_min)
                    fu = F(u)
                end
                # If no break then forget the oldest point and go onto next iteration
                    # This means assigning b->a, c->b, u-> and forgetting about a
                a, b, c, fa, fb, fc = b, c, u, fb, fc, fu
            end
            # Return solution once out of the loop
                # Swap a and c so that a<b<c
                if a>c
                a , c  = c, a
                fa, fc = fc, fa
                end
            # Minimum bracketed in [a,c] with intermediate point b so that fb<fa,fc
            # println("Minimum bracketed in [$a,$c] with intermediate point b=$b \n Function values: F(a)=$fa, F(b)=$fb, F(c)=$fc")
            return a,c,b,fa,fc,fb
        end

## VFI- Continuous

# defining parameters

# Paramters
    # Generate structure for parameters using Parameters module
    # We can set default values for our parameters
    @with_kw struct Par
        # Model Parameters
        z::Float64 = 1    ; # Productivity
        α::Float64 = 1/3  ; # Production function
        β::Float64 = 0.98 ; # Discount factor
        σ::Float64 = 2 ;
        η::Float64 = 1 ;
        δ::Float64 = 0.05 ;
        l::Float64 = 0.4 ;
        M = ((z*α*β)/(1-β+β*δ))^(1/(1-α))
        N = z*M^α-δ*M
        χ = (z*(1-α)*M^α)/(l^(σ+η)*N^σ)
        # VFI Paramters
        max_iter::Int64   = 2000  ; # Maximum number of iterations
        dist_tol::Float64 = 1E-9   ; # Tolerance for distance
        # Howard's Policy Iterations
        H_tol::Float64    = 1E-9  ; # Tolerance for policy function iteration
        N_H::Int64        = 20    ; # Maximum number of policy iterations
        # Minimum consumption for numerical optimization
        c_min::Float64    = 1E-16;
        l_min::Float64    = 1E-16;
        l_max::Float64    = 1;
    end

    p = Par()

    #Steady state constants
    # @unpack z, α, β, σ, η, l, δ = p
    # global M = ((z*α*β)/(1-β+β*δ))^(1/(1-α))
    # global N = z*M^α-δ*M
    # # Getting the right χ for l_ss= 0.4
    # global χ = (z*(1-α)*M^α)/(l^(σ+η)*N^σ)
## Grid
    function Make_K_Grid(n_k,θ_k,p::Par,scale_type="Poly")
        # Get SS
        k_ss,y_ss,c_ss,r_ss,w_ss = SS_values(p)
        # Get k_grid
        if θ_k≠1
            if scale_type=="Poly"
            k_grid = PolyRange(1E-5,2*k_ss;θ=θ_k,N=n_k) ; # Curved grid between 0 and 2*k_ss
            elseif scale_type=="Exp"
            # k_grid = ExpRange(1E-5,2*k_ss;θ=θ_k,N=n_k) ; # Curved grid between 0 and 2*k_ss
            else
            error("scale_type must be either Poly or Exp")
            end
        else
        k_grid = range(1E-5,2*k_ss,length=n_k)
        end
        # Return
        return k_grid
    end

## Model Parameters

# Generate structure of model objects
    @with_kw struct Model
        # Parameters
        p::Par = Par() # Model paramters in their own structure
        # Grids
        θ_k::Float64    = 1     # Curvature of k_grid
        n_k::Int64      = 20    # Size of k_grid
        n_k_fine::Int64 = 1000  # Size of fine grid for interpolation
        k_grid          = Make_K_Grid(n_k,θ_k,p)    # k_grid for model solution
        k_grid_fine     = Make_K_Grid(n_k_fine,1,p) # Fine grid for interpolation
        # Value and policy functions
        V         = Array{Float64}(undef,n_k)       # Value Function
        G_kp      = Array{Float64}(undef,n_k)       # Policy Function
        G_c       = Array{Float64}(undef,n_k)       # Policy Function
        G_l       = Array{Float64}(undef,n_k)       # Policy Function
        V_fine    = Array{Float64}(undef,n_k_fine)  # Value Function on fine grid
        G_kp_fine = Array{Float64}(undef,n_k_fine)  # Policy Function on fine grid
        G_c_fine  = Array{Float64}(undef,n_k_fine)  # Policy Function on fine grid
        G_l_fine  = Array{Float64}(undef,n_k_fine)  # Policy Function on fine grid
        # Anaytical Solutions
        V_a       = Array{Float64}(undef,n_k_fine)  # Analytical Value Function on fine grid
        G_kp_a    = Array{Float64}(undef,n_k_fine)  # Analytical Policy Function on fine grid
        G_c_a     = Array{Float64}(undef,n_k_fine)  # Analytical Policy Function on fine grid
        Euler     = Array{Float64}(undef,n_k_fine)  # Errors in Euler equation
    end

    M = Model()

## Utility function

    function utility(k,kp,l,p::Par)
        @unpack z, α, η, σ, χ, c_min, l_min, δ= p
        y = (a,b) -> z*(a^α)*b^(1-α)+(1-δ)*a
        c = (a,b,d) -> y(a,b)-d
        lb = (a) -> (χ*a^(η+1)/(η+1))
        if l ≤ l_min
            lab = lb(l_min)
            if c(k,kp,l_min) ≤ c_min
                cons = c_min^(1-σ)/(1-σ)
            else
                cons = c(k,kp,l_min)^(1-σ)/(1-σ)
            end
        else
            lab  = lb(l)
            if c(k,kp,l) ≤ c_min
                cons = c_min^(1-σ)/(1-σ)
            else
                cons = c(k,kp,l)^(1-σ)/(1-σ)
            end
        end
        return cons-lab
    end

    # Derivative of utility function wrt labor
    function d_utility_l(k,kp,l,p::Par)
        @unpack α,δ,σ,z,η,c_min,χ = p
        c = z*(k^α)*(l^(1-α))+(1-δ)*k-kp
        d_u = 0
        if c>c_min
            d_u = c^(-σ)
        else
            d_u = c_min^(-σ)
        end
        return d_u*z*(k^α)*(1-α)*(l^(-α))-χ*(l^η)
    end

## Steady State values fn
    # Steady state values (funciton)
    function SS_values(p::Par)
        # This function takes in parameters and provides steady state values
        # Parameters: productivity (z), returns to scale (a) and discount factor (b)
        # Output: values for capital, production, consumption, rental rate, wage
        @unpack z, α, β, l, M, N = p
        l_ss = l
        k_ss = M*l_ss
        c_ss = N*l_ss
        y_ss = z*l_ss*M^α
        r_ss = α*z*M^(α-1)
        w_ss = (1-α)*z*M^α
        return l_ss,k_ss,y_ss,c_ss,r_ss,w_ss
    end

    # Test steady state function
    l_ss,k_ss,y_ss,c_ss,r_ss,w_ss = SS_values(p)
    println(" ")
    println("------------------------")
    println(" Steady State variables")
    println("   Quantities: l = $l_ss; k = $k_ss; y = $y_ss; c = $c_ss;")
    println("   Prices:     r = $r_ss; w = $w_ss;")
    println("------------------------")
    println(" ")

    # Path to steady state
    function path(k,kp,l,M::Model)
        # This function takes in k, k', l and parameters and provides
        # Output: values for capital, production, consumption, rental rate, wage
        @unpack p, n_k = M
        @unpack z, α, β, δ = p
        w = (a,b) -> (1-α)*(a^α)*(b^(-α))
        r = (a,b) -> α*(a^(a-1))*(b^(1-α))
        y = (a,b) -> (a^α)*b^(1-α)+(1-δ)*a
        c = (a,b,d) -> y(a,b)-d
        wage = w.(k,l)
        int = r.(k,l)
        cons = c.(k,kp,l)
        output = y.(k,l)
        labor = get_l.(k,kp,fill(M,length(kp)))
        return wage, int, cons, labor, output
    end

## Optimal labot- Continuous

    # Function to find optimal labor
    function get_l(k,kp,M::Model)::Float64
        @unpack p = M
        L(x) = -utility(k,kp,x,p)
        result = maximize(L,0,1)
        converged(result) || error("Failed to converge in $(iterations(result)) iterations")
        l = maximizer(result)
        return l
    end

## Euler Errors

    # # Euler Equation
    # function Euler_Error(k,kp,kpp,l,lp,p::Par)
    #     # Return percentage error in Euler equation
    #     @unpack z, α, β, σ = p
    #     LHS = 1 ./(z.*k.^α.*l.^(1 .-α).-kp).^σ
    #     RHS = β.*α.*z.*kp.^(α.-1).*lp.^(1 .-α)./(z.*kp.^α.*lp.^(1 .-α).-kpp).^σ
    #     return (RHS/LHS-1)*100
    # end

    # Function that returns the percentage error in the euler equation
    function Euler_Error(k,kp,kpp,l,lp,p::Par)
        # Return percentage error in Euler equation
        @unpack z, α, β, σ, δ = p
        k_in = k.>0
        kp_in = kp.>0
        kpp_in = kpp.>0
        l_in = l.>0
        lp_in = lp.>0
        k = k[k_in.*kp_in.*kpp_in.*l_in.*lp_in]
        kp = kp[k_in.*kp_in.*kpp_in.*l_in.*lp_in]
        kpp = kpp[k_in.*kp_in.*kpp_in.*l_in.*lp_in]
        l = l[k_in.*kp_in.*kpp_in.*l_in.*lp_in]
        lp = lp[k_in.*kp_in.*kpp_in.*l_in.*lp_in]
        LHS = (z.*(k.^α).*(l.^(1-α)).+(1-δ).*k.-kp).^(-σ)
        RHS = β.*(α.*z.*((lp./kp).^(1-α)).+(1-δ)).*((z.*(kp.^α).*(lp.^(1-α)).+(1-δ).*kp.-kpp).^(-σ))
        return (RHS./LHS.-1).*100
    end

    # Euler error by number of points in grid
    function Euler_grid(M::Model)
        @unpack n_k_fine, k_grid_fine, p = M
        # Interpolation of G_kp
        G_kp_ip = ScaledInterpolations(M.k_grid,M.G_kp, BSpline(Cubic(Line(OnGrid()))))
        G_l_ip = ScaledInterpolations(M.k_grid,M.G_l, BSpline(Cubic(Line(OnGrid()))))
        # Euler error of numerical policy function on grid
        Euler = zeros(n_k_fine)
        for i=1:n_k_fine
            k   = k_grid_fine[i]
            kp  = G_kp_ip(k)
            kpp = G_kp_ip(kp)
            l   = G_l_ip(k)
            lp  = G_l_ip(kp)
            Euler[i] = Euler_Error(k,kp,kpp,l,lp,p)
        end
        # Return
        return Euler
    end

## VFI fixed point
    function VFI_Fixed_Point(T::Function,M::Model)
        # Unpack model structure
        @unpack p, n_k, θ_k, k_grid = M
        # VFI paramters
        @unpack max_iter, dist_tol = p
        # Initialize variables for loop
        V_old  = zeros(n_k)     ; # Initialize value function, here I just do 0, not the best option
        # V_old  = utility.(collect(k_grid),zeros(n_k),p) ; # Start at utility with zero savings
        V_dist = 1              ; # Initialize distance
        println(" ")
        println("------------------------")
        println("VFI - n_k=$n_k - θ_k=$θ_k")
        for iter=1:max_iter
            # Update value function
            V_new, G_kp, G_c, G_l = T(Model(M,V=copy(V_old)))
                # println("T(V) = $V_new")
                # println("  V  = $V_old")
            # Update distance and iterations
            V_dist = maximum(abs.(V_new./V_old.-1))
            # Update old function
            V_old  = V_new
            # Report progress
            if mod(iter,100)==0
                println("   VFI Loop: iter=$iter, dist=",100*V_dist,"%")
            end
            # Check convergence and return results
            if V_dist<=dist_tol
                println("VFI - n_k=$n_k - θ_k=$θ_k")
                println("Iterations = $iter and Distance = ",100*V_dist,"%")
                println("------------------------")
                println(" ")
                # Interpolate to fine grid
                V_ip = ScaledInterpolations(M.k_grid,V_new, BSpline(Cubic(Line(OnGrid()))))
                    V_fine = V_ip.(collect(M.k_grid_fine))
                G_kp_ip = ScaledInterpolations(M.k_grid,G_kp, BSpline(Cubic(Line(OnGrid()))))
                    G_kp_fine = G_kp_ip.(collect(M.k_grid_fine))
                G_c_ip = ScaledInterpolations(M.k_grid,G_c, BSpline(Cubic(Line(OnGrid()))))
                    G_c_fine = G_c_ip.(collect(M.k_grid_fine))
                G_l_ip = ScaledInterpolations(M.k_grid,G_l, BSpline(Cubic(Line(OnGrid()))))
                    G_l_fine = G_l_ip.(collect(M.k_grid_fine))
                # Update model
                M = Model(M; V=V_new,G_kp=G_kp,G_c=G_c,G_l=G_l,V_fine=V_fine,G_kp_fine=G_kp_fine,G_c_fine=G_c_fine,G_l_fine=G_l_fine)
                # Euler error
                Euler = Euler_Error(M.k_grid_fine,G_kp_fine,G_kp_ip.(collect(G_kp_fine)),G_l_fine,G_l_ip.(collect(G_kp_fine)),p)
                # Update model
                M = Model(M; Euler=Euler)
                # Return results
                return M::Model
            end
        end
        # If loop ends there was no convergence -> Error!
        error("Error in VFI - Solution not found")
    end

## Bellman Max
# Bellman operator - Continuous Choice- multivariable optimization
function T_cts_max(M::Model)
    @unpack p, n_k, k_grid, V, G_kp, G_c, G_l = M
    @unpack z, α, β, c_min, δ, l_min = p
    l_ss,k_ss,y_ss,c_ss,r_ss,w_ss = SS_values(p)
    # println("Initial V = $V")
    # Define the interpolation of the current value function (V) for values of Vp
    Vp = ScaledInterpolations(k_grid,V, BSpline(Cubic(Line(OnGrid()))))
        # println("Interpolation test: Vp(k_min)=",Vp(k_grid[1])," Vp(k/2)=",Vp((k_grid[1]+k_grid[end])/2))
    # Define the derivative (might not be used, depends on algorithm)
        # For many methods we don't need to use ForwardDiff, we have direct measures,
        # for example for cubic splines. I am using ForwardDiff as an example
    # dVp(x) = ForwardDiff.derivative(Vp,x)

        # println("Interpolation test: Obj_Fun(k_min)=",Obj_Fun(k_grid[1])," Obj_Fun(k/2)=",Obj_Fun((k_grid[1]+k_grid[end])/2))
    get_kp_max(k) = min(z*(k^α)*(1^(1-α)) + (1-δ)*k - c_min , k_grid[end])
    kp_min = k_grid[1]
    l_max = 1.0
    Obj_Fun(k,x1,x2) = -utility(k,x1,x2,p) - β*Vp(x1)
    for i = 1:n_k
        # Min and max kp given current k
        kp_max = get_kp_max(k_grid[i])
        k_grid[i] < k_ss || (kp_min = k_ss)
        k_grid[i] > k_ss || (kp_max = k_ss)
            # Note: we can be smarter here.
            # We know that if k_grid[i]<k_ss then the kp_max is actually k_ss
            # Similarly, if k_grid[i]>k_ss then kp_min = k_ss
        # Check if lower bound binds
        k_0, l_0, k_corner, l_corner = corner(k_grid[i],Vp,Model(n_k=100))
        #x0 = [k_gr, l_0]
            # No need in this problem because without capital there is no production
            # also no consumption. So value would be -Inf.
            # I am leaving two ways of getting the derivative
        # dObj_min = ForwardDiff.derivative(Obj_Fun,kp_min)
        # dObj_min = -d_utility(k_grid[i],kp_min,p) - β*dVp(kp_min)
        # dObj_min = -1
        if k_corner
            # Define objective function: Right hand side of Bellman operator
                # Optimizer will minimize so we need the negative of the function
             V[i]    = Obj_Fun(k_grid[i],k_0,l_0)
             G_kp[i] = k_0
             G_l[i] = l_0
             println("\n corner solution for capital")
             println("k=$k_0 ; l=$l_0; VF=$(V[i])")
        elseif l_corner
            #println("\n corner solution for labor")
            min_result = optimize(x->Obj_Fun(k_grid[i],x,l_0),kp_min,kp_max)
            # Record results
                V[i]     = -min_result.minimum
                G_kp[i]  = min_result.minimizer
                G_l[i]   = l_0
                println("\n corner solution for labor")
                println("k=$(min_result.minimizer); l=$l_0; VF=$(V[i])")
        # # Check if upper bound binds
        #     # This can bind but shouldn't I will leave code to check but comment it
        #     # I am leaving two ways of getting the derivative
        # # dObj_max = ForwardDiff.derivative(Obj_Fun,kp_max)
        # # dObj_max = -d_utility(k_grid[i],kp_max,p) - β*dVp(kp_max)
        # dObj_max = 1
        # if dObj_max<0
        #     V[i]    = -Obj_Fun(kp_max)
        #     G_kp[i] = kp_max
        # else
        # Bracket the solution for k=k_grid[i]
            # In this case this is not necessary, we know the min is bracketed by [kp_min and kp_max]
            # We can still try to use mnbrak to get a better bracket, but in this case the algorithm
            # stops in the initial evaluation, it just verifies that we do have a bracket
            #kp_min, kp_max = mnbrak(kp_min,(kp_max+kp_min)/2,Obj_Fun,kp_min,kp_max)
        else
            # Maximize
            upper = [kp_max,1]
            lower = [kp_min,0]
            init_val = [(kp_min+kp_max)/2,l_ss]
            #inner_optimizer = NelderMead()
            min_result = Optim.optimize(x->Obj_Fun(k_grid[i],x[1],x[2]),lower,upper,init_val,NelderMead(), Optim.Options(iterations=10^6))
                #min_result = optimize(Obj_Fun,kp_min,kp_max)
                #println("maximum = $(-min_result.minimum) with argmin = $(min_result.minimizer) in "*"$(min_result.iterations) iterations")
            # Check result
                converged(min_result) || error("Failed to solve Bellman max in $(iterations(min_result)) iterations")
            # Record results
                V[i]     = -min_result.minimum
                G_kp[i]  = min_result.minimizer[1]
                G_l[i]   = min_result.minimizer[2] #get_l(k_grid[i],G_kp[i],M)
                # println("   Maximization result k_grid[$i] - kp=",min_result.minimizer," V(k)=",min_result.minimum," in $(iterations(min_result)) iterations")
        end # Upper bound if
        # Lower bound if
    end # loop of k_grid[i]
    # Fill in policies
    y = (a,b) -> (a^α)*b^(1-α)+(1-δ)*a
    c = (a,b,d) -> y(a,b)-d
    G_c = c.(k_grid,G_kp,G_l)
    # Return Results
        # println("T(V) = $V")
    return V, G_kp, G_c, G_l
end

## Graphs
# Graphs
function VFI_Graphs(M::Model,VFI_Type)
    gr()
    # Value Function Analytical vs 200
        plot(M.k_grid_fine,M.k_grid_fine,linewidth=4,label = "Analytical Solution",title="Value Function",legend=(0.75,0.2),foreground_color_legend = nothing,background_color_legend = nothing)
        plot(M.k_grid,M.V,linetype=:scatter,marker=(:diamond,4),markercolor=RGB(0.5,0.1,0.1),label="VFI - n_k=$(M.n_k) - θ_k=$(M.θ_k)")
        plot!(M.k_grid_fine,M.V_fine,linewidth=2.5,linestyle=(:dash),linecolor=RGB(0.4,0.4,0.4),label=nothing)
        xlabel!("Capital")
        ylabel!("Value")
        savefig("./graphs/VFI_"*VFI_Type*"_V_$(M.n_k)_$(M.θ_k).png")
    # Capital Policy Function Analytical vs 200
        plot(M.k_grid_fine,M.k_grid_fine,lw=1,linecolor=RGB(0.6,0.6,0.6),label=nothing,title="Policy Function - K",legend=(0.75,0.2),foreground_color_legend = nothing,background_color_legend = nothing)
        #plot!(M.k_grid_fine,M.G_kp_a,linewidth=3,label = "Analytical Solution",title="Policy Function - K",legend=(0.75,0.2),foreground_color_legend = nothing,background_color_legend = nothing)
        plot!(M.k_grid,M.G_kp,linetype=:scatter,marker=(:diamond,4),markercolor=RGB(0.5,0.1,0.1),label="VFI - n_k=$(M.n_k) - θ_k=$(M.θ_k)")
        plot!(M.k_grid_fine,M.G_kp_fine,linewidth=2.5,linestyle=(:dash),linecolor=RGB(0.4,0.4,0.4),label=nothing)
        xlabel!("Capital")
        ylabel!("Capital")
    savefig("./graphs/VFI_"*VFI_Type*"_G_kp_$(M.n_k)_$(M.θ_k).png")
        # Euler Error 200
        plot(M.k_grid_fine,zeros(M.n_k_fine),lw=1,linecolor=RGB(0.6,0.6,0.6),label=nothing,title = "Euler Equation Error (%)",legend=(0.75,0.2),foreground_color_legend = nothing,background_color_legend = nothing)
        plot!(M.k_grid_fine,M.Euler,linetype=:scatter,label="VFI - n_k=$(M.n_k) - θ_k=$(M.θ_k)")
        xlabel!("Capital")
        ylabel!("Percentage Points")
    savefig("./graphs/VFI_"*VFI_Type*"_Euler_$(M.n_k)_$(M.θ_k).png")

    println("\n     Graphs Completed for VFI_$VFI_Type - n_k=$(M.n_k) - θ_k=$(M.θ_k) \n")
end

## Initial Guess

function corner(k,V::AbstractInterpolation,M::Model)
    @unpack n_k, k_grid = M
    l = fill(0.0,n_k)
    V_aux = fill(0.0,n_k)
    for i in 1:n_k
        l[i] = get_l(k,k_grid[i],M)
        V_aux[i] = utility(k,k_grid[i],l[i],p)+V(k_grid[i])
    end
    index = findmax(V_aux)[2]
    k_0 = k_grid[index]
    l_0 = l[index]
    if k_0 == k_grid[end] || k_0 == k_grid[1]
        #print("Corner solution for capital \n")
        k_corner = true
    else
        k_corner = false
    end
    if l_0 == 1 || l_0 == 0
        #print("Corner solution for labor \n")
        l_corner = true
    else
        l_corner = false
    end
    return k_0, l_0, k_corner, l_corner
end

## Bellman operator univariate

# Bellman operator for the continuous choice capital only where labor is solved for with the intratemporal FOC with bisection
function T_univariate_max(M::Model)
    @unpack p, n_k, k_grid, V, G_kp, G_c, G_l = M
    @unpack β,α,z,δ,η,σ,c_min = p
    # Interpolate current iteration of value function v(k') so I can evaluate it at any k'
    Vp = ScaledInterpolations(k_grid,V, BSpline(Cubic(Line(OnGrid())))) # Impose monotonicity because I know it is increasing in capital
    # Define boundaries on capital tomorrow k' and on labor l
    l_min = 1E-16
    l_max = 1.0
    get_kp_max(k,l) = z*(k^α)*(l^(1-α)) + (1-δ)*k - c_min # Function because the max k' depends on the value of capital today k and labor l
    # Define function that finds optimal labor l given (k,k') and returns objective function
    function Obj_fn(k,kp,p::Par)
        min_result = optimize(x->d_utility_l(k,kp,x,p).^2,l_min,l_max,Brent())
        l = min_result.minimizer
        return -utility(k,kp,l,p) - β*Vp.(kp),l
    end
    # optimization of the objective function for each capital level in the grid
    # Initialization
    #V_new = zeros(n_k)
    kp_min = 0.0001
    # Maximize
    for i in 1:n_k
        kp_max = min(get_kp_max(k_grid[i],1.0),0.9999*k_grid[end])
        min_result = optimize(x->Obj_fn(k_grid[i],x,p)[1],kp_min,kp_max,Brent())
        # Check result
        #converged(min_result) || error("Failed to solve Bellman max for capital =" k_grid[i]" in $(iterations(min_result)) iterations")
        # Record results
        V[i] = -min_result.minimum
        G_kp[i] = min_result.minimizer
        G_l[i] = Obj_fn(k_grid[i],G_kp[i],p::Par)[2]
    end
    # Fill in policy for consumption
    G_c = z.*(collect(k_grid).^α).*(G_l.^(1-α)) .- G_kp
    # Return results
    return V, G_kp, G_c, G_l
end

## Solve

    # Execute VFI and plot graphs

    # Execute Numerical VFI - equally spaced grid
        @time M_20  = VFI_Fixed_Point(T_cts_max,Model(n_k=20))
        @time M_50  = VFI_Fixed_Point(T_cts_max,Model(n_k=50))
        @time M_100 = VFI_Fixed_Point(T_cts_max,Model(n_k=100))
        # Graphs
            VFI_Graphs(M_20,"cts_max")
            VFI_Graphs(M_50,"cts_max")
            VFI_Graphs(M_100,"cts_max")
    # Execute Numerical VFI - curved spaced grid
        @time M_20c  = VFI_Fixed_Point(T_cts_max,Model(n_k=20,θ_k=2.5))
        @time M_50c  = VFI_Fixed_Point(T_cts_max,Model(n_k=50,θ_k=2.5))
        @time M_100c = VFI_Fixed_Point(T_cts_max,Model(n_k=100,θ_k=3.5))
        # Graphs
            VFI_Graphs(M_20c,"cts_max")
            VFI_Graphs(M_50c,"cts_max")
            VFI_Graphs(M_100c,"cts_max")

    # Execute Numerical VFI - equally spaced grid
        @time M_20u  = VFI_Fixed_Point(T_univariate_max,Model(n_k=20))
        @time M_50u  = VFI_Fixed_Point(T_univariate_max,Model(n_k=50))
        @time M_100u = VFI_Fixed_Point(T_univariate_max,Model(n_k=100))
        # Graphs
            VFI_Graphs(M_20u,"cts_max")
            VFI_Graphs(M_50u,"cts_max")
            VFI_Graphs(M_100u,"cts_max")
    # Execute Numerical VFI - curved spaced grid
        @time M_20cu  = VFI_Fixed_Point(T_univariate_max,Model(n_k=20,θ_k=2.5))
        @time M_50cu  = VFI_Fixed_Point(T_univariate_max,Model(n_k=50,θ_k=2.5))
        @time M_100cu = VFI_Fixed_Point(TT_univariate_max,Model(n_k=100,θ_k=2.5))
        # Graphs
            VFI_Graphs(M_20cu,"cts_max")
            VFI_Graphs(M_50cu,"cts_max")
            VFI_Graphs(M_100cu,"cts_max")

plot(M_20.k_grid_fine,M_20.G_kp_fine)
plot(M_20.k_grid,M_20.G_l)

max_stars_repo_name: py-r-hans/AdvMacro

max_issues_repo_name: py-r-hans/AdvMacro

max_forks_repo_name: py-r-hans/AdvMacro


text:
# this gives the d and e so that
# [Z A B;
#  0 d e]
#
# is the fixed point
function tail_de(a::AbstractVector{T}; branch=findmax) where T<:Real
    m = length(a)
    C = [view(a,m-1:-1:1) Vcat(-a[end]*Eye(m-2), Zeros{T}(1,m-2))]
    λ, V = eigen(C)
    n2, j = branch(abs2.(λ))
    isreal(λ[j]) || throw(DomainError(a, "Real-valued QL factorization does not exist. Try ql(complex(A)) to see if a complex-valued QL factorization exists."))
    n2 ≥ a[end]^2 || throw(DomainError(a, "QL factorization does not exist. This could indicate that the operator is not Fredholm or that the dimension of the kernel exceeds that of the co-kernel. Try again with the adjoint."))
    c = sqrt((n2 - a[end]^2)/real(V[1,j])^2)
    c*real(V[end:-1:1,j])
end

function tail_de(a::AbstractVector{T}; branch=findmax) where T
    m = length(a)
    C = [view(a,m-1:-1:1) Vcat(-a[end]*Eye(m-2), Zeros{T}(1,m-2))]
    λ, V = eigen(C)::Eigen{float(T),float(T),Matrix{float(T)},Vector{float(T)}}
    n2, j = branch(abs2.(λ))
    n2 ≥ abs2(a[end]) || throw(DomainError(a, "QL factorization does not exist. This could indicate that the operator is not Fredholm or that the dimension of the kernel exceeds that of the co-kernel. Try again with the adjoint."))
    c_abs = sqrt((n2 - abs2(a[end]))/abs2(V[1,j]))
    c_sgn = -sign(λ[j])/sign(V[1,j]*a[end-1] - V[2,j]*a[end])
    c_sgn*c_abs*V[end:-1:1,j]
end

# this calculates the QL decomposition of X and corrects sign
function ql_X!(X)
    s = sign(real(X[2,end]))
    F = ql!(X)
    if s ≠ sign(real(X[1,end-1])) # we need to normalise the sign if ql! flips it
        F.τ[1] = 2 - F.τ[1] # 1-F.τ[1] is the sign so this flips it
        X[1,1:end-1] *= -1
    end
    F
end

function ql(Op::TriToeplitz{T}; kwds...) where T<:Real
    Z,A,B = Op.dl.value, Op.d.value, Op.du.value
    d,e = tail_de([Z,A,B]; kwds...) # fixed point of QL but with two Qs, one that changes sign
    X = [Z A B; zero(T) d e]
    F = ql_X!(X)
    t,ω = F.τ[2],X[1,end]
    QL(_BandedMatrix(Hcat([zero(T), e, X[2,2], X[2,1]], [ω, X[2,3], X[2,2], X[2,1]] * Ones{T}(1,∞)), ∞, 2, 1), Vcat(F.τ[1],Fill(t,∞)))
end

ql(Op::TriToeplitz{T}) where T = ql(InfToeplitz(Op))

# ql for Lower hessenberg InfToeplitz
function ql_hessenberg(A::InfToeplitz{T}; kwds...) where T
    l,u = bandwidths(A)
    @assert u == 1
    a = reverse(A.data.args[1])
    de = tail_de(a; kwds...)
    X = [transpose(a); zero(T) transpose(de)]::Matrix{float(T)}
    F = ql_X!(X) # calculate data for fixed point
    factors = _BandedMatrix(Hcat([zero(T); X[1,end-1]; X[2,end-1:-1:1]], [0; X[2,end:-1:1]] * Ones{float(T)}(1,∞)), ∞, l+u, 1)
    QLHessenberg(factors, Fill(F.Q,∞))
end

# remove one band of A
function ql_pruneband(A; kwds...)
    l,u = bandwidths(A)
    A_hess = A[:,u:end]
    Q,L = ql_hessenberg(A_hess; kwds...)
    p = size(_pertdata(bandeddata(parent(L))),2)+u+1 # pert size
    dat = (UpperHessenbergQ((Q').q[1:(p+l)])) * A[1:p+l+1,1:p]
    pert = Array{eltype(dat)}(undef, l+u+1,size(dat,2)-1)
    for j = 1:u
        pert[u-j+1:end,j] .= view(dat,1:l+j+1,j)
    end
    for j = u+1:size(pert,2)
        pert[:,j] .= view(dat,j-u+1:j+l+1,j)
    end
    H = _BandedMatrix(Hcat(pert, dat[end-l-u:end,end]*Ones{eltype(dat)}(1,∞)), ∞, l+1,u-1)
    Q,H
end

# represent Q as a product of orthogonal operations
struct ProductQ{T,QQ<:Tuple} <: AbstractQ{T}
    Qs::QQ
end

ProductQ(Qs::AbstractMatrix...) = ProductQ{mapreduce(eltype,promote_type,Qs),typeof(Qs)}(Qs)

adjoint(Q::ProductQ) = ProductQ(reverse(map(adjoint,Q.Qs))...)

size(Q::ProductQ, dim::Integer) = size(dim == 1 ? Q.Qs[1] : last(Q.Qs), dim == 2 ? 1 : dim)

function lmul!(Q::ProductQ, v::AbstractVecOrMat)
    for j = length(Q.Qs):-1:1
        lmul!(Q.Qs[j], v)
    end
    v
end

getindex(Q::ProductQ{<:Any,<:Tuple{Vararg{LowerHessenbergQ}}}, i::Integer, j::Integer) = (Q')[j,i]'

function (*)(A::ProductQ{T}, x::AbstractVector{S}) where {T,S}
    TS = promote_op(matprod, T, S)
    lmul!(A, Base.copymutable(convert(AbstractVector{TS},x)))
end

# LQ where Q is a product of orthogonal operations
struct QLProduct{T,QQ<:Tuple,LL} <: Factorization{T}
    Qs::QQ
    L::LL
end

QLProduct(Qs::Tuple, L::AbstractMatrix{T}) where T = QLProduct{T,typeof(Qs),typeof(L)}(Qs, L)
QLProduct(F::QLHessenberg) = QLProduct(tuple(F.Q), F.L)

# iteration for destructuring into components
Base.iterate(S::QLProduct) = (S.Q, Val(:L))
Base.iterate(S::QLProduct, ::Val{:L}) = (S.L, Val(:done))
Base.iterate(S::QLProduct, ::Val{:done}) = nothing

QLProduct{T}(A::QLProduct) where {T} = QLProduct(convert.(AbstractMatrix{T}, A.Qs), convert(AbstractMatrix{T}, A.L))
Factorization{T}(A::QLProduct{T}) where {T} = A
Factorization{T}(A::QLProduct) where {T} = QLProduct{T}(A)
AbstractMatrix(F::QLProduct) = F.Q * F.L
AbstractArray(F::QLProduct) = AbstractMatrix(F)
Matrix(F::QLProduct) = Array(AbstractArray(F))
Array(F::QLProduct) = Matrix(F)

function show(io::IO, mime::MIME{Symbol("text/plain")}, F::QLProduct)
    summary(io, F); println(io)
    println(io, "Q factor:")
    show(io, mime, F.Q)
    println(io, "\nL factor:")
    show(io, mime, F.L)
end

@inline getL(F::QLProduct) = getfield(F, :L)
@inline getQ(F::QLProduct) = ProductQ(F.Qs...)

function getproperty(F::QLProduct, d::Symbol)
    if d == :L
        return getL(F)
    elseif d == :Q
        return getQ(F)
    else
        getfield(F, d)
    end
end

Base.propertynames(F::QLProduct, private::Bool=false) =
    (:L, :Q, (private ? fieldnames(typeof(F)) : ())...)

function _inf_ql(A::AbstractMatrix{T}; kwds...) where T
    _,u = bandwidths(A)
    u ≤ 0 && return QLProduct(tuple(Eye{float(T)}(∞)), A)
    u == 1 && return QLProduct(ql_hessenberg(A; kwds...))
    Q1,H1 = ql_pruneband(A; kwds...)
    F̃ = ql(H1; kwds...)
    QLProduct(tuple(Q1, F̃.Qs...), F̃.L)
end

ql(A::InfToeplitz; kwds...) = _inf_ql(A; kwds...)
ql(A::PertToeplitz; kwds...) = _inf_ql(A; kwds...)

ql(A::Adjoint{<:Any,<:InfToeplitz}) = ql(BandedMatrix(A))
ql(A::Adjoint{<:Any,<:PertToeplitz}) = ql(BandedMatrix(A))


text:
module _IncrEn
export IncrEn
using Statistics: std
    """
        Incr = IncrEn(Sig)

    Returns the increment entropy (`Incr`) estimate of the data sequence
    (`Sig`) using the default parameters:
    embedding dimension = 2, time delay = 1, quantifying resolution = 4,
    logarithm = base 2,

        Incr = IncrEn(Sig::AbstractArray{T,1} where T<:Real; m::Int=2, tau::Int=1, R::Int=4, Logx::Real=2, Norm::Bool=false)

    Returns the increment entropy (`Incr`) estimate of the data sequence
    (`Sig`) using the specified 'keyword' arguments:

    # Arguments:
    `m`     - Embedding Dimension, an integer > 1   \n
    `tau`   - Time Delay, a positive integer    \n
    `R`     - Quantifying resolution, a positive scalar    \n
    `Logx`  - Logarithm base, a positive scalar (enter 0 for natural log) \n
    `Norm`  - Normalisation of IncrEn value: \n
              [false]  no normalisation - default
              [true]   normalises w.r.t embedding dimension (m-1).

    # See also `PermEn`, `SyDyEn`, `MSEn`

    # References:
        [1] Xiaofeng Liu, et al.,
            "Increment entropy as a measure of complexity for time series."
            Entropy
            18.1 (2016): 22.1.

        ***   "Correction on Liu, X.; Jiang, A.; Xu, N.; Xue, J. - Increment
            Entropy as a Measure of Complexity for Time Series,
            Entropy 2016, 18, 22."
            Entropy
            18.4 (2016): 133.

        [2] Xiaofeng Liu, et al.,
            "Appropriate use of the increment entropy for
            electrophysiological time series."
            Computers in biology and medicine
            95 (2018): 13-23.

    """
    function IncrEn(Sig::AbstractArray{T,1} where T<:Real; m::Int=2, tau::Int=1,
        R::Int=4, Logx::Real=2, Norm::Bool=false)

    Logx == 0  ? Logx = exp(1) : nothing

    (size(Sig,1) > 10) ? nothing :  error("Sig:   must be a numeric vector")
    (m > 1) ? nothing :  error("m:     must be an integer > 1")
    (tau>0) ? nothing :  error("tau:   must be an integer > 0")
    (R > 0) ? nothing :  error("R:     must be a positive integer > 0")
    (Logx>0) ? nothing : error("Logx:  must be a positive number > 0")

    Vi = diff(Sig)
    N = size(Vi,1)-((m-1)*tau)
    Vk = zeros(N,m)
    for k = 1:m
        Vk[:,k] = Vi[1+(k-1)*tau:N+(k-1)*tau]
    end

    Sk = sign.(Vk)
    Temp = std(Vk,dims=2)[:]
    Qk = min.(R, floor.((abs.(Vk)*R)./repeat(Temp,outer=(1,m))))
    Qk[any(Temp.==0,dims=2), :] .= 0  #should that be all()
    Wk = Sk.*Qk
    Wk[Wk.==-0] .= 0
    Px = unique(Wk,dims=1)
    Counter = zeros(Int,size(Px,1));
    for k = 1:size(Px,1)
        Counter[k] = sum(all(Wk .- transpose(Px[k,:]) .==0 ,dims=2))
    end
    Ppi = Counter/N

    if size(Px,1) > (2*R + 1)^m
        @warn("Error with probability estimation'")
    elseif round(sum(Ppi),digits=5) != 1
        @warn("Error with probability estimation")
    end
    Incr = -sum(Ppi.*(log.(Logx, Ppi)))
    if Norm
        Incr = Incr/(m-1);
    end

    return Incr
    end

end

"""
Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub
"""


text:
########### Defines the parameters & plot for Fig 3 b,d,f,h,j ###########
c = 4.0; # set the host contact rate as a constant 4.0
# import required variables/modules common to all simulations
include("../Model_scripts/fig3params.jl");
include("../Model_scripts/Functions.jl");
# now that we have all the requirements, let's run the simulation
run_simulation();
Np = N';
########### Figure 3b ###########
# plot the number of infected/uninfected hosts
labels = vcat("Infected", ["" for i in 1:1:n_parasites]);
plot_population_numbers(
    hcat(labels...),
    "c=4.0\n \nNumber of infected and uninfected hosts",
    100,
    size(Np)[1],
    "$figure_directory/graph_3b.png"
)
########### Figure 3d ###########
# calculate the weighted average mortality
ui_w_avg = calculate_mortality()
# calculate the weighted average β
βy_w_avg = calculate_horizontal_transmission()
# calculated the weighted H0
H0_w = calculate_H0()
# calculate the weighted V0
V0_w = calculate_V0()
# plot the average weighted R0
plot_reproductive_rate(5, size(Np)[1], "$figure_directory/graph_3d.png")
########### Figure 3f ###########
# plot the average V0
plot_vertical_reproductive_ratio(1, size(Np)[1], "$figure_directory/graph_3f.png")
########### Figure 3h ###########
# calculate the average virulence
vir_w_avg = calculate_average_virulence()
# plot the average horizontal transmission & weighted virulence
plot_horizontal_and_virulence(1, size(Np)[1], "$figure_directory/graph_3h.png")
########### Figure 3j ###########
# calculate & plot the population evenness through time
plot_evenness(1, size(Np)[1], "$figure_directory/graph_3j.png")


text:
module Shared_funcs
#
  # Waiting for https://github.com/JuliaLang/julia/pull/6884
  using GZip
#
  function format_percent(v)
    return strip(strip(@sprintf("%0.2f", 100.0 * v), '0'), '.') * "%"
  end
  function comp_density(A::AbstractMatrix)
    return countnz(A) / length(A)
  end
  function _read_matrix_mtx(filename, infoonly::Bool=false)
    function skewsymmetric!(M::AbstractMatrix)
        m,n = size(M)
        m == n || throw(DimensionMismatch())
        for i=1:n, j=1:n
            if M[i,j] != 0
                M[j,i] = -M[i,j]
            end
        end
        return M
    end
    function symmetric!(M::AbstractMatrix)
        m,n = size(M)
        m == n || throw(DimensionMismatch())
        for i=1:n, j=1:n
            if M[i,j] != 0
                M[j,i] = M[i,j]
            end
        end
        return M
    end
    function hermitian!(M::AbstractMatrix)
        m,n = size(M)
        m == n || throw(DimensionMismatch())
        for i=1:n, j=1:n
            if M[i,j] != 0
                M[j,i] = conj(M[i,j])
            end
        end
        return M
    end
    #println(pwd())
    mmfile = Nothing()
    if (endswith(filename, ".gz"))
      mmfile = GZip.open(filename,"r")
    else
      mmfile = open(filename,"r")
    end
    # Read first line
    firstline = chomp(readline(mmfile))
    tokens = split(firstline)
    if length(tokens) != 5
      throw(ParseError(string("Not enough words on first line: ", ll)))
    end
    if tokens[1] != "%%MatrixMarket"
      throw(ParseError(string("Not a valid MatrixMarket header:", ll)))
    end
    (head1, rep, field, symm) = map(lowercase, tokens[2:5])
    if head1 != "matrix"
        throw(ParseError("Unknown MatrixMarket data type: $head1 (only \"matrix\" is supported)"))
    end
    eltype = field == "real" ? Float64 :
             field == "complex" ? Complex128 :
             field == "pattern" ? Bool :
             throw(ParseError("Unsupported field $field (only real and complex are supported)"))

    symlabel = symm == "general" ? identity :
               symm == "symmetric" ? symmetric! :
               symm == "hermitian" ? hermitian! :
               symm == "skew-symmetric" ? skewsymmetric! :
               throw(ParseError("Unknown matrix symmetry: $symm (only general, symmetric, skew-symmetric and hermitian are supported)"))
    # Skip all comments and empty lines
    ll   = readline(mmfile)
    while length(chomp(ll))==0 || (length(ll) > 0 && ll[1] == '%')
        ll = readline(mmfile)
    end
    # Read matrix dimensions (and number of entries) from first non-comment line
    dd = map(parseint, split(ll))
    if length(dd) < (rep == "coordinate" ? 3 : 2)
        throw(ParseError(string("Could not read in matrix dimensions from line: ", ll)))
    end
    rows = dd[1]
    cols = dd[2]
    entries = (rep == "coordinate") ? dd[3] : (rows * cols)
    if infoonly
        return (rows, cols, entries, rep, field, symm)
    end
    if rep == "coordinate"
        rr = Array(Int, entries)
        cc = Array(Int, entries)
        xx = Array(eltype, entries)
        for i in 1:entries
            flds = split(readline(mmfile))
            rr[i] = parseint(flds[1])
            cc[i] = parsefloat(flds[2])
            if eltype == Complex128
                xx[i] = Complex128(parsefloat(flds[3]), parsefloat(flds[4]))
            elseif eltype == Float64
                xx[i] = parsefloat(flds[3])
            else
                xx[i] = true
            end
        end
        return symlabel(sparse(rr, cc, xx, rows, cols))
    end
    return symlabel(reshape([parsefloat(readline(mmfile)) for i in 1:entries], (rows,cols)))
  end
  function read_matrix_mtx(filename)
    return _read_matrix_mtx(filename)
  end
  function random_matrix(seed, scale, dense, n, m, show)
    function toindex(x,n) return clamp(int((n)* 0.5*(1.0+x)), 1, n) end
    rng = MersenneTwister(seed)
    B = A = scale * rand(rng, (m, n))
    if dense != 100
      #zeros = rand(rng, 1:(m*n), zero_n) # Needs next version of Julia
      zero_n = int(round((m*n) * (1.0 - dense / 100.0)))
      if (zero_n > 0)
        nz_i = [1:m*n]
        while length(nz_i) + zero_n > (m*n)
          splice!(nz_i, toindex(rand(rng), length(nz_i)))
        end
        B = zeros(Float64, (m,n))
        for i in nz_i B[i] = A[i] end
      end
    end
    if (show)
      println("A", B)
    end
    return B
  end
  function random_matrix(params)
    seed = int(get(params, "seed", int(time_ns()) % 32768))
    scale = int(get(params, "scale", 1.0))
    dense = int(get(params, "dense", 100 ))
    n = int(get(params, "n", 10))
    m = int(get(params, "m", n))
    show = get(params, "show", false)
    println("/seed:", seed)
    return random_matrix(seed, scale, dense, n, m, show)
  end
  function str_range(str, list = [])
  	function single_range(str)
      beg = int(split(str, "-")[1])
      nd = int(split(split(str, "-")[2],":")[1])
      stp = contains(str,":") ? int(split(str, ":")[2]) : 1
      cnt = int((nd-beg+1)/stp)
    	return [range(beg,stp,cnt)]
    end
    if str == "all" return 1:length(list) end
    rngs = Array(Int, 0)
    for x in split(str, ",")
      if length(x) != 0 rngs = vcat(rngs, contains(x, "-") ? single_range(x) : [int(x)]) end
    end
    return rngs
  end
  function print_choose(list, autochoose::Bool = false)
    if length(list) == 0 return [] end
    if length(list) == 1 && autochoose return list end
    cols = [:white, :yellow]
    for i in 1:length(list)
      print_with_color(cols[1+i%2], " $i. $(list[i]) \n")
    end
    print(" > ")
    return str_range(strip(readline()), list)
  end
end

max_stars_repo_name: jadnohra/mlcp-lib

max_issues_repo_name: jadnohra/mlcp-lib

max_forks_repo_name: jadnohra/mlcp-lib


text:
struct Proposal{T<:ProposalStyle, P}
    type :: T
    proposal :: P
end

# Random draws
Base.rand(p::Proposal{<:ProposalStyle, <:Distribution}) = rand(p.proposal)
function Base.rand(p::Proposal{<:ProposalStyle, <:AbstractArray})
    return map(rand, p.proposal)
end

# Densities
function Distributions.logpdf(p::Proposal{<:ProposalStyle, <:UnivariateDistribution}, v)
    return sum(logpdf(p.proposal, v))
end
function Distributions.logpdf(p::Proposal{<:ProposalStyle, <:MultivariateDistribution}, v)
    return sum(logpdf(p.proposal, v))
end
function Distributions.logpdf(p::Proposal{<:ProposalStyle, <:MatrixDistribution}, v)
    return sum(logpdf(p.proposal, v))
end
function Distributions.logpdf(p::Proposal{<:ProposalStyle, <:AbstractArray}, v)
    return sum(map(x -> logpdf(x[1], x[2]), zip(p.proposal, v)))
end
function Distributions.logpdf(p::Proposal{<:ProposalStyle, <:Function}, v)
    return logpdf(p.proposal(v), v)
end

###############
# Random Walk #
###############

function propose(
    proposal::Proposal{RandomWalk, <:Distribution},
    model::DensityModel,
    t
)
    return t + rand(proposal)
end

function propose(
    proposal::Proposal{RandomWalk, <:AbstractArray},
    model::DensityModel,
    t
)
    return t + rand(proposal)
end

function q(
    proposal::Proposal{RandomWalk, <:Distribution},
    t,
    t_cond
)
    return sum(logpdf(proposal, t - t_cond))
end

function q(
    proposal::Proposal{RandomWalk, <:AbstractArray},
    t,
    t_cond
)
    return sum(logpdf(proposal, t - t_cond))
end

##########
# Static #
##########

propose(p::Proposal{RandomWalk}, m::DensityModel) = propose(Proposal(Static(), p.proposal), m)
function propose(
    proposal::Proposal{Static, <:Distribution},
    model::DensityModel,
    t=nothing
)
    return rand(proposal)
end

function propose(
    p::Proposal{Static, <:AbstractArray},
    model::DensityModel,
    t=nothing
)
    props = map(x -> rand(x), p.proposal)
    return props
end

function q(
    proposal::Proposal{Static, <:Distribution},
    t,
    t_cond
)
    return sum(logpdf(proposal, t))
end

function q(
    proposal::Proposal{Static, <:AbstractArray},
    t,
    t_cond
)
    return sum(logpdf(proposal, t))
end

############
# Function #
############

function propose(
    proposal::Proposal{<:ProposalStyle, <:Function},
    model::DensityModel
)
    p = proposal.proposal
    return rand(proposal.proposal())
end

function propose(
    proposal::Proposal{<:ProposalStyle, <:Function},
    model::DensityModel,
    t
)
    p = proposal.proposal
    return rand(proposal.proposal(t))
end

function q(
    proposal::Proposal{<:ProposalStyle, <:Function},
    t,
    t_cond
)
    p = proposal.proposal
    return sum(logpdf.(p(t_cond), t))
end


text:
using Flux, MLDataPattern, Mill, JsonGrinder, JSON, HierarchicalUtils, StatsBase
using JsonGrinder: DictEntry, Entry

# load files in data/documents and parse them
sch = JsonGrinder.schema(readdir("data/documents", join=true), x->open(JSON.parse, x))
# suggest default extractor with some keys as field
extractor = suggestextractor(sch)
# show whole extractor
printtree(extractor)
# we see that there are some dictionaries with lots of keys, let's examine schema
# list_lens lets us iterate over all elements in a way we know their position in schema
# this prints lengths of children of all dict entries.
for i in list_lens(sch)
    e = get(sch, i)
    if e isa DictEntry
        @info i length(e.childs)
    end
end

# that's a lots of numbers, let's see histogram
length_hist = StatsBase.countmap([length(get(sch, i).childs) for i in list_lens(sch) if get(sch, i) isa DictEntry])

# we see highest lengths are 103 and 13, let's set 13 as a threshold
extractor = suggestextractor(sch, (; key_as_field=13))
# show new extractor
printtree(extractor)
# this extractor looks much better
# but still, some values are very sparse,
# let's print all parts of schema where each value is observed only once
for i in list_lens(sch)
    e = get(sch, i)
    if e isa Entry && maximum(values(e.counts)) == 1
        @info i
    end
end

#  we can see lots of leaves under `bib_entries`, which is cased by uniqueness of keys here
# but apart from that, we can see other interesting fields
# [ Info: (@lens _.childs[:metadata].childs[:authors].items.childs[:middle].items)
# [ Info: (@lens _.childs[:metadata].childs[:authors].items.childs[:last])
# [ Info: (@lens _.childs[:metadata].childs[:authors].items.childs[:affiliation].childs[:location].childs[:region])
# [ Info: (@lens _.childs[:paper_id])
# [ Info: (@lens _.childs[:body_text].items.childs[:text])
# [ Info: (@lens _.childs[:body_text].items.childs[:ref_spans].items.childs[:start])
# [ Info: (@lens _.childs[:body_text].items.childs[:ref_spans].items.childs[:end])
# [ Info: (@lens _.childs[:back_matter].items.childs[:text])
# [ Info: (@lens _.childs[:back_matter].items.childs[:cite_spans].items.childs[:ref_id])
# [ Info: (@lens _.childs[:back_matter].items.childs[:cite_spans].items.childs[:start])
# [ Info: (@lens _.childs[:back_matter].items.childs[:cite_spans].items.childs[:text])
# [ Info: (@lens _.childs[:back_matter].items.childs[:cite_spans].items.childs[:end])
# [ Info: (@lens _.childs[:back_matter].items.childs[:ref_spans].items.childs[:start])
# [ Info: (@lens _.childs[:back_matter].items.childs[:ref_spans].items.childs[:text])
# [ Info: (@lens _.childs[:back_matter].items.childs[:ref_spans].items.childs[:end])

# let's remove some of them from extractor
delete!(extractor.dict, :paper_id)
delete!(extractor.dict[:metadata].dict[:authors].item.dict, :last)
delete!(extractor.dict[:metadata].dict[:authors].item.dict, :middle)

# we can also notice, that some long texts are extracted as categorical variables, e.g.
extractor[:body_text].item[:text]
extractor[:body_text].item[:section]
# let's replace them manually by string extractors
# note that we need to use the .dict, as the [] accessor on item is just readonly syntax-sugar
extractor[:body_text].item.dict[:text] = ExtractString()
extractor[:body_text].item.dict[:section] = ExtractString()

# this concludes example about examining schema and modifying extractor accordingly.
using JsonGrinder: is_intable, is_floatable, unify_types, extractscalar
function string_multi_representation_scalar_extractor()
	vcat([
	(e -> unify_types(sch[:paper_id]) <: String,
		(e, uniontypes) -> MultipleRepresentation((
			ExtractCategorical(top_n_keys(e, 20), uniontypes),
			extractscalar(unify_types(e), e, uniontypes)
		)))
	], JsonGrinder.default_scalar_extractor())
end

top_n_keys(e::Entry, n::Int) = map(x->x[1], sort(e.counts |> collect, by=x->x[2], rev=true)[begin:min(n, end)])
suggestextractor(sch, (;
	scalar_extractors=string_multi_representation_scalar_extractor(),
	key_as_field=13,
	)
) |> printtree
unify_types(sch[:paper_id]) <: String


text:
using Documenter, GuidedProposals

makedocs(;
    modules=[GuidedProposals],
    format=Documenter.HTML(
        mathengine = Documenter.MathJax(
            Dict(
                :TeX => Dict(
                    :equationNumbers => Dict(
                        :autoNumber => "AMS"
                    ),
                    :Macros => Dict(
                        :dd => "{\\textrm d}",
                        :RR => "\\mathbb{R}",
                        :wt => ["\\widetilde{#1}", 1]
                    ),
                    :packages => Dict(
                        "[+]" => ["color"]
                    ),
                ),
                :loader => Dict(
                    :load => ["[tex]/color"],
                ),
            ),
        ),
        collapselevel = 1,
    ),
    pages=[
        "Home" => "index.md",
        "Get started" => joinpath("get_started", "overview.md"),
        "User manual" => Any[
            "Guided proposals" => joinpath("manual", "guid_prop.md"),
            "Multiple observations" => joinpath("manual", "multiple_obs.md"),
            "Log-likelihoods" => joinpath("manual", "log_likelihood.md"),
            "Reparameterizations" => joinpath("manual", "reparameterizations.md"),
            "BFFG algorithm" => joinpath("manual", "bffg.md"),
            "Computing path functionals" => joinpath("manual", "path_functionals.md"),
            "(TODO) ODE systems" => joinpath("manual", "ode_type.md"),
        ],
        "How to..." => Any[
            "(TODO) Sample diffusion bridges" => joinpath("how_to_guides", "sample_bridges.md"),
            "(TODO) Sample diffusions in first-passage time setting" => joinpath("how_to_guides", "first_passage_time.md"),
            "...smooth trajectories" => joinpath("how_to_guides", "smoothing.md"),
            "...do blocking" => joinpath("how_to_guides", "blocking.md"),
            "...infer parameters" => joinpath("how_to_guides", "parameter_inference.md")
        ],
        "Tutorials" => Any[
            "(TODO) Understanding in-place vs out-of-place" => joinpath("tutorials", "in_place_vs_out_of_place.md"),
            "Smoothing trajectories of the partially observed FitzHugh–Nagumo diffusion" => joinpath("tutorials", "smoothing.md"),
            "Setting up a blocking scheme for a simple smoothing problem" => joinpath("tutorials", "blocking.md"),
            "Parameter inference" => joinpath("tutorials", "parameter_inference.md"),
        ],
        "Index" => "module_index.md",
    ],
    repo="https://github.com/JuliaDiffusionBayes/GuidedProposals.jl/blob/{commit}{path}#L{line}",
    sitename="GuidedProposals.jl",
    authors="Sebastiano Grazzi, Frank van der Meulen, Marcin Mider, Moritz Schauer",
    assets=String[],
)

deploydocs(;
    repo="github.com/JuliaDiffusionBayes/GuidedProposals.jl",
)


text:
# Full Wrapper over NIVISA, api and constants
# don't make visa a module since it needs to export all constants and error messages into other instrument modules
include("ctypes.jl")
include("visa_common.jl")
include("visa_errors.jl")
include("visa_api.jl")
include("visa_utils.jl")
include("Ivi_common.jl")


text:
export GNhis, getGNhis, updateHis!

"""
mutable struct jInv.InverseSolve.GNhis

stored iteration history of Gauss-Newton methods and variants

Fields:

Jc             - objective function values
dJ             - norm of (projected) gradients
F              - misfits
Dc             - data
Rc             - value of regularizer
alphas         - regularization parameters
Active         - active sets
stepNorm       - norm of step
lsIter         - nuber of line search iterations
timeMisfit     - time to evaluate misfits
timeReg        - time to evaluate regularizer
hisLinSol      - history of linear solver (if applicable)
timeLinSol     - time for linear solves
timeGradMisfit - time for gradient computation

"""
mutable struct GNhis
	Jc::Array
	dJ::Array
	F::Array
	Dc::Array
	Rc::Array
	alphas::Array
	Active::Array
	stepNorm::Array
	lsIter::Array
	timeMisfit::Array
	timeReg::Array
	hisLinSol::Array
	timeLinSol::Array
	timeGradMisfit::Array
end

"""
function jInv.InverseSolve.getGNhis(maxIter,maxIterCG)

constructs GNhis

Input:

maxIter
maxIterCG

"""
function getGNhis(maxIter,maxIterCG)
	Jc = zeros(maxIter+1)
	dJ = zeros(maxIter+1)
	F  = zeros(maxIter+1)
	Dc = []
	Rc = zeros(maxIter+1)
	alphas = zeros(maxIter+1)
	Active = zeros(maxIter+1)
	stepNorm = zeros(maxIter+1)
	lsIter = zeros(Int,maxIter+1)
	timeMisfit = zeros(maxIter+1,4)
	timeReg = zeros(maxIter+1)
	hisLinSol = []
	timeLinSol = zeros(maxIter+1,1)
	timeGradMisfit = zeros(maxIter+1,2)

	return GNhis(Jc,dJ,F,Dc,Rc,alphas,Active,stepNorm,lsIter,timeMisfit,timeReg,hisLinSol,timeLinSol,timeGradMisfit)
end

function updateHis!(iter::Int64,His::GNhis,Jc::Real,dJ::Real,Fc,Dc,Rc::Real,alpha::Real,nActive::Int64,stepNorm::Real,lsIter::Int,timeMisfit::Vector,timeReg::Real)
	His.Jc[iter+1]            = Jc
	His.dJ[iter+1]            = dJ
	His.F[iter+1]             = Fc
	push!(His.Dc,Dc)
	His.Rc[iter+1]            = Rc
	His.alphas[iter+1]        = alpha
	His.Active[iter+1]        = nActive
	His.stepNorm[iter+1]      = stepNorm
	His.lsIter[iter+1]        = lsIter
	His.timeMisfit[iter+1,:] += timeMisfit
	His.timeReg[iter+1]      += timeReg[]
end