text:
"""
    latlong_to_xy(latitude, longitude, utm_zone; isnorth=true, units="m")

Converts arrays of points from latitude and longitude to x and y in meters in a local
coordinate frame based on the point with the lowest magnitude latitude,

# Arguments
- `latitude::Array{Float,N}`
- `longitude::Array{Float,N}`
- `isnorth::Float`: specifies if the point is in the northern hemisphere (defaul: true)
"""
function latlong_to_xy(latitude, longitude, utm_zone; isnorth=true)

    # get number of points
    npoints = length(latitude)

    # set up utm region
    utmregion = gd.UTMfromLLA(utm_zone, isnorth, wgs84)

    # get zero point index
    zp = argmin(latitude)

    # get zero point lat long
    minlla = gd.LLA(latitude[zp], longitude[zp])

    # get zero point utm
    minutm = utmregion(minlla)

    # get x and y for all points
    x = zeros(npoints)
    y = zeros(npoints)
    for i in 1:npoints
        # lat and long for current point
        lla = gd.LLA(latitude[i], longitude[i])

        # convert to utm coordinates
        utm = utmregion(lla)

        # get x location
        x[i] = utm.x - minutm.x

        # get y location
        y[i] = utm.y - minutm.y

    end

    return x, y
end

"""
    hermite_spline(x, x0, x1, y0, dy0, y1, dy1)

Produces the y and (optionally) dy values for a hermite cubic spline
interpolating between two end points with known slopes

# Arguments
- `x::Float`: x position of output y
- `x0::Float`: x position of upwind endpoint of spline
- `x1::Float`: x position of downwind endpoint of spline
- `y0::Float`: y position of upwind endpoint of spline
- `dy0::Float`: slope at upwind endpoint of spline
- `y1::Float`: y position of downwind endpoint of spline
- `dy1::Float`: slope at downwind endpoint of spline
"""
function hermite_spline(x, x0, x1, y0, dy0, y1, dy1; return_deriv=false)

    # initialize coefficients for parametric Hermite cubic spline
    c3 = (2.0*(y1))/(x0^3 - 3.0*x0^2*x1 + 3.0*x0*x1^2 - x1^3) -
        (2.0*(y0))/(x0^3 - 3.0*x0^2*x1 + 3.0*x0*x1^2 - x1^3) +
        (dy0)/(x0^2 - 2.0*x0*x1 + x1^2) +
        (dy1)/(x0^2 - 2.0*x0*x1 + x1^2)

    c2 = (3.0*(y0)*(x0 + x1))/(x0^3 - 3.0*x0^2*x1 + 3.0*x0*x1^2 - x1^3) -
        ((dy1)*(2.0*x0 + x1))/(x0^2 - 2.0*x0*x1 + x1^2) - ((dy0)*(x0 +
        2.0*x1))/(x0^2 - 2.0*x0*x1 + x1^2) - (3.0*(y1)*(x0 + x1))/(x0^3 -
        3.0*x0^2*x1 + 3.0*x0*x1^2 - x1^3)

    c1 = ((dy0)*(x1^2 + 2.0*x0*x1))/(x0^2 - 2.0*x0*x1 + x1^2) + ((dy1)*(x0^2 +
        2.0*x1*x0))/(x0^2 - 2.0*x0*x1 + x1^2) - (6.0*x0*x1*(y0))/(x0^3 -
        3.0*x0^2*x1 + 3.0*x0*x1^2 - x1^3) + (6.0*x0*x1*(y1))/(x0^3 -
        3.0*x0^2*x1 + 3.0*x0*x1^2 - x1^3)

    c0 = ((y0)*(- x1^3 + 3.0*x0*x1^2))/(x0^3 - 3.0*x0^2*x1 + 3.0*x0*x1^2 -
        x1^3) - ((y1)*(- x0^3 + 3.0*x1*x0^2))/(x0^3 - 3.0*x0^2*x1 +
        3.0*x0*x1^2 - x1^3) - (x0*x1^2*(dy0))/(x0^2 - 2.0*x0*x1 + x1^2) -
        (x0^2*x1*(dy1))/(x0^2 - 2.0*x0*x1 + x1^2)

    # Solve for y and dy values at the given point
    y = c3*x^3 + c2*x^2 + c1*x + c0
    dy_dx = c3*3*x^2 + c2*2*x + c1

    if return_deriv
        return y, dy_dx
    else
        return y
    end
end

"""
    overlap_area_func(turbine_y, turbine_z, rotor_diameter, wake_center_y,
    wake_center_z, wake_diameter; tol=1E-6)

Produces the y and (optionally) dy values for a hermite cubic spline
interpolating between two end points with known slopes

All dimensions should be in meters

# Arguments
- `turbine_y::Float`: cross wind location of turbine hub
- `turbine_z::Float`: vertical location of turbine hub
- `rotor_diameter::Float`
- `wake_center_y::Float`: cross wind location of wake center
- `wake_center_z::Float`: vertical location of wake center
- `wake_diameter::Float`
- `tol::Float`: extra distance for comparisons. Default 1E-6
"""
# calculates the overlap area between a given wake and a rotor area
function overlap_area_func(turbine_y, turbine_z, rotor_diameter, wake_center_y,
    wake_center_z, wake_diameter; tol=1E-6)

    # distance between wake center and rotor center
    if (wake_center_z > (turbine_z + tol)) || (wake_center_z < (turbine_z - tol))
        OVdYd = sqrt((wake_center_y-turbine_y)^2 + (wake_center_z - turbine_z)^2)
    elseif (wake_center_y > (turbine_y + tol))
        OVdYd = wake_center_y - turbine_y
    elseif (turbine_y > (wake_center_y + tol))
        OVdYd = turbine_y - wake_center_y
    else
        OVdYd = 0.0
    end

    # find rotor radius
    OVr = rotor_diameter/2.0

    # find wake radius
    OVRR = wake_diameter/2.0

    # make sure the distance from wake center to turbine hub is positive
    OVdYd = abs(OVdYd)

    # determine if there is overlap
    if (OVdYd < (OVr+OVRR)) # if the rotor overlaps the wake zone

        # check that turbine and wake centers are not perfectly aligned
        if (OVdYd > (0.0 + tol))

            # check if the rotor is wholly contained in the wake
            if ((OVdYd + OVr) < OVRR + tol)
                # wake_overlap = pi*OVr*OVr
                wake_overlap = 3.1415926535897*OVr*OVr
            elseif ((OVdYd + OVRR) < OVr + tol)
                # wake_overlap = pi*OVRR*OVRR
                wake_overlap = 3.1415926535897*OVRR*OVRR
            else
                # calculate the distance from the wake center to the chord connecting the lens cusps
                OVL = (-OVr*OVr+OVRR*OVRR+OVdYd*OVdYd)/(2.0*OVdYd)

                OVz = sqrt(OVRR*OVRR-OVL*OVL)
                OVz2 = sqrt(OVr*OVr-(OVdYd-OVL)*(OVdYd-OVL))

                wake_overlap = OVRR*OVRR*acos(OVL/OVRR) + OVr*OVr*acos((OVdYd-OVL)/OVr) - OVL*OVz - (OVdYd-OVL)*OVz2
            end

        # perfect overlap case where the wake is larger than the rotor
        elseif (OVRR > OVr)
            wake_overlap = pi*OVr*OVr
        # perfect overlap case where the rotor is larger than the wake
        else
            wake_overlap = pi*OVRR*OVRR
        end

    # case with no overlap
    else
        wake_overlap = 0.0
    end

end

"""
    smooth_max_ndim(x; s=100.0)

Calculate the smoothmax (a.k.a. softmax or LogSumExponential) of the elements in x.

Based on John D. Cook's writings at
(1) https://www.johndcook.com/blog/2010/01/13/soft-maximum/
and
(2) https://www.johndcook.com/blog/2010/01/20/how-to-compute-the-soft-maximum/

# Arguments
- `x::Float`: first value for comparison
- `y::Float`: second value for comparison
- `s::Float` : controls the level of smoothing used in the smooth max
"""
function smooth_max(x, y; s=10.0)

    # LogSumExponential Method - used this in the past
    # g = (x*exp(s*x)+y*exp(s*y))/(exp(s*x)+exp(s*y))

    # non-overflowing version of Smooth Max function (see ref 2 above)
    max_val = max(x, y)
    min_val = min(x, y)
    r = (log(1.0 + exp(s*(min_val - max_val))) + s*max_val)/s

    return r

end

"""
    smooth_max(x; s=10.0)

Calculate the smoothmax (a.k.a. softmax or LogSumExponential) of the elements in x.

Based on John D. Cook's writings at
(1) https://www.johndcook.com/blog/2010/01/13/soft-maximum/
and
(2) https://www.johndcook.com/blog/2010/01/20/how-to-compute-the-soft-maximum/

And based on article in FeedlyBlog
(3) https://blog.feedly.com/tricks-of-the-trade-logsumexp/

# Arguments
- `x::Array{Float,1}` : vector with all the input values
- `s::Float` : controls the level of smoothing used in the smooth max
"""
function smooth_max(x; s=10.0)

    # non-overflowing version of Smooth Max function (see ref 2 and 3 above)

    # get the maximum value and the index of maximum value
    max_val, max_ind = findmax(x)

    # get the indices of x
    indices = collect(1:length(x))

    # remove the index of the maximum value
    splice!(indices, max_ind)

    # LogSumExp with smoothing factor s
    r = (log(sum([1.0; exp.(s*(x[indices].-max_val))])) + s*max_val)/s

    return r

end

"""
    closeBndryLists(bndryPts_x, bndryPts_y)

Appends the 1st element to the end of each array for a closed boundary.
Note, this will not function properly if there is only one region.
For only one region, use `closeBndryList(bndryPts_x, bndryPts_y)` (note the
singular, not plural 'List' in the function title)

# Arguments
- `bndryPts_x::Array{Float,1}` : N-D array of x-coordinates for the vertices
        around N-many closed boundaries
- `bndryPts_y::Array{Float,1}` : N-D array of y-coordinates for the vertices
        around N-many closed boundaries
"""
function closeBndryLists(region_bndry_x, region_bndry_y)
    # Determine how many regions and points per region were passed
    nRegions = length(region_bndry_x)
    # For every region
    for i in 1:nRegions
        # Append the initial points to the end of that row (if needed)
        region_bndry_x[i], region_bndry_y[i] = closeBndryList(region_bndry_x[i], region_bndry_y[i])
    end

    return region_bndry_x, region_bndry_y
end

"""
    closeBndryList(bndryPts_x, bndryPts_y)

Appends the 1st element to the end of the coordinate arrays if it is not already
repeated. Note, this will only work on 1-D arrays. For an array of 1-D arrays,
use `closeBndryLists(bndryPts_x, bndryPts_y)` (note the
plural, not singular 'Lists' in the function title)

# Arguments
- `bndryPts_x::Array{Float,1}` : 1-D array of x-coordinates for the vertices
        around a singlar closed boundary
- `bndryPts_y::Array{Float,1}` : 1-D array of y-coordinates for the vertices
        around a singlar closed boundary
"""
function closeBndryList(bndry_x, bndry_y)
    # If the end point isn't already a repeat of the staring point
    if !((bndry_x[1] == bndry_x[end]) && (bndry_y[1] == bndry_y[end]))
        # Append the first point to the end of the array
        bndryPts_x_clsd = push!(bndry_x, bndry_x[1])
        bndryPts_y_clsd = push!(bndry_y, bndry_y[1])
    end

    return bndryPts_x_clsd, bndryPts_y_clsd
end

"""
    PointsOnCircum(center_x, center_y, r, n = 100)

Given a circle center, radius, and number of discrete points, returns an array
of discrete points along the circle's circumference

# Arguments
- `center_x::Float64` : cartesian x-coordinate for the center of the circle
- `center_y::Float64` : cartesian y-coordinate for the center of the circle
- `r::Float64` : distance from circle's center to the circumference points
- `n::Float64` : defaults to 100, is the number of discrete evenly-spaced points
        that will be returned along the circle's circumference
"""
function DiscreteCircum(center_x, center_y, r, n = 100)
    bndry_x = zeros(n+1)
    bndry_y = zeros(n+1)
    for i in 0:n
        bndry_x[i+1] = center_x + cos(2*pi/n*i)*r
        bndry_y[i+1] = center_y + sin(2*pi/n*i)*r
    end
    return bndry_x, bndry_y
end

"""
    calcMinorAngle(bndry_x, bndry_y, bndry_z=[0,0,0])

Given three points in space, calculates the magnitude of the non-reflex angle
formed at the center point. Created to be used in VR_bounary_startup()

# Arguments
- `bndry_x::Array{Float,1}` : 1-D array of x-coordinates for the vertices
        around a singlar closed boundary
- `bndry_y::Array{Float,1}` : 1-D array of y-coordinates for the vertices
        around a singlar closed boundary
- `bndry_z::Array{Float,1}` : 1-D array of z-coordinates for the vertices
        around a singlar closed boundary. Default to [0,0,0] for X-Y plane
"""
function calcMinorAngle(bndry_x, bndry_y, bndry_z=[0,0,0])
    # Calculates the magnitude of the non-reflex angle formed at the center point
    ABx = bndry_x[2]-bndry_x[1]
    ABy = bndry_y[2]-bndry_y[1]
    ABz = bndry_z[2]-bndry_z[1]

    BCx = bndry_x[2]-bndry_x[3]
    BCy = bndry_y[2]-bndry_y[3]
    BCz = bndry_z[2]-bndry_z[3]

    Num = (ABx*BCx) + (ABy*BCy) + (ABz*BCz) # Dot Product

    Denom = sqrt(ABx^2 + ABy^2 + ABz^2) * sqrt(BCx^2 + BCy^2 + BCz^2) # Multiplication of magnitudes
    Theta = acosd(Num/Denom) # Get the angle formed

    # If it's greater than 180, get the
    if (Theta > 180.0)
        Theta = 360.0 - Theta    # Get the explementary angle
    end

    return Theta
end

"""
    calcSmallestAngle(bndry_x_clsd, bndry_y_clsd)

Given a 1-D closed array of boundary verticies (with first point repeated at the
end) it determines the smallest non-reflex angle created by any three
consecutive verticies along the boundary. Created to be used in
VR_bounary_startup()

# Arguments
- `bndry_x::Array{Float,1}` : 1-D array of x-coordinates for the vertices
        around a singlar closed boundary
- `bndry_y::Array{Float,1}` : 1-D array of y-coordinates for the vertices
        around a singlar closed boundary
"""
function calcSmallestAngle(bndry_x_clsd, bndry_y_clsd)
    num_angles = length(bndry_x_clsd)-1
    # Loop the second point back on so we get the angle at the beginning
    if !(((bndry_x_clsd[1] == bndry_x_clsd[end-1]) && (bndry_y_clsd[1] == bndry_y_clsd[end-1])) && ((bndry_x_clsd[2] == bndry_x_clsd[end]) && (bndry_y_clsd[2] == bndry_y_clsd[end])))
        bndryPts_x_loopd = vcat(bndry_x_clsd, bndry_x_clsd[2])
        bndryPts_y_loopd = vcat(bndry_y_clsd, bndry_y_clsd[2])
    end

    #- Calculate the smallest angle -#
    smallest_angle = 360
    for i in 1:num_angles
        temp_angle = calcMinorAngle(bndryPts_x_loopd[i:i+2], bndryPts_y_loopd[i:i+2])

        if (temp_angle < smallest_angle)
            smallest_angle = temp_angle
        end
    end

    return smallest_angle
end

"""
    getPerimeterLength(bndry_x_clsd, bndry_y_clsd)

Given a 1-D closed array of boundary verticies (with first point repeated at the
end) returns the length along the perimeter. Created to be used in
VR_bounary_startup()

# Arguments
- `bndry_x::Array{Float,1}` : 1-D array of x-coordinates for the vertices
        around a singlar closed boundary
- `bndry_y::Array{Float,1}` : 1-D array of y-coordinates for the vertices
        around a singlar closed boundary
"""
function getPerimeterLength(bndry_x_clsd, bndry_y_clsd)
    num_bndry_pts = length(bndry_x_clsd)-1
    nLength = zeros(num_bndry_pts)
    for i in 1:num_bndry_pts
        nLength[i] = sqrt(
                 (bndry_x_clsd[i+1]-bndry_x_clsd[i])^2
                +(bndry_y_clsd[i+1]-bndry_y_clsd[i])^2)
    end

    return nLength
end

"""
    coordDist(x1, y1, x2, y2)

Given a two points (x1, y1) and (x2, y2), returns the euclidean distance between
them

# Arguments
- `x1::Float64` : x-coord of the first point
- `y1::Float64` : y-coord of the first point
- `x2::Float64` : x-coord of the second point
- `y2::Float64` : y-coord of the second point
"""
function coordDist(x1, y1, x2, y2)
    xDiff = x1 - x2
    yDiff = y1 - y2
    return sqrt(xDiff^2 + yDiff^2)
end


text:
# Text provided under a Creative Commons Attribution license, CC-BY.  All code is made available under the FSF-approved BSD-3 license.  (c) Lorena A. Barba, Gilbert F. Forsyth 2017. Thanks to NSF for support via CAREER award #1149784.
# [@LorenaABarba](https://twitter.com/LorenaABarba)
# Translation to Julia by [@miguelraz_](https://twitter.com/miguelraz_)

## 12 steps to Navier–Stokes
# ======
# ***

# You should have completed Steps [1](./01_Step_1.jl) and [2](./02_Step_2.jl) before continuing.
# This Jupyter notebook continues the presentation of the **12 steps to Navier–Stokes**,
# the practical module taught in the interactive CFD class of [Prof. Lorena Barba](http://lorenabarba.com).

## Step 3: Diffusion Equation in 1-D
# -----
# ***

# The one-dimensional diffusion equation is:

# $$\frac{\partial u}{\partial t}= \nu \frac{\partial^2 u}{\partial x^2}$$

# The first thing you should notice is that —unlike the previous two simple equations we have studied— this equation has a second-order derivative. We first need to learn what to do with it!

# ### Discretizing $\frac{\partial ^2 u}{\partial x^2}$

# The second-order derivative can be represented geometrically as the line tangent to the curve given by the first derivative.  We will discretize the second-order derivative with a Central Difference scheme: a combination of Forward Difference and Backward Difference of the first derivative.  Consider the Taylor expansion of $u_{i+1}$ and $u_{i-1}$ around $u_i$:

# $u_{i+1} = u_i + \Delta x \frac{\partial u}{\partial x}\bigg|_i + \frac{\Delta x^2}{2} \frac{\partial ^2 u}{\partial x^2}\bigg|_i + \frac{\Delta x^3}{3!} \frac{\partial ^3 u}{\partial x^3}\bigg|_i + O(\Delta x^4)$

# $u_{i-1} = u_i - \Delta x \frac{\partial u}{\partial x}\bigg|_i + \frac{\Delta x^2}{2} \frac{\partial ^2 u}{\partial x^2}\bigg|_i - \frac{\Delta x^3}{3!} \frac{\partial ^3 u}{\partial x^3}\bigg|_i + O(\Delta x^4)$

# If we add these two expansions, you can see that the odd-numbered derivative terms will cancel each other out.  If we neglect any terms of $O(\Delta x^4)$ or higher (and really, those are very small), then we can rearrange the sum of these two expansions to solve for our second-derivative.

# $u_{i+1} + u_{i-1} = 2u_i+\Delta x^2 \frac{\partial ^2 u}{\partial x^2}\bigg|_i + O(\Delta x^4)$

# Then rearrange to solve for $\frac{\partial ^2 u}{\partial x^2}\bigg|_i$ and the result is:

# $$\frac{\partial ^2 u}{\partial x^2}=\frac{u_{i+1}-2u_{i}+u_{i-1}}{\Delta x^2} + O(\Delta x^2)$$

# ### Back to Step 3

# We can now write the discretized version of the diffusion equation in 1D:

# $$\frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t}=\nu\frac{u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}}{\Delta x^2}$$

# As before, we notice that once we have an initial condition, the only unknown is $u_{i}^{n+1}$, so we re-arrange the equation solving for our unknown:

# $$u_{i}^{n+1}=u_{i}^{n}+\frac{\nu\Delta t}{\Delta x^2}(u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n})$$

# The above discrete equation allows us to write a program to advance a solution in time. But we need an initial condition. Let's continue using our favorite: the hat function. So, at $t=0$, $u=2$ in the interval $0.5\le x\le 1$ and $u=1$ everywhere else. We are ready to number-crunch!

nₓ = 41 # n\_<TAB>x
Δx = 2 / (nₓ - 1)
nₜ = 20    #nₜ is the number of timesteps we want to calculate
nᵤ = 0.3
σ = 0.2
Δt = σ * Δx^2 / nᵤ  #Δt is the amount of time each timestep covers (delta t)

u = ones(nₓ)      #as before, we initialize u with every value equal to 1.
indexes = Int.((.5/Δx):(1/Δx))
u[indexes] .= 2 #then set u = 2 between 0.5 and 1 as per our I.C.s

uₙ = ones(nₓ) #initialize our placeholder array un, to hold the time-stepped solution

u = numpy.ones(nx)      #a numpy array with nx elements all equal to 1.
u[int(.5 / dx):int(1 / dx + 1)] = 2  #setting u = 2 between 0.5 and 1 as per our I.C.s

for n in nₜ  #iterate through time
    uₙ = copy(u) ##copy the existing values of u into uₙ
    for i in 1:nₓ  ##now we'll iterate through the u array

     ###This is the line from Step 1, copied exactly.  Edit it for our new equation.
     ###then uncomment it and run the cell to evaluate Step 2
           ###u[i] = uₙ[i] - c * dt / dx * (uₙ[i] - uₙ[i-1])
    u[i] = uₙ[i] + nu * Δt / Δx^2 * (uₙ[i+1] - 2 * uₙ[i] + uₙ[i-1])
    end
end

plot(range(start = 0, stop = 2, length = nₓ), u);

## Learn More

# For a careful walk-through of the discretization of the diffusion equation with finite differences (and all steps from 1 to 4), watch **Video Lesson 4** by Prof. Barba on YouTube.

# TODO
# YouTubeVideo('y2WaK7_iMRI')

ext: jl

max_stars_repo_name: mkitti/CFDJulia

max_issues_repo_name: mkitti/CFDJulia

max_forks_repo_name: mkitti/CFDJulia


text:
using ITensors
using ITensorNetworks
using ITensorsVisualization

using ITensorNetworks: inds_network, project_boundary, truncation_projectors
using ITensorNetworks: Models

model = Models.Model"ising"()
βc = Models.critical_point(model)
β = 1.001 * βc
@show β / βc

N = (4, 4)
ndims = length(N)

space = [QN(0) => 2]
#space = 2
tn_inds = inds_network(N...; linkdims=space, addtags="S=1/2")
A = Models.local_boltzmann_weight("ising", Val(ndims); β)

tn = itensor.((A,), tn_inds)

# Project periodic boundary indices
# onto the 1 state
state = 1
tn = project_boundary(tn, state)

cutoff = 1e-15
_maxdim = 100

tn_split, U, Ud = truncation_projectors(tn; cutoff=cutoff, maxdim=_maxdim)

tn_projected = [tn_split, U, Ud]
tn_projected_flat = vcat(vec.(tn_projected)...)

@show noncommoninds(tn_projected_flat...)

Z_approx = @visualize contract(tn_projected_flat...)
Z_exact = contract(vec(tn)...)

@disable_warn_order begin
  @show Z_approx[] / Z_exact[]
end


text:
module day14
using Test
using Debugger

function readinput(filename)
    input = readlines(filename)
    inputs = []
    mask_re = r"^mask = (.{31})$"
    mem_re = r"^mem\[([0-9]*)] = ([0-9]*)$"
    for i in input
        i[1] == "#" && continue
        if occursin("mask", i)
            maskliteral = i[8:end]
            mask = mask_from_literal(maskliteral)
            push!(inputs, mask)
        end
        m = match(mem_re, i)
        if !(m isa Nothing)
            push!(inputs, [parse(Int, m[1]), parse(Int, m[2])])
        end
    end
    inputs
end

mask_from_literal(x) = [(idx, c == 'X' ? -1 : parse(Int, c)) for (idx, c) in enumerate(x)]
all_perms(n) = Iterators.product(ntuple(_ -> [0,1], n)...)
int2mask(x) = digits(x, base=2, pad=36) |> reverse
bitarray2int(x) = sum(2^(i-1)*σ for (i, σ) in enumerate(reverse(x)))
mask2int(x) = [i[2] for i in x] |> bitarray2int

function maskit(x, mask)
    b = int2mask(x)
    for (idx, val) in mask
        b[idx] = val
    end
    bitarray2int(b)
end

function masks_from_mem_addr(mem_addr, current_mask)
    cm = deepcopy(current_mask)
    mem_addr_to_mask = int2mask(mem_addr)
    mem_addr_to_mask = zip(1:length(mem_addr_to_mask), mem_addr_to_mask) |> collect
    floating_bits = [(idx, bit) for (idx, bit) in cm if bit == -1]
    n_bits = length(floating_bits)
    perms = all_perms(n_bits)
    masks = []
    for perm in perms
        p = collect(perm)
        change = [(j[1], p[i]) for (i, j) in enumerate(floating_bits)]
        m = deepcopy(mem_addr_to_mask)
        for (i, v) in change
            m[i] = (i, v)
        end
        push!(masks, m)
    end
    masks
end

@testset "masks_from_mem_addr" begin
    mem_addr = 42
    mask_literal = "000000000000000000000000000000X1001X"
    current_mask = mask_from_literal(mask_literal)

    masks = masks_from_mem_addr(mem_addr, current_mask)
    @show masks
end

function solve(filename)
    input = readinput(filename)
    mem = Dict{Int64,Int64}()
    current_mask = []
    for i in input
        if isa(i[1], Int)
            mem_addr = i[1]
            val = i[2]
#            mem_addr_mask_literal = int2mask(mem_addr)
#            mem_addr_mask = zip(1:length(mem_addr_mask_literal), mem_addr_mask_literal) |> collect
            masks = masks_from_mem_addr(mem_addr, current_mask)
            @show mem_addr, current_mask
            for mask in masks
                @show mask
                new_mem_addr = mask2int(mask) #maskit(mem_addr, mask)
                @show new_mem_addr
#                @show mem_addr, new_mem_addr
            end
            #mem[mem_addr] = maskit(i[2], current_mask)
        else
            current_mask = i
            continue
        end
    end
    sum(values(mem))
end

@testset "solve" begin
    arg = "input/day14-test2.txt"
    solve(arg)
end

@testset "int2mask" begin
   @test int2mask(11)[end-4:end] == [0, 1, 0, 1, 1]
end

@testset "readinput" begin
    arg = "input/day14-test2.txt"
    res = readinput(arg)
#    @test res[1][end] == (36, -1)
#    @test res[2] == [42, 100]
#    @test res[4] == [26, 1]
end

end


text:
using StochasticDiffEq, Test, Random
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
using DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_linear
Random.seed!(100)
prob = prob_sde_linear

integrator = init(prob,EM(),dt=1//2^(4),tstops = [0.33])

for (u,t) in tuples(integrator)
  @show u,t
end

sol = solve(prob,EM(),dt=1//2^(4),tstops = [0.33])

@test 0.33 ∈ sol.t

sol = solve(prob,EM(),tstops = [0.33,0.80,1.0])

@test sol.t == [0.0,0.33,0.80,1.0]

sol = solve(prob,SRIW1(),tstops = [0.33])

@test 0.33 ∈ sol.t

# check reverse time and negative start times
for (i, tdir) in enumerate([-1.; 1.])
  @info i
  prob2 = remake(prob_sde_linear, tspan=(tdir*1.0, 0.0))
  integrator = init(prob2,SRIW1())
  tstops = tdir .* [0,0.33,0.80,1]
  for tstop in tstops
    add_tstop!(integrator, tstop)
  end
  solve!(integrator)
  for tstop in tstops
    @test tstop ∈ integrator.sol.t
  end
end


text:
include("../src/Model.jl")
include("../src/Regression.jl")


text:
using Flux, Test, CUDA
using Zygote: pullback

@info "Testing GPU Support"
CUDA.allowscalar(false)

function gpu_gradtest(f, args...)
  args_gpu = gpu.(args)

  l_cpu, back_cpu = pullback((x...) -> f(x...), args...)
  g_cpu = back_cpu(1f0)[1]

  l_gpu, back_gpu = pullback((x...) -> f(x...), args_gpu...)
  g_gpu = back_gpu(1f0)[1]

  @test l_cpu ≈ l_gpu   rtol=1e-4 atol=1e-4
  @test g_gpu isa CuArray
  @test g_cpu ≈ collect(g_gpu)   rtol=1e-4 atol=1e-4
end

include("cuda.jl")
include("losses.jl")
include("layers.jl")

if CUDA.has_cudnn()
  @info "Testing Flux/CUDNN"
  include("cudnn.jl")
  include("curnn.jl")
else
  @warn "CUDNN unavailable, not testing GPU DNN support"
end


text:
# Check if two arrays of edges are the same, regardless of directions
function isSame(primsEdges, primsMST)
    for edge in primsEdges
        if(!(in(Edge(edge.dst, edge.src), primsMST)) && !(in(edge, primsMST)))
            return false
        end
    end
    true
end

# Testing PrimsI

function testPrimsI(vertices, connections, attempts)
	sames = 0
	for i in 1 : attempts
		g = randomWeightedGraph(vertices, connections)
		baseLibraryPrims = LightGraphs.prim_mst(g);
		primI = parallelPrims(g, vertices)
		if(isSame(baseLibraryPrims, primI) && isSame(primI, baseLibraryPrims))
			sames += 1
		end
	end
	if (sames == attempts)
		println("Correct")
	else
		println("Rip")
	end
end

# Testing PrimsII
function testPrimsII(vertices, connections, numberOfThreads, attempts)
	sames = 0
	for i in 1 : attempts
		g = randomWeightedGraph(vertices, connections)
		baseLibraryPrims = LightGraphs.prim_mst(g);
		primII = prims(g, numberOfThreads)
		if(isSame(baseLibraryPrims, primII) && isSame(primII, baseLibraryPrims))
			sames += 1
		end
	end
	if (sames == attempts)
		println("Correct")
	else
		println("Rip")
	end
end

# Testing A*
function testAStar(graphWidth, attempts)
	sames = 0
	for i in 1 : attempts
		startVertex = rand(1 : graphWidth * graphWidth)
		endVertex = rand(1 : graphWidth * graphWidth)
		graph = makeMaze(makeGrid(graphWidth));
		baseLibraryAStar = LightGraphs.a_star(graph, startVertex, endVertex)
		ourAStar = astar(graph, graphWidth, startVertex, endVertex)
		if(isSame(baseLibraryAStar, ourAStar) && isSame(ourAStar, baseLibraryAStar))
			sames += 1
		end
	end

	if (sames == attempts)
		println("Correct")
	else
		println("Rip")
	end
end


text:
module RegistryTests

using Pkg, UUIDs, LibGit2, Test
using Pkg: depots1
using Pkg.REPLMode: pkgstr
using Pkg.Types: PkgError

include("utils.jl")

function setup_test_registries(dir = pwd())
    # Set up two registries with the same name, with different uuid
    pkg_uuids = ["c5f1542f-b8aa-45da-ab42-05303d706c66", "d7897d3a-8e65-4b65-bdc8-28ce4e859565"]
    reg_uuids = ["e9fceed0-5623-4384-aff0-6db4c442647a", "a8e078ad-b4bd-4e09-a52f-c464826eef9d"]
    for i in 1:2
        regpath = joinpath(dir, "RegistryFoo$(i)")
        mkpath(joinpath(regpath, "Example"))
        write(joinpath(regpath, "Registry.toml"), """
            name = "RegistryFoo"
            uuid = "$(reg_uuids[i])"
            repo = "https://github.com"
            [packages]
            $(pkg_uuids[i]) = { name = "Example$(i)", path = "Example" }
            """)
        write(joinpath(regpath, "Example", "Package.toml"), """
            name = "Example$(i)"
            uuid = "$(pkg_uuids[i])"
            repo = "https://github.com/JuliaLang/Example.jl.git"
            """)
        write(joinpath(regpath, "Example", "Versions.toml"), """
            ["0.5.1"]
            git-tree-sha1 = "8eb7b4d4ca487caade9ba3e85932e28ce6d6e1f8"
            """)
        write(joinpath(regpath, "Example", "Deps.toml"), """
            ["0.5"]
            julia = "0.6-1.0"
            """)
        write(joinpath(regpath, "Example", "Compat.toml"), """
            ["0.5"]
            julia = "0.6-1.0"
            """)
        LibGit2.with(LibGit2.init(regpath)) do repo
            LibGit2.add!(repo, "*")
            LibGit2.commit(repo, "initial commit"; author=TEST_SIG, committer=TEST_SIG)
        end
    end
end

function test_installed(registries)
    @test setdiff(
        UUID[r.uuid for r in registries],
        UUID[r.uuid for r in Pkg.Types.collect_registries()]
        ) == UUID[]
end

function is_pkg_available(pkg::PackageSpec)
    uuids = UUID[]
    for registry in Pkg.Types.collect_registries()
        reg_dict = Pkg.Types.read_registry(joinpath(registry.path, "Registry.toml"))
        for (uuid, pkginfo) in reg_dict["packages"]
            push!(uuids, UUID(uuid))
        end
    end
    return in(pkg.uuid, uuids)
end

function with_depot2(f)
    Base.DEPOT_PATH[1:2] .= Base.DEPOT_PATH[2:-1:1]
    f()
    Base.DEPOT_PATH[1:2] .= Base.DEPOT_PATH[2:-1:1]
end

@testset "registries" begin
    temp_pkg_dir() do depot; mktempdir() do depot2
        insert!(Base.DEPOT_PATH, 2, depot2)
        # set up registries
        regdir = mktempdir()
        setup_test_registries(regdir)
        generalurl = Pkg.Types.DEFAULT_REGISTRIES[1].url # hehe
        General = RegistrySpec(name = "General", uuid = "23338594-aafe-5451-b93e-139f81909106",
            url = generalurl)
        Foo1 = RegistrySpec(name = "RegistryFoo", uuid = "e9fceed0-5623-4384-aff0-6db4c442647a",
            url = joinpath(regdir, "RegistryFoo1"))
        Foo2 = RegistrySpec(name = "RegistryFoo", uuid = "a8e078ad-b4bd-4e09-a52f-c464826eef9d",
            url = joinpath(regdir, "RegistryFoo2"))

        # Packages in registries
        Example  = PackageSpec(name = "Example",  uuid = "7876af07-990d-54b4-ab0e-23690620f79a")
        Example1 = PackageSpec(name = "Example1", uuid = "c5f1542f-b8aa-45da-ab42-05303d706c66")
        Example2 = PackageSpec(name = "Example2", uuid = "d7897d3a-8e65-4b65-bdc8-28ce4e859565")

        # Add General registry
        ## Pkg REPL
        for reg in ("General",
                    "23338594-aafe-5451-b93e-139f81909106",
                    "General=23338594-aafe-5451-b93e-139f81909106")
            pkgstr("registry add $(reg)")
            test_installed([General])

            pkgstr("registry up $(reg)")
            test_installed([General])
            pkgstr("registry rm $(reg)")
            test_installed([])
        end
        ## Registry API
        for reg in ("General",
                    RegistrySpec("General"),
                    RegistrySpec(name = "General"),
                    RegistrySpec(name = "General", url = generalurl),
                    RegistrySpec(uuid = "23338594-aafe-5451-b93e-139f81909106"),
                    RegistrySpec(name = "General", uuid = "23338594-aafe-5451-b93e-139f81909106"))
            Pkg.Registry.add(reg)
            test_installed([General])
            @test is_pkg_available(Example)
            Pkg.Registry.update(reg)
            test_installed([General])
            Pkg.Registry.rm(reg)
            test_installed([])
            @test !is_pkg_available(Example)
        end

        # Add registry from URL/local path.
        pkgstr("registry add $(Foo1.url)")
        test_installed([Foo1])
        @test is_pkg_available(Example1)
        @test !is_pkg_available(Example2)
        with_depot2(() -> pkgstr("registry add $(Foo2.url)"))
        test_installed([Foo1, Foo2])
        @test is_pkg_available(Example1)
        @test is_pkg_available(Example2)

        # reset installed registries
        rm.(joinpath.(Base.DEPOT_PATH[1:2], "registries"); force=true, recursive=true)

        Registry.add(RegistrySpec(url = Foo1.url))
        test_installed([Foo1])
        @test is_pkg_available(Example1)
        @test !is_pkg_available(Example2)
        with_depot2(() -> Registry.add(RegistrySpec(url = Foo2.url)))
        test_installed([Foo1, Foo2])
        @test is_pkg_available(Example1)
        @test is_pkg_available(Example2)

        # Behaviour with conflicting registry names
        @test_throws PkgError pkgstr("registry up RegistryFoo")
        @test_throws PkgError Registry.update("RegistryFoo")
        @test_throws PkgError Registry.update(RegistrySpec("RegistryFoo"))
        @test_throws PkgError Registry.update(RegistrySpec(name = "RegistryFoo"))
        @test_throws PkgError pkgstr("registry remove RegistryFoo")
        @test_throws PkgError Registry.rm("RegistryFoo")
        @test_throws PkgError Registry.rm(RegistrySpec("RegistryFoo"))
        @test_throws PkgError Registry.rm(RegistrySpec(name = "RegistryFoo"))

        pkgstr("registry up $(Foo1.uuid)")
        pkgstr("registry update $(Foo1.name)=$(Foo1.uuid)")
        Registry.update(RegistrySpec(uuid = Foo1.uuid))
        Registry.update(RegistrySpec(name = Foo1.name, uuid = Foo1.uuid))

        test_installed([Foo1, Foo2])
        pkgstr("registry rm $(Foo1.uuid)")
        test_installed([Foo2])
        @test !is_pkg_available(Example1)
        @test is_pkg_available(Example2)
        Registry.add(RegistrySpec(url = Foo1.url))
        test_installed([Foo1, Foo2])
        @test is_pkg_available(Example1)
        @test is_pkg_available(Example2)
        pkgstr("registry rm $(Foo1.name)=$(Foo1.uuid)")
        test_installed([Foo2])
        @test !is_pkg_available(Example1)
        @test is_pkg_available(Example2)
        pkgstr("registry rm $(Foo2.name)")
        test_installed([])
        @test !is_pkg_available(Example1)
        @test !is_pkg_available(Example2)

        Registry.add(RegistrySpec(url = Foo1.url))
        with_depot2(() -> Registry.add(RegistrySpec(url = Foo2.url)))
        test_installed([Foo1, Foo2])
        @test is_pkg_available(Example1)
        @test is_pkg_available(Example2)
        Registry.rm(RegistrySpec(uuid = Foo1.uuid))
        test_installed([Foo2])
        @test !is_pkg_available(Example1)
        @test is_pkg_available(Example2)
        Registry.add(RegistrySpec(url = Foo1.url))
        test_installed([Foo1, Foo2])
        @test is_pkg_available(Example1)
        @test is_pkg_available(Example2)
        Registry.rm(RegistrySpec(name = Foo1.name, uuid = Foo1.uuid))
        test_installed([Foo2])
        @test !is_pkg_available(Example1)
        @test is_pkg_available(Example2)
        Registry.rm(RegistrySpec(Foo2.name))
        test_installed([])
        @test !is_pkg_available(Example1)
        @test !is_pkg_available(Example2)

        # multiple registries on the same time
        pkgstr("registry add General $(Foo1.url)")
        with_depot2(() -> pkgstr("registry add $(Foo2.url)"))
        test_installed([General, Foo1, Foo2])
        @test is_pkg_available(Example)
        @test is_pkg_available(Example1)
        @test is_pkg_available(Example2)
        pkgstr("registry up General $(Foo1.uuid) $(Foo2.name)=$(Foo2.uuid)")
        pkgstr("registry rm General $(Foo1.uuid) $(Foo2.name)=$(Foo2.uuid)")
        test_installed([])
        @test !is_pkg_available(Example)
        @test !is_pkg_available(Example1)
        @test !is_pkg_available(Example2)

        Registry.add([RegistrySpec("General"),
                      RegistrySpec(url = Foo1.url)])
        with_depot2(() -> Registry.add([RegistrySpec(url = Foo2.url)]))
        test_installed([General, Foo1, Foo2])
        @test is_pkg_available(Example)
        @test is_pkg_available(Example1)
        @test is_pkg_available(Example2)
        Registry.update([RegistrySpec("General"),
                         RegistrySpec(uuid = Foo1.uuid),
                         RegistrySpec(name = Foo2.name, uuid = Foo2.uuid)])
        Registry.rm([RegistrySpec("General"),
                     RegistrySpec(uuid = Foo1.uuid),
                     RegistrySpec(name = Foo2.name, uuid = Foo2.uuid)])
        test_installed([])
        @test !is_pkg_available(Example)
        @test !is_pkg_available(Example1)
        @test !is_pkg_available(Example2)

        # Trying to add a registry with the same name as existing one
        pkgstr("registry add $(Foo1.url)")
        @test_throws PkgError pkgstr("registry add $(Foo2.url)")
        @test_throws PkgError Registry.add([RegistrySpec(url = Foo2.url)])

    end end

    # issue #711
    temp_pkg_dir() do depot; mktempdir() do depot2
        insert!(Base.DEPOT_PATH, 2, depot2)
        Registry.add("General")
        with_depot2(() -> Registry.add("General"))
        # This add should not error because depot/Example and depot2/Example have the same uuid
        Pkg.add("Example")
        @test isinstalled((name = "Example", uuid = UUID("7876af07-990d-54b4-ab0e-23690620f79a")))
    end end
end

end # module


text:
# NOTE: it's simpler to just manually define the overloads of ChainRules
# so let's not generate them automatically for now, until we can handle
# the function signatures better, e.g wrap tracked-types over supported
# types of the rrule signatures automatically
Base.getindex(A::TrackedArrayType, indices::Int...) = trace(getindex, A, indices...)
Base.reshape(A::TrackedArrayType, shape::Int...) = trace(reshape, A, shape...)
Base.reshape(A::TrackedArrayType, shape::NTuple{N, Int}) where N = trace(reshape, A, shape)
Base.vect(X::Vararg{<:TrackedArrayType}) = trace(Base.vect, X)
Base.hcat(Xs::TrackedArrayType...) = trace(Base.hcat, Xs...)

Base.:(+)(A::TrackedType, B::TrackedType) = trace(+, A, B)
Base.:(-)(A::TrackedType, B::TrackedType) = trace(-, A, B)
Base.:(*)(A::TrackedType, B::TrackedType) = trace(*, A, B)
Base.:(-)(A::TrackedType) = trace(-, A)

Base.prod(A::TrackedArrayType) = trace(prod, A)
Base.sum(A::TrackedArrayType) = trace(sum, A)
Base.abs(X::TrackedType) = trace(abs, X)
Base.abs2(X::TrackedType) = trace(abs2, X)


text:
using Documenter, MinHash

DocMeta.setdocmeta!(MinHash, :DocTestSetup, :(using MinHash); recursive=true)

makedocs(;
    modules=[MinHash],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/jakobnissen/MinHash.jl/blob/{commit}{path}#L{line}",
    sitename="MinHash.jl",
    authors="Jakob Nybo Nissen",
    assets=String[],
)

deploydocs(;
    repo="github.com/jakobnissen/MinHash.jl.git",
)


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

include("vonmises.jl")

# Elasticity problems
abstract ElasticPlasticProblem            <: AbstractProblem
abstract PlaneStressElasticPlasticProblem <: ElasticPlasticProblem

function get_unknown_field_name{P<:ElasticPlasticProblem}(::Type{P})
    return "displacement"
end

function get_unknown_field_type{P<:ElasticPlasticProblem}(::Type{P})
    return Vector{Float64}
end

# 3D Elasticity problems
function ElasticPlasticProblem(dim::Int=3, elements=[])
    return Problem{ElasticPlasticProblem}("elasticplastic problem", dim, elements)
end

# 2D Plane stress elasticity problems
function PlaneStressElasticPlasticProblem(dim::Int=2, elements=[])
    return Problem{PlaneStressElasticPlasticProblem}("plane stress elasticplastic problem", dim, elements)
end

function get_residual_vector{P<:PlaneStressElasticPlasticProblem}(problem::Problem{P}, element::Element, ip::IntegrationPoint, time::Number; variation=nothing)
    r = zeros(Float64, problem.dim, length(element))

    J = get_jacobian(element, ip, time)

    # internal forces
    if haskey(element, "youngs modulus") && haskey(element, "poissons ratio")
        if !haskey(element, "integration points")
            if P == PlaneStressElasticPlasticProblem
                last_stress = zeros(2,2)
                last_strain = zeros(2,2)
            else
                last_stress = zeros(3,3)
                last_strain = zeros(3,3)
            end
        else
            for each_ip in element("integration points", time)
                if isapprox(each_ip.xi, ip.xi)
                    last_stress = ip("stress", time)
                    last_strain = ip("stress", time)
                    break
                end
            end
        end
        u = element("displacement", time, variation)
        grad = element(ip, time, Val{:grad})
        gradu = grad*u

        # deformation gradient
        F = I + gradu
        E = 1/2*(F'*F - I)

        #E = 1/2*(gradu + gradu') # finite strain (total)

        # material
        young = element("youngs modulus", ip, time)
        poisson = element("poissons ratio", ip, time)
        stress_y = element("yield stress", time).data
        dstrain = E - last_strain
        de_v = [dstrain[1,1], dstrain[2,2], dstrain[1,2]]
        material_model = element("material model", time)
        s = last_stress
        de = copy(ForwardDiff.get_value(dstrain))

        if P == PlaneStressElasticPlasticProblem
            C = stiffnessTensorPlaneStress(young, poisson)
            s_v = [s[1,1], s[2,2], s[1,2]]
            de_ = [de[1,1], de[2,2], de[1,2]]
            problem_stress_type = :PlaneStressElasticPlasticProblem
        else
            C = stiffnessTensor(young, poisson)
            s_v = [s[1,1], s[2,2], s[3,3], s[2,3], s[1,3], s[1,2]]
            de_ = [de[1,1], de[2,2], de[3,3], de[2,3], de[1,3], de[1,2]]
            problem_stress_type = :ElasticPlasticProblem
        end
        dep = zeros(3)
        stress_inc, dep = calculate_stress(de_,
                                           s_v,
                                           C,
                                           stress_y,
                                           Val{:vonMises},
                                           Val{problem_stress_type})
        info("%% ", dep)
        s_v += C * (de_v - dep)
        info("--: ", ForwardDiff.get_value(s_v))
        # stress
        if P == PlaneStressElasticPlasticProblem
            S = [s_v[1] s_v[3];
                 s_v[3] s_v[2]]
        else
        S = [s_v[1] s_v[6] s_v[5];
             s_v[6] s_v[2] s_v[4];
             s_v[5] s_v[4] s_v[3]]
        end
        r += F*S*grad*det(J)
    end

    # external forces - volume load
    if haskey(element, "displacement load")
        basis = element(ip, time)
        b = element("displacement load", ip, time)
        r -= b*basis*det(J)
    end

    # external forces - surface traction force
    if haskey(element, "displacement traction force")
        basis = element(ip, time)
        T = element("displacement traction force", ip, time)
        JT = transpose(J)
        s = size(JT, 2) == 1 ? JT : cross(JT[:,1], JT[:,2])
        r -= T*basis*norm(s)
    end

    return vec(r)
end

#=
function get_residual_vector{P<:ElasticPlasticProblem}(problem::Problem{P}, element::Element, ip::IntegrationPoint, time::Number; variation=nothing)
    r = zeros(Float64, problem.dim, length(element))

    J = get_jacobian(element, ip, time)

    info("_____________________")
    # internal forces
    if haskey(element, "youngs modulus") && haskey(element, "poissons ratio")

        if !haskey(element, "integration points")
            if P == PlaneStressElasticPlasticProblem
                last_stress = zeros(2,2)
                last_strain = zeros(2,2)
            else
                last_stress = zeros(3,3)
                last_strain = zeros(3,3)
            end
        else
            for each_ip in element("integration points", time)
                if isapprox(each_ip.xi, ip.xi)
                    last_stress = ip("stress", time)
                    last_strain = ip("stress", time)
                    break
                end
            end
        end
        u = element("displacement", time, variation)
        grad = element(ip, time, Val{:grad})
        gradu = grad*u

        # deformation gradient
        F = I + gradu

        # material
        young = element("youngs modulus", ip, time)
        poisson = element("poissons ratio", ip, time)
        mu = young/(2*(1+poisson))
        lambda = young*poisson/((1+poisson)*(1-2*poisson))
        if P == PlaneStressElasticityProblem
            lambda = 2*lambda*mu/(lambda + 2*mu)  # <- correction for 2d problems
        end

        # strain
        E = 1/2*(F'*F - I)
        #E = 1/2*(gradu + gradu') # finite strain (total)

        young = element("youngs modulus", ip, time)
        poisson = element("poissons ratio", ip, time)
        stress_y = element("yield stress", time).data
        dstrain = E - last_strain
        material_model = element("material model", time)
        s = last_stress
        de = ForwardDiff.get_value(dstrain)

        if P == PlaneStressElasticPlasticProblem
            C = stiffnessTensorPlaneStress(young, poisson)
            s_v = [s[1,1], s[2,2], s[1,2]]
            de_ = [de[1,1], de[2,2], de[1,2]]
            problem_stress_type = :PlaneStressElasticPlasticProblem
        else
            C = stiffnessTensor(young, poisson)
            s_v = [s[1,1], s[2,2], s[3,3], s[2,3], s[1,3], s[1,2]]
            de_ = [de[1,1], de[2,2], de[3,3], de[2,3], de[1,3], de[1,2]]
            problem_stress_type = :ElasticPlasticProblem
        end

        stress_inc, lambda = plastic_multiplier = calculate_stress(de_,
                                                           s_v,
                                                           C,
                                                           stress_y,
                                                           Val{:vonMises},
                                                           Val{problem_stress_type})

        # dep = lambda * dfds(s)
        # upate_material_parameters!(...)
        s_new = s_v + stress_inc
        #S = [s_v[1] s_v[6] s_v[5];
    #         s_v[6] s_v[2] s_v[4];
#             s_v[5] s_v[4] s_v[3]]
        S = [s_new[1] s_new[3];
             s_new[3] s_new[2]]
        # S = C * (E - dep)

        info("Stress: ", vec(ForwardDiff.get_value(S)))
        # stress
        #S = lambda*trace(E)*I + 2*mu*E

        r += F*S*grad*det(J)

    end

    # external forces - volume load
    if haskey(element, "displacement load")
        basis = element(ip, time)
        b = element("displacement load", ip, time)
        r -= b*basis*det(J)
    end

    # external forces - surface traction force
    if haskey(element, "displacement traction force")
        basis = element(ip, time)
        T = element("displacement traction force", ip, time)
        JT = transpose(J)
        s = size(JT, 2) == 1 ? JT : cross(JT[:,1], JT[:,2])
        r -= T*basis*norm(s)
    end

    return vec(r)
end
=# #fff


text:
# If you follow the equations in the PDF in the /doc directory, you will see all functions here straightforward.

# TODO: I cound that "utils.jl" has been `include`ed in the main `VectorSphericalWaves` module. Is there a more neat way? when I `include("utils.jl")` here, I get lots of warnings.

using StaticArrays

# exporting main functions, separating real and imaginary parts
export M_mn_wave_SeparateRealImag
export N_mn_wave_SeparateRealImag

export M_mn_wave_SeparateRealImag_SMatrix
export N_mn_wave_SeparateRealImag_SMatrix

#############################################################################################
# The following functions are the same as before, but are desinged to be automatically differentiable (e.g., with Zygote).
# The main difference is that we don't use any complex numbers here. Complex numbers and vectors are represented as 1x2 and nx2 Matrices, respectively. The first and second columns represent the real and imaginary parts.
#############################################################################################
#############################################################################################
# calculate B(θ), C(θ), P(θ)
"""
    equation C.19, returns an array
It returns a 3x2 Matrix, the first and second columns are the real and imaginary parts, respectively.
"""
function B_mn_of_θ_SeparateRealImag(m::I, n::I, θ::R) where {R <: Real, I <: Integer}
    # equation C.19
    return hcat(
        vcat(zero(R), τₘₙ(m, n, θ), zero(R)),
        vcat(zero(R), zero(R), πₘₙ(m, n, θ))
    )
end

"""
    Same as B_mn_of_θ_SeparateRealImag, but return SMatrix
It returns a 3x2 SMatrix, the first and second columns are the real and imaginary parts, respectively.
"""
function B_mn_of_θ_SeparateRealImag_SMatrix(m::I, n::I, θ::R) where {R <: Real, I <: Integer}
    # equation C.19
    return hcat(
        SVector(zero(R), τₘₙ(m, n, θ), zero(R)),
        SVector(zero(R), zero(R), πₘₙ(m, n, θ))
    )
end

"""
    equation C.20, returns an array
It returns a 3x2 Matrix, the first and second columns are the real and imaginary parts, respectively.
"""
function C_mn_of_θ_SeparateRealImag(m::I, n::I, θ::R) where {R <: Real, I <: Integer}
    # equation C.20
    return hcat(
        vcat(zero(R), zero(R), -1 * τₘₙ(m, n, θ)),
        vcat(zero(R), πₘₙ(m, n, θ), zero(R))
    )
end

"""
    Same as C_mn_of_θ_SeparateRealImag, but return SMatrix
It returns a 3x2 SMatrix, the first and second columns are the real and imaginary parts, respectively.
"""
function C_mn_of_θ_SeparateRealImag_SMatrix(m::I, n::I, θ::R) where {R <: Real, I <: Integer}
    # equation C.20
    return hcat(
        SVector(zero(R), zero(R), -1 * τₘₙ(m, n, θ)),
        SVector(zero(R), πₘₙ(m, n, θ), zero(R))
)
end

"""
    equation C.21, returns an array
It returns a 3x2 Matrix, the first and second columns are the real and imaginary parts, respectively.
"""
function P_mn_of_θ_SeparateRealImag(m::I, n::I, θ::R) where {R <: Real, I <: Integer}
    # equation C.21
    return hcat(
        P_mn_of_θ(m, n, θ),
        vcat(zero(R), zero(R), zero(R))
)
end

"""
    Same as P_mn_of_θ_SeparateRealImag, returns SMatrix
It returns a 3x2 SMatrix, the first and second columns are the real and imaginary parts, respectively.
"""
function P_mn_of_θ_SeparateRealImag_SMatrix(m::I, n::I, θ::R) where {R <: Real, I <: Integer}
    # equation C.21
    return hcat(
        P_mn_of_θ_SVector(m, n, θ),
        SVector(zero(θ), zero(θ), zero(θ))
    )
end

function convert_from_fun_of_θ_to_fun_of_θ_ϕ(fun_tobe_converted::Function, m::I, n::I, θ::R, ϕ::R) where {R <: Real, I <: Integer}
    # equation C.16, C.17, and C.18
    coeff = Zygote.dropgrad((-one(R))^m * sqrt_factorial_n_plus_m_over_factorial_n_minus_m(m, n))
    B_of_θ_coef = convert(R, coeff) .* fun_tobe_converted(m, n, θ)
    B_of_θ_coef_real = B_of_θ_coef[:,1]
    B_of_θ_coef_imag = B_of_θ_coef[:,2]
    exp_imϕ_real = cos(Zygote.dropgrad(m) * ϕ)
    exp_imϕ_imag = sin(Zygote.dropgrad(m) * ϕ)

    return hcat(
        B_of_θ_coef_real .* exp_imϕ_real - B_of_θ_coef_imag .* exp_imϕ_imag,
        B_of_θ_coef_real .* exp_imϕ_imag + B_of_θ_coef_imag .* exp_imϕ_real
    )
end

#############################################################################################
# calculate B(θ,ϕ), C(θ,ϕ), P(θ,ϕ)
function B_mn_of_θ_ϕ_SeparateRealImag(m::I, n::I, θ::R, ϕ::R) where {R <: Real, I <: Integer}
    # equation C.16
    return convert_from_fun_of_θ_to_fun_of_θ_ϕ(B_mn_of_θ_SeparateRealImag, m, n, θ, ϕ)
end

function C_mn_of_θ_ϕ_SeparateRealImag(m::I, n::I, θ::R, ϕ::R) where {R <: Real, I <: Integer}
    # equation C.17
    return convert_from_fun_of_θ_to_fun_of_θ_ϕ(C_mn_of_θ_SeparateRealImag, m, n, θ, ϕ)
end

function P_mn_of_θ_ϕ_SeparateRealImag(m::I, n::I, θ::R, ϕ::R) where {R <: Real, I <: Integer}
    # equation C.18
    return convert_from_fun_of_θ_to_fun_of_θ_ϕ(P_mn_of_θ_SeparateRealImag, m, n, θ, ϕ)
end

# same as above, but returns *_SMatrix
function B_mn_of_θ_ϕ_SeparateRealImag_SMatrix(m::I, n::I, θ::R, ϕ::R) where {R <: Real, I <: Integer}
    # equation C.16
    return convert_from_fun_of_θ_to_fun_of_θ_ϕ(B_mn_of_θ_SeparateRealImag_SMatrix, m, n, θ, ϕ)
end

function C_mn_of_θ_ϕ_SeparateRealImag_SMatrix(m::I, n::I, θ::R, ϕ::R) where {R <: Real, I <: Integer}
    # equation C.17
    return convert_from_fun_of_θ_to_fun_of_θ_ϕ(C_mn_of_θ_SeparateRealImag_SMatrix, m, n, θ, ϕ)
end

function P_mn_of_θ_ϕ_SeparateRealImag_SMatrix(m::I, n::I, θ::R, ϕ::R) where {R <: Real, I <: Integer}
    # equation C.18
    return convert_from_fun_of_θ_to_fun_of_θ_ϕ(P_mn_of_θ_SeparateRealImag_SMatrix, m, n, θ, ϕ)
end

#############################################################################################
# calculate spherical Bessel and Hankel functions and their derivarive
function SeparateRealImag_for_Bessel_Hankel_and_derivatives(fun_tobe_converted::Function, n::I, x_r::R, x_i::R) where {R <: Real, I <: Integer}
    bessel_complex = fun_tobe_converted(n, complex(x_r, x_i))
    return hcat(real(bessel_complex), imag(bessel_complex))
end

function spherical_Bessel_j_n_SeparateRealImag(n::I, x_r::R, x_i::R) where {R <: Real, I <: Integer}
    return SeparateRealImag_for_Bessel_Hankel_and_derivatives(spherical_Bessel_j_n, n, x_r, x_i)
end

function spherical_Bessel_y_n_SeparateRealImag(n::I, x_r::R, x_i::R) where {R <: Real, I <: Integer}
    return SeparateRealImag_for_Bessel_Hankel_and_derivatives(spherical_Bessel_y_n, n, x_r, x_i)
end

function spherical_Hankel_h1_n_SeparateRealImag(n::I, x_r::R, x_i::R) where {R <: Real, I <: Integer}
    return SeparateRealImag_for_Bessel_Hankel_and_derivatives(spherical_Hankel_h1_n, n, x_r, x_i)
end

function one_over_x_by_∂_x_j_n_by_∂x_SeparateRealImag(n::I, x_r::R, x_i::R) where {R <: Real, I <: Integer}
    # n_over_x = complex_divide(n * ones(size(x_r)), zeros(size(x_r)), x_r, x_i) # TODO: find a way to make sure all inputs have the same size
    n_over_x = complex_divide(Zygote.dropgrad(n), zero(x_r), x_r, x_i)
    return (spherical_Bessel_j_n_SeparateRealImag(n - 1, x_r, x_i) - complex_multiply(n_over_x, spherical_Bessel_j_n_SeparateRealImag(n, x_r, x_i)))
end

function one_over_x_by_∂_x_y_n_by_∂x_SeparateRealImag(n::I, x_r::R, x_i::R) where {R <: Real, I <: Integer}
    # n_over_x = complex_divide(n * ones(size(x_r)), zeros(size(x_r)), x_r, x_i) # TODO: find a way to make sure all inputs have the same size
    n_over_x = complex_divide(Zygote.dropgrad(n), zero(x_r), x_r, x_i)
    return (spherical_Bessel_y_n_SeparateRealImag(n - 1, x_r, x_i) - complex_multiply(n_over_x, spherical_Bessel_y_n_SeparateRealImag(n, x_r, x_i)))
end

function one_over_x_by_∂_x_h_n_by_∂x_SeparateRealImag(n::I, x_r::R, x_i::R) where {R <: Real, I <: Integer}
    der_j = one_over_x_by_∂_x_j_n_by_∂x_SeparateRealImag(n, x_r, x_i)
    der_y = one_over_x_by_∂_x_y_n_by_∂x_SeparateRealImag(n, x_r, x_i)
    return der_j + complex_multiply(zero(R), one(R), der_y[1], der_y[2])
    # return (one_over_x_by_∂_x_j_n_by_∂x_SeparateRealImag(n, x_r, x_i) + complex_multiply(0, 1, one_over_x_by_∂_x_y_n_by_∂x_SeparateRealImag(n, x_r, x_i)...))
end

function get_radial_function_and_special_derivative_given_kind_SeparateRealImag(kind)
    """
    return the radial function that is appropriate with the type of VSWF
    """

    if kind in ["regular", "incoming", 1]
        radial_function = spherical_Bessel_j_n_SeparateRealImag
        radial_function_special_derivative = one_over_x_by_∂_x_j_n_by_∂x_SeparateRealImag
    elseif kind in ["irregular", "outgoing", 2]
        radial_function = spherical_Hankel_h1_n_SeparateRealImag
        radial_function_special_derivative = one_over_x_by_∂_x_h_n_by_∂x_SeparateRealImag
    else
        throw(DomainError(""" 'kind' has to be one of the following: ["regular", "incoming", "irregular", "outgoing"] """))
    end

    return radial_function, radial_function_special_derivative
end

"""Calculate vector spherical wave M
This version takes kr_r and kr_i as inputs (i.e., wavevector multiplied by r)
"""
function M_mn_wave_SeparateRealImag(m::I, n::I, kr_r::R, kr_i::R, θ::R, ϕ::R, kind="regular") where {R <: Real, I <: Integer}
    radial_function, _ = get_radial_function_and_special_derivative_given_kind_SeparateRealImag(kind)
    gamma_by_radial = convert(R, Zygote.dropgrad(γ_mn(m, n))) .* radial_function(Zygote.dropgrad(n), kr_r, kr_i)
    return complex_multiply(vcat((gamma_by_radial), (gamma_by_radial), (gamma_by_radial)), C_mn_of_θ_ϕ_SeparateRealImag(Zygote.dropgrad(m), Zygote.dropgrad(n), θ, ϕ))

    # TODO: write this in a more elegant way: vcat(gamma_by_radial, gamma_by_radial, gamma_by_radial)
end

"""Calculate vector spherical wave M
This version takes k_r, k_i and r as inputs.
"""
function M_mn_wave_SeparateRealImag(m::I, n::I, k_r::R, k_i::R, r::R, θ::R, ϕ::R, kind="regular") where {R <: Real, I <: Integer}
    return M_mn_wave_SeparateRealImag(m, n, k_r*r, k_i*r, θ, ϕ, kind)
end

"""
    Same as `M_mn_wave_SeparateRealImag`, but returns SMatrix
"""
function M_mn_wave_SeparateRealImag_SMatrix(m::I, n::I, kr_r::R, kr_i::R, θ::R, ϕ::R, kind="regular") where {R <: Real, I <: Integer}
    radial_function, _ = get_radial_function_and_special_derivative_given_kind_SeparateRealImag(kind)
    gamma_by_radial = convert(R, Zygote.dropgrad(γ_mn(m, n))) .* radial_function(Zygote.dropgrad(n), kr_r, kr_i)
    return complex_multiply_SMatrix(gamma_by_radial .* SMatrix{3,2}(ones(3, 2)), C_mn_of_θ_ϕ_SeparateRealImag_SMatrix(Zygote.dropgrad(m), Zygote.dropgrad(n), θ, ϕ))
    # TODO: write this in a more elegant way: vcat(gamma_by_radial, gamma_by_radial, gamma_by_radial)
end

function M_mn_wave_SeparateRealImag_SMatrix(m::I, n::I, k_r::R, k_i::R, r::R, θ::R, ϕ::R, kind="regular") where {R <: Real, I <: Integer}
    return M_mn_wave_SeparateRealImag_SMatrix(m, n, k_r*r, k_i*r, θ, ϕ, kind)
end

function N_mn_wave_SeparateRealImag(m::I, n::I, kr_r::R, kr_i::R, θ::R, ϕ::R, kind="regular") where {R <: Real, I <: Integer}
    radial_function, radial_function_special_derivative  = get_radial_function_and_special_derivative_given_kind_SeparateRealImag(kind)
    ceoff = complex_multiply(complex_divide(Zygote.dropgrad(n * (n + 1)), zero(kr_r), kr_r, kr_i), radial_function(Zygote.dropgrad(n), kr_r, kr_i))
    rad_fun_der = radial_function_special_derivative(Zygote.dropgrad(n), kr_r, kr_i)
    return convert(R, Zygote.dropgrad(γ_mn(m, n))) .* (
        complex_multiply([ceoff; ceoff; ceoff], P_mn_of_θ_ϕ_SeparateRealImag(Zygote.dropgrad(m), n, θ, ϕ))
        + complex_multiply([rad_fun_der; rad_fun_der; rad_fun_der], B_mn_of_θ_ϕ_SeparateRealImag(Zygote.dropgrad(m), n, θ, ϕ))
    )
    # TODO: write this in a more elegant way: [ceoff; ceoff; ceoff]
end

function N_mn_wave_SeparateRealImag(m::I, n::I, k_r::R, k_i::R, r::R, θ::R, ϕ::R, kind="regular") where {R <: Real, I <: Integer}
    return N_mn_wave_SeparateRealImag(m, n, k_r*r, k_i*r, θ, ϕ, kind)
end

"""
    Same as `N_mn_wave_SeparateRealImag`, but returns SMatrix
"""
function N_mn_wave_SeparateRealImag_SMatrix(m::I, n::I, kr_r::R, kr_i::R, θ::R, ϕ::R, kind="regular") where {R <: Real, I <: Integer}
    radial_function, radial_function_special_derivative  = get_radial_function_and_special_derivative_given_kind_SeparateRealImag(kind)
    ceoff = complex_multiply(complex_divide(n * (n + 1), zero(kr_r), kr_r, kr_i), radial_function(n, kr_r, kr_i))
    rad_fun_der = radial_function_special_derivative(n, kr_r, kr_i)
    return convert(R, γ_mn(m, n)) .* (
        complex_multiply_SMatrix(ceoff .* SMatrix{3,2}(ones(3, 2)), P_mn_of_θ_ϕ_SeparateRealImag_SMatrix(m, n, θ, ϕ))
        + complex_multiply_SMatrix(rad_fun_der .* SMatrix{3,2}(ones(3, 2)), B_mn_of_θ_ϕ_SeparateRealImag_SMatrix(m, n, θ, ϕ))
    )
    # TODO: write this in a more elegant way: [ceoff; ceoff; ceoff]
end

function N_mn_wave_SeparateRealImag_SMatrix(m::I, n::I, k_r::R, k_i::R, r::R, θ::R, ϕ::R, kind="regular") where {R <: Real, I <: Integer}
    return N_mn_wave_SeparateRealImag_SMatrix(m, n, k_r*r, k_i*r, θ, ϕ, kind)
end


text:
# include("AlgAss/AlgAssTypes.jl")
include("AlgAss/AlgAss.jl")
include("AlgAss/AlgAssElem.jl")
include("AlgAss/Map.jl")
include("AlgAss/CSAMaxOrd.jl")


text:
abstract type AbstractAGCFormulation <: AbstractServiceFormulation end
struct PIDSmoothACE <: AbstractAGCFormulation end

"""
Steady State deviation of the frequency
"""
function add_variables!(psi_container::PSIContainer, ::Type{SteadyStateFrequencyDeviation})
    variable_name = make_variable_name(SteadyStateFrequencyDeviation)
    time_steps = model_time_steps(psi_container)
    variable = add_var_container!(psi_container, variable_name, time_steps)
    for t in time_steps
        variable[t] = JuMP.@variable(psi_container.JuMPmodel, base_name = "ΔF_{$(t)}")
    end
end

"""
This function add the upwards scheduled regulation variables for power generation output to the model
"""
function AddVariableSpec(
    ::Type{ActivePowerVariable},
    ::Type{U},
    psi_container::PSIContainer,
) where {U <: PSY.Area}
    return AddVariableSpec(;
        variable_name = make_variable_name(ActivePowerVariable, U),
        binary = false,
    )
end

"""
This function adds the smooth ACE Variable
"""
function AddVariableSpec(
    ::Type{SmoothACE},
    ::Type{U},
    psi_container::PSIContainer,
) where {U <: PSY.AggregationTopology}
    return AddVariableSpec(;
        variable_name = make_variable_name(SmoothACE, U),
        binary = false,
    )
end

"""
This function add the upwards scheduled regulation variables for power generation output to the model
"""
function AddVariableSpec(
    ::Type{DeltaActivePowerUpVariable},
    ::Type{U},
    psi_container::PSIContainer,
) where {U <: PSY.Area}
    return AddVariableSpec(;
        variable_name = make_variable_name(DeltaActivePowerUpVariable, U),
        binary = false,
        lb_value_func = x -> 0.0,
    )
end

"""
This function add the downwards scheduled regulation variables for power generation output to the model
"""
function AddVariableSpec(
    ::Type{DeltaActivePowerDownVariable},
    ::Type{U},
    psi_container::PSIContainer,
) where {U <: PSY.Area}
    return AddVariableSpec(;
        variable_name = make_variable_name(DeltaActivePowerDownVariable, U),
        binary = false,
        lb_value_func = x -> 0.0,
    )
end

#= Commented out since not in use. These functions will substitute balancing_auxiliary_variables!
"""
This function add the upwards scheduled regulation variables for power generation output to the model
"""
function AddVariableSpec(
    ::Type{AdditionalDeltaActivePowerUpVariable},
    ::Type{U},
    psi_container::PSIContainer,
) where {U <: PSY.Area}
    return AddVariableSpec(;
        variable_name = make_variable_name(AdditionalDeltaActivePowerUpVariable, U),
        binary = false,
        lb_value_func = x -> 0.0,
        expression_name = :emergency_up,
    )
end

"""
This function add the variables for power generation output to the model
"""
function AddVariableSpec(
    ::Type{AdditionalDeltaActivePowerDownVariable},
    ::Type{U},
    psi_container::PSIContainer,
) where {U <: PSY.Area}
    return AddVariableSpec(;
        variable_name = make_variable_name(AdditionalDeltaActivePowerDownVariable, U),
        binary = false,
        lb_value_func = x -> 0.0,
        expression_name = :emergency_dn,
    )
end
=#

function AddVariableSpec(
    ::Type{T},
    ::Type{PSY.Area},
    ::PSIContainer,
) where {T <: AreaMismatchVariable}
    return AddVariableSpec(;
        variable_name = make_variable_name(AreaMismatchVariable),
        binary = false,
    )
end

function AddVariableSpec(
    ::Type{T},
    ::Type{PSY.Area},
    ::PSIContainer,
) where {T <: LiftVariable}
    AddVariableSpec(;
        variable_name = make_variable_name(LiftVariable),
        binary = false,
        lb_value_func = x -> 0.0,
    )
end

function balancing_auxiliary_variables!(psi_container, sys)
    area_names = [PSY.get_name(a) for a in PSY.get_components(PSY.Area, sys)]
    time_steps = model_time_steps(psi_container)
    R_up_emergency = JuMPVariableArray(undef, area_names, time_steps)
    R_dn_emergency = JuMPVariableArray(undef, area_names, time_steps)
    assign_variable!(
        psi_container,
        make_variable_name(AdditionalDeltaActivePowerUpVariable, PSY.Area),
        R_up_emergency,
    )
    assign_variable!(
        psi_container,
        make_variable_name(AdditionalDeltaActivePowerDownVariable, PSY.Area),
        R_dn_emergency,
    )
    emergency_up =
        add_expression_container!(psi_container, :emergency_up, area_names, time_steps)
    emergency_dn =
        add_expression_container!(psi_container, :emergency_dn, area_names, time_steps)
    for t in time_steps, a in area_names
        R_up_emergency[a, t] = JuMP.@variable(
            psi_container.JuMPmodel,
            base_name = "Re_up_{$(a),$(t)}",
            lower_bound = 0.0
        )
        emergency_up[a, t] = R_up_emergency[a, t] + 0.0
        R_dn_emergency[a, t] = JuMP.@variable(
            psi_container.JuMPmodel,
            base_name = "Re_dn_{$(a),$(t)}",
            lower_bound = 0.0
        )
        emergency_dn[a, t] = R_dn_emergency[a, t] + 0.0
    end

    return
end

function absolute_value_lift(psi_container::PSIContainer, areas)
    time_steps = model_time_steps(psi_container)
    area_names = [PSY.get_name(a) for a in areas]
    container_lb = JuMPConstraintArray(undef, area_names, time_steps)
    assign_constraint!(psi_container, "absolute_value_lb", container_lb)
    container_ub = JuMPConstraintArray(undef, area_names, time_steps)
    assign_constraint!(psi_container, "absolute_value_ub", container_ub)
    mismatch = get_variable(psi_container, :area_mismatch)
    z = get_variable(psi_container, :lift)

    for t in time_steps, a in area_names
        container_lb[a, t] =
            JuMP.@constraint(psi_container.JuMPmodel, mismatch[a, t] <= z[a, t])
        container_ub[a, t] =
            JuMP.@constraint(psi_container.JuMPmodel, -1 * mismatch[a, t] <= z[a, t])
    end

    JuMP.add_to_expression!(
        psi_container.cost_function,
        sum(z[a, t] for t in time_steps, a in area_names) * SERVICES_SLACK_COST,
    )
    return
end

"""
Expression for the power deviation given deviation in the frequency. This expression allows updating the response of the frequency depending on commitment decisions
"""
function frequency_response_constraint!(psi_container::PSIContainer, sys::PSY.System)
    time_steps = model_time_steps(psi_container)
    services = PSY.get_components(PSY.AGC, sys)
    area_names = [PSY.get_name(PSY.get_area(s)) for s in services]
    frequency_response = 0.0
    for area in PSY.get_components(PSY.Area, sys)
        frequency_response += PSY.get_load_response(area)
    end
    for g in PSY.get_components(PSY.RegulationDevice, sys, x -> PSY.get_available(x))
        d = PSY.get_droop(g)
        response = 1 / d
        frequency_response += response
    end

    @assert frequency_response >= 0.0
    # This value is the one updated later in simulation based on the UC result
    inv_frequency_reponse = 1 / frequency_response
    area_balance = get_variable(psi_container, ActivePowerVariable, PSY.Area)
    frequency = get_variable(psi_container, "Δf")
    R_up = get_variable(psi_container, DeltaActivePowerUpVariable, PSY.Area)
    R_dn = get_variable(psi_container, DeltaActivePowerDownVariable, PSY.Area)
    R_up_emergency =
        get_variable(psi_container, AdditionalDeltaActivePowerUpVariable, PSY.Area)
    R_dn_emergency =
        get_variable(psi_container, AdditionalDeltaActivePowerUpVariable, PSY.Area)

    container = JuMPConstraintArray(undef, time_steps)
    assign_constraint!(psi_container, "frequency_response", container)

    for s in services, t in time_steps
        system_balance = sum(area_balance.data[:, t])
        total_reg = JuMP.AffExpr(0.0)
        for a in area_names
            JuMP.add_to_expression!(total_reg, R_up[a, t])
            JuMP.add_to_expression!(total_reg, -1 * R_dn[a, t])
            JuMP.add_to_expression!(total_reg, R_up_emergency[a, t])
            JuMP.add_to_expression!(total_reg, -1 * R_dn_emergency[a, t])
        end
        container[t] = JuMP.@constraint(
            psi_container.JuMPmodel,
            frequency[t] == -inv_frequency_reponse * (system_balance + total_reg)
        )
    end
    return
end

function smooth_ace_pid!(psi_container::PSIContainer, services::Vector{PSY.AGC})
    time_steps = model_time_steps(psi_container)
    area_names = [PSY.get_name(PSY.get_area(s)) for s in services]
    RAW_ACE = add_expression_container!(psi_container, :RAW_ACE, area_names, time_steps)
    SACE = get_variable(psi_container, SmoothACE, PSY.Area)
    SACE_pid = JuMPConstraintArray(undef, area_names, time_steps)
    assign_constraint!(psi_container, "SACE_pid", SACE_pid)

    Δf = get_variable(psi_container, make_variable_name("Δf"))

    for (ix, service) in enumerate(services)
        kp = PSY.get_K_p(service)
        ki = PSY.get_K_i(service)
        kd = PSY.get_K_d(service)
        B = PSY.get_bias(service)
        Δt = convert(Dates.Second, psi_container.resolution).value
        a = PSY.get_name(PSY.get_area(service))
        for t in time_steps
            # Todo: Add initial Frequency Deviation
            RAW_ACE[a, t] = -10 * B * Δf[t] + 0.0
            SACE[a, t] =
                JuMP.@variable(psi_container.JuMPmodel, base_name = "SACE_{$(a),$(t)}")
            if t == 1
                SACE_ini =
                    get_initial_conditions(psi_container, ICKey(AreaControlError, PSY.AGC))[ix]
                sace_exp = SACE_ini.value + kp * ((1 + Δt / (kp / ki)) * (RAW_ACE[a, t]))
                SACE_pid[a, t] =
                    JuMP.@constraint(psi_container.JuMPmodel, SACE[a, t] == sace_exp)
                continue
            end
            SACE_pid[a, t] = JuMP.@constraint(
                psi_container.JuMPmodel,
                SACE[a, t] ==
                SACE[a, t - 1] +
                kp * (
                    (1 + Δt / (kp / ki) + (kd / kp) / Δt) * (RAW_ACE[a, t]) +
                    (-1 - 2 * (kd / kp) / Δt) * (RAW_ACE[a, t - 1])
                )
            )
        end
    end
    return
end

function aux_constraints!(psi_container::PSIContainer, sys::PSY.System)
    time_steps = model_time_steps(psi_container)
    area_names = [PSY.get_name(a) for a in PSY.get_components(PSY.Area, sys)]
    aux_equation = JuMPConstraintArray(undef, area_names, time_steps)
    assign_constraint!(psi_container, "balance_aux", aux_equation)
    area_mismatch = get_variable(psi_container, :area_mismatch)
    SACE = get_variable(psi_container, SmoothACE, PSY.Area)
    R_up = get_variable(psi_container, DeltaActivePowerUpVariable, PSY.Area)
    R_dn = get_variable(psi_container, DeltaActivePowerDownVariable, PSY.Area)
    R_up_emergency =
        get_variable(psi_container, AdditionalDeltaActivePowerUpVariable, PSY.Area)
    R_dn_emergency =
        get_variable(psi_container, AdditionalDeltaActivePowerUpVariable, PSY.Area)

    for t in time_steps, a in area_names
        aux_equation[a, t] = JuMP.@constraint(
            psi_container.JuMPmodel,
            -1 * SACE[a, t] ==
            (R_up[a, t] - R_dn[a, t]) +
            (R_up_emergency[a, t] - R_dn_emergency[a, t]) +
            area_mismatch[a, t]
        )
    end
    return
end


text:
function atmosphere(mtime)

  global Θv,TA,TS,ρA,ustar,LH,SH,E,LW,TAU,Qnet,SW
  global KAm,KAt
  global UA1,qA1,ΘA1
  global qA,ΘA,UA

  ########################################################################
  # Atmospheric explicit time stepping
  ########################################################################
  TAUS=sign(UA[1]);

  UAst[1]=UA[1]-TAU*TAUS/(ΔZA*ρA[1])*ΔT-(1-Aimp)*(
      KAm[1].*(UA[1]-UA[2])/ΔZA^2*ΔT);
  UAst[end]=UA[end]-(1-Aimp)*(
      KAm[end]*(UA[end]-UAup)/ΔZA^2*ΔT);
  UAst[2:end-1]=UA[2:end-1]+(1-Aimp)*(
      KAm[1:end-2].*(UA[1:end-2]-UA[2:end-1])/ΔZA^2*ΔT-
      KAm[2:end-1].*(UA[2:end-1]-UA[3:end  ])/ΔZA^2*ΔT);

  qAst[1]=qA[1]+E/(ΔZA*ρO)*ΔT-
      (1-Aimp)*(KAt[1]*(qA[1]-qA[2])/ΔZA^2*ΔT+
      (-γcq[1].*KAt[1])/ΔZA*ΔT);
  qAst[end]=qA[end]-(1-Aimp)*(KAt[end]*(qA[end]-q0)/ΔZA^2*ΔT+
      (γcq[end-1].*KAt[end-1]-γcq[end].*KAt[end])/ΔZA*ΔT);
  qAst[2:end-1]=qA[2:end-1]+(1-Aimp)*(
      KAt[1:end-2].*(qA[1:end-2]-qA[2:end-1])/ΔZA^2*ΔT-
      KAt[2:end-1].*(qA[2:end-1]-qA[3:end  ])/ΔZA^2*ΔT+
      (γcq[1:end-2].*KAt[1:end-2]-γcq[2:end-1].*KAt[2:end-1])/ΔZA*ΔT);

  ΘAst[1]=ΘA[1]+(SH)/cpa/(ΔZA*ΘA[1])*ΔT-
        (1-Aimp)*(KAt[1]*(ΘA[1]-ΘA[2])/ΔZA^2*ΔT+
        (-γct[1].*KAt[1])/ΔZA*ΔT);
  ΘAst[end]=ΘA[end]-
        (1-Aimp)*(KAt[end]*(ΘA[end]-ΘA0)/ΔZA^2*ΔT+
        (γct[end-1].*KAt[end-1]-γct[end].*KAt[end])/ΔZA*ΔT);
  ΘAst[2:end-1]=ΘA[2:end-1]+(1-Aimp)*(
        KAt[1:end-2].*(ΘA[1:end-2]-ΘA[2:end-1])/ΔZA^2*ΔT-
        KAt[2:end-1].*(ΘA[2:end-1]-ΘA[3:end  ])/ΔZA^2*ΔT+
        (γct[1:end-2].*KAt[1:end-2]-γct[2:end-1].*KAt[2:end-1])/ΔZA*ΔT);

  if Aimp==0
      UA1=UAst; qA1=qAst; ΘA1=ΘAst;
  else
      du=-Aimp*KAm[1:end-1].*ΔT./ΔZA^2; dl=-Aimp*KAm[2:end].*ΔT./ΔZA^2;
      d=(1 .+Aimp.*ΔT./ΔZA^2 .*([KAm[1:end-1]; 0]+[0; KAm[2:end]]));
      A=inv(Tridiagonal(dl,d,du));
      UA1=A*UAst;
      du=-Aimp*KAt[1:end-1].*ΔT./ΔZA^2; dl=-Aimp*KAt[2:end].*ΔT./ΔZA^2;
      d=(1 .+Aimp.*ΔT./ΔZA^2 .*([KAt[1:end-1]; 0]+[0; KAt[2:end]]));
      A=inv(Tridiagonal(dl,d,du));
      qA1=A*qAst; ΘA1=A*ΘAst;
  end

  TA=(ΘA1).*(PA./PA[1]).^(2/7).-d2k;
  ρA=PA./(Rgas.*(TA.+d2k));

  if moist==1
      qsat=saltsat*cvapor_fac*exp.(-cvapor_exp./(TA.+d2k))./ρA;
      dq=(qA.-qsat).*((qA1.-qsat).>0);
      dq=min.(dq,qA1);
      # change of temperature due to condensation
      dΘA1=Av./cpa.*(PA[1]./PA).^(2/7).*dq;
      # change of temperature due to change in staurated vapor pressure
      dΘA2=Av.^2 ./Rv ./cpa.*(PA[1]./PA).^(2*2/7).*qsat./ΘA1.^2 .*dq;
      #println("dΘA")
      #println(dΘA1)
      #println("============")
      dΘA=dΘA1+dΘA2; #[J/kg][j-1 kg K]
      #println([i maximum(dTHETA) maximum(dq)])
      ind=qA1.>qsat;
      qA1[ind]=qA1[ind]-dq[ind];
      ΘA1[ind]=ΘA1[ind]+dΘA[ind];
      TA=(ΘA1).*(PA./PA[1]).^(2/7).-d2k;
  end

  qA=qA1; ΘA=ΘA1; UA=UA1;
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualNetworkPropertiesFormat <: SwaggerModel
    addressSpace::Any # spec type: Union{ Nothing, AddressSpace } # spec name: addressSpace
    dhcpOptions::Any # spec type: Union{ Nothing, DhcpOptions } # spec name: dhcpOptions
    subnets::Any # spec type: Union{ Nothing, Vector{Subnet} } # spec name: subnets
    virtualNetworkPeerings::Any # spec type: Union{ Nothing, Vector{VirtualNetworkPeering} } # spec name: virtualNetworkPeerings
    resourceGuid::Any # spec type: Union{ Nothing, String } # spec name: resourceGuid
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState
    enableDdosProtection::Any # spec type: Union{ Nothing, Bool } # spec name: enableDdosProtection
    enableVmProtection::Any # spec type: Union{ Nothing, Bool } # spec name: enableVmProtection
    ddosProtectionPlan::Any # spec type: Union{ Nothing, SubResource } # spec name: ddosProtectionPlan
    bgpCommunities::Any # spec type: Union{ Nothing, VirtualNetworkBgpCommunities } # spec name: bgpCommunities
    ipAllocations::Any # spec type: Union{ Nothing, Vector{SubResource} } # spec name: ipAllocations

    function VirtualNetworkPropertiesFormat(;addressSpace=nothing, dhcpOptions=nothing, subnets=nothing, virtualNetworkPeerings=nothing, resourceGuid=nothing, provisioningState=nothing, enableDdosProtection=false, enableVmProtection=false, ddosProtectionPlan=nothing, bgpCommunities=nothing, ipAllocations=nothing)
        o = new()
        validate_property(VirtualNetworkPropertiesFormat, Symbol("addressSpace"), addressSpace)
        setfield!(o, Symbol("addressSpace"), addressSpace)
        validate_property(VirtualNetworkPropertiesFormat, Symbol("dhcpOptions"), dhcpOptions)
        setfield!(o, Symbol("dhcpOptions"), dhcpOptions)
        validate_property(VirtualNetworkPropertiesFormat, Symbol("subnets"), subnets)
        setfield!(o, Symbol("subnets"), subnets)
        validate_property(VirtualNetworkPropertiesFormat, Symbol("virtualNetworkPeerings"), virtualNetworkPeerings)
        setfield!(o, Symbol("virtualNetworkPeerings"), virtualNetworkPeerings)
        validate_property(VirtualNetworkPropertiesFormat, Symbol("resourceGuid"), resourceGuid)
        setfield!(o, Symbol("resourceGuid"), resourceGuid)
        validate_property(VirtualNetworkPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(VirtualNetworkPropertiesFormat, Symbol("enableDdosProtection"), enableDdosProtection)
        setfield!(o, Symbol("enableDdosProtection"), enableDdosProtection)
        validate_property(VirtualNetworkPropertiesFormat, Symbol("enableVmProtection"), enableVmProtection)
        setfield!(o, Symbol("enableVmProtection"), enableVmProtection)
        validate_property(VirtualNetworkPropertiesFormat, Symbol("ddosProtectionPlan"), ddosProtectionPlan)
        setfield!(o, Symbol("ddosProtectionPlan"), ddosProtectionPlan)
        validate_property(VirtualNetworkPropertiesFormat, Symbol("bgpCommunities"), bgpCommunities)
        setfield!(o, Symbol("bgpCommunities"), bgpCommunities)
        validate_property(VirtualNetworkPropertiesFormat, Symbol("ipAllocations"), ipAllocations)
        setfield!(o, Symbol("ipAllocations"), ipAllocations)
        o
    end
end # type VirtualNetworkPropertiesFormat

const _property_map_VirtualNetworkPropertiesFormat = Dict{Symbol,Symbol}(Symbol("addressSpace")=>Symbol("addressSpace"), Symbol("dhcpOptions")=>Symbol("dhcpOptions"), Symbol("subnets")=>Symbol("subnets"), Symbol("virtualNetworkPeerings")=>Symbol("virtualNetworkPeerings"), Symbol("resourceGuid")=>Symbol("resourceGuid"), Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("enableDdosProtection")=>Symbol("enableDdosProtection"), Symbol("enableVmProtection")=>Symbol("enableVmProtection"), Symbol("ddosProtectionPlan")=>Symbol("ddosProtectionPlan"), Symbol("bgpCommunities")=>Symbol("bgpCommunities"), Symbol("ipAllocations")=>Symbol("ipAllocations"))
const _property_types_VirtualNetworkPropertiesFormat = Dict{Symbol,String}(Symbol("addressSpace")=>"AddressSpace", Symbol("dhcpOptions")=>"DhcpOptions", Symbol("subnets")=>"Vector{Subnet}", Symbol("virtualNetworkPeerings")=>"Vector{VirtualNetworkPeering}", Symbol("resourceGuid")=>"String", Symbol("provisioningState")=>"ProvisioningState", Symbol("enableDdosProtection")=>"Bool", Symbol("enableVmProtection")=>"Bool", Symbol("ddosProtectionPlan")=>"SubResource", Symbol("bgpCommunities")=>"VirtualNetworkBgpCommunities", Symbol("ipAllocations")=>"Vector{SubResource}")
Base.propertynames(::Type{ VirtualNetworkPropertiesFormat }) = collect(keys(_property_map_VirtualNetworkPropertiesFormat))
Swagger.property_type(::Type{ VirtualNetworkPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualNetworkPropertiesFormat[name]))}
Swagger.field_name(::Type{ VirtualNetworkPropertiesFormat }, property_name::Symbol) =  _property_map_VirtualNetworkPropertiesFormat[property_name]

function check_required(o::VirtualNetworkPropertiesFormat)
    true
end

function validate_property(::Type{ VirtualNetworkPropertiesFormat }, name::Symbol, val)
end


text:
using Logging

# define safe loggers that use raw streams,
# since we can't use regular streams that lock upon use
for level in [:trace, :debug, :info, :warn, :error, :fatal]
    @eval begin
        macro $(Symbol("safe_$level"))(ex...)
            macrocall = :(@placeholder $(ex...))
            # NOTE: `@placeholder` in order to avoid hard-coding @__LINE__ etc
            macrocall.args[1] = Symbol($"@$level")
            quote
                old_logger = global_logger()
                global_logger(Logging.ConsoleLogger(Core.stderr, old_logger.min_level))
                ret = $(esc(macrocall))
                global_logger(old_logger)
                ret
            end
        end
    end
end


text:
# using Ipopt
include(joinpath(pwd(),"src/Ipopt.jl"))
using Test, ForwardDiff, LinearAlgebra

function linear_interp(x0,xf,T)
    n = length(x0)
    X = [copy(Array(x0)) for t = 1:T]

    for t = 1:T
        for i = 1:n
            X[t][i] = (xf[i]-x0[i])/(T-1)*(t-1) + x0[i]
        end
    end

    return X
end

mutable struct Hopper{T,S}
    mb::T
    ml::T
    Jb::T
    Jl::T
    r::T
    μ::T
    g::T
    k::S
    Δt::T
end

# Dimensions
nq = 5 # configuration dim
nu = 2 # control dim
nc = 1 # number of contact points
nf = 1 # number of faces for friction cone pyramid
nβ = nc*nf
ns = 2

nx = nq+nu+nc+nβ+nc+ns
np = nq+nβ+4nc
T = 10 # number of time steps to optimize

# Parameters
g = 9.81 # gravity
Δt = 0.1 # time step
μ = 0.5  # coefficient of friction
mb = 10. # body mass
ml = 1.  # leg mass
Jb = 2.5 # body inertia
Jl = 0.25 # leg inertia

# Kinematics
r = 0.7
p1(q) = [q[1] + q[3]*sin(q[5]), q[2] - q[3]*cos(q[5])]

# Methods
M(h::Hopper,q) = Diagonal([h.mb+h.ml, h.mb+h.ml, h.ml, h.Jb, h.Jl])
∇V(h::Hopper,q) = [0., (h.mb+h.ml)*h.g, 0., 0., 0.]

C(h::Hopper,qk,qn) = zeros(nq)

function ϕ(::Hopper,q)
    q[2] - q[3]*cos(q[5])
end

N(::Hopper,q) = ([0., 1., -cos(q[5]), 0., q[3]*sin(q[5])])'

function P(::Hopper,q)
        ([1., 0., sin(q[5]), 0., q[3]*cos(q[5])])'
end

B(::Hopper,q) = [0. 0. 0. 1. -1.;
                 0. 0. 1. 0. 0.]

model = Hopper(mb,ml,Jb,Jl,r,μ,g,p1,Δt)

function unpack(x)
    q = x[1:nq]
    u = x[nq .+ (1:nu)]
    y = nc == 1 ? x[nq+nu+nc] : x[nq+nu .+ (1:nc)]
    β = nβ == 1 ? x[nq+nu+nc+nβ] : x[nq+nu+nc .+ (1:nβ)]
    ψ = nc == 1 ? x[nq+nu+nc+nβ+nc] : x[nq+nu+nc+nβ .+ (1:nc)]
    s = x[nq+nu+nc+nβ+nc .+ (1:ns)]

    return q,u,y,β,ψ,s
end

W = Diagonal([1.0,1.0,1.0,1.0e-1,1.0e-1])
R = Diagonal([1.0e-1,1.0e-3])
Wf = Diagonal(5.0*ones(nq))
q0 = [0., r, r, 0., 0.]
qf = [2.0, r, r, 0., 0.]
q0_ref = [0., model.r, 0.5*model.r, 0., 0.]
qf_ref = [2.0, model.r, 0.5*model.r, 0., 0.]
Q0_ref = [q0,linear_interp(q0_ref,qf_ref,T+1)[1:end-1]...,qf]
uf = zeros(nu)
w = -W*qf
wf = -Wf*qf
rr = -R*uf
obj_c = [0.5*(Q0_ref[t+2]'*W*Q0_ref[t+2] + uf'*R*uf) for t = 1:T]
obj_cf = [0.5*(Q0_ref[t+2]'*Wf*Q0_ref[t+2] + uf'*R*uf) for t = 1:T]

Q0 = linear_interp(q0,qf,T+2)

qpp = Q0[2]
qp = Q0[2]

function obj(z)
    _sum = 0.
    for t = 1:T
        q,u,y,β,ψ,s = unpack(z[(t-1)*nx .+ (1:nx)])

        if t != T
            _sum += 0.5*q'*W*q + w'*q + 0.5*u'*R*u + rr'*u + obj_c[t] + λ_global[(t-1)*ns .+ (1:ns)]'*s + 0.5*mu_prev*s'*s
        else
            _sum += 0.5*q'*Wf*q + wf'*q + 0.5*u'*R*u + rr'*u + obj_cf[t] + λ_global[(t-1)*ns .+ (1:ns)]'*s + 0.5*mu_prev*s'*s
        end
    end
    return _sum
end

function con!(c,z)

    for t = 1:T
        q,u,y,β,ψ,s = unpack(z[(t-1)*nx .+ (1:nx)])

        if t == 1
            _qpp = qpp
            _qp = qp
        elseif t == 2
            _qpp = qp
            _qp = z[(t-2)*nx .+ (1:nq)]
        else
            _qpp = z[(t-3)*nx .+ (1:nq)]
            _qp = z[(t-2)*nx .+ (1:nq)]
        end

        c[(t-1)*np .+ (1:np)] .= [(ϕ(model,q));
                                  ((model.μ*y)^2 - β'*β);
                                  (1/model.Δt*(M(model,_qpp)*(_qp - _qpp) - M(model,_qp)*(q - _qp)) - model.Δt*∇V(model,_qp) + B(model,q)'*u +  N(model,q)'*y + P(model,q)'*β);
                                  (P(model,q)*(q-_qp)/model.Δt + 2.0*β*ψ);
                                  s[1] - ϕ(model,q)*y;
                                  s[2] - ((model.μ*y)^2 - β'*β)*ψ]
     end
     return c
end

n = T*nx
m = T*np
x_L = zeros(T*nx)
x_U = 1.0e8*ones(T*nx)
for t = 1:T
    x_L[(t-1)*nx .+ (1:nq+nu)] .= -1.0e8
    x_L[(t-1)*nx + 3] = model.r/2.
    x_U[(t-1)*nx + 3] = model.r
    x_L[(t-1)*nx+nq+nu+nc .+ (1:nβ)] .= -1.0e8
    # x_L[(t-1)*nx+nq+nu+nc+nβ+nc .+ (1:ns)] .= -1.0e8
end
# cI_idx_t = zeros(Bool,np)
# cI_idx_t[1:nc+nc] .= 1
# cI_idx = zeros(Bool,m)
# for t = 1:T
#     cI_idx[(t-1)*np .+ (1:np)] .= cI_idx_t
# end

g_L = zeros(m)
g_U = zeros(m)
for t = 1:T
    g_U[(t-1)*np .+ (1:2nc)] .= 1.0e8
    # g_U[(t-1)*np + nq+nβ+2nc .+ (1:2nc)] .= 1.0e8
end

u0 = 1.0e-5*ones(nu)
y0 = 1.0e-5*ones(1)[1]
β0 = 1.0e-5*ones(nβ)[1]
ψ0 = 1.0e-5*ones(1)[1]
s0 = 1.0e-5*ones(ns)

x0 = zeros(T*nx)
for t = 1:T
    x0[(t-1)*nx .+ (1:nx)] .= [Q0[t+2];u0;y0;β0;ψ0;s0]
end

function row_col!(row,col,r,c)
    for cc in c
        for rr in r
            push!(row,convert(Int,rr))
            push!(col,convert(Int,cc))
        end
    end
    return row, col
end

function eval_f(x::Vector{Float64})
  x_global .= x
  return obj(x)
end

function eval_g(x::Vector{Float64}, g::Vector{Float64})
	con!(g,x)

	return g
end

function eval_grad_f(x::Vector{Float64}, grad_f::Vector{Float64})
	grad_f .= ForwardDiff.gradient(obj,x)
end

function eval_jac_g(x::Vector{Float64}, mode, rows::Vector{Int32}, cols::Vector{Int32}, values::Vector{Float64})
  if mode == :Structure
		row = []
		col = []
    row_col!(row,col,(1:m),(1:n))
		rows .= vec(row)
		cols .= vec(col)
  else

  	values .= vec(ForwardDiff.jacobian(con!,zeros(m),x))
  end
end

function convergence(x,tol)
	s_tmp = zeros(ns*T)
	for t = 1:T
		q,u,y,β,ψ,s = unpack(x[(t-1)*nx .+ (1:nx)])
		s_tmp[(t-1)*ns .+ (1:ns)] = copy(s)
	end

	cond1 = norm(s_tmp,Inf) < tol

    c_tmp = zeros(m)
	con!(c_tmp,x_global)
	for t = 1:T
		c_tmp[(t-1)*np .+ (1:2nc)] .= min.(0.0,c_tmp[(t-1)*np .+ (1:2:nc)])
	end
	cond2 = norm(c_tmp,Inf) < tol

	return cond1 && cond2
end
convergence(x0,1.0)

# This tests callbacks.
function intermediate(alg_mod::Int, iter_count::Int,
  obj_value::Float64, inf_pr::Float64, inf_du::Float64, mu::Float64,
  d_norm::Float64, regularization_size::Float64, alpha_du::Float64, alpha_pr::Float64,
  ls_trials::Int)

	# println("CALLBACK")
  # println("   mu_prev: $(mu_prev)")
  # println("   mu: $(mu)")
  # println("   x: $(x_global)")

  if mu != mu_prev
    println("-mu difference-")
    println("-AL update-")
	for t = 1:T
		q,u,y,β,ψ,s = unpack(x_global[(t-1)*nx .+ (1:nx)])

    	λ_global[(t-1)*ns .+ (1:ns)] .+= 1/mu_prev*s
	end
    global mu_prev = mu
  end

  return true#!convergence(x_global,1.0e-3) #iter_count < 1  # Interrupts after one iteration.
end

function eval_h(x::Vector{Float64}, mode, rows::Vector{Int32}, cols::Vector{Int32}, obj_factor::Float64, lambda::Vector{Float64}, values::Vector{Float64})
  if mode == :Structure

  else

  end
end

prob = createProblem(n, x_L, x_U, m, g_L, g_U, m*n, n*n,
                     eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h)
addOption(prob,"hessian_approximation","limited-memory")
addOption(prob,"mu_strategy","monotone")

addOption(prob,"constr_viol_tol",1.0e-2)
addOption(prob,"tol",1.0e-2)

# addOption(prob,"")
setIntermediateCallback(prob, intermediate)

global x_global = zeros(n)
global λ_global = zeros(ns*T)
global mu_prev = 0.1

prob.x .= copy(x0)
x_global .= prob.x
@time solvestat = solveProblem(prob)

λ_global
mu_prev
# @test Ipopt.ApplicationReturnStatus[solvestat] == :User_Requested_Stop

convergence(prob.x,1.0e-3)
for t = 1:T
	for t = 1:T
		q,u,y,β,ψ,s = unpack(x_global[(t-1)*nx .+ (1:nx)])
		println("s: $(s)")
	end
end


text:
module KrylovMethods

	import Base.BLAS

	include("cg.jl")
	include("blockCG.jl")
	include("cgls.jl")
	include("bicgstb.jl")
	include("blockBiCGSTB.jl")
	include("gmres.jl")
	include("lanczosBidiag.jl")
	include("ssor.jl")
	include("lsqr.jl")
	include("lanczosTridiag.jl")
	include("lanczos.jl")
	include("minres.jl")
	include("gs.jl")
end


text:
xfgridsections = Dict([("xf-comment", 0), ("xf-header", 1), ("xf-dimension", 2),
("xf-node", 10), ("xf-periodic-face", 18), ("xf-cell", 12), ("xf-face", 13), ("xf-face-tree", 59),
  ("xf-cell-tree", 58), ("xf-face-parents", 61), ("xf-zone", 45)])

celltype_nodes_per_cell = Dict([(0,nothing),(1,3),(2,4),(3,4),(4,6),(5,5),(6,5)])
celltype_faces_per_cell = Dict([(0,nothing),(1,3),(2,4),(3,4),(4,8),(5,5),(6,6)])
facetype_nodes_per_face = Dict([(0,nothing),(2,2),(3,3),(4,4)])

struct NamedZone
    ID::Int
    type::String
    name::String
end

struct NodeSection
    zoneID::Int
    firstIdx::Int
    lastIdx::Int
    coords::Union{Nothing, Array{Real, 2}}
end

struct CellSection
    zoneID::Int
    firstIdx::Int
    lastIdx::Int
    type::Int
    cellType::Int
    mixedCellTypes::Union{Nothing, Array{Int, 1}}
end

struct Face
    type::Int
    nodeIDs::Array{Int,1}
    cr::Int
    cl::Int
end

struct FaceSection
    zoneID::Int
    firstIdx::Int
    lastIdx::Int
    type::Int
    faceType::Int
    FaceTypes::Union{Nothing, Array{Face, 1}}
end

# TODO: implent and test other xf-types ...

struct gridSectionStruct
    gridIDx::UInt8
    content::Union{String, Nothing}
end

struct TGrid
    header::String
    dimensions::Int
    cellSection::Array{CellSection,1}
    faceSection::Array{FaceSection,1}
    nodeSection::Array{NodeSection,1}
    zone::Array{NamedZone,1}
end

max_stars_repo_name: c-schubert/linearMeshIO

max_issues_repo_name: c-schubert/linearMeshIO

max_forks_repo_name: c-schubert/linearMeshIO


text:
include("utils.jl")
include("../Applications/KP.jl")

n = 20
k = 2
inst_id = 1
density = 0.5
M = 4

const heuristicmode = "HEURCG"
#const heuristicmode = "DUAL"
const formulation = "AGG"
#const formulation = "CONF"
const constraint_uncertainty = false

instance = "$(inst_id)_$(n)_$(M)_$(density).txt"
data = read_data(instance,constraint_uncertainty)
K = 1:k
N = 1:data.n

# The following executions avoid warming-up issues
const STARTTIME = time()
@info "Executing functions to warm up ..."
oldstd = stdout
redirect_stdout(open("/dev/null", "w"))
scenario_generation(false,false)
redirect_stdout(oldstd) # recover original stdout
@info "warm-up finished!"

# Now comes the instance that is tested, passed by the arguments
n = parse(Int, ARGS[1])
k = parse(Int, ARGS[2])
inst_id = parse(Int, ARGS[3])
constraint_uncertainty = parse(Bool, ARGS[4])
formulation = ARGS[5]
M = 4

instance = "$(inst_id)_$(n)_$(M)_$(density).txt"
if constraint_uncertainty
    suffix = "_WU_"
else
    suffix = "_"
end
instance_being_solved = "Solving instance $inst_id with k=$k and M=$M and density=$density and ϵ=$ϵ "
if constraint_uncertainty
    instance_being_solved *="and constraint uncertainty ..."
else
    instance_being_solved *="and deterministic constraint ..."
end
@info instance_being_solved
data = read_data(instance,constraint_uncertainty)
K = 1:k
N = 1:data.n

outputfilename = "/../results/KP$(suffix)scenario_generation_$(formulation)_$(n)_$(k).csv"
reset_timer!(to::TimerOutput)
STARTTIME = time()
runtime = @elapsed lb, ub, solution, gap, NUMBER_GENSC, NUMITER, NUMBER_GENSOL = scenario_generation(false,true)
outputfile = open(string(dirname(@__FILE__),outputfilename), "a")
output = "$n $inst_id $k $NUMBER_GENSC $NUMBER_GENSOL $ϵ $runtime $NUMITER $lb $ub $gap $(TimerOutputs.time(to["heuristic"])/10^9) $(TimerOutputs.time(to["gen_sol"])/10^9)"
println(output)
println()
println(outputfile, output)
close(outputfile)
show(to)

# outputfilename = "/../results/KP$(suffix)rc_generation_$(n)_$(k)_$(density).csv"
# reset_timer!(to::TimerOutput)
# STARTTIME = time()
# runtime = @elapsed lb, ub, incumbent, gap, it = rcg()
# outputfile = open(string(dirname(@__FILE__),outputfilename), "a")
# output = "$n $inst_id $k $M $runtime $it $lb $ub $gap"
# println(output)
# println()
# println(outputfile, output)
# close(outputfile)
# show(to)

max_stars_repo_name: mjposs/min-max-min

max_issues_repo_name: mjposs/min-max-min

max_forks_repo_name: mjposs/min-max-min


text:
import Base: ==, checkbounds, collect, eltype, filter, first, firstindex, getindex,
    in, isempty, iterate, keys, last, lastindex, length, map, pairs, show,
    IteratorEltype, IteratorSize

export NamingContext,
    freenames,
    freshnames,
    pushfirst,
    split

"""
    UniqueIndices(start, log_skip)

Immutable sequence of indices which can be kept linear in a context by splitting:

````
[1, 2, 3, ...]
-> ([1, 3, 5, ...], [2, 4, 6, ...])
-> (([1, 5, 9, ...], [3, 7, 11, ...]), [2, 4, 6, ...])
-> etc.
````
"""
struct UniqueIndices
    start::Int
    log_skip::Int
end

nextindex(u::UniqueIndices) = UniqueIndices(u.start + 2^u.log_skip, u.log_skip)

function splitindex(u::UniqueIndices)
    left = UniqueIndices(u.start, u.log_skip + 1)
    right = UniqueIndices(u.start + 2^u.log_skip, u.log_skip + 1)
    (left, right)
end

"""
    NamingContext([freenames]; namehint = :x)

Keeps track of free variables and their indices, and can generate fresh names based on them.
Fresh names are generated based on `namehint` and different from all `freenames`.

This behaves mostly like an immutable array of `Name`s.
"""
struct NamingContext
    freenames::Vector{Name}
    namehint::Name
    indices::UniqueIndices
end

NamingContext(freenames = Name[]; namehint = :x) =
    NamingContext(reverse!(collect(freenames)), namehint, UniqueIndices(1, 0))

freenames(Γ::NamingContext) = Γ.freenames

==(Γ₁::NamingContext, Γ₂::NamingContext) = Γ₁.freenames == Γ₂.freenames
checkbounds(::Type{Bool}, Γ::NamingContext, i::Int) = firstindex(Γ) ≤ i ≤ lastindex(Γ)
# cat(Γ₁::NamingContext, Γs::NamingContext...) =
#     NamingContext(cat(Γ₁.freenames, freenames.(Γs)...))
collect(Γ::NamingContext) = Γ.freenames
# copy(Γ::NamingContext) = NamingContext(copy(Γ.freenames))
eltype(::Type{NamingContext}) = Symbol
first(Γ::NamingContext) = firstindex(Γ.freenames)
firstindex(Γ::NamingContext) = firstindex(Γ.freenames)
filter(f, Γ::NamingContext) = NamingContext(filter(f, Γ.freenames))
getindex(Γ::NamingContext, key) = getindex(Γ.freenames, reverseind(Γ.freenames, key))
in(name, Γ::NamingContext) = name in Γ.freenames
isempty(Γ::NamingContext) = isempty(Γ.freenames)
keys(Γ::NamingContext) = keys(Γ.freenames)
last(Γ::NamingContext) = last(Γ.freenames)
lastindex(Γ::NamingContext) = lastindex(Γ.freenames)
length(Γ::NamingContext) = length(Γ.freenames)
map(f, Γs::NamingContext...) = NamingContext(map(f, freenames.(Γs)...))
pairs(Γ::NamingContext) = zip(keys(Γ.freenames), Iterators.reverse(Γ.freenames))

function show(io::IO, Γ::NamingContext)
    print(io, "NamingContext([")
    join(io, reverse(sprint.(show, Γ.freenames)), ", ")
    print(io, "])")
end

struct FreshNames
    Γ::NamingContext
end

iterate(iter::FreshNames) = iterate(iter, iter.Γ.indices)
function iterate(iter::FreshNames, state)
    while true
        newname = Symbol(iter.Γ.namehint, state.start)
        if newname ∉ iter.Γ.freenames
            return newname, nextindex(state)
        end
        state = nextindex(state)
    end
end

IteratorSize(::Type{FreshNames}) = Base.IsInfinite()
IteratorEltype(::Type{FreshNames}) = Base.HasEltype()
eltype(::Type{FreshNames}) = Symbol

"""
    freshnames(Γ)

An (infinite) iterator for the fresh names that can be generated from `Γ`.
"""
freshnames(Γ::NamingContext) = FreshNames(Γ)

"""
    pushfirst(Γ, names...)

Construct a new `NamingContext` by adding `names` as free names to `Γ`.  The new names will be added
to the beginning, thus having the lowest index in the new context.
"""
pushfirst(Γ::NamingContext, names...) =
    NamingContext([Γ.freenames; names...], Γ.namehint, Γ.indices)

"""
    freshname(Γ)

Generate a new name fresh in `Γ`, and a new context `Γ′`, which is `Γ` updated by this name.
Names are constructed based on `Γ.namehint`.
"""
function freshname(Γ::NamingContext)
    (newname, newindex) = iterate(FreshNames(Γ))
    Γ′ = NamingContext([Γ.freenames; newname], Γ.namehint, newindex)
    (newname, Γ′)
end

"""
    split(Γ)

Separate `Γ` into two contexts containing the same free variables, but ensuring that different fresh
names are generated by them.
"""
function split(Γ::NamingContext)
    l, r = splitindex(Γ.indices)
    NamingContext(Γ.freenames, Γ.namehint, l), NamingContext(Γ.freenames, Γ.namehint, r)
end


text:
using SparseArrays: SparseMatrixCSC
using OSQP

# undo https://github.com/oxfordcontrol/OSQP.jl/blob/2e1dba7be28d7e25bac0427ac34a5f65651438bf/src/types.jl#L33

function unsafe_convert(::Type{SparseMatrixCSC}, c::OSQP.Ccsc)
  m = c.m
  n = c.n
  nzmax = c.nzmax
  nzval = [unsafe_load(c.x, i) for i=1:nzmax]
  rowval = [unsafe_load(c.i, i) for i=1:nzmax] .+ 1
  colptr = [unsafe_load(c.p, i) for i=1:(n+1)] .+ 1
  SparseMatrixCSC(m, n, colptr, rowval, nzval)
end

function unsafe_load_data(d::OSQP.Data)
  A = unsafe_convert(SparseMatrixCSC, unsafe_load(d.A))
  P = unsafe_convert(SparseMatrixCSC, unsafe_load(d.P))
  l = [unsafe_load(d.l, i) for i=1:d.m]
  u = [unsafe_load(d.u, i) for i=1:d.m]
  q = [unsafe_load(d.u, i) for i=1:d.n]

  (A=A, P=P, l=l, u=u, q=q)
end

unsafe_load_data(m::OSQP.Model) = unsafe_load_data(unsafe_load(unsafe_load(r_m.workspace).data))

using GeometryReferenceGenericSolver
using GeometryTypes

xy_plane = 1e6 .* (Point(-1,0,0), Point(1,0,0), Point(0, 1, 0))

p = Point(15099.411397006435, 63329.17689456945, 489.2398433903766)
p = Point(56270.48118496456, 53948.45819773468, 34.889982642960504)
p = Point(37517.76357686068, 84362.38878880453, 233.50717418757495)

c1 = xy_plane
c2 = (p, )

using GeometryReferenceGenericSolver
#r_p = GeometryReferenceGenericSolver.find_minimum_distance(c1, c2).model.optimizer.inner
#r_m = find_minimum_distance2(c1, c2).model

max_stars_repo_name: goretkin/GeometryReferenceGenericSolver

max_issues_repo_name: goretkin/GeometryReferenceGenericSolver

max_forks_repo_name: goretkin/GeometryReferenceGenericSolver


text:
abstract type HiddenLayer{T <: Number} end

mutable struct MutableHiddenLayer{T <: Number} <: HiddenLayer{T}
    n_neurons::Int
    neurons::Vector{Neurons}
    HH::Union{Matrix{T}, Nothing}
    TH::Union{Matrix{T}, Nothing}
    initialized::Bool

    function MutableHiddenLayer(::Type{T}) where {T <: Number}
        n_neurons = 0
        neurons = Vector{Neurons}()
        HH = nothing
        TH = nothing
        initialized = false
        new{T}(n_neurons, neurons, HH, TH, initialized)
    end
end

struct ImmutableHiddenLayer{T <: Number} <: HiddenLayer{T}
    n_neurons::Int
    neurons::Vector{Neurons}
end

function is_initialized(hidden_layer::T) where {T <: MutableHiddenLayer}
    hidden_layer.initialized
end

function initialize!(hidden_layer::T₁,
                     n_outputs::Int) where {T₁ <: MutableHiddenLayer}

    param(::MutableHiddenLayer{T}) where {T} = T

    L = hidden_layer.n_neurons
    Q = n_outputs

    hidden_layer.HH = zeros(param(hidden_layer), L + 1, L + 1)
    hidden_layer.TH = zeros(param(hidden_layer), Q, L + 1)
    hidden_layer.initialized = true

    hidden_layer
end

function make_immutable(hidden_layer::MutableHiddenLayer)
    param(::HiddenLayer{T}) where {T} = T

    n_neurons = hidden_layer.n_neurons
    neurons = deepcopy(hidden_layer.neurons)
    ImmutableHiddenLayer{param(hidden_layer)}(n_neurons, neurons)
end

# Get H, the nonlinear representation of the sample array X. First, we project
# X onto R random hyperplanes (where R is the number of neurons in the hidden
# layer) by taking the product of the input weights W and X. Then, we apply
# some nonlinear operation f (the activation function for each neuron) to each
# element of the projection.
function project(hidden_layer::T₁,
                 samples::T₂) where {T₁ <: HiddenLayer,
                                     T₂ <: AbstractMatrix}
    param(::HiddenLayer{T}) where {T} = T

    X = samples
    L = hidden_layer.n_neurons
    N = last(size(X))

    H = ones(param(hidden_layer), L + 1, N)
    i₀ = 1

    for (i, neurons) in enumerate(hidden_layer.neurons)
        L₀ = neurons.n_neurons
        f = neurons.activation_function
        W = neurons.weights

        i₁ = i₀ + L₀ - 1

        H[i₀:i₁,:] .= f.(W * X)

        i₀ = i₁ + 1
    end

    H
end