text:
module poly2D

using QHull, Statistics, LinearAlgebra, Clipper, Devices, PyPlot, LazySets

"""
"""

function  checkConvex(V)
# Given a set of points determine if they form a convex polygon

    px=V[:,1];
    py=V[:,2];

    isConvex = false;

    numPoints = length(px);
    if numPoints < 4
        isConvex = true;
        return isConvex
    end

    # can determine if the polygon is convex based on the direction the angles
    # turn.  If all angles are the same direction, it is convex.
    v1 = [px[1] - px[end], py[1] - py[end]];
    v2 = [px[2] - px[1], py[2] - py[1]];
    signPoly = sign(det([v1'; v2']));

    # check subsequent vertices
    for k = 2:numPoints-1
        v1 = v2;
        v2 = [px[k+1] - px[k], py[k+1] - py[k]];
        curr_signPoly = sign(det([v1'; v2']));
        # check that the signs match
        if curr_signPoly != signPoly
            isConvex = false;
            return isConvex
        end
    end

    # check the last vectors
    v1 = v2;
    v2 = [px[1] - px[end], py[1] - py[end]];
    curr_signPoly = sign(det([v1'; v2']));
    if curr_signPoly != signPoly
        isConvex = false;
    else
        isConvex = true;
    end

    return isConvex
end

function createLine(p1,p2)

    # first input parameter is first point, and second input is the
    # second point.
    line = [p1[1] p1[2] p2[1]-p1[1] p2[2]-p1[2]];

end

function distanceMat(points1,points2)

    numPoints1=size(points1,1)
    numPoints2=size(points2,1)

    distMat=zeros(numPoints1,numPoints2)
    for i=1:numPoints1
        point_i=points1[i,:]
        for j=1:numPoints2
            point_j=points2[j,:]
            distMat[i,j]=sqrt(sum((point_i-point_j).^2))
        end
    end

    return distMat

end

function expandPolygon(V, dist)
    numVertices = size(V,1)
    poly = Devices.Polygon([Devices.Point(V[i,1],V[i,2]) for i=1:numVertices])
    poly_ = offset(poly, dist)

    poly__ = poly_[1].p
    numVertices_ = size(poly__,1)
    V_ = zeros(numVertices_,2)
    for i=1:numVertices_
        V_[i,1] = poly__[i][1]
        V_[i,2] = poly__[i][2]
    end

    V__ = poly2D.expandPolygonSides(V, dist, collect(1:numVertices))

    V_out = zeros(size(V__,1),2)
    posVec = zeros(size(V__,1),1)
    for i = 1:size(V__,1)
        for j = 1:size(V_,1)
            if sqrt(sum((V__[i,:] - V_[j,:]).^2)) <= 1
                V_out[i,:] =  V_[j,:]
                posVec[i] = j
            end
        end
    end
    pos = setdiff(collect(1:size(V_,1)), unique(posVec))

    if length(pos) >= 0
        V_out[V_out[:,1] .== 0,: ] .= V_[pos,:]
    end

    return V_out

end

function expandPolygonSide(poly_in, dist, edge)
    poly=copy(poly_in);

    # eventually copy first point at the end to ensure closed polygon
    if sum(poly[end, :] == poly[1,:]) != 2
        poly = [poly; poly[1,:]'];
    end

    # number of vertices of the polygon
    N = size(poly, 1)-1;

    # find lines parallel to polygon edges located at distance DIST
    lines = zeros(N, 4);
    for i = 1:N
        side = createLine(poly[i,:], poly[i+1,:]);
        if i==edge
            lines[i, 1:4] = parallelLine(side, dist);
        else
            lines[i, 1:4] = side;
        end
    end

    # compute intersection points of consecutive lines
    lines = [lines;lines[1,:]'];
    poly_out = zeros(N, 2);
    for i = 1:N
        pa1=lines[i,:];
        pa2=lines[i+1,:];
        poly_aux = intersectLines(pa1, pa2);
        poly_out[i,1:2] = poly_aux;
    end

    poly_out = [poly_out[end,:]'; poly_out[1:N-1,:]];

end

function expandPolygonSides(poly, sep, edges)

    numEdgesTotales = size(poly,1)
    conjuntoEdgesTotales = 1:numEdgesTotales
    conjuntoEdgesNoSeleccionados = setdiff(conjuntoEdgesTotales,edges)
    numEdges = length(edges);

    poly_ = copy(poly);

    if numEdges >= 1
        for i in edges
            poly_aux = expandPolygonSide(poly_, sep, i);
            poly_ = copy(poly_aux);
        end
        poly_out = poly_;
    end

    return poly_out
end

function inpoly_mat(points, poly, edge)

    nvrt = size(points,1);
    nnod = size(poly,1);
    nedg = size(edge,1);

    stat = falses(nvrt,1);
    bnds = falses(nvrt,1);

#----------------------------------- loop over polygon edges
    for epos = 1 : size(edge,1)
        inod = edge[epos,1] ;
        jnod = edge[epos,2] ;

    #------------------------------- calc. edge bounding-box
        yone = poly[inod,2];
        ytwo = poly[jnod,2];
        xone = poly[inod,1];
        xtwo = poly[jnod,1];

        xmin = min(xone, xtwo);
        xmax = max(xone, xtwo);

        ymin = yone;
        ymax = ytwo;

        ydel = ytwo - yone;
        xdel = xtwo - xone;

    #------------------------------- find top VERT(:,2)<YONE
        ilow = 1;
        iupp = nvrt;

        while (ilow < iupp - 1)    # binary search
            imid = ilow + Int(floor((iupp-ilow) / 2));

            if (points[imid,2] < ymin)
                ilow = imid;
            else
                iupp = imid ;
            end
        end

        if (points[ilow,2] >= ymin)
            ilow = ilow - 1 ;
        end

    #------------------------------- calc. edge-intersection
        for jpos = ilow+1:nvrt

            if (bnds[jpos])
                 continue;
            end

            xpos = points[jpos,1];
            ypos = points[jpos,2];

            if (ypos <= ymax)
                if (xpos >= xmin)
                    if (xpos <= xmax)

                #------------------- compute crossing number
                    mul1 = ydel * (xpos - xone) ;
                    mul2 = xdel * (ypos - yone) ;

                        if (ypos == yone &&  xpos == xone )

                    #------------------- BNDS -- match about ONE
                            bnds[jpos]= true ;
                            stat[jpos]= true ;

                        elseif (ypos == ytwo &&  xpos == xtwo )

                    #------------------- BNDS -- match about TWO
                            bnds[jpos]= true ;
                            stat[jpos]= true ;

                        elseif (mul1 < mul2)

                            if (ypos >= yone && ypos <  ytwo)

                        #------------------- advance crossing number
                                stat[jpos] = ~stat[jpos] ;
                            end
                        end

                    end
                else

                    if (ypos >= yone && ypos <  ytwo)
                #------------------- advance crossing number
                        stat[jpos] = !stat[jpos] ;
                    end

                end
            else

                break ;            # done -- due to the sort

            end

        end

    end

    return stat, bnds
end

function inpoly(points, poly)

    nnod = size(poly,1);
    nvrt = size(points,1);

    edge = [(1:nnod-1) (2:nnod); [nnod 1]];
#----------------------------------- prune points using bbox
    mask =  (points[:,1] .>= minimum(poly[:,1])) .* (points[:,1] .<= maximum(poly[:,1])) .* (points[:,2] .>= minimum(poly[:,2])) .* (points[:,2] .<= maximum(poly[:,2]));

    vert = points[mask, :];

#----------------------------------- sort points via y-value
    swap = poly[edge[:,2],2] .< poly[edge[:,1],2];

    edge[swap,[1,2]] = edge[swap,[2,1]];

    ivec = sortperm(points[:,2]);
    points = points[ivec,:];

    #-- call the native m-code version

    stat, bnds = inpoly_mat(points, poly, edge);

    stat[ivec] = stat;
    bnds[ivec] = bnds;

    #STAT[mask] = stat;
    #BNDS[mask] = bnds;

    STAT = stat;
    BNDS = bnds;

    return STAT, BNDS
end

function intersectEdges(edge1, edge2)

    x1_ini  = edge1[1];
    y1_ini  = edge1[2];
    x1_fin  = edge1[3];
    y1_fin  = edge1[4];
    dx1 = x1_fin - x1_ini;
    dy1 = y1_fin - y1_ini;
    m1 = dy1 / dx1;

    x2_ini  = edge2[1];
    y2_ini  = edge2[2];
    x2_fin  = edge2[3];
    y2_fin  = edge2[4];
    dx2 = x2_fin - x2_ini;
    dy2 = y2_fin - y2_ini;
    m2 = dy2 / dx2;

    min_e1_x = min(x1_ini,x1_fin)
    min_e2_x = min(x2_ini,x2_fin)
    max_e1_x = max(x1_ini,x1_fin)
    max_e2_x = max(x2_ini,x2_fin)

    min_e1_y = min(y1_ini,y1_fin)
    min_e2_y = min(y2_ini,y2_fin)
    max_e1_y = max(y1_ini,y1_fin)
    max_e2_y = max(y2_ini,y2_fin)

# tolerance for precision
    tol = 1e-3; #1e-14;

# initialize result array
    x0  = 0;
    y0  = 0;

# indices of parallel edges
    par = abs(m1 - m2) < tol;

# Parallel edges have no intersection -> return [NaN NaN]
    if par
        x0 = NaN;
        y0 = NaN;
    end

# Process non parallel cases

    # compute intersection points of supporting lines
    delta = dx2 * dy1 - dx1 * dy2;
    x0 = ((y2_ini - y1_ini) * dx1 * dx2 + x1_ini * dy1 * dx2 - x2_ini * dy2 * dx1) / delta;
    y0 = ((x2_ini - x1_ini) * dy1 * dy2 + y1_ini * dx1 * dy2 - y2_ini * dx2 * dy1) / -delta;

    if  x0 < min_e1_x-tol || x0 < min_e2_x-tol || x0 > max_e1_x+tol || x0 > max_e2_x+tol ||
        y0 < min_e1_y-tol || y0 < min_e2_y-tol || y0 > max_e1_y+tol || y0 > max_e2_y+tol
        point = [NaN NaN];
#    elseif (x0-x1_ini)^2 + (y0-y1_ini)^2 <= 1  || (x0-x1_fin)^2 + (y0-y1_fin)^2 <= 1 ||  (x0-x2_ini)^2 + (y0-y2_ini)^2 <= 1  || (x0-x2_fin)^2 + (y0-y2_fin)^2  <= 1
#        point = [NaN NaN];
    else
        point = [x0 y0];
    end

end

function expantLine(line, f)
    #intersectLines([x_0, y_0, x_1, y_1],2)
    x_0 = line[1]
    y_0 = line[2]
    x_1 = line[3]
    y_1 = line[4]
    m = (y_1 - y_0) / (x_1 - x_0)
    d = sqrt((y_1-y_0)^2+(x_1-x_0)^2)
    df = d*f
    y=x*m + a
end

function intersectLines(line1, line2)
    # intersectLines([x1_0, y1_0, x1_1, y1_1], [x2_0, y2_0, x2_1, y2_1])

    # extract tolerance
    tol = 1e-14;

    # Check parallel and colinear lines

    # coordinate differences of origin points
    dx = line2[1] - line1[1];
    dy = line2[2] - line1[2];

    # indices of parallel lines
    denom = line1[3] .* line2[4] - line2[3] .* line1[4];
    par = abs(denom) < tol;

    # initialize result array
    x0 = 0;
    y0 = 0;

    # initialize result for parallel lines
    if par
        x0 = Inf;
        y0 = Inf;
        point = [x0 y0];
        return;
    end

    # Extract coordinates of itnersecting lines

    # extract base coordinates of first lines
    x1 =  line1[1];
    y1 =  line1[2];
    dx1 = line1[3];
    dy1 = line1[4];

    # extract base coordinates of second lines
    x2 =  line2[1];
    y2 =  line2[2];
    dx2 = line2[3];
    dy2 = line2[4];

    # re-compute coordinate differences of origin points
    dx = line2[1] - line1[1];
    dy = line2[2] - line1[2];

    # Compute intersection points

    x0 = (x2 .* dy2 .* dx1 - dy .* dx1 .* dx2 - x1 .* dy1 .* dx2) ./ denom ;
    y0 = (dx .* dy1 .* dy2 + y1 .* dx1 .* dy2 - y2 .* dx2 .* dy1) ./ denom ;

    # concatenate result
    point = [x0 y0];
end

function intersectPoly2d(V1, V2)

    V1 = [V1; V1[1,:]'];
    V2 = [V2; V2[1,:]'];
    N1 = size(V1, 1)
    N2 = size(V2, 1)

    # Loop over segments of V1
    i=1
    X=[NaN NaN NaN NaN]
    for n1 = 1:N1-1
        for n2 = 1:N2-1
            e1 = [V1[n1,:]' V1[n1+1,:]']
            e2 = [V2[n2,:]' V2[n2+1,:]']
            pInt = intersectEdges(e1, e2)
            if abs(pInt[1])>0
                if i==1
                    X[:] = [n1 n2 pInt]
                else
                    X = [X; [n1 n2 pInt]]
                end
                i=i+1
            end
        end
    end

    return X

end

function isPointInPolygon(points, V_)

    V=copy(V_);

    numPointsToCheck=size(points,1);
    A, b, c = vert2con(V);

    flagInPoly=[];
    for i=1:numPointsToCheck
        p_i=points[i,:];
        flagInPoly_=A*p_i.<b;
        flagInPoly_i=prod(flagInPoly_);
        push!(flagInPoly,flagInPoly_i);

    end

    return flagInPoly
end

function findNonConvexVert(Vin)

    V = copy(Vin);

    V_ = poly2D.expandPolygon(V, -.5)
    V_ch = poly2D.convHull(V_)

    flagNonConvex = isPointInPolygon(V, V_ch)

    idNonConv = collect(1:size(V,1))
    idNonConv = idNonConv[flagNonConvex .== true]

    return idNonConv
end

function isRectInPoly(rect,poly)

    pos_x = rect[1];
    pos_y = rect[2];
    largo_x = rect[3];
    largo_y = rect[4];
    theta = rect[5];

    Rmat = rotationMatrix(theta);
    cr = [pos_x; pos_y];
    p1 = (Rmat * ([pos_x; pos_y] - cr) + cr)';
    p2 = (Rmat * ([pos_x + largo_x; pos_y] - cr) + cr)';
    p3 = (Rmat * ([pos_x + largo_x; pos_y + largo_y] - cr) + cr)';
    p4 = (Rmat * ([pos_x; pos_y + largo_y] - cr) + cr)';

    if collect(inpoly(p1, poly)[1])[1] && collect(inpoly(p2, poly)[1])[1] && collect(inpoly(p3, poly)[1])[1] &&
        collect(inpoly(p4, poly)[1])[1] && !(intersectPoly2d([p1; p2; p3; p4], poly)[1,1]>0)
        return true
    else
        return false
    end

end

function lineAngle(line)

    # angle of one line with horizontal
    theta = mod(atan(line[4], line[3]) + 2*pi, 2*pi);

end

function parallelLine(line, dist)

    # use a distance. Compute position of point located at distance DIST on
    # the line orthogonal to the first one.
    point = pointOnLine([line[1] line[2] line[4] -line[3]], dist);

    res = [point' (line[3:4])'];
end

function pointOnLine(line, pos)

    angle = lineAngle(line);
    point = [line[1] + pos .* cos(angle), line[2] + pos .* sin(angle)];

end

function polyArea(V)

    if V != []
      X=V[:,1];
      Y=V[:,2];
      numPoints=size(X,1);

      area = 0;   # Accumulates area
      j = numPoints;

      for i = 1:numPoints
        area = area + (X[j]+X[i])*(Y[j]-Y[i]);
        j = i;  #j is previous vertex to i
      end

      return out=abs(area/2);
    else
      return out=0;
    end
end

function randomPointInPolygon(poly, nPts)

    # polygon extreme coordinates
    xmin = minimum(poly[:,1]);  xmax = maximum(poly[:,1]);
    ymin = minimum(poly[:,2]);  ymax = maximum(poly[:,2]);

    # compute size of box
    boxSizeX = xmax - xmin;
    boxSizeY = ymax - ymin;

    # allocate memory for result
    points = zeros(nPts, 2);

    # contains indices of remaining points to process
    ind = 1;

    # iterate until all points have been sampled within the polygon
    while ind<=nPts
        x = rand() .* boxSizeX .+ xmin;
        y = rand() .* boxSizeY .+ ymin;
        points[ind, 1] = x;
        points[ind, 2] = y;
        #isPointInPolygon
        if collect(inpoly(points[ind, :]', poly)[1])[1]
            ind=ind+1;
        end
    end

    return points

end

function randomPointInExtrude(poly, alt, nPts)

    # polygon extreme coordinates
    xmin = minimum(poly[:,1]);  xmax = maximum(poly[:,1]);
    ymin = minimum(poly[:,2]);  ymax = maximum(poly[:,2]);
    zmin = minimum(0);  zmax = maximum(alt);

    # compute size of box
    boxSizeX = xmax - xmin;
    boxSizeY = ymax - ymin;
    boxSizeZ = zmax - zmin;

    # allocate memory for result
    points = zeros(nPts, 3);

    # contains indices of remaining points to process
    ind = 1;

    # iterate until all points have been sampled within the polygon
    while ind<=nPts
        x = rand() .* boxSizeX .+ xmin;
        y = rand() .* boxSizeY .+ ymin;
        z = rand() .* boxSizeZ .+ zmin;
        points[ind, 1] = x;
        points[ind, 2] = y;
        points[ind, 3] = z;
        #isPointInPolygon
        if collect(inpoly(points[ind, 1:2]', poly)[1])[1]
            ind=ind+1;
        end
    end

    return points

end

function rotationMatrix(theta)
    R=[cos(theta) -sin(theta); sin(theta) cos(theta)];
end

function convHull(V)

    n = size(V,1)
    v = [[V[i,1], V[i,2]] for i=1:n]
    ch_v = LazySets.convex_hull(v)

    V_out = [0 0]; for i=1:length(ch_v); V_out=[V_out; ch_v[i]']; end; V_out = V_out[2:end,:]

    return V_out
end

function simplicesConvHull(V_)
    V = copy(V_);
    k_ = QHull.chull(V);
    k = k_.simplices;
end

function verticesConvHull(V_)
    V = copy(V_);
    k_ = QHull.chull(V);
    k = k_.vertices;
    n = length(k);
    kk = vcat(k,k);
    K = [kk[n:2*n-1] kk[1:n]];
end

function vert2con(V_)

    V = copy(V_);

    k = verticesConvHull(V);
    c = mean!([1. 1.], V[unique(k),:]);
    V = V-repeat(c,size(V,1),1);
    A = NaN*zeros(size(k,1),size(V,2));
    rc=0;
    for ix = 1:size(k,1)
        F = V[k[ix,:],:];
        #if rank(F, 1e-5) == size(F,1)
            rc=rc+1;
            A[rc,:]=F\ones(size(F,1),1);
        #end
    end
    A=A[1:rc,:];
    b=ones(size(A,1),1);
    b=b+A*c';

    return A,b,c

end

function plotScatter3d(x,y,z, fig=nothing, ax=nothing, ax_mat=nothing, color="red", marker="o", alpha=.3)

    plt = pyimport("matplotlib.pyplot")
    pyimport("mpl_toolkits.mplot3d")

    if fig === nothing
        pygui(true)

        fig = plt.figure()
        ax = fig[:add_subplot](111, projection = "3d")

        min_ax_x = minimum(x)
        min_ax_y = minimum(y)
        min_ax_z = minimum(z)
        max_ax_x = maximum(x)
        max_ax_y = maximum(y)
        max_ax_z = maximum(z)
    else
        min_ax_x = ax_mat[1,1]
        min_ax_y = ax_mat[2,1]
        min_ax_z = ax_mat[3,1]
        max_ax_x = ax_mat[1,2]
        max_ax_y = ax_mat[2,2]
        max_ax_z = ax_mat[3,2]
        min_ax_x = minimum([min_ax_x minimum(x)])
        min_ax_y = minimum([min_ax_y minimum(y)])
        min_ax_z = minimum([min_ax_z minimum(z)])
        max_ax_x = maximum([max_ax_x maximum(x)])
        max_ax_y = maximum([max_ax_y maximum(y)])
        max_ax_z = maximum([max_ax_z maximum(z)])

    end

    ax.set_xlim(min_ax_x, max_ax_x)
    ax.set_ylim(min_ax_y, max_ax_y)
    ax.set_zlim(min_ax_z, max_ax_z)

    ax[:scatter3D](x, y, z, color=color, marker=marker, alpha=alpha)

    ax.set_xlabel("X")
    ax.set_ylabel("Y")
    ax.set_zlabel("Z")

    plt.show()

    ax_mat = [min_ax_x max_ax_x; min_ax_y max_ax_y; min_ax_z max_ax_z]

    return fig, ax, ax_mat

end

function angulosLados(V)
    x = V[:,1]
    x = [x;x[1]]
    y = V[:,2]
    y = [y;y[1]]
    dx = x[2:end] .- x[1:end - 1]
    dy = y[2:end] .- y[1:end - 1]

    theta = zeros(size(V, 1), 1)
    for i = 1:size(V, 1)
        if dy[i] > 0 && dx[i] > 00
            theta[i] = atan(dy[i] / dx[i])
        elseif dy[i] > 0 && dx[i] < 0
            theta[i] = (90 + 180 / pi * atan(abs(dx[i]) / dy[i]) ) * pi / 180
        elseif dy[i] < 0 && dx[i] < 0
            theta[i] = (180 + 180 / pi * atan(abs(dy[i]) / abs(dx[i])) ) * pi / 180
        elseif dy[i] < 0 && dx[i] > 0
            theta[i] = (270 + 180 / pi * atan(dx[i] / abs(dy[i])) ) * pi / 180
        else
            theta[i] = 0
        end
    end
    return theta
end

function rectangle(p, w, h, theta)
    R = poly2D.rotationMatrix(theta);
    x0 = p[1]
    y0 = p[2]
    cr = [x0; y0]
    p1 = R*([x0; y0]-cr) + cr
    p2 = R*([x0+w; y0]-cr) + cr
    p3 = R*([x0+w; y0+h]-cr) + cr
    p4 = R*([x0; y0+h]-cr) + cr

    V = [p1';p2';p3';p4']

    return V
end

function distPointLine(q, p1, p2)
    dist = sqrt(sum( ((q-p1) - ((q-p1)'*(p2-p1)) / ((p2-p1)'*(p2-p1)) * (p2-p1) ).^2))
    return dist
end

export  checkConvex, verticesConvHull, convHull, createLine, distanceMat, expandPolygonSide, expandPolygonSides, inpoly, inpoly_mat,
        intersectEdges, intersectLines, intersectPoly2d, isPointInPolygon, isRectInPoly, lineAngle, parallelLine,
        pointOnLine, polyArea, randomPointInPolygon, rotationMatrix, vert2con, plotScatter3d, findNonConvexVert,
        angulosLados, rectangle, distPointLine
end


text:
function laplace_filter(x::Matrix{Tv}) where {Tv<:AbstractFloat}

    (n1, n2) = size(x)
    x1 = zeros(Tv, n1, n2)
    for i2 = 2 : n2-1
        for i1 = 2 : n1-1
            x1[i1,i2] = (8.0 * x[i1,i2] - x[i1-1,i2-1] - x[i1,i2-1] - x[i1+1,i2-1]
                                        - x[i1-1,i2  ]              - x[i1+1,i2  ]
                                        - x[i1-1,i2+1] - x[i1,i2+1] - x[i1+1,i2+1])
        end
    end
    return x1
end

# function laplace_filter(img::Matrix{Tv}, iflag::Int64) where {Tv<:AbstractFloat}
#     (m, n) = size(img)
#     # d = zeros(eltype(img), m, n)
#     # for ix = 2 : n-1
#     #     for iz = 2 : m-1
#     #         d[iz, ix] = 4*img[iz,ix]-img[iz-1,ix]-img[iz+1,ix]-img[iz,ix-1]-img[iz,ix+1]
#     #     end
#     # end
#     a = -1.0 * ones(m-1); a[end] = 0.0
#     b =    2 * ones(m  ); b[1]   = 1.0; b[end] = 1.0
#     c = -1.0 * ones(m-1); c[1]   = 0.0;
#     V = spdiagm((a,b,c), (-1,0,1), m, m)
#
#     a = -1.0 * ones(n-1); a[end] = 0.0
#     b =    2 * ones(n  ); b[1]   = 1.0; b[end] = 1.0
#     c = -1.0 * ones(n-1); c[1]   = 0.0;
#     H = spdiagm((a,b,c), (-1,0,1), n, n)
#
#     if iflag == 1
#        r = V * img + img * H'
#     elseif iflag == 2
#        r = V'* img + img * H
#     end
#
#     return r
# end

ext: jl


text:
function solution()
    sum = 0
    for x in 3:999
        if  mod(x, 3) == 0 || mod (x, 5) == 0
            sum += x
        end
    end
    println(sum)
end

max_stars_repo_name: patrickbryan/Project-Euler-Solutions

max_issues_repo_name: patrickbryan/Project-Euler-Solutions

max_forks_repo_name: patrickbryan/Project-Euler-Solutions


text:
mutable struct Person
    name:: String
    height::Float64
end

fieldnames(Person) #> (:name, :height)
Person.types  #> svec(String, Float64)

code_lowered(+, (Int, Int))
# 1-element Array{Core.CodeInfo,1}:
#  CodeInfo(
# 53 1 ─ %1 = Base.add_int(%%x, %%y)::Any                                     │
#    └──      return %1                                                       │
# )

code_typed(+, (Int, Int))
# 1-element Array{Any,1}:
#  CodeInfo(
# 53 1 ─ %1 = Base.add_int(%%x, %%y)::Int64                                   │
#    └──      return %1                                                       │
# ) => Int64


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mc01od(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    KMAX = 10
    REZ = Array{Float64,1}(undef, KMAX)
    IMZ = Array{Float64,1}(undef, KMAX)
    REP = Array{Float64,1}(undef, KMAX+1)
    IMP = Array{Float64,1}(undef, KMAX+1)
    DWORK = Array{Float64,1}(undef, 2*KMAX+2)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    K = parse(BlasInt, vs[1])
    if ( K<0 || K>KMAX )
        @error "Illegal K=$K"
    end
    vs = String[]
    _isz,_jsz = (K,2)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
        REZ[i],IMZ[i] = parsex.(Float64, vs[2i-1:2i])
    end
    close(f)

# interp call 1

    INFO = SLICOT.mc01od!(K, REZ, IMZ, REP, IMP)
    @test INFO == 0
    INFO == 0 || return

# unable to translate write loop:
#  write , REP(I+1), IMP(I+1), I = 0,K
# interp output 1
    println(io, "P:")
    show(io, "text/plain", REP[1:K+1]+im*IMP[1:K+1])
    println(io)

end # run_mc01od()


text:
# overwrite similar from indices, because similar on `Diagonal` tends to create
# `SparseArray` for some horrible reason

function similar_from_indices(T::Type, ind::IndexTuple, A::Diagonal, CA::Symbol)
    sz = similarstructure_from_indices(T, ind, A, CA)
    return similar(A.diag, T, sz)
end
function similar_from_indices(T::Type, poA::IndexTuple, poB::IndexTuple,
                                p1::IndexTuple, p2::IndexTuple,
                                A::Diagonal, B::AbstractArray, CA::Symbol, CB::Symbol)
    sz = similarstructure_from_indices(T, poA, poB, p1, p2, A, B, CA, CB)
    similar(A.diag, T, sz)
end
function similar_from_indices(T::Type, poA::IndexTuple, poB::IndexTuple,
                                p1::IndexTuple, p2::IndexTuple,
                                A::Diagonal, B::Diagonal, CA::Symbol, CB::Symbol)
    sz = similarstructure_from_indices(T, poA, poB, p1, p2, A, B, CA, CB)
    similar(A.diag, T, sz)
end

function contract!(α, A::AbstractArray, CA::Symbol, B::Diagonal, CB::Symbol,
        β, C::AbstractArray,
        oindA::IndexTuple, cindA::IndexTuple, oindB::IndexTuple, cindB::IndexTuple,
        indCinoAB::IndexTuple, syms::Union{Nothing, NTuple{3,Symbol}} = nothing)

    pA = (oindA...,cindA...)
    (length(pA) == ndims(A) && TupleTools.isperm(pA)) ||
        throw(IndexError("invalid permutation of length $(ndims(A)): $pA"))
    pB = (oindB...,cindB...)
    (length(pB) == ndims(B) && TupleTools.isperm(pB)) ||
        throw(IndexError("invalid permutation of length $(ndims(B)): $pB"))
    (length(oindA) + length(oindB) == ndims(C)) ||
        throw(IndexError("non-matching output indices in contraction"))
    (ndims(C) == length(indCinoAB) && isperm(indCinoAB)) ||
        throw(IndexError("invalid permutation of length $(ndims(C)): $indCinoAB"))

    sizeA = i->size(A, i)
    sizeB = i->size(B, i)
    sizeC = i->size(C, i)

    csizeA = sizeA.(cindA)
    csizeB = sizeB.(cindB)
    osizeA = sizeA.(oindA)
    osizeB = sizeB.(oindB)

    csizeA == csizeB ||
        throw(DimensionMismatch("non-matching sizes in contracted dimensions"))
    sizeAB = let osize = (osizeA..., osizeB...)
        i->osize[i]
    end
    sizeAB.(indCinoAB) == size(C) ||
        throw(DimensionMismatch("non-matching sizes in uncontracted dimensions"))

    Bd = B.diag
    if CA == :N && CB == :N
        @unsafe_strided A Bd C begin
            _contract!(α, A, Bd, β, C, oindA, cindA, oindB, cindB, indCinoAB)
        end
    elseif CA == :C && CB == :N
        @unsafe_strided A Bd C begin
            _contract!(α, conj(A), Bd, β, C, oindA, cindA, oindB, cindB, indCinoAB)
        end
    elseif CA == :N && CB == :C
        @unsafe_strided A Bd C begin
            _contract!(α, A, conj(Bd), β, C, oindA, cindA, oindB, cindB, indCinoAB)
        end
    elseif CA == :C && CB == :C
        @unsafe_strided A Bd C begin
            _contract!(α, conj(A), conj(Bd), β, C, oindA, cindA, oindB, cindB, indCinoAB)
        end
    else
        throw(ArgumentError("unknown conjugation flag $CA and $CB"))
    end
    return C
end
function _contract!(α, A::UnsafeStridedView, Bd::UnsafeStridedView,
        β, C::UnsafeStridedView, oindA, cindA, oindB, cindB, indCinoAB)

    sizeA = i->size(A, i)
    csizeA = sizeA.(cindA)
    osizeA = sizeA.(oindA)

    if length(oindB) == 1 # length(cindA) == length(cindB) == 1
        A2 = permutedims(A, (oindA..., cindA...))
        C2 = permutedims(C, TupleTools.invperm(indCinoAB))
        B2 = sreshape(Bd, ((one.(osizeA))..., csizeA...))
        totsize = (osizeA..., csizeA...)
        if α != 1
            if β == 0
                Strided._mapreducedim!((x,y)->α*x*y, +, zero, totsize, (C2, A2, B2))
            elseif β == 1
                Strided._mapreducedim!((x,y)->α*x*y, +, nothing, totsize, (C2, A2, B2))
            else
                Strided._mapreducedim!((x,y)->α*x*y, +, y->β*y, totsize, (C2, A2, B2))
            end
        else
            if β == 0
                return Strided._mapreducedim!(*, +, zero, totsize, (C2, A2, B2))
            elseif β == 1
                Strided._mapreducedim!(*, +, nothing, totsize, (C2, A2, B2))
            else
                Strided._mapreducedim!(*, +, y->β*y, totsize, (C2, A2, B2))
            end
        end
    elseif length(oindB) == 0
        strideA = i->stride(A, i)
        newstrides = (strideA.(oindA)..., strideA(cindA[1]) + strideA(cindA[2]))
        totsize = (osizeA..., csizeA[1])
        A2 = UnsafeStridedView(A.ptr, totsize, newstrides, A.offset, A.op)
        B2 = sreshape(Bd, ((one.(osizeA))..., csizeA[1]))
        C2 = permutedims(C, TupleTools.invperm(indCinoAB))

        if α != 1
            if β == 0
                Strided._mapreducedim!((x,y)->α*x*y, +, zero, totsize, (C2, A2, B2))
            elseif β == 1
                Strided._mapreducedim!((x,y)->α*x*y, +, nothing, totsize, (C2, A2, B2))
            else
                Strided._mapreducedim!((x,y)->α*x*y, +, y->β*y, totsize, (C2, A2, B2))
            end
        else
            if β == 0
                return Strided._mapreducedim!(*, +, zero, totsize, (C2, A2, B2))
            elseif β == 1
                Strided._mapreducedim!(*, +, nothing, totsize, (C2, A2, B2))
            else
                Strided._mapreducedim!(*, +, y->β*y, totsize, (C2, A2, B2))
            end
        end
    else # length(oindB) == 2
        if β != 1
            rmul!(C, β)
        end
        A2 = sreshape(permutedims(A, (oindA..., cindA...)), (osizeA..., 1))
        C2 = permutedims(C, TupleTools.invperm(indCinoAB))
        B2 = sreshape(Bd, ((one.(osizeA))..., length(Bd)))

        sC = strides(C2)
        newstrides = (Base.front(Base.front(sC))..., sC[end-1] + sC[end])
        totsize = (osizeA..., length(Bd))

        C3 = UnsafeStridedView(C2.ptr, totsize, newstrides, C2.offset, C2.op)
        if α != 1
            Strided._mapreducedim!((x,y)->α*x*y, +, nothing, totsize, (C3, A2, B2))
        else
            Strided._mapreducedim!(*, +, nothing, totsize, (C3, A2, B2))
        end
    end
    return C
end

function contract!(α, A::Diagonal, CA::Symbol, B::AbstractArray, CB::Symbol,
        β, C::AbstractArray,
        oindA::IndexTuple, cindA::IndexTuple, oindB::IndexTuple, cindB::IndexTuple,
        indCinoAB::IndexTuple, syms::Union{Nothing, NTuple{3,Symbol}} = nothing)

    indCinoAB′ = map(i->i <= length(oindA) ? length(oindB)+i : i-length(oindA), indCinoAB)
    contract!(α, B, CB, A, CA, β, C, oindB, cindB, oindA, cindA, indCinoAB′)
    return C
end


text:
import POMDPs: MDP, Solver

mutable struct DQNPolicySolver <: Solver
    epochs::Int64
    nnet_spec::Vector{Int64}
    epsilon::Float64
    epsilon_discount::Float64
    learning_rate::Float64
    discount::Float64
    max_frame_iterations::Int64
    print_every_n::Int64
end

DQNPolicySolver(epochs; nnet_spec::Vector{Int64} = [10, 5], epsilon = 0.5, epsilon_discount = 0.999, learning_rate = 0.005, discount = 0.99, max_frame_iterations = 100, print_every_n = 25) =
    DQNPolicySolver(epochs, nnet_spec, epsilon, epsilon_discount, learning_rate, discount, max_frame_iterations, print_every_n)

function solve(solver::DQNPolicySolver, pomdp::MDP; verbose = true)

    policy = DQNPolicy(pomdp; nnet_spec = solver.nnet_spec, nnet_optimizer = mx.ADAM(lr = solver.learning_rate), epsilon = solver.epsilon, epsilon_discount = solver.epsilon_discount)
    last_rewards = zeros(Float64, solver.epochs)

    for i = 1:solver.epochs

        if verbose && (i == 1 || i % solver.print_every_n == 0)
            println("Epoch: $i, Epsilon: $(round(policy.epsilon, 2))")
        end

        previous_state::Vector{Float64} = initial_state(pomdp)
        rewards = 0.0

        for j=1:solver.max_frame_iterations

            current_action_i, current_values = action(policy, previous_state)
            current_action = action_index(pomdp, current_action_i)
            current_state, current_reward = generate_sr(pomdp, previous_state, current_action)
            current_state_terminal = isterminal(pomdp, current_state)

            if j % 1 == 0
                println("Episode $i step $j reward $rewards, Current action:  $current_action, current reward: $current_reward")
            end

            next_value = if current_state_terminal 0 else maximum(values(policy, current_state)) end
            new_value = current_reward + solver.discount * next_value
            current_values[current_action_i + 1] = new_value

            update(policy, previous_state, current_values)

            previous_state = current_state
            rewards += current_reward

            if current_state_terminal
                println("Episode $i step $j reward $rewards, Current action:  $current_action, current reward: $current_reward")
                break
            end
        end

        if 0 < solver.epsilon_discount < 1
            policy.epsilon *= solver.epsilon_discount
        end

        last_rewards[i] = rewards

        if verbose && (i > 100 && i % solver.print_every_n == 0)
            cum_sum_reward = sum(last_rewards[i-101:i])
            cum_avg_reward = mean(last_rewards[i-101:i])
            println("Epoch: $i, Cummulative reward: $(round(cum_avg_reward, 2))/$(round(cum_sum_reward, 2)), Epsilon: $(policy.epsilon)")
        end
    end

    return policy

end

max_stars_repo_name: dmitrijsc/practical-rl

max_issues_repo_name: dmitrijsc/practical-rl

max_forks_repo_name: dmitrijsc/practical-rl


text:
module BioSimulatorPetriNets

using LightGraphs, TikzGraphs, BioSimulator

struct BioSimPetriNet
    g::SimpleDiGraph{Int}
    # node metadata
    species_nodes::Vector{Int}
    reaction_nodes::Vector{Int}
    node_labels::Vector{String}
    species_styles::Dict{Int,String}
    reaction_styles::Dict{Int,String}
    # edge metadata
    edge_labels::Dict{Tuple{Int,Int},String}
    edge_styles::Dict{Tuple{Int,Int},String}
end

# style for species nodes
const SPC_STYLE = "draw, rounded corners, fill=blue!10, font=\\large"

# style for reaction nodes
const RXN_STYLE = "draw, rounded corners, thick, fill=red!10, font=\\large"

# style for species ---> reaction edges
const SPCRXN_STYLE = "-stealth, draw, rounded corners=5pt, line width=0.5mm, solid, xshift=-2"

# style for reaction ---> species edges
const RXNSPC_STYLE = "-stealth, draw, rounded corners=5pt, line width=0.5mm, dashed, xshift=2, red, swap"

function construct_petri_net(model::Network)
    # initialize directed graph
    s = number_species(model)
    r = number_reactions(model)
    g = SimpleDiGraph(s + r)

    # define mapping between index and species
    species = species_list(model)
    idx2spc = Dict(id => i for (i, id) in enumerate(keys(species)))

    # initialize fields for storing edge metadata
    reactions = reaction_list(model)
    edge_set = Tuple{Int,Int}[]                 # pairs, i ---> j
    stoc_set = Dict{Tuple{Int,Int},Int}()       # stoichiometry for (i,j)
    edge_styles = Dict{Tuple{Int,Int},String}() # display style for (i,j)

    # construct edges; reactions numbered s+1 thru s+r
    for (k, reaction) in enumerate(values(reactions))
        j = k + s # node index in directed graph
        reactants = reaction.reactants
        products  = reaction.products

        # iterate over reactants to add species ---> reaction edges
        for (reactant, v) in reactants
            i = idx2spc[reactant] # retrieve node index for reactant
            e = (i, j)            # species i ---> reaction j

            edge_styles[e] = SPCRXN_STYLE # set edge style
            v > 1 && (stoc_set[e] = v)    # display stoichiometry > 1

            push!(edge_set, e)
        end

        # iterate over products to add reaction ---> species edges
        for (product, v) in products
            i = idx2spc[product] # retrieve node index for product
            e = (j, i)          # reaction j ---> species i

            edge_styles[e] = RXNSPC_STYLE # set edge style
            v > 1 && (stoc_set[e] = v)    # display stoichiometry > 1

            push!(edge_set, e)
        end
    end

    # add edges to directed graph
    for edge in edge_set
        add_edge!(g, edge[1], edge[2])
    end

    # set labels and styles for species nodes, 1:s
    species_nodes  = collect(1:s)
    species_labels = [string(v) for v in keys(species)]
    species_styles = Dict( i => SPC_STYLE for i in species_nodes )

    # set labels and styles for reaction nodes, s+1:s+r
    reaction_nodes = collect(s+1:s+r)
    reaction_labels = [string(v) for v in keys(reactions)]
    reaction_styles = Dict( i => RXN_STYLE for i in reaction_nodes)

    # put labels together
    node_labels = [species_labels; reaction_labels]
    edge_labels = Dict(e => string(v) for (e, v) in stoc_set)

    return BioSimPetriNet(
        g,
        species_nodes,
        reaction_nodes,
        node_labels,
        species_styles,
        reaction_styles,
        edge_labels,
        edge_styles
    )
end

function visualize_petri_net(x::BioSimPetriNet; layout = Layouts.Layered(), kwargs...)
    # unpack graph and metadata
    graph           = x.g
    labels          = x.node_labels
    species_styles  = x.species_styles
    reaction_styles = x.reaction_styles
    edge_labels     = x.edge_labels
    edge_styles     = x.edge_styles

    # replace underscores (_) with dots (⋅) to separate complexes
    labels = map(x -> replace(x, "_" => "\$\\cdot\$"), labels)

    # pass info to TikzGraphs
    return TikzGraphs.plot(graph, layout,
        labels,
        node_styles = merge(species_styles, reaction_styles),
        edge_labels = edge_labels,
        edge_styles = edge_styles;
        kwargs...
    )
end

function visualize_petri_net(x::Network; kwargs...)
    petri_net = construct_petri_net(x)

    return visualize_petri_net(x; kwargs...)
end

export BioSimPetriNet, construct_petri_net, visualize_petri_net

end # module


text:
function genw_sparse(N,m,sd,p)

N    = Int(N);
w    = zeros((N,N));

# recurrent connections
for i=1:N
  for j=1:N
    d = Normal(m,sd)
    b = Bernoulli(p)

    w[i,j] = rand(d)*rand(b);
  end
end

# No autapse
for i=1:N
    w[i,i] = 0.;
end

return w

end


text:
# load packages
using RPOMDPs, RPOMDPModels, RPOMDPToolbox, RobustValueIteration
using DataFrames, Query, CSV, ProgressMeter
const rpbvi = RobustValueIteration

# setup results structures
# solutions: holds policies and full sim results
# data: holds solution values and copmutation time
factors = DataFrame(ID = Int[], Problem = String[], Short_Name = String[],
    Solution = String[], Reward = String[], Information_Function = String[],
    Uncertainty_Size = Float64[], Dynamics = String[])
probnames = ["Crying Baby", "Tiger", "SimpleBaby", "SimpleTiger"]
shortnames = ["baby", "tiger", "simplebaby", "simpletiger"]
soltypes = ["Standard", "Robust"]
rewardtypes = ["Standard", "Information"]
infofuncs = ["Simple", "Approximate Entropy"]
uncsizes = [0.001, 0.01, 0.1, 0.2]
dyntypes = ["Nominal", "Ambiguous"]
respsols = ["Policy", "Simulation Values"]
headerfacts = ["ID","Problem", "Short Name", "Solution", "Reward",
    "Information Function", "Uncertainty Size", "Dynamics"]
headersols = vcat(factors, respsols)

ind = 0
for pname in probnames, sol in soltypes, r in rewardtypes, i in infofuncs, u in uncsizes, d in dyntypes
    ind += 1
    sname = shortnames[findin(probnames,[pname])[1]]
    push!(factors, [ind, pname, sname, sol, r, i, u, d])
end

for i in 1:size(factors,1)
    if factors[:Reward][i] == "Standard"
        factors[:Information_Function][i] = "NA"
    end
    if (factors[:Solution][i] == "Standard") && (factors[:Dynamics][i] == "Nominal")
        factors[:Uncertainty_Size][i] = 0.0
    end
end
new_factors = unique(factors[:, 2:8])

# simple info reward function
const ra11 = [1.0, -1.0] # belief-reward alpha vector
const ra12 = [-1.0, 1.0] # belief-reward alpha vector
const ralphas1 = [ra11, ra12]

# complicated info reward function
const ra21 = [1.0, 0.0]
const ra22 = [0.9, 0.1]
const ra23 = [0.8, 0.2]
const ra24 = [0.7, 0.3]
const ralphas2 = [ra21, ra22, ra23, ra24, 1-ra24, 1-ra23, 1-ra22, 1-ra21]

# set up problems - WIP
function build(sname, robust, info, err, rewardfunc)
    prob = nothing
    if rewardfunc == "Simple"
        r = ralphas1
    elseif rewardfunc == "Approximate Entropy"
        r = ralphas2
    else
        r = nothing
    end
    if sname == "baby"
        if robust == "Robust"
            prob = (info == "Standard") ? BabyRPOMDP(err) : BabyInfoRPOMDP(err, r)
        else
            prob = (info == "Standard") ? BabyPOMDP() : BabyInfoPOMDP(r)
        end
    elseif sname == "maze1x3"
        if robust == "Robust"
            prob = info == "Standard" ? Maze13RPOMDP(err) : Maze13RIPOMDP(err, r)
        else
            prob = info == "Standard" ? Maze13POMDP() : Maze13IPOMDP(r)
        end
    elseif sname == "maze4x3"
        if robust == "Robust"
            prob = info == "Standard" ? Maze43RPOMDP(err) : Maze43RIPOMDP(err, r)
        else
            prob = info == "Standard" ? Maze43POMDP() : Maze43IPOMDP(r)
        end
    elseif sname == "aircraftid"
        if robust == "Robust"
            prob = info == "Standard" ? AircraftRPOMDP(err) : AircraftRIPOMDP(err, r)
        else
            prob = info == "Standard" ? AircraftPOMDP() : AircraftIPOMDP(r)
        end
    elseif sname == "paint"
        if robust == "Robust"
            prob = info == "Standard" ? PaintRPOMDP(err) :PaintRIPOMDP(err, r)
        else
            prob = info == "Standard" ? PaintPOMDP() : PaintIPOMDP(r)
        end
    elseif sname == "shuttle"
        if robust == "Robust"
            prob = info == "Standard" ? ShuttleRPOMDP(err) : ShuttleRIPOMDP(err, r)
        else
            prob = info == "Standard" ? ShuttlePOMDP() : ShuttleIPOMDP(r)
        end
    elseif sname == "tiger"
        if robust == "Robust"
            prob = info == "Standard" ? TigerRPOMDP(0.95, err) : TigerInfoRPOMDP(0.95, err, r)
        else
            prob = info == "Standard" ? TigerPOMDP(0.95) : TigerInfoPOMDP(0.95, r)
        end
    elseif sname == "simplebaby"
        if robust == "Robust"
            prob = info == "Standard" ? SimpleBaby2RPOMDP(-15.0, -10.0, 0.9, err) : BabyInfoRPOMDP(err, r)
        else
            prob = info == "Standard" ? Baby2POMDP(-15.0, -10.0, 0.9) : BabyInfoPOMDP(r)
        end
    elseif sname == "simpletiger"
        if robust == "Robust"
            prob = info == "Standard" ? SimpleTigerRPOMDP(0.95, err) : TigerInfoRPOMDP(0.95, err, r)
        else
            prob = info == "Standard" ? TigerPOMDP(0.95) : TigerInfoPOMDP(0.95, r)
        end
    end
    prob
end

##################################################

sname = "simplebaby"
sversion = "5.1"
nreps = 500

# all sname runs
# dfexp = @from run in new_factors begin
#             @where run.Short_Name == sname
#             @select run
#             @collect DataFrame
# end

# non-info sname runs
dfexp = @from run in new_factors begin
            @where run.Short_Name == sname && run.Reward == "Standard"
            @select run
            @collect DataFrame
end

nrows = size(dfexp,1)
probs = Array{Union{POMDP,IPOMDP,RPOMDP,RIPOMDP}}(nrows)
simprobs = Array{Union{POMDP,IPOMDP,RPOMDP,RIPOMDP}}(nrows)
# problems and sim problems
for i = 1:nrows
    sname = dfexp[:Short_Name][i]
    robust = dfexp[:Solution][i]
    info = dfexp[:Reward][i]
    err = dfexp[:Uncertainty_Size][i]
    rewardfunc = dfexp[:Information_Function][i]
    probs[i] = build(sname, robust, info , err, rewardfunc)
    simrobust = (dfexp[:Dynamics][i] == "Nominal") ? "Standard" : "Robust"
    simprobs[i] = build(sname, simrobust, info, err, rewardfunc)
end

bs = [[b, 1-b] for b in 0.0:0.05:1.0]
nsteps = 100
solver = RPBVISolver(beliefpoints = bs, max_iterations = nsteps)
policies = Array{AlphaVectorPolicy}(nrows)
soldynamics = Array{AlphaVectorPolicy}(nrows) # worst-case dynamics
policies[1] = RPBVI.solve(solver, probs[1])
(dfexp[:Dynamics][1] == "Ambiguous") && (soldynamics[1] = RPBVI.solve(solver, simprobs[1]))

for i = 2:nrows
    println("\rPolicy $i")
    if probs[i] == probs[i-1]
        policies[i] = policies[i-1]
    else
        policies[i] = RPBVI.solve(solver, probs[i])
    end
end

function findpolicy(simprob::Union{POMDP,IPOMDP,RPOMDP,RIPOMDP}, probs, policies)
    ind = 0
    counter = 1
    while ind < 1
        (simprob == probs[counter]) && (ind = counter)
        counter += 1
    end
    policies[ind]
end

for i = 1:nrows
    soldynamics[i] = findpolicy(simprobs[i], probs, policies)
end

function meanci(data::Vector{Float64})
    n = length(data)
    m = mean(data)
    s = std(data)
    tstar = 1.962
    hw = tstar * s / sqrt(n)
    (m - hw, m + hw)
end

psim = RolloutSimulator(max_steps = nsteps)
simvals = [Vector{Float64}(nreps) for _ in 1:nrows] # simulated values
simps = [Vector{Float64}(nreps) for _ in 1:nrows] # simulated percent correct
ves = Array{Float64}(nrows) # expected values
vms = Array{Float64}(nrows) # mean of sim values
vss = Array{Float64}(nrows) # std dev of sim values
vci = Array{Tuple{Float64,Float64}}(nrows) # 95% ci of mean of sim values
pms = Array{Float64}(nrows) # mean percent correct
pss = Array{Float64}(nrows) # std percent correct
pci = Array{Tuple{Float64,Float64}}(nrows) # 95% ci of mean of sim percent correct

# rollout generates state/obs/rewards based on the sim problem, and beliefs based on the solution policy
srand(923475)
for i in 1:nrows
    println("Run $i")
    println("Solution Type: ", new_factors[:Solution][i], " Dynamics, ", new_factors[:Reward][i], " Reward")
    prob = probs[i]
    simprob = simprobs[i]
    sol = policies[i]
    ves[i] = value(sol, initial_belief(prob))
    bu = updater(sol)
    println("Simulation Type: ", new_factors[:Dynamics][i], " ", typeof(simprob))
    @showprogress 1 "Simulating value..." for j = 1:nreps
        if typeof(simprob) <: Union{RPOMDP,RIPOMDP}
            sv, sp = simulate_worst(psim, simprob, sol, bu, soldynamics[i].alphas)
        else
            sv, sp = simulate(psim, simprob, sol, bu)
        end
        simvals[i][j] = sv
        simps[i][j] = sp
    end
    vms[i] = mean(simvals[i])
    vss[i] = std(simvals[i])
    vci[i] = meanci(simvals[i])
    pms[i] = mean(simps[i])
    pss[i] = std(simps[i])
    pci[i] = meanci(simps[i])
end

ids = collect(1:size(dfexp,1))
dfexp[:ID] = ids
rdata = DataFrame(ID = ids, ExpectedValue = ves, SimMean = vms, SimStd = vss,
    SimCI = vci, CorrectMean = pms, CorrectStd = pss, CorrectCI = pci)
ndec = 3
rdata_round = DataFrame(ID = ids, ExpectedValue = round.(ves, ndec),
    SimMean = round.(vms, ndec),
    SimStd = round.(vss, ndec),
    SimCI = [round.(vci[i], ndec) for i in 1:length(vci)],
    CorrectMean = round.(pms, ndec),
    CorrectStd = round.(pss, ndec),
    CorrectCI = [round.(pci[i], ndec) for i in 1:length(pci)])
df = join(dfexp, rdata_round, on = :ID)
simdata = hcat(ids, hcat(simvals...)') |> DataFrame
@show rdata

path = joinpath(homedir(),".julia\\v0.6\\RobustInfoPOMDP\\data")
fnresults = string("exp_results_", sname, "_", sversion, ".csv")
fnsim = string("exp_sim_values_", sname, "_", sversion, ".csv")
CSV.write(joinpath(path, fnresults), df)
CSV.write(joinpath(path, fnsim), simdata)

using Plots; gr()
rnum = 1
valfunc1 = [maximum(dot(policies[rnum].alphas[i], [b, 1-b]) for i = 1:length(policies[rnum].alphas))
        for b = 0.0:0.01:1.0]
valfunc13 = [maximum(dot(policies[13].alphas[i], [b, 1-b]) for i = 1:length(policies[13].alphas))
        for b = 0.0:0.01:1.0]
x = collect(0:0.01:1.0)
vplot = plot(x, valfunc1, xticks = 0:0.1:1, label = "Nominal",
        # title = "Tiger POMDP Value Function",
        xlab = "Belief, P(State = Tiger Left)",
        # xlab = "Belief, P(State = Hungry)",
        ylab = "Expected Total Discounted Reward",
        legend = :bottomleft,
        # legend = :topright,
        line = :dash,
        linealpha = 0.9)
plot!(x, valfunc13, color = :red, linealpha = 0.8, label = "Robust: 0.2")
fn = string("value_function_", sname, "_", sversion, ".pdf")
path = joinpath(homedir(),".julia\\v0.6\\RobustInfoPOMDP\\data\\figures\\",fn)
savefig(vplot, path)

max_stars_repo_name: ajkeith/RobustInfoPOMDP

max_issues_repo_name: ajkeith/RobustInfoPOMDP

max_forks_repo_name: ajkeith/RobustInfoPOMDP


text:
##############
# RTN | LVLH #
##############

export rRTNtoECI
export rECItoRTN

"""
Compute the radial, along-track, cross-track (RTN) to Earth-centered inertial
rotation matrix. If applied to a position vector in the RTN frame, it will
transform that vector to into the equivalent position vector in the ECI frame.

The RTN frame is also commonly refered to as the local-vertical, local-horizontal (LVLH) frame.

Arguments:
- `x::AbstractVector{<:Real}`: Inertial state (position and velocity) of primary (observing) satellite

Returns:
- `R_rtn2eci::SMatrix{3,3}`: Rotation matrix transforming _from_ the RTN frame _to_ the ECI frame
"""
function rRTNtoECI(x::AbstractVector{<:Real})
    r = x[idx1t3]
    v = x[idx4t6]
    n = cross(r, v)

    R = normalize(r)
    N = normalize(n)
    T = cross(N, R)

    R_rtn2eci = hcat(R, T, N)

    return R_rtn2eci
end

"""
Compute the Earth-centered inertial to radial, along-track, cross-track (RTN)
rotation matrix. If applied to a position vector in the ECI frame, it will
transform that vector into the equivalent position vector in the RTN frame.

The RTN frame is also commonly refered to as the local-vertical, local-horizontal (LVLH) frame.

Arguments:
- `x::AbstractVector{<:Real}`: Inertial state (position and velocity) of primary (observing) satellite

Returns:
- `R_eci2rtn::SMatrix{3,3}`: Rotation matrix transforming _from_ the ECI frame _to_ the RTN frame
"""
function rECItoRTN(x::AbstractVector{<:Real})
    r = x[idx1t3]
    v = x[idx4t6]
    n = cross(r, v)

    R = normalize(r)
    N = normalize(n)
    T = cross(N, R)

    R_eci2rtn = hcat(R, T, N)'

    return R_eci2rtn
end

export sECItoRTN
"""
Compute the radial, along-track, cross-track (RTN) coordinates of a target satellite in the primary satellites RTN frame.

The RTN frame is also commonly refered to as the local-vertical, local-horizontal (LVLH) frame.

Arguments:
- `x::AbstractVector{<:Real}`: Inertial state (position and velocity) of primary (observing) satellite
- `xt::AbstractVector{<:Real}`: Inertial state (position and velocity) of the target satellite

Returns:
- `xrtn::AbstractVector{<:Real}`: Position and velocity of the target relative of the observing satellite in the RTN.
"""
function sECItoRTN(x::AbstractVector{<:Real}, xt::AbstractVector{<:Real})
    if length(xt) >= 6
        return sECItoRTN(x, xt[idx1t6])
    else
        return sECItoRTN(x, xt[idx1t3])
    end
end

"""
Arguments:
- `x::AbstractVector{<:Real}`: Inertial state (position and velocity) of primary (observing) satellite
- `xt::SVector{3,<:Real}`: Inertial state (position) of the target satellite

Returns:
- `xrtn::SVector{3,<:Real}`: Position of the target relative of the observing satellite in the RTN.
"""
function sECItoRTN(x::AbstractVector{<:Real}, xt::SVector{3,<:Real})
    # Create RTN rotation matrix
    R_eci2rtn = rECItoRTN(x)

    # Transform Position
    r     = x[idx1t3]
    rho   = xt[idx1t3] - r
    r_rtn = R_eci2rtn*rho

    return r_rtn
end

"""
Arguments:
- `x::AbstractVector{<:Real}`: Inertial state (position and velocity) of primary (observing) satellite
- `xt::SVector{6,<:Real}`: Inertial state (position and velocity) of the target satellite

Returns:
- `xrtn::SVector{6,<:Real}`: Position and velocity of the target relative of the observing satellite in the RTN.
"""
function sECItoRTN(x::AbstractVector{<:Real}, xt::SVector{6,<:Real})
    # Create RTN rotation matrix
    R_eci2rtn = rECItoRTN(x)

    # Transform Position
    r     = x[idx1t3]
    rho   = xt[idx1t3] - r
    r_rtn = R_eci2rtn*rho

    v       = x[idx4t6]
    f_dot   = norm(cross(r, v)) / norm(r)^2
    omega   = SVector(0, 0, f_dot)
    rho_dot = xt[idx4t6] - v
    v_rtn   = R_eci2rtn*rho_dot - cross(omega, r_rtn)

    return vcat(r_rtn, v_rtn)
end

export sRTNtoECI
"""
Compute the Earth-center

The RTN frame is also commonly refered to as the local-vertical, local-horizontal (LVLH) frame.

Arguments:
- `x::AbstractVector{<:Real}`: Inertial state (position and velocity) of primary (observing) satellite
- `xrtn::AbstractVector{<:Real}`: Inertial state (position and velocity) of the target satellite

Returns:
- `xt::AbstractVector{<:Real}`: Position and velocity of the target relative of the observing satellite in the RTN.
"""
function sRTNtoECI(x::AbstractVector{<:Real}, xrtn::AbstractVector{<:Real})
    if length(xrtn) >= 6
        return sRTNtoECI(x, xrtn[idx1t6])
    else
        return sRTNtoECI(x, xrtn[idx1t3])
    end
end

"""
Arguments:
- `x::AbstractVector{<:Real}`: Inertial state (position and velocity) of primary (observing) satellite
- `xrtn::SVector{3,<:Real}`: Inertial state (position) of the target satellite

Returns:
- `xt::SVector{3,<:Real}`: Position of the target relative of the observing satellite in the RTN.
"""
function sRTNtoECI(x::AbstractVector{<:Real}, xrtn::SVector{3,<:Real})
    # Create RTN rotation matrix
    R_rtn2eci = rRTNtoECI(x)

    # Transform position
    r   = x[idx1t3]
    r_t = R_rtn2eci*xrtn + r

    return r_t
end

"""
Arguments:
- `x::AbstractVector{<:Real}`: Inertial state (position and velocity) of primary (observing) satellite
- `xrtn::SVector{6,<:Real}`: Inertial state (position and velocity) of the target satellite

Returns:
- `xt::SVector{6,<:Real}`: Position and velocity of the target relative of the observing satellite in the RTN.
"""
function sRTNtoECI(x::AbstractVector{<:Real}, xrtn::SVector{6,<:Real})
    # Create RTN rotation matrix
    R_rtn2eci = rRTNtoECI(x)

    # Transform position
    r     = x[idx1t3]
    r_rtn = xrtn[idx1t3]
    r_t   = R_rtn2eci*r_rtn + r

    # Transform velocity
    v     = x[idx4t6]
    v_rtn = xrtn[idx4t6]
    f_dot = norm(cross(r, v)) / norm(r)^2
    omega = SVector(0, 0, f_dot)
    v_t   = R_rtn2eci*(v_rtn + cross(omega, r_rtn)) + v

    return vcat(r_t, v_t)
end

#######
# ROE #
#######

"""
Convert Keplerian relative orbital element state to ROE relative state
"""
function sOSCtoROE(oe_c::SVector{6,<:Real}, oe_d::SVector{6,<:Real})
end

"""
Convert ROE relative state and chief relative orbital elements and return
osculating orbital elements
"""
function sROEtoOSC(oe_c::SVector{6,<:Real}, roe_d::SVector{6,<:Real})
end

"""
Convert inertial states to relative orbital element
"""
function sECItoROE(x_c::SVector{6,<:Real}, x_d::SVector{6,<:Real})
end

"""
Convert ROE relative state and chief absolute inertial state and return
osculating orbital elements
"""
function sROEtoECI(x_c::SVector{6,<:Real}, roe_d::SVector{6,<:Real})
end

"""
Convert RTN state and chief absolute inertial state and return relative
orbital element state
"""
function sRTNtoROE(x_c::SVector{6,<:Real}, rtn_d::SVector{6,<:Real})
end

"""
Covert relative orbital element state and chief inertial state and return
RTN state.
"""
function sROEtoRTN(x_c::SVector{6,<:Real}, roe_d::SVector{6,<:Real})
end


text:
using StoppingTutorials
using Test

file_path(file_name) = joinpath(@__DIR__,"..", "src", file_name)

@testset "StoppingTutorials.jl" begin
  # Should we test something here?
end

for (i,(title,filename)) in enumerate(StoppingTutorials.files)
  @testset "StoppingTutorials: test $(title)" begin
    include(file_path(filename))
  end
end


text:
T_reflecting = [-1.1 1.1; 1.0 -1.0]
c_reflecting = [1.0; -1.0]
b_reflecting = 10.0

P_upr_reflecting = [0.0 1.0]
P_lwr_reflecting = [1.0 0.0] # both reflecting boundaries

model = BoundedFluidQueue(T_reflecting,c_reflecting,P_lwr_reflecting,P_upr_reflecting,b_reflecting)


text:
using BinaryBuilder

compiler_target = triplet(platform_key(ARGS[end]))
if compiler_target == "unknown-unknown-unknown"
    error("This is not a typical build_tarballs.jl!  Must provide exactly one platform as the last argument!")
end
deleteat!(ARGS, length(ARGS))

# Encode the compiler target into the name
name = "Binutils-$(compiler_target)"
version = v"2.24"

# Collection of sources required to build Binutils
sources = [
    "https://ftp.gnu.org/gnu/binutils/binutils-2.24.tar.bz2" =>
	"e5e8c5be9664e7f7f96e0d09919110ab5ad597794f5b1809871177a0f0f14137",
    "https://github.com/tpoechtrager/apple-libtapi.git" =>
    "e56673694db395e25b31808b4fbb9a7005e6875f",
    "https://github.com/tpoechtrager/cctools-port.git" =>
    "ecb84d757b6f011543504967193375305ffa3b2f",
    "./bundled",
]

# Bash recipe for building across all platforms
script = """
# FreeBSD build system for binutils apparently requires that uname sit in /usr/bin/
ln -sf \$(which uname) /usr/bin/uname

# On MacOS, we don't actually install Binutils.  We install libtapi and cctools.  Le sigh.
if [[ $(compiler_target) == *apple* ]]; then
    mkdir -p \${WORKSPACE}/srcdir/apple-libtapi/build
    cd \${WORKSPACE}/srcdir/apple-libtapi/build

    # Install libtapi
    cmake ../src/apple-llvm/src \\
        -DLLVM_INCLUDE_TESTS=OFF \\
        -DCMAKE_BUILD_TYPE=RELEASE \\
        -DCMAKE_INSTALL_PREFIX=\${prefix}
    make -j\${nproc} VERBOSE=1
    make install

    # Install cctools.  Someday, try this without disabling the clang assembler!
    cd \${WORKSPACE}/srcdir/cctools-port/cctools
    ./configure --prefix=\${prefix} \\
        --target=$(compiler_target) \\
        --host=\${MACHTYPE} \\
        --with-libtapi=\${prefix}
    make -j\${nproc}
    make install
else
    cd \${WORKSPACE}/srcdir/binutils-*/

    #atomic_patch -p1 "\${WORKSPACE}/srcdir/patches/binutils_COFF_bfd.patch"
    update_configure_scripts

    ./configure --prefix=\${prefix} \\
        --target=$(compiler_target) \\
        --host=\${MACHTYPE} \\
        --with-sysroot="\${prefix}/$(compiler_target)/sys-root" \\
        --enable-multilib \\
        --disable-werror

    make -j\${nproc}
    make install
fi

# Finally, create a bunch of symlinks stripping out the target so that
# things like `nm` "just work", as long as we've got our path set properly.
for f in \${prefix}/bin/$(compiler_target)-*; do
    fbase=\$(basename \$f)
    ln -s \$fbase \${prefix}/bin/\${fbase#$(compiler_target)-}
done
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [Platform("x86_64", "linux")]

# The products that we will ensure are always built
products = prefix -> [
    ExecutableProduct(prefix, "$(compiler_target)-as", :as),
    ExecutableProduct(prefix, "$(compiler_target)-nm", :nm),
    ExecutableProduct(prefix, "$(compiler_target)-ld", :ld),
    ExecutableProduct(prefix, "$(compiler_target)-strip", :strip),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/staticfloat/KernelHeadersBuilder/releases/download/v4.12.0-0/build_KernelHeaders.v4.12.0.jl",
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; skip_audit=true)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using Statistics,Serialization,StatsPlots

path_to_data = "../data/mutation/beta_l50/"
filename_prefix = "../figures/"

mkpath(filename_prefix)

include("$(path_to_data)covid_par_ini.jl")
worker_results = deserialize(open("$(path_to_data)batchdata.dat"))

datapoint = fld(T,datat)

include("../emp_traj.jl")

empn = 200 # number of periods for empirical plots

# Total number of infected: model vs. empirical data for Germany
totinfmean = mean(hcat((results -> results[:totinftraj]).(worker_results)...), dims=2)
totinfstd = std(hcat((results -> results[:totinftraj]).(worker_results)...), dims=2)
fac = nhh/50000
nzer= (Int(ceil(virustime/datat))-1)*datat+2*corlatent
emptotinftraj = vcat(zeros(nzer),fac*emptotinf)
nne = size(emptotinf)[1]

emptotinftraj = vcat(zeros(nzer),fac*emptotinfnew1)
nnen = size(emptotinfnew1)[1]
nnen = empn
totinfmeansc = mean(hcat(detfrac*(results -> results[:totinftraj]).(worker_results)...), dims=2)
totinfstdsc = std(hcat(detfrac*(results -> results[:totinftraj]).(worker_results)...), dims=2)
pl5b = plot(collect(1:nnen),[totinfmeansc[nzer+1:nzer+nnen], totinfmeansc[nzer+1:nzer+nnen].-totinfstdsc[nzer+1:nzer+nnen],totinfmeansc[nzer+1:nzer+nnen].+totinfstdsc[nzer+1:nzer+nnen]], linestyle= [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:blue :black :black],label = false, xlabel="days", ylabel = "reported tot. inf.")
plot!(pl5b, detfrac*emptotinftraj[nzer+1:nzer+nnen], linecolor = [:green], label=false, linewidth = [2])
savefig(pl5b,"$(filename_prefix)baseline_empirical_total_infected.pdf")

# R-value: model vs. empirical data for Germany
RKIR0mean = mean(hcat((results -> results[:RKIR0traj]).(worker_results)...), dims=2)
RKIR0std = std(hcat((results -> results[:RKIR0traj]).(worker_results)...), dims=2)
nn1 = size(RKIR0mean)[1]

nzer=(Int(ceil(virustime/datat))-1)*datat+2*corlatent
empR0smtraj = vcat(zeros(nzer),empR0sm)
RKIR0smmean = mean(hcat((results -> results[:RKIR0smtraj]).(worker_results)...), dims=2)
RKIR0smstd = std(hcat((results -> results[:RKIR0smtraj]).(worker_results)...), dims=2)
nnhc = nn1 - size(RKIR0smmean)[1]
nne = size(empR0sm)[1]
RKIR0smmean = vcat(zeros(nnhc),RKIR0smmean)
RKIR0smstd = vcat(zeros(nnhc),RKIR0smstd)

nzer=(Int(ceil(virustime/datat))-1)*datat+2*corlatent
empR0smtraj = vcat(zeros(nzer),empR0new1)
RKIR0smmean = mean(hcat((results -> results[:RKIR0smtraj]).(worker_results)...), dims=2)
RKIR0smstd = std(hcat((results -> results[:RKIR0smtraj]).(worker_results)...), dims=2)
nnhc = nn1 - size(RKIR0smmean)[1]
nne = size(empR0new1)[1]
nne = empn
RKIR0smmean = vcat(zeros(nnhc),RKIR0smmean)
RKIR0smstd = vcat(zeros(nnhc),RKIR0smstd)
for i=1:7
    RKIR0smmean[nzer+i] = 1/0
    empR0smtraj[nzer+i] = 1/0
end

for i in 1:length(RKIR0smmean)
    if RKIR0smmean[i] == Inf
        RKIR0smmean[i] = RKIR0smmean[i-1]
        RKIR0smstd[i] = 0.0
    end
end
for i=1:7
    RKIR0smmean[nzer+i] = 1/0
    empR0smtraj[nzer+i] = 1/0
end

# exclude Inf values
RKIR0s = hcat((results -> results[:RKIR0smtraj]).(worker_results)...)
x = 0
for i in 1:length(RKIR0s)
    if RKIR0s[i] == Inf
        RKIR0s[i] = NaN
        global x = i
    end
end

nanmean1(x) = mean(filter(!isnan,x))
nanstd1(x) = std(filter(!isnan,x))

RKIR0smmean= []
RKIR0smstd= []
for i in 1:length(RKIR0s[:,1])
    RKIR0s_row = RKIR0s[i,:]
    push!(RKIR0smmean,nanmean1(RKIR0s_row))
    push!(RKIR0smstd,nanstd1(RKIR0s_row))
end
nnhc = nn1 - size(RKIR0smmean)[1]
RKIR0smmean = vcat(zeros(nnhc),RKIR0smmean)
RKIR0smstd = vcat(zeros(nnhc),RKIR0smstd)
pl16b = plot(collect(1:nne),[RKIR0smmean[nzer+1:nzer+nne], RKIR0smmean[nzer+1:nzer+nne].-RKIR0smstd[nzer+1:nzer+nne],RKIR0smmean[nzer+1:nzer+nne].+RKIR0smstd[nzer+1:nzer+nne]], linestyle= [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:blue :black :black],label = false, xlabel = "days", ylabel = "R0")
plot!(pl16b, empR0smtraj[nzer+1:nzer+nne], linecolor = [:green], label=false, linewidth = [2])
savefig(pl16b,"$(filename_prefix)baseline_empirical_R-value.pdf")

# Casualites: model vs. empirical data for Germany
totcasmean = mean(cat(100*(results -> results[:totcastraj]).(worker_results)...,dims=3), dims=3)
totcasstd = std(cat(100*(results -> results[:totcastraj]).(worker_results)...,dims=3), dims=3)
nzer= (Int(ceil(virustime/datat+2*corlatent/datat))-1)
nzerw= (Int(ceil(virustime/datat+2*corlatent/datat))-1)
empcastraj = vcat(zeros(datat*nzer),fac*empcasnew1)
nne = Int(ceil(size(empcasnew1)[1]/datat))
nne = Int(ceil(empn / datat))
pl7b = plot(datat*collect(1:nne),[totcasmean[nzer+1:nzer+nne,3], totcasmean[nzer+1:nzer+nne,3].-totcasstd[nzer+1:nzer+nne,3],totcasmean[nzer+1:nzer+nne,3].+totcasstd[nzer+1:nzer+nne,3]], linestyle = [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:blue :black :black] )
plot!(pl7b, empcastraj[datat*nzer+1:empn+datat*nzer], linecolor = [:green], linewidth = [2], legend = false, xlabel = "days", ylabel = "casualties [%]")
savefig(pl7b,"$(filename_prefix)baseline_empirical_casualties.pdf")

# GDP loss: model vs. emprical data for Germany
gdplossmean = mean(hcat((results -> results[:gdplosstraj]).(worker_results)...), dims=2)
gdplossstd = std(hcat((results -> results[:gdplosstraj]).(worker_results)...), dims=2)
nnew = Int(ceil(empn/datat))
pl10a = plot(datat*collect(1:nnew),[gdplossmean[nzerw+1:nzerw+nnew], gdplossmean[nzerw+1:nzerw+nnew].-gdplossstd[nzerw+1:nzerw+nnew],gdplossmean[nzerw+1:nzerw+nnew].+gdplossstd[nzerw+1:nzerw+nnew]], linestyle= [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:blue :black :black], label=false, xlabel = "days", ylabel = "GDP loss[%]")
# calculate average GDP loss in quarters
q1loss = (10*0+sum(gdplossmean[4:6]))/13
q2loss = sum(gdplossmean[7:19])/13
q3loss = sum(gdplossmean[20:32])/13
scatter!(pl10a, [21], [0.018], markersize = 5, markercolor = [:green], label = false)
scatter!(pl10a, [112], [0.113], markersize = 5, markercolor = [:green], label = false)
scatter!(pl10a, [201], [0.037], markersize = 5, markercolor = [:green], label = false)
savefig(pl10a,"$(filename_prefix)baseline_empirical_gdploss.pdf")

# Workers on short time: model vs. empirical data for Germany
shorttimemean = mean(hcat((results -> results[:shorttimetraj]).(worker_results)...), dims=2)
shorttimestd = std(hcat((results -> results[:shorttimetraj]).(worker_results)...), dims=2)

pl26b1 = plot(datat*collect(1:nnew),[shorttimemean[nzerw+1:nzerw+nnew], shorttimemean[nzerw+1:nzerw+nnew].-shorttimestd[nzerw+1:nzerw+nnew],shorttimemean[nzerw+1:nzerw+nnew].+shorttimestd[nzerw+1:nzerw+nnew]], linestyle= [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:blue :black :black], label=false, xlabel = "days", ylabel = "short time [%]")
# tatsächlich ausgezahlte kurzarbeit from april to july with workforce 33 mio
scatter!(pl26b1, [23], [0.078], markersize = 5, markercolor = [:green], label = false)
scatter!(pl26b1, [53], [0.180], markersize = 5, markercolor = [:green], label = false)
scatter!(pl26b1, [83], [0.179], markersize = 5, markercolor = [:green], label = false)
scatter!(pl26b1, [114], [0.140], markersize = 5, markercolor = [:green], label = false)
scatter!(pl26b1, [145], [0.128], markersize = 5, markercolor = [:green], label = false)
scatter!(pl26b1, [176], [0.078], markersize = 5, markercolor = [:green], label = false)
scatter!(pl26b1, [201], [0.067], markersize = 5, markercolor = [:green], label = false)

savefig(pl26b1,"$(filename_prefix)baseline_empirical_shorttime.pdf")

# Excess household saving per GDP: model vs. empirical data for Germany
gdpmean = mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
gdpstd = std(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
totalsavmean = mean(hcat((results -> results[:totalsavtraj]).(worker_results)...), dims=2)
totalsavstd = std(hcat((results -> results[:totalsavtraj]).(worker_results)...), dims=2)
totalsavmean_0 = (totalsavmean/nhh .- totalsavmean[2]/nhh)./(gdpmean[2]*52)
totalsavstd_0 =totalsavstd/nhh - totalsavstd/nhh
pl28a = plot(datat*collect(1:nnew),[totalsavmean_0[nzerw+1:nzerw+nnew], (totalsavmean_0[nzerw+1:nzerw+nnew].-totalsavstd_0[nzerw+1:nzerw+nnew]/(gdpmean[2]*52)),(totalsavmean_0[nzerw+1:nzerw+nnew].+totalsavstd_0[nzerw+1:nzerw+nnew]/(gdpmean[2]*52))], linestyle= [:solid :dot :dot], linewidth = [2 1 1], linecolor = [:blue :black :black], label=false, ylabel = "Excess HH Savings per GDP per Capita", xlabel = "days")
scatter!(pl28a, [21], [0.004], markersize = 5, markercolor = [:green], label = false)
scatter!(pl28a, [112], [0.019], markersize = 5, markercolor = [:green], label = false)
scatter!(pl28a, [201], [0.026], markersize = 5, markercolor = [:green], label = false)
savefig(pl28a,"$(filename_prefix)baseline_emprical_exccess_savings.pdf")


text:
"""
    TimeMachine.Callbacks

A suite of callback functions to be used with the TimeMachine.jl ODE solvers.
"""
module Callbacks

import MPI, DiffEqBase

"""
    TimeMachine.Callbacks.initialize!(f!::F, integrator)

Initialize a callback event for callbacks of type `F`. By default this does nothing, but
can be extended for new callback events.
"""
function initialize!(f!, integrator)
end

"""
    TimeMachine.Callbacks.finalize!(f!::F, integrator)

Finalize a callback event for callbacks of type `F`. By default this does nothing, but
can be extended for new callback events.
"""
function finalize!(f!, integrator)
end

export EveryXWallTimeSeconds, EveryXSimulationTime, EveryXSimulationSteps

"""
    EveryXWallTimeSeconds(f!, Δwt, comm::MPI.Comm;
                          atinit=false)

Trigger `f!(integrator)` every `Δwt` wallclock seconds.

An MPI communicator `comm` must be provided to synchronize timing across all ranks.

[`Callbacks.initialize!`](@ref) and [`Callbacks.finalize!`](@ref) can be defined for `f!`.

If `atinit=true`, then `f!(integrator)` will additionally be triggered at initialization,
otherwise the first trigger will be after `Δwt` seconds.
"""
function EveryXWallTimeSeconds(f!, Δwt, comm::MPI.Comm;
                               atinit=false)
    wt_next = 0.0

    function _initialize(c, u, t, integrator)
        wt = MPI.Allreduce(time(), max, comm)
        wt_next = wt + Δwt
        initialize!(c.affect!, integrator)
        if atinit
            c.affect!(integrator)
        end
    end

    function _finalize(c, u, t, integrator)
        finalize!(c.affect!, integrator)
    end

    function condition(u, t, integrator)
        wt = MPI.Allreduce(time(), max, comm)
        if wt >= wt_next
            while wt >= wt_next
                wt_next += Δwt
            end
            return true
        else
            return false
        end
    end

    if isdefined(DiffEqBase, :finalize!)
        DiffEqBase.DiscreteCallback(condition, f!; initialize=_initialize, finalize=_finalize)
    else
        DiffEqBase.DiscreteCallback(condition, f!; initialize=_initialize)
    end
end

"""
    EveryXSimulationTime(f!, Δt;
                         atinit=false)

Trigger `f!(integrator)` every `Δt` simulation time.

[`Callbacks.initialize!`](@ref) and [`Callbacks.finalize!`](@ref) can be defined for `f!`.

If `atinit=true`, then `f!` will additionally be triggered at initialization. Otherwise
the first trigger will be after `Δt` simulation time.
"""
function EveryXSimulationTime(f!, Δt;
                              atinit=false)
    t_next = zero(Δt)

    function _initialize(c, u, t, integrator)
        t_next = Δt
        initialize!(c.affect!, integrator)
        if atinit
            c.affect!(integrator)
        end
    end

    function _finalize(c, u, t, integrator)
        finalize!(c.affect!, integrator)
    end

    function condition(u, t, integrator)
        if t >= t_next
            while t >= t_next
                t_next += Δt
            end
            return true
        else
            return false
        end
    end
    if isdefined(DiffEqBase, :finalize!)
        DiffEqBase.DiscreteCallback(condition, f!; initialize=_initialize, finalize=_finalize)
    else
        DiffEqBase.DiscreteCallback(condition, f!; initialize=_initialize)
    end
end

"""
    EveryXSimulationSteps(f!, Δsteps;
                         atinit=false)

Trigger `f!(integrator)` every `Δsteps` simulation steps.

[`Callbacks.initialize!`](@ref) and [`Callbacks.finalize!`](@ref) can be defined for `f!`.

If `atinit==true`, then `f!` will additionally be triggered at initialization. Otherwise
the first trigger will be after `Δsteps`.
"""
function EveryXSimulationSteps(f!, Δsteps;
                               atinit=false)
    steps = 0
    steps_next = 0

    function _initialize(c, u, t, integrator)
        steps = 0
        steps_next = Δsteps
        initialize!(c.affect!, integrator)
        if atinit
            c.affect!(integrator)
        end
    end

    function _finalize(c, u, t, integrator)
        finalize!(c.affect!, integrator)
    end

    function condition(u, t, integrator)
        steps += 1
        if steps >= steps_next
            steps_next += Δsteps
            return true
        else
            return false
        end
    end

    if isdefined(DiffEqBase, :finalize!)
        DiffEqBase.DiscreteCallback(condition, f!; initialize=_initialize, finalize=_finalize)
    else
        DiffEqBase.DiscreteCallback(condition, f!; initialize=_initialize)
    end
end

end # module


text:
# Apply a log-transform on a Gamma-distributed variable to sample from it using a normal proposal
# Use forward-mode automatic differentiation to compute the gradient of the log-target

using Klara

plogtarget(p::Real, v::Vector) = (v[1]-1)*p-exp(p)/v[2]+p

p = BasicContUnvParameter(:p, logtarget=plogtarget, nkeys=3, diffopts=DiffOptions(mode=:forward))

model = likelihood_model([Constant(:k), Constant(:θ), p], isindexed=false)

mcsampler = MALA()

mcrange = BasicMCRange(nsteps=100000, burnin=10000)

v0 = Dict(:k=>2., :θ=>1., :p=>log(10.))

job = BasicMCJob(model, mcsampler, mcrange, v0)

run(job)

chain = output(job)

mean(exp.(chain.value))

acceptance(chain, diagnostics=false)