text:
#=~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~=#
# Problem set 7 solutions
# Written by Tyler Ransom
# Commented by Giuseppe Grasso
#=~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~=#

using Random
using LinearAlgebra
using Statistics
using Distributions
using Optim
using DataFrames
using CSV
using HTTP
using GLM
using FreqTables
using ForwardDiff
using LineSearches
using SMM
cd("/Users/peppegrass/Documents/GitHub/fall-2020/ProblemSets/PS6-ccp/") # GG: sets working directory
pwd() ## GG: prints working directory
readdir() # GG: equivalent to -ls- to see elements of working directory

function allwrap()

    #=
    GG: In this problem set, we will practice estimating models by
    Generalized Method of Moments (GMM) and Simulated Method of Moments (SMM).
    =#

    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    # Question 1 - GMM for simple regression model
    # data & model are SAME AS PS2, QUESTION 2
    #
    # 1. Estimate the linear regression model from Question 2 of Problem Set 2 by GMM.
    # Write down the moment function as in slide #8 of the Lecture 9 slide deck and use Optim for estimation.
    # Use the N × N Identity matrix as your weighting matrix. Check your answer using the closed-form matrix formula for the OLS estimator.
    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS1-julia-intro/nlsw88.csv"
    df = CSV.read(HTTP.get(url).body)
    X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
    y = log.(df.wage)

    function ols_gmm(β, X, y) # GG: GMM function
        K = size(X,2) # GG: don't really needed for GMM, unless for computing RMSE in step below
        N = size(y,1) # GG: don't really needed for GMM, unless for computing RMSE in step below
        ŷ = X*β
        g = y .- ŷ # GG: force y and y_hat to be as close as possible
        J = g'*I*g # GG: GMM objective function
        return J
    end

    # we can also estimate the RMSE of the regression by ADDING ANOTHER MOMENT (and assuming homoskedasticity)
    function ols_gmm_with_σ(θ, X, y)
        K = size(X,2)
        N = size(y,1)
        β = θ[1:end-1]
        σ = θ[end]
        g = y .- X*β
        g = vcat(g,( (N-1)/(N-K) )*var(g) .- σ^2) # this builds on the assumption that MSE: SSR/(N-K) = σ² # GG: performing DOFs corrections
        J = g'*I*g
        return J
    end

    β_hat_gmm = optimize(b -> ols_gmm(b, X, y), rand(size(X,2)), LBFGS(), Optim.Options(g_tol=1e-6, iterations=100_000, show_trace=true))
    println(β_hat_gmm.minimizer)
    println("ans: ",X\y) # GG: X\y gives the least squares; cool uh!?
    β_hat_gmm = optimize(b -> ols_gmm_with_σ(b, X, y), rand(size(X,2)+1), LBFGS(), Optim.Options(g_tol=1e-6, iterations=100_000, show_trace=true))
    println(β_hat_gmm.minimizer)
    println("ans: ",vcat(X\y,sqrt( sum( (y .- X*β_hat_gmm.minimizer[1:end-1]).^2 )/( size(y,1)-size(X,2) ) )))

    # GG: of course ols_gmm(b, X, y) and ols_gmm_with_σ(b, X, y) yield same results but latter with MSE

    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    # Question 2 - GMM for multinomial logit model
    # data & model are SAME AS PS2, QUESTION 5
    #
    # 2. Estimate the multinomial logit model from Question 5 of Problem Set 2 by the following means:
    #   (a) Maximum likelihood (i.e. re-run your code [or mine] from Question 5 of Problem Set 2)
    #   (b) GMM with the MLE estimates as starting values. Your g object should be a vector of dimension N ×J where
    #   N is the number of rows of the X matrix and J is the dimension of the choice set. Each element, g should equal d − P,
    #   where d and P are “stacked” vectors of dimension N×J
    #   (c) GMM with random starting values
    # Compare your estimates from part (b) and (c). Is the objective function globally concave?
    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    df = dropmissing(df, :occupation)
    df[df.occupation.==8 ,:occupation] .= 7
    df[df.occupation.==9 ,:occupation] .= 7
    df[df.occupation.==10,:occupation] .= 7
    df[df.occupation.==11,:occupation] .= 7
    df[df.occupation.==12,:occupation] .= 7
    df[df.occupation.==13,:occupation] .= 7

    #df = @transform(df, white = race.==1)
    df[:,:white] = df.race .== 1
    X = [ones(size(df,1),1) df.age df.white df.collgrad]
    y = df.occupation

    # get starting values for mlogit from series of binary logits
    svals = zeros(size(X,2),7)
    for j=1:7
        tempname = Symbol(string("occ",j))
        df[:,tempname] = df.occupation.==j
    end
    # GG: this is to derive some "reasonable starting values" for a multinomial logit model
    # GG: what do we mean for "reasonable starting values". You don't wanna have neither a vector of 0 nor a vector of random numbers
    # GG: So what starting values for a Multinomial logit. Estimate a bunch of binary logits and then...
    svals[:,1] = coef(lm(@formula(occ1 ~ age + white + collgrad), df))
    svals[:,2] = coef(lm(@formula(occ2 ~ age + white + collgrad), df))
    svals[:,3] = coef(lm(@formula(occ3 ~ age + white + collgrad), df))
    svals[:,4] = coef(lm(@formula(occ4 ~ age + white + collgrad), df))
    svals[:,5] = coef(lm(@formula(occ5 ~ age + white + collgrad), df))
    svals[:,6] = coef(lm(@formula(occ6 ~ age + white + collgrad), df))
    svals[:,7] = coef(lm(@formula(occ7 ~ age + white + collgrad), df))
    svals = svals[:,1:6] .- svals[:,7] # GG: ... impose the identification normalization (i.e. you subtract for the first 6 estimates the one of the 7th)
    svals = svals[:] # GG: then vectorize it and use it as a starting value

    function mlogit_mle(α, X, y)

        K = size(X,2)
        J = length(unique(y))
        N = length(y)
        bigY = zeros(N,J)
        for j=1:J
            bigY[:,j] = y.==j
        end
        bigα = [reshape(α,K,J-1) zeros(K)]

        P = exp.(X*bigα) ./ sum.(eachrow(exp.(X*bigα)))

        loglike = -sum( bigY.*log.(P) ) # GG: MLE objective function

        return loglike
    end

    α_zero = zeros(6*size(X,2))
    α_rand = rand(6*size(X,2))
    # α_true = [-.5120867,-.6950393,.3165616,-1.646016,-2.300697,-.7656569]
    α_true = [.1910213,-.0335262,.5963968,.4165052,-.1698368,-.0359784,1.30684,-.430997,.6894727,-.0104578,.5231634,-1.492475,-2.26748,-.0053001,1.391402,-.9849661,-1.398468,-.0142969,-.0176531,-1.495123,.2454891,-.0067267,-.5382892,-3.78975]
    α_start = α_true.*rand(size(α_true))
    println(size(α_true))
    α_hat_optim = optimize(a -> mlogit_mle(a, X, y), svals, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
    α_hat_mle = α_hat_optim.minimizer
    println("mlogit MLE estimates: ",α_hat_mle)

    function mlogit_gmm(α, X, y)
        K = size(X,2)
        J = length(unique(y))
        N = length(y)
        bigY = zeros(N,J)
        for j=1:J
            bigY[:,j] = y.==j
        end
        bigα = [reshape(α,K,J-1) zeros(K)]

        P = exp.(X*bigα) ./ sum.(eachrow(exp.(X*bigα)))

        g = zeros((J-1)*K)
        for j=1:(J-1)
            for k = 1:K
                g[(j-1)*K+k] = mean( (bigY[:,j] .- P[:,j]) .* X[:,k] ) # GG: GMM (exactly-identified) objective function
                # GG: (j-1)*K moment conditions; (j-1) equations and K betas; same number of elements as number non-zero elements in bigαlpha
            end
        end
        J = N*g⋅g #GG: exactly identified; no weighting matrix
        return J
    end

    function mlogit_gmm_overid(α, X, y)
        K = size(X,2)
        J = length(unique(y))
        N = length(y)
        bigY = zeros(N,J)
        for j=1:J
            bigY[:,j] = y.==j
        end
        bigα = [reshape(α,K,J-1) zeros(K)]

        P = exp.(X*bigα) ./ sum.(eachrow(exp.(X*bigα))) # GG: same as before P = exp.(X*bigα) ./ sum.(eachrow(exp.(X*bigα)))

        g = bigY[:] .- P[:] # GG: now in over-identified case; we vectorize big Y and P with the [:] operator and subtract the two
        # GG: this will be equivalent to NONLINEAR LEAST SQUARES

        J = g'*I*g # GG: now we have weighting matrix
        return J
    end
    #α_zero = zeros(6*size(X,2))
    #α_rand = rand(6*size(X,2))
    α_true = [.1910213,-.0335262,.5963968,.4165052,-.1698368,-.0359784,1.30684,-.430997,.6894727,-.0104578,.5231634,-1.492475,-2.26748,-.0053001,1.391402,-.9849661,-1.398468,-.0142969,-.0176531,-1.495123,.2454891,-.0067267,-.5382892,-3.78975]
    #α_start = α_true.*rand(size(α_true))
    #println("starting values: ",α_rand)

    # GG: here we're playing around with different optimization algorithms, starting values, and later using autodiff

    # LBFGS
    α_hat_optim = optimize(a -> mlogit_gmm_overid(a, X, y), α_true .+ .0001*rand(size(α_true)), LBFGS(), Optim.Options(g_tol = 1e-5, iterations=1_000, show_trace=true, show_every=100))
    α_hat_gmm = α_hat_optim.minimizer
    println("LBFGS GMM minimum: ",α_hat_optim.minimum)
    println("LBFGS GMM estimates: ",α_hat_gmm)

    # use autodiff with LBFGS
    td = TwiceDifferentiable(b -> mlogit_gmm_overid(b, X, y), svals; autodiff = :forward)
    α_hat_optim_ad = optimize(td, α_true .+ .0001*rand(size(α_true)), LBFGS(), Optim.Options(g_tol = 1e-5, iterations=10_000, show_trace=true, show_every=100))
    α_hat_gmm_ad = α_hat_optim_ad.minimizer
    println("LBFGS GMM autodiff minimum: ",α_hat_optim_ad.minimum)
    println("LBFGS GMM autodiff estimates: ",α_hat_gmm_ad)

    # use autodiff with Newton
    α_hat_optim_ad_newton = optimize(td, α_true .+ .0001*rand(size(α_true)), Newton(linesearch = BackTracking()), Optim.Options(g_tol = 1e-5, iterations=10_000, show_trace=true, show_every=100))
    α_hat_gmm_ad_newton = α_hat_optim_ad_newton.minimizer
    println("Newton GMM autodiff minimum: ",α_hat_optim_ad_newton.minimum)
    println("Newton GMM autodiff estimates: ",α_hat_gmm_ad_newton)

    # autodiff from bad starting values
    α_hat_optim_ad_fr_bad = optimize(td, α_rand.*sign.(α_true), LBFGS(), Optim.Options(g_tol = 1e-10, f_tol = 1e-10, x_tol = 1e-10, iterations=10_000, show_trace=true, show_every=50))
    α_hat_gmm_ad_fr_bad = α_hat_optim_ad_fr_bad.minimizer
    println("LBFGS GMM autodiff minimum: ",α_hat_optim_ad_fr_bad.minimum)
    println("LBFGS GMM autodiff estimates: ",α_hat_gmm_ad_fr_bad)

    # GG: putting different answers in a DataFrame and comparing them
    # GG: found out close answers; the one using bad starting value didn't get close to other solutions; some times GMM tolerance values needed to be played around with
    compare = DataFrame(mle = α_hat_mle, gmm = α_hat_gmm, gmm_ad = α_hat_gmm_ad, gmm_ad_newton = α_hat_gmm_ad_newton, gmm_ad_fr_bad = α_hat_gmm_ad_fr_bad)
    println(compare)
    # GG: objective function globally concave, but couldn't get to global maximum from random starting values, although pretty close

    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    # Question 3 - SIMULATE DATA from a logit model
    #
    # 3. SIMULATE A DATASET FROM A MULTINOMIAL LOGIT MODEL, and then estimate its parameter values and verify that the estimates are close to the
    # parameter values you set. That is, for a given sample size N, choice set dimension J and parameter vector β, write a function that outputs
    # data X and Y. I will let you choose N, J, β and the number of covariates in X (K), but J should be larger than 2 and K should be larger than 1.
    # If you haven’t done this before, you may want to follow these steps:
    #   (a) Generate X using a random number generator—rand() or randn().
    #   (b) Set values for β such that conformability with X and J is satisﬁed
    #   (c) Generate the N×J matrix of choice probabilities P
    #   (d) Draw the preference shocks ε as a N×1 vector of U[0,1] random numbers
    #   (e) Generate Y as follows:
    #     • Initialize Y as an N ×1 vector of 0s
    #     • Update Y i = ∑ j=1 J 1 P ik > εi  ∑k=j J  [see text of problem set]
    #   (f) An alternative way to generate choices would be to draw a N × J matrix of ε’s from a T1EV distribution.
    #   This distribution is already deﬁned in the Distributions package. Then Y i = argmax j X i β j + ε i j .
    #   I’ll show you an example of how to do that in the solutions code for this problem set.
    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    function sim_logit(N=100_000,J=4) # putting "=" after the arguments here sets default values for these arguments
        # generate X, β and P
        X = hcat(ones(N),randn(N),2 .+ 2 .* randn(N))
        if J==4
            β = hcat([1, -1, 0.5],[-2, 0.5, 0.3],[0, -0.5, 2],zeros(3)) # GG: numbers that square nicely with 4
        else
            β = 2 .* rand(size(X,2),J) .- 1
        end
        P = exp.(X*β) ./ sum.(eachrow(exp.(X*β))) # GG: usual one-liner to generate multinomial ogit choice probabilities; eachrow is summation
        # draw choices
        draw = rand(N)
        Y = zeros(N) # Initializing Y as an Nx1 vector of 0s
        for j=1:J
            Ytemp = sum(P[:,j:J];dims=2) .> draw # GG: then updating it (summing across j) based on whether choice probability is above or below draw [see formula of problem set]
            Y += Ytemp # GG: for each iteration of the loop adding Ytemp to Y
        end
        return Y,X
    end
    ySim,XSim = sim_logit() # GG: it will use "default values" for N and J set in line 244 through s"et to the default operator =""
    α_hat_sim_optim = optimize(a -> mlogit_mle(a, XSim, ySim), rand(9), LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
    α_hat_sim_mle = α_hat_sim_optim.minimizer
    println("simulated mlogit MLE estimates: ",α_hat_sim_mle) # GG: and here are our simulated estimates; values pretty close because sample size big and SEs low

    # now with drawing from Gumbel instead (#GG: i.e T1EV)
    # GG: python does not have the nice Distributions package as in Julia
    function sim_logit_w_gumbel(N=100_000,J=4) # putting = here sets default values for these arguments
        # generate X and β (No P needed)
        X = hcat(ones(N),randn(N),2 .+ 2 .* randn(N))
        if J==4
            β = hcat([1, -1, 0.5],[-2, 0.5, 0.3],[0, -0.5, 2],zeros(3))
        else
            β = 2 .* rand(size(X,2),J) .- 1
        end
        # draw choices
        ϵ = rand(Gumbel(0,1),N,J) # GG: using Julia's Distributions package; Gumbel(0,1) is "standard" T1EV
        Y = argmax.(eachrow(X*β .+ ϵ)) # GG: getting argmax of eachrow; eachrow allows avoiding to have to write a loop; tells Julia to do argmax computation row-wise and find column j maximizer
        # GG: that's just another way of deriving the Y; the loop in the previous function is just as fine
        return Y,X
    end
    ySim,XSim = sim_logit_w_gumbel()
    α_hat_sim_optim = optimize(a -> mlogit_mle(a, XSim, ySim), rand(9), LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
    α_hat_sim_mle = α_hat_sim_optim.minimizer
    println("simulated mlogit MLE estimates: ",α_hat_sim_mle)

    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    # Question 4 - run SMM.parallelNormal()
    #
    # 4. Use SMM.jl to run the example code on slide #21 of the Lecture 9 slide deck.
    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    # GG: SMM.jl has a nice count-down that Optim does not have
    MA = SMM.parallelNormal()
    dc = SMM.history(MA.chains[1])
    dc = dc[dc[:accepted].==true, :]
    println(describe(dc)) # GG: can verify you get the mean of the MVN distribution

    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    # Question 5 - estimate mlogit by SMM
    # GG: Wanna use SIMULATED METHOD OF MOMENTS (using code of Q3) to ESTIMATE RESULT OF Q2
    #
    # 5. Use your code from Question 3 to estimate the multinomial logit model from Question 2 using SMM and the code example from slide #18 of
    # the Lecture 9 slide deck.
    #:::::::::::::::::::::::::::::::::::::::::::::::::::
    # GG: X and y of Q2
    X = [ones(size(df,1),1) df.age df.white df.collgrad]
    y = df.occupation
    # GG: chunk below was commented out by Tyler beacuse the set of moment conditions used were just not working out for him
    #function mlogit_smm(α, X, y, D)
    #    K = size(X,2)
    #    J = length(unique(y))
    #    N = length(y)
    #    bigY = zeros(N,J)
    #    bigỸ = zeros(N,J)
    #    for j=1:J
    #        bigY[:,j] = y.==j
    #    end
    #    bigα = [reshape(α,K,J-1) zeros(K)]

    #    T = promote_type(eltype(α),eltype(X))
    #    gmodel = zeros(T,(J-1)*K,D)
    #    gdata  = zeros(T,(J-1)*K)

    #    # data moments
    #    for j=1:(J-1)
    #        gdata[(j-1)*K+1] = mean( (bigY[:,j]) )
    #        for k = 2:K
    #            gdata[(j-1)*K+k] = cov( bigY[:,j], X[:,k] )
    #        end
    #    end

    #    Random.seed!(1234)
    #    # simulated model moments
    #    for d=1:D
    #        # draw choices
    #        ε = rand(Gumbel(0,1),N,J)
    #        ỹ = argmax.(eachrow(X*bigα .+ ε))
    #        for j=1:J
    #            bigỸ[:,j] = ỹ.==j
    #        end
    #        # data moments
    #        for j=1:(J-1)
    #            gmodel[(j-1)*K+1,d] = mean( (bigỸ[:,j]) )
    #            for k = 2:K
    #                gmodel[(j-1)*K+k,d] = cov( bigỸ[:,j], X[:,k] )
    #            end
    #        end
    #    end

    #    # minimize squared difference between data and moments
    #    g = mean(gmodel .- gdata; dims=2)
    #    J = N*g⋅g
    #    return J
    #

    # GG: the function below for the SIMULATED METHOD OF MOMENTS is the one from the SLIDES OF LECTURE 9 (mind the SCROLLER; gotta use HTML and not PDF version)
    function mlogit_smm_overid(α, X, y, D) # GG: D is the number of draws
        K = size(X,2)
        J = length(unique(y))
        N = length(y)
        bigY = zeros(N,J)
        bigỸ = zeros(N,J)
        for j=1:J
            bigY[:,j] = y.==j # GG: bigY is our real data that we're trying to MATCH
        end
        bigα = [reshape(α,K,J-1) zeros(K)] # GG: bigα is our parameter vector that we're gonna try to guess value of, so that we get values so that simulated data match actual ones as closely as possible

        Random.seed!(1234) # GG: cardinal rule, set up the seed; don't forget that!
        for d=1:D
            # draw choices
            ε = rand(Gumbel(0,1),N,J) # GG: if we change the distribution of ε to a multivariate normal, we can do MNP without having to do any Monte Carlo integration. We don't need to have a closed form solution for this.
            # GG: and if we want to do GEV models, we need to specify nest-specific and alternative-specific component; just two errors but same principle.
            ỹ = argmax.(eachrow(X*bigα .+ ε)) # GG: we're goona draw ỹ through ε
            for j=1:J
                bigỸ[:,j] .+= (ỹ.==j)*(1/D) # GG: bigỹ is our SIMULATED data
                # GG: we set the jth column of bigỸ to a dummy if ỹ is equal to j
            end
        end

        g = bigY[:] .- bigỸ[:] # GG: exactly like g in Question 2; very same principle

        J = g'*I*g
        return J
    end
    println(size(X,2))
    println(length(unique(y))-1)
    println("size of svals: ",size(α_true .+ .0001*rand(size(α_true))))
    bigα = [reshape(α_true,size(X,2),length(unique(y))-1) zeros(size(X,2))]
    println("size of bigα: ",size(bigα))
    td = TwiceDifferentiable(th -> mlogit_smm_overid(th, X, y, 2_000), α_true .+ .0001*rand(size(α_true)); autodiff = :forward)
    α_hat_smm_overid = optimize(td, α_true .+ .0001*rand(size(α_true)), LBFGS(), Optim.Options(g_tol=1e-8, x_tol=1e-8, f_tol=1e-8, iterations=100_000, show_trace=true, show_every=5))
    println(α_hat_smm_overid.minimizer)
    return nothing
end

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 6
#:::::::::::::::::::::::::::::::::::::::::::::::::::
@time allwrap()


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
function add_storage_size_constraints(m, p, b; _n="")
    # TODO add formal types for storage (i.e. "b")

	# Constraint (4b)-1: Lower bound on Storage Energy Capacity
	@constraint(m,
        m[Symbol("dvStorageEnergy"*_n)][b] >= p.s.storage.attr[b].min_kwh
    )

	# Constraint (4b)-2: Upper bound on Storage Energy Capacity
	@constraint(m,
        m[Symbol("dvStorageEnergy"*_n)][b] <= p.s.storage.attr[b].max_kwh
    )

	# Constraint (4c)-1: Lower bound on Storage Power Capacity
	@constraint(m,
        m[Symbol("dvStoragePower"*_n)][b] >= p.s.storage.attr[b].min_kw
    )

	# Constraint (4c)-2: Upper bound on Storage Power Capacity
	@constraint(m,
        m[Symbol("dvStoragePower"*_n)][b] <= p.s.storage.attr[b].max_kw
    )
end

function add_general_storage_dispatch_constraints(m, p, b; _n="")
    # Constraint (4a): initial state of charge
	@constraint(m,
        m[Symbol("dvStoredEnergy"*_n)][b, 0] == p.s.storage.attr[b].soc_init_pct * m[Symbol("dvStorageEnergy"*_n)][b]
    )

    #Constraint (4n): State of charge upper bound is storage system size
    @constraint(m, [ts in p.time_steps],
        m[Symbol("dvStoredEnergy"*_n)][b,ts] <= m[Symbol("dvStorageEnergy"*_n)][b]
    )

    # Constraint (4j): Minimum state of charge
	@constraint(m, [ts in p.time_steps],
        m[Symbol("dvStoredEnergy"*_n)][b, ts] >= p.s.storage.attr[b].soc_min_pct * m[Symbol("dvStorageEnergy"*_n)][b]
    )

    #Constraint (4j): Dispatch from storage is no greater than power capacity
	@constraint(m, [ts in p.time_steps],
        m[Symbol("dvStoragePower"*_n)][b] >= m[Symbol("dvDischargeFromStorage"*_n)][b, ts]
    )

end

function add_elec_storage_dispatch_constraints(m, p, b; _n="")

	# Constraint (4g): state-of-charge for electrical storage - with grid
	@constraint(m, [ts in p.time_steps_with_grid],
        m[Symbol("dvStoredEnergy"*_n)][b, ts] == m[Symbol("dvStoredEnergy"*_n)][b, ts-1] + p.hours_per_time_step * (
            sum(p.s.storage.attr[b].charge_efficiency * m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for t in p.techs.elec)
            + p.s.storage.attr[b].grid_charge_efficiency * m[Symbol("dvGridToStorage"*_n)][b, ts]
            - m[Symbol("dvDischargeFromStorage"*_n)][b,ts] / p.s.storage.attr[b].discharge_efficiency
        )
	)

	# Constraint (4h): state-of-charge for electrical storage - no grid
	@constraint(m, [ts in p.time_steps_without_grid],
        m[Symbol("dvStoredEnergy"*_n)][b, ts] == m[Symbol("dvStoredEnergy"*_n)][b, ts-1] + p.hours_per_time_step * (
            sum(p.s.storage.attr[b].charge_efficiency * m[Symbol("dvProductionToStorage"*_n)][b,t,ts] for t in p.techs.elec)
            - m[Symbol("dvDischargeFromStorage"*_n)][b, ts] / p.s.storage.attr[b].discharge_efficiency
        )
    )

	# Constraint (4i)-1: Dispatch to electrical storage is no greater than power capacity
	@constraint(m, [ts in p.time_steps],
        m[Symbol("dvStoragePower"*_n)][b] >=
            sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for t in p.techs.elec) + m[Symbol("dvGridToStorage"*_n)][b, ts]
    )

	#Constraint (4k)-alt: Dispatch to and from electrical storage is no greater than power capacity
	@constraint(m, [ts in p.time_steps_with_grid],
        m[Symbol("dvStoragePower"*_n)][b] >= m[Symbol("dvDischargeFromStorage"*_n)][b, ts] +
            sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for t in p.techs.elec) + m[Symbol("dvGridToStorage"*_n)][b, ts]
    )

	#Constraint (4l)-alt: Dispatch from electrical storage is no greater than power capacity (no grid connection)
	@constraint(m, [ts in p.time_steps_without_grid],
        m[Symbol("dvStoragePower"*_n)][b] >= m[Symbol("dvDischargeFromStorage"*_n)][b,ts] +
            sum(m[Symbol("dvProductionToStorage"*_n)][b, t, ts] for t in p.techs.elec)
    )

    # Remove grid-to-storage as an option if option to grid charge is turned off
    if !(p.s.storage.attr[b].can_grid_charge)
        for ts in p.time_steps_with_grid
            fix(m[Symbol("dvGridToStorage"*_n)][b, ts], 0.0, force=true)
        end
	end

    if p.s.storage.attr[b].minimum_avg_soc_fraction > 0
        avg_soc = sum(m[Symbol("dvStoredEnergy"*_n)][b, ts] for ts in p.time_steps) /
                   (8760. / p.hours_per_time_step)
        @constraint(m, avg_soc >= p.s.storage.attr[b].minimum_avg_soc_fraction *
            sum(m[Symbol("dvStorageEnergy"*_n)][b])
        )
    end
end

function add_hot_thermal_storage_dispatch_constraints(m, p, b; _n="")

    # # Constraint (4f)-1: (Hot) Thermal production sent to storage or grid must be less than technology's rated production
	# # Constraint (4f)-1a: BoilerTechs
	if !isempty(p.techs.boiler)
	# 	if !isempty(p.SteamTurbineTechs)
    #         @constraint(m, BoilerTechProductionFlowCon[b in p.s.storage.types.hot, t in p.techs.boiler, ts in p.time_steps],
    #                 m[Symbol("dvProductionToStorage"*_n)][b,t,ts] + m[Symbol("dvThermalToSteamTurbine"*_n)][t,ts]  <=
    #                 m[Symbol("dvThermalProduction"*_n)][t,ts]
    #                 )
    #     else
            @constraint(m, BoilerTechProductionFlowCon[b in p.s.storage.types.hot, t in p.techs.boiler, ts in p.time_steps],
                    m[Symbol("dvProductionToStorage"*_n)][b,t,ts]  <=
                    m[Symbol("dvThermalProduction"*_n)][t,ts]
                    )
    #     end
    end

    # # Constraint (4f)-1b: SteamTurbineTechs
	# if !isempty(p.SteamTurbineTechs)
	# 	@constraint(m, SteamTurbineTechProductionFlowCon[b in p.s.storage.types.hot, t in p.SteamTurbineTechs, ts in p.time_steps],
	# 		m[Symbol("dvProductionToStorage"*_n)][b,t,ts] <=
	# 		p.production_factor[t,ts] * m[Symbol("dvThermalProduction"*_n)][t,ts]
	# 		)
	# end

    # # Constraint (4g): CHP Thermal production sent to storage or grid must be less than technology's rated production
	if !isempty(p.techs.chp)
	# 	if !isempty(p.SteamTurbineTechs)
    #         @constraint(m, CHPTechProductionFlowCon[b in p.s.storage.types.hot, t in p.techs.chp, ts in p.time_steps],
    #                 m[Symbol("dvProductionToStorage"*_n)][b,t,ts] + m[Symbol("dvProductionToWaste"*_n)][t,ts] + m[Symbol("dvThermalToSteamTurbine"*_n)][t,ts] <=
    #                 m[Symbol("dvThermalProduction"*_n)][t,ts]
    #                 )
    #     else
            @constraint(m, CHPTechProductionFlowCon[b in p.s.storage.types.hot, t in p.techs.chp, ts in p.time_steps],
                    m[Symbol("dvProductionToStorage"*_n)][b,t,ts] + m[Symbol("dvProductionToWaste"*_n)][t,ts] <=
                    m[Symbol("dvThermalProduction"*_n)][t,ts]
                    )
    #     end
	end

    # Constraint (4j)-1: Reconcile state-of-charge for (hot) thermal storage
	@constraint(m, [b in p.s.storage.types.hot, ts in p.time_steps],
    m[Symbol("dvStoredEnergy"*_n)][b,ts] == m[Symbol("dvStoredEnergy"*_n)][b,ts-1] + (1/p.s.settings.time_steps_per_hour) * (
        sum( p.s.storage.attr[b].charge_efficiency * m[Symbol("dvProductionToStorage"*_n)][b,t,ts] for t in p.techs.heating) -
        m[Symbol("dvDischargeFromStorage"*_n)][b,ts] / p.s.storage.attr[b].discharge_efficiency -
        p.s.storage.attr[b].thermal_decay_rate_fraction * m[Symbol("dvStorageEnergy"*_n)][b]
        )
    )

    #Constraint (4n)-1: Dispatch to and from thermal storage is no greater than power capacity
	@constraint(m, [b in p.s.storage.types.hot, ts in p.time_steps],
        m[Symbol("dvStoragePower"*_n)][b] >=
        m[Symbol("dvDischargeFromStorage"*_n)][b,ts] +
        sum(m[Symbol("dvProductionToStorage"*_n)][b,t,ts] for t in p.techs.heating)
    )
    # TODO missing thermal storage constraints from API ???

end

function add_cold_thermal_storage_dispatch_constraints(m, p, b; _n="")

    # Constraint (4f)-2: (Cold) Thermal production sent to storage or grid must be less than technology's rated production
	if !isempty(p.techs.cooling)
		@constraint(m, CoolingTechProductionFlowCon[b in p.s.storage.types.cold, t in p.techs.cooling, ts in p.time_steps],
    	        m[Symbol("dvProductionToStorage"*_n)][b,t,ts]  <=
				m[Symbol("dvThermalProduction"*_n)][t,ts]
				)
	end

    # Constraint (4j)-2: Reconcile state-of-charge for (cold) thermal storage
	@constraint(m, ColdTESInventoryCon[b in p.s.storage.types.cold, ts in p.time_steps],
    m[Symbol("dvStoredEnergy"*_n)][b,ts] == m[Symbol("dvStoredEnergy"*_n)][b,ts-1] + (1/p.s.settings.time_steps_per_hour) * (
        sum(p.s.storage.attr[b].charge_efficiency * m[Symbol("dvProductionToStorage"*_n)][b,t,ts] for t in p.techs.cooling) -
        m[Symbol("dvDischargeFromStorage"*_n)][b,ts]/p.s.storage.attr[b].discharge_efficiency -
        p.s.storage.attr[b].thermal_decay_rate_fraction * m[Symbol("dvStorageEnergy"*_n)][b]
        )
    )

    #Constraint (4n)-2: Dispatch to and from thermal storage is no greater than power capacity
    @constraint(m, [b in p.s.storage.types.cold, ts in p.time_steps],
        m[Symbol("dvStoragePower"*_n)][b] >= m[Symbol("dvDischargeFromStorage"*_n)][b,ts] +
        sum(m[Symbol("dvProductionToStorage"*_n)][b,t,ts] for t in p.techs.cooling)
    )
end

function add_storage_sum_constraints(m, p; _n="")

	##Constraint (8c): Grid-to-storage no greater than grid purchases
	@constraint(m, [ts in p.time_steps_with_grid],
      sum(m[Symbol("dvGridPurchase"*_n)][ts, tier] for tier in 1:p.s.electric_tariff.n_energy_tiers) >=
      sum(m[Symbol("dvGridToStorage"*_n)][b, ts] for b in p.s.storage.types.elec)
    )
end

ext: jl

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
# https://pythonhosted.org/rasterstats/manual.html#raster-data-sources

println("Welcome... computing land use by region and year [1992-2019]...")

cd(@__DIR__)
using Pkg
Pkg.activate(".")
Pkg.instantiate()

#Pkg.add("PyCall")
#Pkg.add("DataFrames")
#Pkg.add("CSV")
#Pkg.add("Tables")
using PyCall, Tables, DataFrames, CSV

vectorBoundariesFile     = "./adminAreas/byLevels/gadm36_1.shp"
rasterDataFile_1992_2015 = "./luse/ESACCI-LC-L4-LCCS-Map-300m-P1Y-1992_2015-v2.0.7.tif"

sqKmPerCell          = (300*300)/(1000*1000)

const rs  = pyimport("rasterstats")
const gpd = pyimport("geopandas")

println("Modules loaded, reading the regions definition shapefile..")

regdf     = gpd.read_file(vectorBoundariesFile)
gid_0     = convert(Vector{String},regdf.GID_0)
name_0    = convert(Vector{String},regdf.NAME_0)
gid_1     = convert(Vector{String},regdf.GID_1)
name_1    = convert(Vector{String},regdf.NAME_1)
nRegions  = length(gid_0)

# see https://datastore.copernicus-climate.eu/documents/satellite-land-cover/D3.3.12-v1.3_PUGS_ICDR_LC_v2.1.x_PRODUCTS_v1.3.pdf
# Attention that the classification is hierarchical, e.g. "60" includes also "61" and "62"
luseIds = [ 10,  11,  12,  20,  30,  40,  50,  60,  61,  62,  70,  71,  72,
            80,  81,  82,  90, 100, 110, 120, 121, 122, 130, 140, 150, 151,
           152, 153, 160, 170, 180, 190, 200, 201, 202, 210, 220]
headersForExport = vcat("gid_0","name_0","gid_1","name_1",string.(luseIds)...)

years = 1992:2019
years = 2015:2019

for year in years
    println("Starting processing year $(year)...")
    if year <= 2015
        band = year-1992+1
        stats   = rs.zonal_stats(vectorBoundariesFile, rasterDataFile_1992_2015, band=band, categorical=true)  # stats=['count', 'mean', 'count','sum']
        stats   = convert(Vector{Dict{Int64,Int64}},stats)
        length(stats) == nRegions || @error "Oug, we have different number of regions in these computed stats!"
    else
        yearlyRasterFile = "./luse/C3S-LC-L4-LCCS-Map-300m-P1Y-$(year)-v2.1.1.tif"
        stats   = rs.zonal_stats(vectorBoundariesFile, yearlyRasterFile, band=1, categorical=true)  # stats=['count', 'mean', 'count','sum']
        stats   = convert(Vector{Dict{Int64,Int64}},stats)
        length(stats) == nRegions || @error "Oug, we have different number of regions in these computed stats!"
    end
    luseValues = [get(stat,luc,0) for stat in stats, luc in luseIds]
    yearlUseValue = hcat(gid_0,name_0,gid_1,name_1,luseValues)
    CSV.write("./out-computedLUseStatsByRegionAndYear/lUseByLev1Regions_$(year).csv",Tables.table(yearlUseValue), delim=';', header=headersForExport)
end

println("** Done!")


text:
# In this file, doctests which differ in the printed Float32 values won't fail
```@meta
DocTestFilters = r"[0-9\.]+f0"
```

"""
    mae(ŷ, y; agg = mean)

Return the loss corresponding to mean absolute error:

    agg(abs.(ŷ .- y))

# Example
```jldoctest
julia> y_model = [1.1, 1.9, 3.1];

julia> Flux.mae(y_model, 1:3)
0.10000000000000009
```
"""
function mae(ŷ, y; agg = mean)
  _check_sizes(ŷ, y)
  agg(abs.(ŷ .- y))
end

"""
    mse(ŷ, y; agg = mean)

Return the loss corresponding to mean square error:

    agg((ŷ .- y) .^ 2)

See also: [`mae`](@ref), [`msle`](@ref), [`crossentropy`](@ref).

# Example
```jldoctest
julia> y_model = [1.1, 1.9, 3.1];

julia> y_true = 1:3;

julia> Flux.mse(y_model, y_true)
0.010000000000000018
```
"""
function mse(ŷ, y; agg = mean)
  _check_sizes(ŷ, y)
  agg(abs2.(ŷ .- y))
end

"""
    msle(ŷ, y; agg = mean, ϵ = eps(ŷ))

The loss corresponding to mean squared logarithmic errors, calculated as

    agg((log.(ŷ .+ ϵ) .- log.(y .+ ϵ)) .^ 2)

The `ϵ` term provides numerical stability.
Penalizes an under-estimation more than an over-estimatation.

# Example
```jldoctest
julia> Flux.msle(Float32[1.1, 2.2, 3.3], 1:3)
0.009084041f0

julia> Flux.msle(Float32[0.9, 1.8, 2.7], 1:3)
0.011100831f0
```
"""
function msle(ŷ, y; agg = mean, ϵ = epseltype(ŷ))
  _check_sizes(ŷ, y)
  agg((log.((ŷ .+ ϵ) ./ (y .+ ϵ))) .^2 )
end

"""
    huber_loss(ŷ, y; δ = 1, agg = mean)

Return the mean of the [Huber loss](https://en.wikipedia.org/wiki/Huber_loss)
given the prediction `ŷ` and true values `y`.

                 | 0.5 * |ŷ - y|^2,            for |ŷ - y| <= δ
    Huber loss = |
                 |  δ * (|ŷ - y| - 0.5 * δ), otherwise
"""
function huber_loss(ŷ, y; agg = mean, δ = ofeltype(ŷ, 1))
   _check_sizes(ŷ, y)
   abs_error = abs.(ŷ .- y)
   #TODO: remove dropgrad when Zygote can handle this function with CuArrays
   temp = Zygote.dropgrad(abs_error .<  δ)
   x = ofeltype(ŷ, 0.5)
   agg(((abs_error .^ 2) .* temp) .* x .+ δ * (abs_error .- x * δ) .* (1 .- temp))
end

"""
    label_smoothing(y::Union{Number, AbstractArray}, α; dims::Int=1)

Returns smoothed labels, meaning the confidence on label values are relaxed.

When `y` is given as one-hot vector or batch of one-hot, its calculated as

    y .* (1 - α) .+ α / size(y, dims)

when `y` is given as a number or batch of numbers for binary classification,
its calculated as

    y .* (1 - α) .+ α / 2

in which case the labels are squeezed towards `0.5`.

α is a number in interval (0, 1) called the smoothing factor. Higher the
value of α larger the smoothing of `y`.

`dims` denotes the one-hot dimension, unless `dims=0` which denotes the application
of label smoothing to binary distributions encoded in a single number.

# Example
```jldoctest
julia> y = Flux.onehotbatch([1, 1, 1, 0, 1, 0], 0:1)
2×6 OneHotMatrix(::Vector{UInt32}) with eltype Bool:
 ⋅  ⋅  ⋅  1  ⋅  1
 1  1  1  ⋅  1  ⋅

julia> y_smoothed = Flux.label_smoothing(y, 0.2f0)
2×6 Matrix{Float32}:
 0.1  0.1  0.1  0.9  0.1  0.9
 0.9  0.9  0.9  0.1  0.9  0.1

julia> y_sim = softmax(y .* log(2f0))
2×6 Matrix{Float32}:
 0.333333  0.333333  0.333333  0.666667  0.333333  0.666667
 0.666667  0.666667  0.666667  0.333333  0.666667  0.333333

julia> y_dis = vcat(y_sim[2,:]', y_sim[1,:]')
2×6 Matrix{Float32}:
 0.666667  0.666667  0.666667  0.333333  0.666667  0.333333
 0.333333  0.333333  0.333333  0.666667  0.333333  0.666667

julia> Flux.crossentropy(y_sim, y) < Flux.crossentropy(y_sim, y_smoothed)
true

julia> Flux.crossentropy(y_dis, y) > Flux.crossentropy(y_dis, y_smoothed)
true
```
"""
function label_smoothing(y::Union{AbstractArray,Number}, α::Number; dims::Int = 1)
    if !(0 < α < 1)
        throw(ArgumentError("α must be between 0 and 1"))
    end
    if dims == 0
        y_smoothed = y .* (1 - α) .+ α*1//2
    elseif dims == 1
        y_smoothed = y .* (1 - α) .+ α* 1 // size(y, 1)
    else
        throw(ArgumentError("`dims` should be either 0 or 1"))
    end
    return y_smoothed
end

"""
    crossentropy(ŷ, y; dims = 1, ϵ = eps(ŷ), agg = mean)

Return the cross entropy between the given probability distributions;
calculated as

    agg(-sum(y .* log.(ŷ .+ ϵ); dims))

Cross entropy is typically used as a loss in multi-class classification,
in which case the labels `y` are given in a one-hot format.
`dims` specifies the dimension (or the dimensions) containing the class probabilities.
The prediction `ŷ` is supposed to sum to one across `dims`,
as would be the case with the output of a [`softmax`](@ref) operation.

For numerical stability, it is recommended to use [`logitcrossentropy`](@ref)
rather than `softmax` followed by `crossentropy` .

Use [`label_smoothing`](@ref) to smooth the true labels as preprocessing before
computing the loss.

See also: [`logitcrossentropy`](@ref), [`binarycrossentropy`](@ref), [`logitbinarycrossentropy`](@ref).

# Example
```jldoctest
julia> y_label = Flux.onehotbatch([0, 1, 2, 1, 0], 0:2)
3×5 OneHotMatrix(::Vector{UInt32}) with eltype Bool:
 1  ⋅  ⋅  ⋅  1
 ⋅  1  ⋅  1  ⋅
 ⋅  ⋅  1  ⋅  ⋅

julia> y_model = softmax(reshape(-7:7, 3, 5) .* 1f0)
3×5 Matrix{Float32}:
 0.0900306  0.0900306  0.0900306  0.0900306  0.0900306
 0.244728   0.244728   0.244728   0.244728   0.244728
 0.665241   0.665241   0.665241   0.665241   0.665241

julia> sum(y_model; dims=1)
1×5 Matrix{Float32}:
 1.0  1.0  1.0  1.0  1.0

julia> Flux.crossentropy(y_model, y_label)
1.6076053f0

julia> 5 * ans ≈ Flux.crossentropy(y_model, y_label; agg=sum)
true

julia> y_smooth = Flux.label_smoothing(y_label, 0.15f0)
3×5 Matrix{Float32}:
 0.9   0.05  0.05  0.05  0.9
 0.05  0.9   0.05  0.9   0.05
 0.05  0.05  0.9   0.05  0.05

julia> Flux.crossentropy(y_model, y_smooth)
1.5776052f0
```
"""
function crossentropy(ŷ, y; dims = 1, agg = mean, ϵ = epseltype(ŷ))
  _check_sizes(ŷ, y)
  agg(.-sum(xlogy.(y, ŷ .+ ϵ); dims = dims))
end

"""
    logitcrossentropy(ŷ, y; dims = 1, agg = mean)

Return the cross entropy calculated by

    agg(-sum(y .* logsoftmax(ŷ; dims); dims))

This is mathematically equivalent to `crossentropy(softmax(ŷ), y)`,
but is more numerically stable than using functions [`crossentropy`](@ref)
and [`softmax`](@ref) separately.

See also: [`binarycrossentropy`](@ref), [`logitbinarycrossentropy`](@ref), [`label_smoothing`](@ref).

# Example
```jldoctest
julia> y_label = Flux.onehotbatch(collect("abcabaa"), 'a':'c')
3×7 OneHotMatrix(::Vector{UInt32}) with eltype Bool:
 1  ⋅  ⋅  1  ⋅  1  1
 ⋅  1  ⋅  ⋅  1  ⋅  ⋅
 ⋅  ⋅  1  ⋅  ⋅  ⋅  ⋅

julia> y_model = reshape(vcat(-9:0, 0:9, 7.5f0), 3, 7)
3×7 Matrix{Float32}:
 -9.0  -6.0  -3.0  0.0  2.0  5.0  8.0
 -8.0  -5.0  -2.0  0.0  3.0  6.0  9.0
 -7.0  -4.0  -1.0  1.0  4.0  7.0  7.5

julia> Flux.logitcrossentropy(y_model, y_label)
1.5791205f0

julia> Flux.crossentropy(softmax(y_model), y_label)
1.5791197f0
```
"""
function logitcrossentropy(ŷ, y; dims = 1, agg = mean)
  _check_sizes(ŷ, y)
  agg(.-sum(y .* logsoftmax(ŷ; dims = dims); dims = dims))
end

"""
    binarycrossentropy(ŷ, y; agg = mean, ϵ = eps(ŷ))

Return the binary cross-entropy loss, computed as

    agg(@.(-y * log(ŷ + ϵ) - (1 - y) * log(1 - ŷ + ϵ)))

Where typically, the prediction `ŷ` is given by the output of a [`sigmoid`](@ref) activation.
The `ϵ` term is included to avoid infinity. Using [`logitbinarycrossentropy`](@ref) is recomended
over `binarycrossentropy` for numerical stability.

Use [`label_smoothing`](@ref) to smooth the `y` value as preprocessing before
computing the loss.

See also: [`crossentropy`](@ref), [`logitcrossentropy`](@ref).

# Examples
```jldoctest
julia> y_bin = Bool[1,0,1]
3-element Vector{Bool}:
 1
 0
 1

julia> y_prob = softmax(reshape(vcat(1:3, 3:5), 2, 3) .* 1f0)
2×3 Matrix{Float32}:
 0.268941  0.5  0.268941
 0.731059  0.5  0.731059

julia> Flux.binarycrossentropy(y_prob[2,:], y_bin)
0.43989f0

julia> all(p -> 0 < p < 1, y_prob[2,:])  # else DomainError
true

julia> y_hot = Flux.onehotbatch(y_bin, 0:1)
2×3 OneHotMatrix(::Vector{UInt32}) with eltype Bool:
 ⋅  1  ⋅
 1  ⋅  1

julia> Flux.crossentropy(y_prob, y_hot)
0.43989f0
```
"""
function binarycrossentropy(ŷ, y; agg = mean, ϵ = epseltype(ŷ))
  _check_sizes(ŷ, y)
  agg(@.(-xlogy(y, ŷ + ϵ) - xlogy(1 - y, 1 - ŷ + ϵ)))
end

"""
    logitbinarycrossentropy(ŷ, y; agg = mean)

Mathematically equivalent to
[`binarycrossentropy(σ(ŷ), y)`](@ref) but is more numerically stable.

See also: [`crossentropy`](@ref), [`logitcrossentropy`](@ref).

# Examples
```jldoctest
julia> y_bin = Bool[1,0,1];

julia> y_model = Float32[2, -1, pi]
3-element Vector{Float32}:
  2.0
 -1.0
  3.1415927

julia> Flux.logitbinarycrossentropy(y_model, y_bin)
0.160832f0

julia> Flux.binarycrossentropy(sigmoid.(y_model), y_bin)
0.16083185f0
```
"""
function logitbinarycrossentropy(ŷ, y; agg = mean)
  _check_sizes(ŷ, y)
  agg(@.((1 - y) * ŷ - logσ(ŷ)))
end

"""
    kldivergence(ŷ, y; agg = mean, ϵ = eps(ŷ))

Return the
[Kullback-Leibler divergence](https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence)
between the given probability distributions.

The KL divergence is a measure of how much one probability distribution is different
from the other. It is always non-negative, and zero only when both the distributions are equal.

# Example
```jldoctest
julia> p1 = [1 0; 0 1]
2×2 Matrix{Int64}:
 1  0
 0  1

julia> p2 = fill(0.5, 2, 2)
2×2 Matrix{Float64}:
 0.5  0.5
 0.5  0.5

julia> Flux.kldivergence(p2, p1) ≈ log(2)
true

julia> Flux.kldivergence(p2, p1; agg = sum) ≈ 2log(2)
true

julia> Flux.kldivergence(p2, p2; ϵ = 0)  # about -2e-16 with the regulator
0.0

julia> Flux.kldivergence(p1, p2; ϵ = 0)  # about 17.3 with the regulator
Inf
```
"""
function kldivergence(ŷ, y; dims = 1, agg = mean, ϵ = epseltype(ŷ))
  _check_sizes(ŷ, y)
  entropy = agg(sum(xlogx.(y), dims = dims))
  cross_entropy = crossentropy(ŷ, y; dims = dims, agg = agg, ϵ = ϵ)
  return entropy + cross_entropy
end

"""
    poisson_loss(ŷ, y)

# Return how much the predicted distribution `ŷ` diverges from the expected Poisson
# distribution `y`; calculated as `sum(ŷ .- y .* log.(ŷ)) / size(y, 2)`.

[More information.](https://peltarion.com/knowledge-center/documentation/modeling-view/build-an-ai-model/loss-functions/poisson).
"""
function poisson_loss(ŷ, y; agg = mean)
  _check_sizes(ŷ, y)
  agg(ŷ .- xlogy.(y, ŷ))
end

"""
    hinge_loss(ŷ, y; agg = mean)

Return the [hinge_loss loss](https://en.wikipedia.org/wiki/Hinge_loss) given the
prediction `ŷ` and true labels `y` (containing 1 or -1); calculated as
`sum(max.(0, 1 .- ŷ .* y)) / size(y, 2)`.

See also: [`squared_hinge_loss`](@ref)
"""
function hinge_loss(ŷ, y; agg = mean)
  _check_sizes(ŷ, y)
  agg(max.(0, 1 .- ŷ .* y))
end

"""
    squared_hinge_loss(ŷ, y)

Return the squared hinge_loss loss given the prediction `ŷ` and true labels `y`
(containing 1 or -1); calculated as `sum((max.(0, 1 .- ŷ .* y)).^2) / size(y, 2)`.

See also: [`hinge_loss`](@ref)
"""
function squared_hinge_loss(ŷ, y; agg = mean)
  _check_sizes(ŷ, y)
  agg((max.(0, 1 .- ŷ .* y)) .^ 2)
end

"""
    dice_coeff_loss(ŷ, y; smooth = 1)

Return a loss based on the dice coefficient.
Used in the [V-Net](https://arxiv.org/abs/1606.04797) image segmentation
architecture.
Similar to the F1_score. Calculated as:

    1 - 2*sum(|ŷ .* y| + smooth) / (sum(ŷ.^2) + sum(y.^2) + smooth)
"""
function dice_coeff_loss(ŷ, y; smooth = ofeltype(ŷ, 1.0))
  _check_sizes(ŷ, y)
  1 - (2 * sum(y .* ŷ) + smooth) / (sum(y .^ 2) + sum(ŷ .^ 2) + smooth) #TODO agg
end

"""
    tversky_loss(ŷ, y; β = 0.7)

Return the [Tversky loss](https://arxiv.org/abs/1706.05721).
Used with imbalanced data to give more weight to false negatives.
Larger β weigh recall more than precision (by placing more emphasis on false negatives)
Calculated as:
    1 - sum(|y .* ŷ| + 1) / (sum(y .* ŷ + β*(1 .- y) .* ŷ + (1 - β)*y .* (1 .- ŷ)) + 1)
"""
function tversky_loss(ŷ, y; β = ofeltype(ŷ, 0.7))
    _check_sizes(ŷ, y)
    #TODO add agg
    num = sum(y .* ŷ) + 1
    den = sum(y .* ŷ + β * (1 .- y) .* ŷ + (1 - β) * y .* (1 .- ŷ)) + 1
    1 - num / den
end

"""
    binary_focal_loss(ŷ, y; agg=mean, γ=2, ϵ=eps(ŷ))

Return the [binary_focal_loss](https://arxiv.org/pdf/1708.02002.pdf)
The input, 'ŷ', is expected to be normalized (i.e. [`softmax`](@ref) output).

For `γ == 0`, the loss is mathematically equivalent to [`Losses.binarycrossentropy`](@ref).

# Example
```jldoctest
julia> y = [0  1  0
            1  0  1]
2×3 Matrix{Int64}:
 0  1  0
 1  0  1

julia> ŷ = [0.268941  0.5  0.268941
            0.731059  0.5  0.731059]
2×3 Matrix{Float64}:
 0.268941  0.5  0.268941
 0.731059  0.5  0.731059

julia> Flux.binary_focal_loss(ŷ, y) ≈ 0.0728675615927385
true
```

See also: [`Losses.focal_loss`](@ref) for multi-class setting

"""
function binary_focal_loss(ŷ, y; agg=mean, γ=2, ϵ=epseltype(ŷ))
    _check_sizes(ŷ, y)
    ŷ = ŷ .+ ϵ
    p_t = y .* ŷ  + (1 .- y) .* (1 .- ŷ)
    ce = -log.(p_t)
    weight = (1 .- p_t) .^ γ
    loss = weight .* ce
    agg(loss)
end

"""
    focal_loss(ŷ, y; dims=1, agg=mean, γ=2, ϵ=eps(ŷ))

Return the [focal_loss](https://arxiv.org/pdf/1708.02002.pdf)
which can be used in classification tasks with highly imbalanced classes.
It down-weights well-classified examples and focuses on hard examples.
The input, 'ŷ', is expected to be normalized (i.e. [`softmax`](@ref) output).

The modulating factor, `γ`, controls the down-weighting strength.
For `γ == 0`, the loss is mathematically equivalent to [`Losses.crossentropy`](@ref).

# Example
```jldoctest
julia> y = [1  0  0  0  1
            0  1  0  1  0
            0  0  1  0  0]
3×5 Matrix{Int64}:
 1  0  0  0  1
 0  1  0  1  0
 0  0  1  0  0

julia> ŷ = softmax(reshape(-7:7, 3, 5) .* 1f0)
3×5 Matrix{Float32}:
 0.0900306  0.0900306  0.0900306  0.0900306  0.0900306
 0.244728   0.244728   0.244728   0.244728   0.244728
 0.665241   0.665241   0.665241   0.665241   0.665241

julia> Flux.focal_loss(ŷ, y) ≈ 1.1277571935622628
true
```

See also: [`Losses.binary_focal_loss`](@ref) for binary (not one-hot) labels

"""
function focal_loss(ŷ, y; dims=1, agg=mean, γ=2, ϵ=epseltype(ŷ))
    _check_sizes(ŷ, y)
    ŷ = ŷ .+ ϵ
    agg(sum(@. -y * (1 - ŷ)^γ * log(ŷ); dims=dims))
end

"""
    siamese_contrastive_loss(ŷ, y; margin = 1, agg = mean)

Return the [contrastive loss](http://yann.lecun.com/exdb/publis/pdf/hadsell-chopra-lecun-06.pdf)
which can be useful for training Siamese Networks. It is given by

    agg(@. (1 - y) * ŷ^2 + y * max(0, margin - ŷ)^2)

Specify `margin` to set the baseline for distance at which pairs are dissimilar.

"""
function siamese_contrastive_loss(ŷ, y; agg = mean, margin::Real = 1)
    _check_sizes(ŷ, y)
    margin < 0 && throw(DomainError(margin, "Margin must be non-negative"))
    return agg(@. (1 - y) * ŷ^2 + y * max(0, margin - ŷ)^2)
end

```@meta
DocTestFilters = nothing
```


text:
using SpecialFunctions

#    _|_|_|    _|_|    _|_|_|    _|      _|
#  _|        _|    _|  _|    _|  _|      _|
#  _|        _|    _|  _|    _|    _|  _|
#    _|_|_|    _|_|    _|    _|      _|

@makepdftype myBW(x;p) = let xth = 2.961477
    AlgebraPDF.amplitudeBWsq(x+xth, p.Δm+xth, p.Γ)
end
mybw = myBW(Ext(Δm=0.1, Γ=3e-3), (0,0.22))

# convolute with average resolution
const σ_mean = 5e-3
mybw_conv = convGauss(mybw, σ_mean)

@testset "convGauss: update pars works" begin
    @test pars(updatepars(mybw_conv, (Δm=0.15,))).Δm == 0.15
end

function updategetmaxdensity(d, Δm)
    xv = range(0,0.22,length=100)
    maximum(updatepars(d,(Δm=Δm,))(xv))
end

# resolution is linear function of energy
resolution_model = FunctionWithParameters((x;p)->σ_mean+(x-0.1)*3e-2; p=∅)
mybw_conv_dep = convGauss(mybw, resolution_model)

@testset "convGauss energy-dependent gets wider" begin
    # if the peak located below 0.1 the resolution the effective resolution  < σ_mean => the peak max is higher
    @test updategetmaxdensity(mybw_conv_dep, 0.05) > updategetmaxdensity(mybw_conv, 0.05)
    # if the peak located above 0.1 the resolution the effective resolution  > σ_mean => the peak max is smaller
    @test updategetmaxdensity(mybw_conv_dep, 0.15) < updategetmaxdensity(mybw_conv, 0.15)
end


text:
module ConstraintTransforms

using ..Model
using ..SensitiveFloats
using ..DeterministicVI: VariationalParams
using Celeste: Const, @aliasscope
using Compat

####################
# Constraint Types #
####################

abstract type Constraint end

struct BoxConstraint <: Constraint
    lower::Float64
    upper::Float64
    scale::Float64
    function BoxConstraint(lower, upper, scale)
        @assert lower > -Inf
        @assert lower < upper
        @assert upper < Inf
        @assert 0.0 < scale < Inf
        return new(lower, upper, scale)
    end
end

struct SimplexConstraint <: Constraint
    lower::Float64
    scale::Float64
    n::Int # length of the simplex in the free parameterization
    function SimplexConstraint(lower, scale, n)
        @assert n >= 2
        @assert 0.0 <= lower < 1/n
        @assert 0.0 < scale < Inf
        return new(lower, scale, n)
    end
end

struct ParameterConstraint{C<:Constraint} <: Constraint
    constraint::C
    indices::UnitRange{Int} # indices in the bound parameterization
end

ParameterConstraint(c::Constraint, i::Vector) = ParameterConstraint(c, first(i):last(i))
ParameterConstraint(c::Constraint, i::Int) = ParameterConstraint(c, i:i)

struct ConstraintBatch
    boxes::Vector{Vector{ParameterConstraint{BoxConstraint}}}
    simplexes::Vector{Vector{ParameterConstraint{SimplexConstraint}}}
    function ConstraintBatch(boxes, simplexes)
        @assert length(boxes) == length(simplexes)
        return new(boxes, simplexes)
    end
end

#################################################
# `to_bound`/`to_bound!` & `to_free`/`to_free!` #
#################################################

# BoxConstraint #
#---------------#

# to_bound

logit(x) = 1.0 / (1.0 + exp(-x))

function to_bound(free::Real, c::BoxConstraint)
    @assert isfinite(free)
    return logit(free / c.scale) * (c.upper - c.lower) + c.lower
end

# to_free

inv_logit(x) = -log(1.0 / x - 1)

function to_free(bound::Real, c::BoxConstraint)
    @assert c.lower < bound < c.upper
    return inv_logit((bound - c.lower) / (c.upper - c.lower)) * c.scale
end

# SimplexConstraint #
#-------------------#

simplexify(x::Real, c::SimplexConstraint) = (1 - c.n * c.lower) * x + c.lower

unsimplexify(x::Real, c::SimplexConstraint) = (x - c.lower) / (1 - c.n * c.lower)

function to_bound!(bound::Vector{T}, bound_range::UnitRange,
                   free::Vector{T}, free_range::UnitRange,
                   c::SimplexConstraint) where {T<:Real}
    for i in 1:length(free_range)
        x = free[free_range[i]]
        @assert isfinite(x)
        bound[bound_range[i]] = x / c.scale
    end
    m = bound[bound_range[1]]
    for i in 1:length(free_range)
        m = max(bound[bound_range[i]], m)
    end
    exp_neg_m = exp(-m)
    bound_sum = exp_neg_m
    for i in bound_range[1:end-1]
        x = exp(bound[i] - m)
        bound[i] = x
        bound_sum += x
    end
    for i in bound_range[1:end-1]
        bound[i] = simplexify(bound[i] / bound_sum, c)
    end
    bound[bound_range[end]] = simplexify(inv(bound_sum) * exp_neg_m, c)
    return bound
end

function to_free!(free::Vector{T}, free_range::UnitRange,
                  bound::Vector{T}, bound_range::UnitRange,
                  c::SimplexConstraint) where {T<:Real}
    log_last = log(unsimplexify(bound[bound_range[end]], c))
    for i in 1:length(free_range)
        ith_term = log(unsimplexify(bound[bound_range[i]], c)) - log_last
        free[free_range[i]] = c.scale * ith_term
    end
    return free
end

# ParameterConstraint #
#---------------------#

function to_bound!(bound::Vector{T}, free::Vector{T},
                   c::ParameterConstraint{BoxConstraint},
                   free_index::Integer) where {T<:Real}
    for i in c.indices
        bound[i] = to_bound(free[free_index], c.constraint)
        free_index += 1
    end
    return free_index
end

function to_bound!(bound::Vector{T}, free::Vector{T},
                   c::ParameterConstraint{SimplexConstraint},
                   free_index::Integer) where {T<:Real}
    free_length = (c.constraint.n - 1)
    free_range = (1:free_length) + free_index - 1
    to_bound!(bound, c.indices, free, free_range, c.constraint)
    return free_index + free_length
end

function to_free!(free::Vector{T}, bound::Vector{T},
                  c::ParameterConstraint{BoxConstraint},
                  free_index::Integer) where {T<:Real}
    for i in c.indices
        free[free_index] = to_free(bound[i], c.constraint)
        free_index += 1
    end
    return free_index
end

function to_free!(free::Vector{T}, bound::Vector{T},
                  c::ParameterConstraint{SimplexConstraint},
                  free_index::Integer) where {T<:Real}
    free_length = (c.constraint.n - 1)
    free_range = (1:free_length) + free_index - 1
    to_free!(free, free_range, bound, c.indices, c.constraint)
    return free_index + free_length
end

# ConstraintBatch #
#-----------------#

function to_bound!(bound::Vector{T}, free::Vector{T},
                   constraints::ConstraintBatch, src::Integer) where {T<:Real}
    free_index = 1
    for constraint in constraints.boxes[src]
        free_index = to_bound!(bound, free, constraint, free_index)
    end
    for constraint in constraints.simplexes[src]
        free_index = to_bound!(bound, free, constraint, free_index)
    end
end

function to_bound!(bound::VariationalParams{T}, free::VariationalParams{T},
                   constraints::ConstraintBatch) where {T}
    @assert length(free) == length(bound)
    for src in 1:length(free)
        to_bound!(bound[src], free[src], constraints, src)
    end
end

function to_free!(free::Vector{T}, bound::Vector{T},
                  constraints::ConstraintBatch, src::Integer) where {T<:Real}
    free_index = 1
    for constraint in constraints.boxes[src]
        free_index = to_free!(free, bound, constraint, free_index)
    end
    for constraint in constraints.simplexes[src]
        free_index = to_free!(free, bound, constraint, free_index)
    end
end

function to_free!(free::VariationalParams{T}, bound::VariationalParams{T},
                  constraints::ConstraintBatch) where {T}
    @assert length(free) == length(bound)
    for src in 1:length(free)
        to_free!(free[src], bound[src], constraints, src)
    end
end

########################
# `enforce`/`enforce!` #
########################

# BoxConstraint #
#---------------#

function enforce(bound::Number, c::BoxConstraint)
    if !(c.lower < bound < c.upper)
        # don't set the value to exact bounds to avoid Infs during transforms
        return max(min(bound, prevfloat(c.upper)), nextfloat(c.lower))
    end
    return bound
end

# SimplexConstraint #
#-------------------#

function enforce!(bound, c::SimplexConstraint)
    bound_sum = zero(eltype(bound))
    for i in eachindex(bound)
        if !(c.lower < bound[i] < 1.0)
            # don't set the value to exact bounds to avoid Infs during transforms
            bound[i] = max(min(bound[i], prevfloat(1.0)), nextfloat(c.lower))
        end
        bound_sum += bound[i]
    end
    if !(isapprox(bound_sum, one(eltype(bound))))
        # normalize in a way that maintains the lower bounds
        rescale = (1 - c.n * c.lower) / (bound_sum - c.n * c.lower)
        for i in eachindex(bound)
            bound[i] = nextfloat(c.lower) + rescale * (bound[i] - c.lower)
        end
    end
    return bound
end

# ParameterConstraint #
#---------------------#

function enforce!(bound::Vector{T}, c::ParameterConstraint{BoxConstraint}) where {T<:Real}
    for i in c.indices
        bound[i] = enforce(bound[i], c.constraint)
    end
    return bound
end

function enforce!(bound::Vector{T}, c::ParameterConstraint{SimplexConstraint}) where {T<:Real}
    enforce!(view(bound, c.indices), c.constraint)
    return bound
end

# ConstraintBatch #
#-----------------#

function enforce!(bound::Vector{T},
                  constraints::ConstraintBatch,
                  src::Integer) where {T<:Real}
    for constraint in constraints.boxes[src]
        enforce!(bound, constraint)
    end
    for constraint in constraints.simplexes[src]
        enforce!(bound, constraint)
    end
end

function enforce!(bound::VariationalParams{T},
                  constraints::ConstraintBatch) where {T}
    for src in 1:length(bound)
        enforce!(bound[src], constraints, src)
    end
end

##########################
# `allocate_free_params` #
##########################

function free_length(boxes::Vector{ParameterConstraint{BoxConstraint}})
    return sum(length(box.indices) for box in boxes)
end

function free_length(simplexes::Vector{ParameterConstraint{SimplexConstraint}})
    return sum(simplex.constraint.n - 1 for simplex in simplexes)
end

function allocate_free_params(bound::VariationalParams{T}, constraints::ConstraintBatch) where {T}
    free = similar(bound)
    for src in 1:length(bound)
        n_simplex_params = free_length(constraints.simplexes[src])
        n_box_params = free_length(constraints.boxes[src])
        free[src] = zeros(T, n_simplex_params + n_box_params)
    end
    return free
end

###################
# Differentiation #
###################

using ForwardDiff
using ForwardDiff: Dual, jacobian!, JacobianConfig

const NUM_ELBO_BOUND_PARAMS = length(CanonicalParams)

const TransformJacobianConfig{M,N,T} = JacobianConfig{N,T,Tuple{Array{Dual{N,T},M},Vector{Dual{N,T}}}}

struct TransformDerivatives{N,T}
    jacobian::Matrix{T}
    hessian::Array{T,3}
    output_dual::Vector{Dual{N,T}}
    inner_cfg::TransformJacobianConfig{1,N,Dual{N,T}}
    outer_cfg::TransformJacobianConfig{2,N,T}
end

# this is a good chunk size because it divides evenly into `length(CanonicalParams)`
const DEFAULT_CHUNK = 11

function TransformDerivatives(
    output::Vector{T}, input::Vector{T},
    ::Type{Val{N}} = Val{DEFAULT_CHUNK}) where {T<:Real,N}

    jacobian = zeros(T, length(output), length(input))
    hessian = zeros(T, length(output), length(input), length(input))
    output_dual = copy!(similar(output, Dual{N,T}), output)
    inner_cfg = JacobianConfig{N}(output_dual, similar(input, Dual{N,T}))
    outer_cfg = JacobianConfig{N}(jacobian, input)
    return TransformDerivatives{N,T}(jacobian, hessian, output_dual, inner_cfg, outer_cfg)
end

function TransformDerivatives(
    output::VariationalParams{T},
    input::VariationalParams{T},
    ::Type{Val{N}} = Val{DEFAULT_CHUNK}) where {T<:Real,N}

    @assert length(output) == length(input)
    @assert all(length(src) == length(output[1]) for src in output)
    @assert all(length(src) == length(input[1]) for src in input)
    return TransformDerivatives(output[1], input[1], Val{N})
end

function differentiate!(transform!::F,
                        derivs::TransformDerivatives,
                        input::Vector{T}) where {F,T<:Number}
    jacobian!(reshape(derivs.hessian, length(derivs.jacobian), length(input)),
              (out, x) -> jacobian!(out, transform!, derivs.output_dual, x, derivs.inner_cfg),
              derivs.jacobian, input, derivs.outer_cfg)
    return derivs
end

# Propagate derivatives of `transform!` back from the output-parameterized SensitiveFloat
# (`sf_output`) to the input-parameterized SensitiveFloat (`sf_input`). Note that the memory
# for the raw output parameters is supplied via `derivs`, whose constructor is
# `TransformDerivatives(output_params, input_params)`.
function propagate_derivatives!(transform!::F,
                                sf_output::SensitiveFloat,
                                sf_input::SensitiveFloat,
                                input_sources::VariationalParams{T},
                                constraints::ConstraintBatch,
                                derivs::TransformDerivatives) where {F,T}
    sf_input.v[] = sf_output.v[]
    n_sources = length(input_sources)
    n_input_params = size(sf_input.d, 1)
    n_output_params = size(sf_output.d, 1)
    input_hessian = reshape(sf_input.h, n_sources, n_input_params, n_sources, n_input_params)
    output_hessian = reshape(sf_output.h, n_sources, n_output_params, n_sources, n_output_params)
    # Julia v0.5 erroneously fails to type-infer `src` for some reason, so we type-annotate
    # it explicitly here. It still unnecessarily boxes/unboxes `src` for some of the uses
    # in the body of the loop, but I'm not sure what to do for that.
    for src::Int64 in 1:n_sources
        differentiate!((y, x) -> transform!(y, x, constraints, src), derivs, input_sources[src])
        backprop_jacobian!(sf_input.d, sf_output.d, derivs.jacobian, src)
        backprop_hessian!(input_hessian, output_hessian, derivs.hessian, derivs.jacobian, sf_output.d, src)
    end
    symmetrize!(sf_input.h)
    return sf_input
end

# equivalent to At_mul_B!(input[:, src], jacobian, output[:, src])
function backprop_jacobian!(input, output, jacobian, src)
    for i in 1:size(input, 1)
        x = zero(eltype(input))
        for j in 1:size(output, 1)
            @inbounds x += jacobian[j, i] * output[j, src]
        end
        @inbounds input[i, src] = x
    end
    return input
end

function backprop_hessian!(input, output, hessian, jacobian, gradient, src)
    first_quad_form!(input, jacobian, output, src)
    for j in 1:size(hessian, 2)
        for k in 1:size(hessian, 3)
            x = zero(eltype(input))
            for i in 1:size(hessian, 1)
                @inbounds x += hessian[i, j, k] * gradient[i, src]
            end
            @inbounds input[src, j, src, k] += x
        end
    end
    return input
end

# equivalent to `At_mul_B!(view(C, src, :, src, :), A, view(B, src, :, src, :) * A)`
# this could be further optimized by assuming the symmetry of `view(B, src, :, src, :)`
function first_quad_form!(C, A, B, src)
    const m = NUM_ELBO_BOUND_PARAMS
    @assert m == size(A, 1)
    n = size(A, 2)
    scratch = Array{Float64, 2}(m, m)
    @aliasscope begin
        for i in 1:m, j in 1:m
            scratch[i, j] = Const(B)[src, j, src, i]
        end
    end
    @aliasscope begin
        for i in 1:n, j in 1:n
            x = zero(eltype(C))
            for k in 1:m
                y = zero(eltype(C))
                @inbounds for l in 1:m
                    @fastmath y += Const(scratch)[l, k] * Const(A)[l, j]
                end
                @inbounds x += Const(A)[k, i] * y
            end
            @inbounds C[src, i, src, j] = x
        end
    end
    return C
end

# ensure exact symmetry (necessary for some numerical linear algebra routines)
function symmetrize!(A, c = 0.5)
    for i in 1:size(A, 1), j in 1:i
        @inbounds A[i, j] = A[j, i] = c * (A[i, j] + A[j, i])
    end
    return A
end

end # module


text:
using ITensors.NDTensors,
      Test

@testset "DenseTensor basic functionality" begin

  A = Tensor(3,4)

  for I in eachindex(A)
    @test A[I] == 0
  end

  randn!(A)

  for I in eachindex(A)
    @test A[I] != 0
  end

  for I in eachindex(A)
    @test A[I] != 0
  end

  @test ndims(A) == 2
  @test dims(A) == (3,4)
  @test inds(A) == (3,4)

  A[1,1] = 11

  @test A[1,1] == 11

  Aview = A[2:3,2:3]

  @test dims(Aview) == (2,2)
  @test A[2,2] == Aview[1,1]

  @test A * 2.0 == 2.0 * A

  Asim = similar(data(A), 10)
  @test eltype(Asim) == Float64
  @test length(Asim) == 10

  B = Tensor(undef,3,4)
  randn!(B)

  C = A+B

  for I in eachindex(C)
    @test C[I] == A[I] + B[I]
  end

  Ap = permutedims(A,(2,1))

  for I in eachindex(A)
    @test A[I] == Ap[NDTensors.permute(I,(2,1))]
  end

  t = Tensor(ComplexF64,100,100)
  randn!(t)
  @test conj(data(store(t))) == data(store(conj(t)))
  @test typeof(conj(t)) <: DenseTensor

  @test Dense(ComplexF64) == Dense{ComplexF64}()
  @test Dense(ComplexF64) == complex(Dense(Float64))

  D = Tensor(ComplexF64,(100,100))
  @test eltype(D) == ComplexF64
  @test ndims(D) == 2
  @test dim(D) == 100^2

  E = Tensor(ComplexF64,undef, (100,100))
  @test eltype(E) == ComplexF64
  @test ndims(E) == 2
  @test dim(E) == 100^2

  F = Tensor((100,100))
  @test eltype(F) == Float64
  @test ndims(F) == 2
  @test dim(F) == 100^2

  G = Tensor(undef, (100,100))
  @test eltype(G) == Float64
  @test ndims(G) == 2
  @test dim(G) == 100^2

  H = Tensor(ComplexF64,undef, 100,100)
  @test eltype(H) == ComplexF64
  @test ndims(H) == 2
  @test dim(H) == 100^2

  I_arr = rand(10,10,10)
  I = Tensor(I_arr, (10,10,10))
  @test eltype(I) == Float64
  @test dim(I) == 1000
  @test Array(I) == I_arr

  J = Tensor(2,2)
  K = Tensor(2,2)
  @test Array(J * K) ≈ Array(J) * Array(K)
end


text:
module RotationsCurrence

include("./constants.jl")
using JuMP
using .Constants

export cons_rotations_currence

#Constraints ensuring that the only allowed rotations are those of Currence Robotics
function cons_rotations_currence(model)

    # lx + ly = 1
    @constraint(model,
    lxly[i=1:length(boxes)], (model[:lpx_length_of_box_parallel_to_x][i] +
    model[:lpy_length_of_box_parallel_to_y][i] == 1)
    )

    #wpx + wpy =1
    @constraint(model,
    wxwy[i=1:length(boxes)], (model[:wpx_width_of_box_parallel_to_x][i] +
    model[:wpy_width_of_box_parallel_to_y][i] == 1)
    )

    #lx + wx = 1
    @constraint(model,
    lxwx[i=1:length(boxes)], (model[:lpx_length_of_box_parallel_to_x][i] +
    model[:wpx_width_of_box_parallel_to_x][i] == 1)
    )

    #ly + wy = 1
    @constraint(model,
    lywy[i=1:length(boxes)], (model[:lpy_length_of_box_parallel_to_y][i] +
    model[:wpy_width_of_box_parallel_to_y][i] == 1)
    )

    #hz = 1, hz, hy, lz, wz = 0
    @constraints(model, begin
    [i=1:length(boxes)], model[:hpz_height_of_box_parallel_to_z][i] == 1
    [i=1:length(boxes)], model[:hpx_height_of_box_parallel_to_x][i] == 0
    [i=1:length(boxes)], model[:hpy_height_of_box_parallel_to_y][i] == 0
    [i=1:length(boxes)], model[:lpz_length_of_box_parallel_to_z][i] == 0
    [i=1:length(boxes)], model[:wpz_width_of_box_parallel_to_z][i] == 0
    end)

    return model
end

end


text:
module ProjC

const VERSION = v"0.1"

greet() = print("Hello World!")

end # module


text:
# for GaussianNode:
using Distributions

################################################################
# STATE OF A NODE
################################################################

"""
A node can be in one of three states:
* unmarked
* temporary
* permanently
This is necessary for graph algorithms, especially for finding
a topological order of the nodes in an SPN.
"""
@enum State unmarked temporary permanently

################################################################
# ABSTRACT NODES
################################################################

"""
Every node of an SPN has to be of a subtype of `Node`.

Every subtype of `Node` needs to have the fields (even though Julia does not enforce it):

* `parents`:  a list of parent nodes
* `scope`: a list of variables which influence the node
* `logval`: the current value of the node
* `logdrv`: the current value of the log derivative
* `state`: a symbol representing the state of the node for graph traversals
"""
abstract type Node end

"""
Inner nodes are either sum nodes or product nodes.

They must have an additional field `children::Vector{Node}` that is a list of their child nodes.
"""
abstract type InnerNode <: Node end

"""
Leaf nodes are nodes without children that represent some univariate probability distribution.
"""
abstract type LeafNode <: Node end

################################################################
# INNER NODES
################################################################

"""
A product node computes the product of the values of its child nodes.
"""
mutable struct ProdNode <: InnerNode
	"The log-likelihood value of this node."
	logval::Float64
    "The log-derivative value of this node."
    logdrv::Float64
    "The node can have a state. Necessary for graph traversals."
    state::State

	"The parent nodes of this node."
	parents::Vector{InnerNode}
	"The child nodes of this node."
	children::Vector{Node}
	"The scope of this node (a list of variables which influence the node)."
	scope::Vector{Int}
end

"""
Creates a new product node.
"""
function ProdNode()
	logval = -Inf  # The default logval is log(0)=-Inf
    logdrv = -Inf  # The default logdrv is log(0)=-Inf
    state = unmarked
	parents = InnerNode[]
	children = Node[]
	scope = Int[]
	ProdNode(logval, logdrv, state, parents, children, scope)
end

"""
A sum node computes a weighted sum of the values of its child nodes.
"""
mutable struct SumNode <: InnerNode
	"The log-likelihood value of this node."
	logval::Float64
    "The log-derivative value of this node."
    logdrv::Float64
    "The index of the child with highest weighted value. Necessary for MPE inference."
    maxidx::Int
    "The node can have a state. Necessary for graph traversals."
    state::State

	"The parent nodes of this node."
	parents::Vector{InnerNode}
	"The child nodes of this node."
	children::Vector{Node}
	"The scope of this node (a list of variables which influence the node)."
	scope::Vector{Int}

	"The weights of the edges connecting the sum node with its child nodes."
	weights::Vector{Float64}
    # counts for each child. The EM algorithm uses these to compute new weights.
    counts::Vector{Float64}
end

"""
Creates a new sum node.
"""
function SumNode()
	logval = -Inf  # The default logval is log(0)=-Inf
    logdrv = -Inf  # The default logdrv is log(0)=-Inf
    maxidx = -1
    state = unmarked
	parents = InnerNode[]
	children = Node[]
	scope = Int[]

	weights = Float64[]
    counts = Float64[]

	SumNode(logval, logdrv, maxidx, state, parents, children, scope, weights, counts)
end

################################################################
# LEAF NODES
# A leaf node represents some univariate probability
# distribution.
################################################################

"""
An indicator node for a random variable X and value k is 1 if X=k
for the current input example from the dataset.
Otherwise it has a value of 0. The corresponding log-values are
log(1) = 0 and log(0)=-Inf.
"""
mutable struct IndicatorNode <: LeafNode
    "the log-likelihood value of this node."
    logval::Float64
    "The log-derivative value of this node."
    logdrv::Float64
    "The node can have a state. Necessary for graph traversals."
    state::State

    "The parent nodes of this node."
    parents::Vector{InnerNode}
    "the scope of this node."
    scope::Vector{Int}
    "The value of the (discrete) variable that the node indicates."
    indicates::Float64
end

"""
	IndicatorNode(varidx::Int, indicates::Int) -> IndicatorNode

Creates a new indicator node for a random variable with floating point values.
"""
function IndicatorNode(varidx::Int, indicates::Float64)
	@assert varidx > 0
    logval = -Inf
    logdrv = -Inf
    state = unmarked
    parents = InnerNode[]
    scope = Int[varidx]
    IndicatorNode(logval, logdrv, state, parents, scope, indicates)
end

"""
Creates a new indicator node for a random variable with integer values.
"""
IndicatorNode(varidx::Int, indicates::Int) = IndicatorNode(varidx, float(indicates))

"""
Creates a new indicator node for a random variable with boolean values.
"""
IndicatorNode(varidx::Int, indicates::Bool) = IndicatorNode(varidx, float(indicates))

"""
A `GaussianNode` represents a univariate Gaussian distribution.
"""
mutable struct GaussianNode <: LeafNode
    "the log-likelihood value of this node."
    logval::Float64
    "The log-derivative value of this node."
    logdrv::Float64
    "The node can have a state. Necessary for graph traversals."
    state::State

    "The parent nodes of this node."
    parents::Vector{InnerNode}
    "the scope of this node."
    scope::Vector{Int}

    "The mean of the Gaussian."
    μ::Float64
    "The standard deviation of the Gaussian"
    σ::Float64
    "The distribution object that defines the Gaussian."
    distr::Distributions.Normal{Float64}
end

"""
    GaussianNode(varidx::Int, μ::Float64, σ::Float64) -> GaussianNode

Creates a new Gaussian node.
"""
function GaussianNode(varidx::Int, μ::Float64, σ::Float64)
    @assert varidx > 0

    logval = -Inf
    logdrv = -Inf
    state = unmarked
    parents = InnerNode[]
    scope = Int[varidx]
    distr = Normal(μ,σ)

    GaussianNode(logval,logdrv,state,parents,scope,μ,σ,distr)
end

################################################################
# CONNECTING NODES
# We can connect nodes by adding the child node to the list of
# children of the parent and the parent to the list of parents
# of the child. Note that the parent has to be an inner node.
################################################################

"""
    connect!(parent::InnerNode, child::Node)

Connects two nodes with an edge.
"""
function connect!(parent::InnerNode, child::Node)
    push!(parent.children, child)
    push!(child.parents, parent)
end

"""
    connect!(parent::SumNode, child::Node; weight=rand())

Connects two nodes with an edge, where the parent node is a
sum node. Because edges emanating from a sum node have weights
assigned to them, we can specify them with a keyword argument.
If no weight is specified, a weight is chosen at random.
"""
function connect!(parent::SumNode, child::Node; weight=rand())
    push!(parent.children, child)
    push!(child.parents, parent)
    push!(parent.weights, weight)
    push!(parent.counts, 0.0)
end

################################################################
# DISPLAYING NODES
################################################################

"""
Display an inner node.
"""
function Base.show(io::IO, n::InnerNode)
   print(io, "$(typeof(n))(parents=$(length(n.parents)), children=$(length(n.children)), logval=$(n.logval))")
end

"""
Display an indicator node.
"""
function Base.show(io::IO, i::IndicatorNode)
   print(io, "$(typeof(i))(parents=$(length(i.parents)), scope=$(i.scope[1]), indicates=$(i.indicates), logval=$(i.logval))")
end

"""
Display a Gaussian node.
"""
function Base.show(io::IO, g::GaussianNode)
   print(io, "$(typeof(g))(parents=$(length(g.parents)), scope=$(g.scope[1]), μ=$(g.μ), σ=$(g.σ), logval=$(g.logval))")
end

################################################################
# EVALUATING NODES
# An SPN is evaluated by an upward pass through the network,
# evaluating one node after another. Evaluating a node means
# setting its value according to the current input.
# Working in logspace is less prone to numerical problems,
# such as underflow resulting from multiplying together several
# very small probabilities.
################################################################

"""
    setInput!(i::IndicatorNode, x::AbstractVector, e::BitVector)

Sets the logval of the indicator node `i` according to the input `x`.
The bitvector `e` represents which values are in the evidence
and which are unknown.

Let X be the variable of the indicator and k the indicated value.
The value is log(1)=0 if X==k or if the value of X is unknown.
Otherwise the value of the indicator is log(0)=-Inf.
"""
function setInput!(i::IndicatorNode, x::AbstractVector, e::BitVector=trues(length(x)))
    idx = i.scope[1]  # get the column index of the variable i indicates
    @assert length(x) >= idx
    if e[idx] == false  # variable not in evidence
        i.logval = 0.0
    elseif x[idx] ≈ i.indicates
        i.logval = 0.0
    else
        i.logval = -Inf
    end
end

"""
    setInput!(g::GaussianNode, x::AbstractVector, e::BitVector)

Sets the logval of the Gaussian node `g` according to the input `x`.
The bitvector `e` represents which values are in the evidence
and which are unknown.
"""
function setInput!(g::GaussianNode, x::AbstractVector, e::BitVector=trues(length(x)))
    idx = g.scope[1]  # get the column index of the variable g represents
    @assert length(x) >= idx
    if e[idx] == false  # variable not in evidence
        g.logval = 0.0
    else
        g.logval = logpdf(g.distr, x[idx])  # evaluate Gaussian
    end
end

"""
    eval!(l::LeafNode; max::Bool=false; max::Bool=false) -> Float64

Returns the log value of the leaf node `l` on the current input.
The logval is already computed by `setInput!`.
"""
function eval!(l::LeafNode; max::Bool=false)
    return l.logval
end

"""
    eval!(p::ProdNode; max::Bool=false) -> Float64

Computes the log value of the product node `p` on the current input.

The value of a product node is the product of the values of its children.
Therefore, the log value is the sum of the log values of its children.
"""
function eval!(p::ProdNode; max::Bool=false)
    childvalues = [child.logval for child in p.children]
    p.logval = sum(childvalues)
    return p.logval
end

"""
    eval!(s::SumNode; max::Bool=false) -> Float64

Compute the log value of the sum node `s` on the current input.

The value of a sum node is the sum of the values of its children.
In log-space this looks a lot more ugly:
log(S_i) = log(sum_j w_ij S_j) = log(sum_j exp(log(w_ij)) * exp(log(S_j)))
= log(sum_j exp(log(w_ij) + log(S_j)).
If the keyword argument `max` is set to true the value of the sum node is its maximum weighted child value.
In log-space this means max(log(w_ij) + s_j).
"""
function eval!(s::SumNode; max::Bool=false)
    childvalues = [child.logval for child in s.children]
    numChildren = length(childvalues)

    sum_val = 0.0
    max_val = -Inf  # the maximum weighted childvalue
    max_idx = -1
    for (w, cval, idx) in zip(s.weights, childvalues, 1:numChildren)
        weighted_cval = cval + log(w)
        sum_val += exp(weighted_cval)
        if weighted_cval > max_val
            max_val = weighted_cval
            max_idx = idx
        end
    end

    if max == false
        s.logval = log(sum_val)
    else
        #s.logval = max_val
        s.logval = childvalues[max_idx]   # TODO: is this correct?
    end
    s.maxidx = max_idx

    return s.logval
end

################################################################
# COMPUTING NODE DERIVATIVES
# We require the derivatives of the likelihood with respect
# to the nodes of the network: ∂S/∂S_i. Because values are
# treated in logspace, derivatives are also in logspace.
# That means, we compute log(∂S/∂S_i).
# The derivatives can be computed using backpropagation.
################################################################

"""
Computes the derivative
"""
function passDerivative!(n::InnerNode)
    for (i, child) in enumerate(n.children)
        # compute l value: l=∂S_n/∂S_i
        if typeof(n) == SumNode
            l = log(n.weights[i])
        else  # n is ProdNode
            l = n.logval - child.logval
        end

        if child.logdrv == -Inf
            child.logdrv = n.logdrv + l
        else
            child.logdrv = addExp(child.logdrv, n.logdrv + l)
        end
    end
end

################################################################
# MISCELLANEOUS FUNCTIONS
################################################################

"""
    normalize!(s::SumNode)

Normalizes the weights of the edges emanating from a sum node,
such that they sum up to 1.
"""
function Base.normalize!(s::SumNode)
    ε::Float64 = 1e-10
    s.weights = s.weights / (sum(s.weights)+ε)
end

"""
    resetCounts!(s::SumNode)

Sets the counts for each child of a `SumNode` to 0.
"""
function resetCounts!(s::SumNode)
    s.counts = zeros(Float64, length(s.weights))
end


text:
using GZip

import Base: eof, close, position
export BamReader, close, value, eof, advance!, eachposition

type BamReader
    bamStream
    done::Bool
    tid_to_refname::Dict{Int32,ASCIIString}
end

type Cigar
    length::UInt32
    operation::UInt8
end

type Read
    data::Array{Int8}
end

function BamReader(bamFileName::UTF8String)
    BamReader(ascii(bamFileName))
end

function BamReader(bamFileName::ASCIIString)
    f = GZip.open(bamFileName)

    # make sure this is a BAM file
    code = read(f, UInt8, 4)
    @assert code == b"BAM\1"

    # get through the header data
    l_text = read(f, Int32)
    skip(f, l_text)

    # make sure the contigs match our reference
    n_ref = read(f, Int32)

    tid_to_refname = Dict{Int32,ASCIIString}()

    for j in 1:n_ref
        l_name = read(f, Int32)
        refName = convert(ASCIIString, read(f, UInt8, l_name)[1:end-1]) # ignore the null terminator
        l_ref = read(f, Int32)
        tid_to_refname[j] = refName
    end

    BamReader(f, false, tid_to_refname)
end

close(read::Read) = GZip.close(read.bamStream)
refid(read::Read) = reinterpret(Int32,read.data[1:4])[1] + 1
position(read::Read) = reinterpret(Int32,read.data[5:8])[1] + 1
l_read_name(read::Read) = read.data[9]
mq(read::Read) = read.data[10]
bin(read::Read) = reinterpret(Int16,read.data[11:12])[1]
n_cigar_op(read::Read) = reinterpret(Int16,read.data[13:14])[1]
flag(read::Read) = reinterpret(Int16,read.data[15:16])[1]
l_seq(read::Read) = reinterpret(Int32,read.data[17:20])[1]
next_refid(read::Read) = reinterpret(Int32,read.data[21:24])[1] + 1
next_position(read::Read) = reinterpret(Int32,read.data[25:28])[1] + 1
tlen(read::Read) = reinterpret(Int32,read.data[29:32])[1] + 1
read_name(read::Read) = ascii(reinterpret(UInt8, read.data[33:31+l_read_name(read)]))
is_paired(read::Read) = (flag(read::Read) & 1) > 0
is_proper_paired(read::Read) = (flag(read::Read) & 2) > 0
is_unmapped(read::Read) = (flag(read::Read) & 4) > 0
mate_unmapped(read::Read) = (flag(read::Read) & 8) > 0
is_reverse(read::Read) = (flag(read::Read) & 16) > 0
mate_reverse(read::Read) = (flag(read::Read) & 32) > 0
is_read1(read::Read) = (flag(read::Read) & 64) > 0
is_read2(read::Read) = (flag(read::Read) & 128) > 0
not_primary(read::Read) = (flag(read::Read) & 256) > 0
fail_qc(read::Read) = (flag(read::Read) & 512) > 0
is_duplicate(read::Read) = (flag(read::Read) & 1024) > 0
supplementary_alignment(read::Read) = (flag(read::Read) & 2048) > 0

function cigar(read::Read)
    start_pos = 33 + l_read_name(read)
    stop_pos = start_pos + n_cigar_op(read) * 4
    cigar_array = reinterpret(UInt32,read.data[start_pos:stop_pos])
    [(c & 15, c >> 4) for c in cigar_array]
end

function seq(read::Read, trans = b"=ACMGRSVTWYHKDBN")
    seq_length = l_seq(read)
    start_pos = 33 + l_read_name(read) + n_cigar_op(read) * 4
    stop_pos = start_pos + div(seq_length + 1, 2)
    seq = reinterpret(UInt8,read.data[start_pos:stop_pos])

    ret = Array{UInt8}(seq_length)
    for i = 1:seq_length
        s = seq[div(i + 1,2)]
        shift = (i % 2) * 4
        mask = 15 << shift
        c = ((s & mask) >> shift) + 1
        ret[i] = trans[c]
    end

    ascii(ret)
end

function qualities(read::Read)
    seq_length = l_seq(read)
    start_pos = 33 + l_read_name(read) + n_cigar_op(read) * 4 + div(seq_length + 1, 2)
    stop_pos = start_pos + seq_length - 1
    a = reinterpret(UInt8,read.data[start_pos:stop_pos])
    for i = 1:length(a)
        a[i] = a[i] + 33
    end

    ascii(a)
end

# reader functions
reference_name(reader::BamReader, tid::Int64) = reader.tid_to_refname[tid]
eof(reader::BamReader) = reader.position == -1

function nextread(r::BamReader)
    f = r.bamStream
    if peek(f) == -1 # eof does not work on the BAM files either in C++ or here (BAM vs. gzip issue?)
        r.done = true
        r.position = -1
        return
    end

    buf = Array(Int32, 1)
    gzread(f, pointer(buf), 4)
    block_size = buf[1]

    read = Array(Int8, block_size) #the raw bytes
    gzread(f, pointer(read), block_size)

    Read(read)
end

## here we want to update the reader
#eachposition(r::BamReader) = BamReaderIterator(r)
#immutable BamReaderIterator
#	reader::BamReader
#end
#Base.start(it::BamReaderIterator) = it.reader.position
#Base.done(it::BamReaderIterator, position) = position == -1
#function Base.next(it::BamReaderIterator, position)
#	pos = it.reader.position
#	advance!(it.reader)
#	pos,it.reader.position
#end


text:
function ueberpruefeeindeutigeLoesung3x2(A1,b1,A2,b2)

if size(A2,1) == 3 && size(A2,2) == 2 && size(b2,1) == 3 && size(b2,2) == 1
    y1 = inv(transpose(A1)*A1)*transpose(A1)*b1;
    y2 = inv(transpose(A2)*A2)*transpose(A2)*b2;
    if almostequal(y1[1], y2[1]) && almostequal(y1[2], y2[2])
        println("Super, ihr habt ein lineares Gleichungssystem aufgestellt, das sich eindeutig lösen lässt!");
        println("Eure Lösung lautet:");
        println("x1 = ",round(y2[1];digits=2));
        println("x2 = ",round(y2[2];digits=2));
    else
        println("Die Größe eures Gleichungssystems stimmt zwar, allerdings lässt es sich nicht eindeutig lösen.")
    end
else
    println("Ihr habt noch einen Fehler bei der Anzahl der Einträge in der Matrix oder dem Vektor. Stellt noch einmal sicher, dass ihr euer lineares Gleichungssystem von gerade um eine Zeile erweitert habt. Das heißt, dass sowohl der Vektor, als auch die Matrix, um eins in die Länge wachsen muss.")
end

end


text:
#=
Copyright (c) 2018-2019: David Sanders and Luis Benet.
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
=#

import IntervalArithmetic
using TaylorSeries
import TaylorIntegration
using TaylorModels: validated_step!, TaylorModel1, TaylorModelN, evaluate

import IntervalArithmetic
const IA = IntervalArithmetic
const TI = TaylorIntegration

function validated_integ(f!, qq0::AbstractVector{T}, δq0::IntervalBox{N,T},
                         t0::T, tmax::T, orderQ::Int, orderT::Int, abstol::T,
                         params=nothing; # not currently used
                         maxsteps::Int=500, parse_eqs::Bool=true,
                         check_property::Function=(t, x)->true) where {N, T<:Real}

    # Set proper parameters for jet transport
    @assert N == get_numvars()
    dof = N

    # Some variables
    zI = zero(IA.Interval{T})
    q0 = IntervalBox(qq0)
    t = t0 + Taylor1(orderT)
    tI = t0 + Taylor1(orderT+1)
    δq_norm = IntervalBox(IA.Interval{T}(-1,1), N)
    q0box = q0 .+ δq_norm

    # Allocation of vectors
    # Output
    tv = Vector{T}(undef, maxsteps+1)
    xv = Vector{IntervalBox{N,T}}(undef, maxsteps+1)
    xTM1v = Matrix{TaylorModel1{TaylorN{T}, T}}(undef, dof, maxsteps+1)

    # Internals: jet transport integration
    x = Vector{Taylor1{TaylorN{T}}}(undef, dof)
    dx = Vector{Taylor1{TaylorN{T}}}(undef, dof)
    xaux = Vector{Taylor1{TaylorN{T}}}(undef, dof)
    xTMN = Vector{TaylorModelN{N, T, T}}(undef, dof)

    # Internals: Taylor1{Interval{T}} integration
    xI = Vector{Taylor1{IA.Interval{T}}}(undef, dof)
    dxI = Vector{Taylor1{IA.Interval{T}}}(undef, dof)
    xauxI = Vector{Taylor1{IA.Interval{T}}}(undef, dof)

    # Set initial conditions
    rem = Vector{IA.Interval{T}}(undef, dof)

    @inbounds for i in eachindex(x)
        qaux = normalize_taylor(qq0[i] + TaylorN(i, order=orderQ), δq0, true)
        x[i] = Taylor1(qaux, orderT)
        dx[i] = x[i]
        xTMN[i] = TaylorModelN(qaux, zI, q0, q0box)

        xI[i] = Taylor1(q0[i] + δq0[i], orderT+1)
        dxI[i] = xI[i]
        rem[i] = zI

        xTM1v[i, 1] = TaylorModel1(deepcopy(x[i]), zI, zI, zI)
    end

    # Output vectors
    @inbounds tv[1] = t0
    @inbounds xv[1] = evaluate(xTMN, δq_norm)

    # Determine if specialized jetcoeffs! method exists (built by @taylorize)
    parse_eqs = parse_eqs && (length(methods(TI.jetcoeffs!)) > 2)
    if parse_eqs
        try
            TI.jetcoeffs!(Val(f!), t, x, dx, params)
        catch
            parse_eqs = false
        end
    end

    # Integration
    nsteps = 1
    while t0 < tmax

        # Validated step of the integration
        δt = validated_step!(f!, t, x, dx, xaux, tI, xI, dxI, xauxI,
            t0, tmax, xTMN, xv, rem, δq_norm, q0, q0box,
            nsteps, orderT, abstol, params, parse_eqs, check_property)

        # New initial conditions and time
        nsteps += 1
        t0 += δt
        @inbounds begin
            t[0] = t0
            tI[0] = t0
            tv[nsteps] = t0
            for i in eachindex(x)
                δtI = IA.Interval(0, δt)
                xTM1v[i, nsteps] = TaylorModel1(deepcopy(x[i]), rem[i], zI, δtI)
                aux = x[i](δt)
                x[i]  = Taylor1(aux, orderT)
                dx[i] = Taylor1(zero(aux), orderT)
                auxI = xTMN[i](δq_norm)
                xI[i] = Taylor1(auxI, orderT+1)
                dxI[i] = xI[i]
            end
        end

        if nsteps > maxsteps
            warn("Maximum number of integration steps $maxsteps reached; exiting. Try increasing `max_steps`.")
            break
        end

    end

    return view(tv, 1:nsteps), view(xv, 1:nsteps), view(xTM1v, :, 1:nsteps)
end


text:
using Test
using kNN
using QuadGK

@testset "Kernels" begin
    ks = kNN.kernels
    for k in keys(ks)
        kernel = ks[k]
        for u in [0.0, 0.00001, 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0]
            @test isapprox(kernel(-u), kernel(u))
        end
        integral, approxerr = quadgk(kernel, -Inf, Inf)
        @test abs(integral - 1.0) < approxerr
    end
end


text:
include("../MuZero.jl")
using .MuZero

using Logging
using TensorBoardLogger
using Dates: now, format

env = MuZero.TTTexample.env
benchmark = MuZero.TTTexample.benchmark

path = "results/tictactoe/" * format(now(),"yyyy-mm-dd-THHMM") * "/"
tblogger=TBLogger(path, min_level=Logging.Info)
@info "Run `tensorboard --logdir results` to see results in real-time"

# TRAIN
with_logger(tblogger) do
  @info "params" params=env.params
  # @info "alphalike params" αnns.f.net.hyper
  # @timed learning_step!(env)
  @info "network params" MuZero.TTTexample.simμNetworkHP
  # @info "Memory Analysis" memory_analysis(env.memory)... "generated by MinMax player"
  @info "Benchmark" run_duel(env, benchmark)...
  mutrain!(env, benchmark=benchmark, path=path)
end


text:
import REPL
using REPL.TerminalMenus
using Plots

options=["sin","cos"]
menu = RadioMenu(options)

choice=request("どの関数を描画したいですか？", menu)

if choice != -1
    plot(eval(Meta.parse(options[choice])),-π:0.01:π)
else
    println("Menu canceled")
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
@maintain_type struct TrajectoryNode{T} <: FieldVector{12,T} t::T; s::T; V::T; A::T; E::T; N::T; ψ::T; κ::T; θ::T; ϕ::T; edge_L::T; edge_R::T end
@maintain_type struct TimeInterpolants{T} <: FieldVector{3,T} s::T; V::T; A::T end
@maintain_type struct SpatialInterpolants{T} <: FieldVector{8,T} E::T; N::T; ψ::T; κ::T; θ::T; ϕ::T; edge_L::T; edge_R::T end

TrajectoryNode(t, ti::TimeInterpolants, si::SpatialInterpolants) = TrajectoryNode(t, ti.s, ti.V, ti.A, si.E, si.N, si.ψ, si.κ, si.θ, si.ϕ, si.edge_L, si.edge_R)
LocalRoadGeometry(tn::TrajectoryNode) = LocalRoadGeometry(tn.ψ, tn.κ, tn.θ, tn.ϕ)

struct TrajectoryTube{T}
    t::Vector{T}         # time (s)
    s::Vector{T}         # longitudinal distance along path (m)
    V::Vector{T}         # longitudinal velocity (m/s)
    A::Vector{T}         # longitudinal acceleration (m/s^2)
    E::Vector{T}         # E coordinate (m)
    N::Vector{T}         # N coordinate (m)
    ψ::Vector{T}         # heading (rad)
    κ::Vector{T}         # curvature (1/m)
    θ::Vector{T}         # grade (rad)
    ϕ::Vector{T}         # bank (rad)
    edge_L::Vector{T}    # left lateral deviation bound (m)
    edge_R::Vector{T}    # right lateral deviation bound (m)

    interp_by_s::Extrapolation{SpatialInterpolants{T},1,
                               GriddedInterpolation{SpatialInterpolants{T},1,SpatialInterpolants{T},Gridded{Linear},Tuple{Vector{T}}},
                               Gridded{Linear},Line{Nothing}}

    function TrajectoryTube{T}(t::Vector{T}, s::Vector{T}, V::Vector{T}, A::Vector{T},
                               E::Vector{T}, N::Vector{T}, ψ::Vector{T}, κ::Vector{T},
                               θ::Vector{T}, ϕ::Vector{T}, edge_L::Vector{T}, edge_R::Vector{T}) where {T}
        @assert length(t) == length(s) == length(V) == length(A) ==
                length(E) == length(N) == length(ψ) == length(κ) ==
                length(θ) == length(ϕ) == length(edge_L) == length(edge_R)
        interp_by_s = extrapolate(
            interpolate((s,), vec(reinterpret(SpatialInterpolants{T}, [E N ψ κ θ ϕ edge_L edge_R]')), Gridded(Linear())),
            Line()
        )
        new(t, s, V, A, E, N, ψ, κ, θ, ϕ, edge_L, edge_R, interp_by_s)
    end
end
function TrajectoryTube(t::Vector{T}, s::Vector{T}, V::Vector{T}, A::Vector{T},
                        E::Vector{T}, N::Vector{T}, ψ::Vector{T}, κ::Vector{T},
                        θ::Vector{T}=zero(t), ϕ::Vector{T}=zero(t),
                        edge_L::Vector{T}=fill(T(4), size(t)), edge_R::Vector{T}=fill(-T(4), size(t))) where {T}
    TrajectoryTube{T}(t, s, V, A, E, N, ψ, κ, θ, ϕ, edge_L, edge_R)
end
Base.length(traj::TrajectoryTube) = length(traj.t)
DifferentialDynamicsModels.duration(traj::TrajectoryTube) = traj.t[end] - traj.t[1]
function (traj::TrajectoryTube)(t)
    i = clamp(searchsortedfirst(traj.t, t, 1, length(traj), Base.Order.ForwardOrdering()) - 1, 1, length(traj)-1)
    A = (traj.V[i+1] - traj.V[i])/(traj.t[i+1] - traj.t[i])    # potentially different from traj.A[i]
    dt = t - traj.t[i]
    ti = TimeInterpolants(traj.s[i] + traj.V[i]*dt + A*dt*dt/2, traj.V[i] + A*dt, A)
    si = traj.interp_by_s(ti.s)
    TrajectoryNode(t, ti, si)
end
function Base.getindex(traj::TrajectoryTube, s)
    i = clamp(searchsortedfirst(traj.s, s, 1, length(traj), Base.Order.ForwardOrdering()) - 1, 1, length(traj)-1)
    A = (traj.V[i+1] - traj.V[i])/(traj.t[i+1] - traj.t[i])    # potentially different from traj.A[i]
    ds = s - traj.s[i]
    if abs(A) < 1e-3 || s > traj.s[end]
        dt = ds/traj.V[i]
    else
        dt = (sqrt(2*A*ds + traj.V[i]^2) - traj.V[i])/A
    end
    t = traj.t[i] + dt
    ti = TimeInterpolants(s, traj.V[i] + A*dt, A)
    si = traj.interp_by_s(s)
    TrajectoryNode(t, ti, si)
end

Base.getindex(traj::TrajectoryTube, i::Integer) = SVector(traj.E[i], traj.N[i])    # this overloading is a bit hacky/brittle
function path_coordinates(traj::TrajectoryTube, x)
    d2min, imin = eltype(x)(Inf), 0
    for i in 1:length(traj)-1
        d2 = distance2((traj[i], traj[i+1]), x)
        if d2 < d2min
            d2min, imin = d2, i
        end
    end
    i = imin
    v = traj[i+1] - traj[i]
    w = x - traj[i]
    ds = sqrt(w⋅w - d2min)
    s = traj.s[i] + ds
    e = sqrt(d2min)*sign(cross(v, w))
    A = (traj.V[i+1] - traj.V[i])/(traj.t[i+1] - traj.t[i])    # potentially different from traj.A[i]
    if abs(A) < 1e-3
        dt = ds/traj.V[i]
    else
        dt = (sqrt(2*A*ds + traj.V[i]^2) - traj.V[i])/A
    end
    t = traj.t[i] + dt
    s, e, t
end
path_coordinates(traj::TrajectoryTube, q::BicycleState) = path_coordinates(traj, SVector(q.E, q.N))

function straight_trajectory(len, vel)
    TrajectoryTube([0., len/vel], # t
                   [0., len],     # s
                   [vel, vel],    # V
                   [0., 0.],      # A
                   [0., 0.],      # E
                   [0., len],     # N
                   [0., 0.],      # ψ
                   [0., 0.])      # κ
end

max_stars_repo_name: Michelle-NYX/CollisionAvoidance

max_issues_repo_name: Michelle-NYX/CollisionAvoidance

max_forks_repo_name: Michelle-NYX/CollisionAvoidance