text:
module Lights
"""
Module for the modelling of light pollution.
"""

using ArchGDAL

include(".\\Utils\\Functions.jl")

export run_light

const agd = ArchGDAL

#= Based on
    https://github.com/thingless/viewshed/blob/master/server/algo.py

function iter_to_runs( visibles, pixels )
    cur_val = Inf
    start_idx = nothing
    out = []
    for (i, val) in enumerate(vcat(visibles, [nothing]))
        if cur_val != val
            if cur_val == true
                # we just ended a run of "True" values
                append!( out, (pixels[start_idx], pixels[i - 1]) )
            cur_val = val
            start_idx = i
    return out
end

function generate_line_segments( radius::Int64, center::Tuple{Int64, Int64} )
    """Generate radii of a circle that are a fixed width apart on the circle.
    Args:
      radius: radius of the circle, in pixels
      center: center of the circle (x, y) as tuple
    Returns:
      iterator of points (center, point on circle)
    """
    ang_step = SPACING / radius  # angle step in radians
    for ang in 0:ang_step:2π
        ang += ang_step
        yield (center, (center[0] + radius * math.cos(ang), center[1] + radius * math.sin(ang)))
    end
end

function generate_visible( src_height::Float32, heightmap::AbstractArray{Float32} )
    """Trace a ray and determine if a region is viewable.
    Args:
      tower_height: the elevation in meters above sea level of your antenna
      heightmap: an enumerable of heights in a given direction
    Returns:
      an enumerable of True/False for visibility
    """

    min_angle = -Inf
    for (i, height) in enumerate(iterate(heightmap))
        if tower_height - height == 0
            angle_to_point = 0
        elseif tower_height > height
            angle_to_point = math.atan(i / (tower_height - height))
        else
            angle_to_point = atan( (height - tower_height) / i ) + π/2
        end
        if angle_to_point >= min_angle
            min_angle = angle_to_point
            yield True
        else
            yield False
        end
end
=#

function viewshed()
end

function run_light( dem_file::AbstractString, source_file::AbstractString, resolution::Integer, intensity::Real, source_height::Real=0.0, observer_height::Real=1.75, rarefraction::Real=0.14286;
                    output_path::AbstractString=".\\light_intensity.tiff"  )
 # *VERSIONE CON INTENSITA'/FONTI MULTIPLE
 #= *
    if any( i -> i < 0, intensity )
        throw(DomainError(intenisty, "`intenisty` must be greater than 0."))
    end
 =#
    if intensity < 0
        throw(DomainError(intenisty, "`intenisty` must be greater than 0."))
    end

    src_layer = agd.getlayer(agd.read(source_file), 0)
 # * src_geoms = agd.getgeom.(collect(src_layer))
    src_geom = agd.getgeom(collect(src_layer)[1])

 #= *
    if any( geom -> agd.geomdim(geom) != 0, src_geoms )
        throw(DomainError(source_file, "The shapefile must contain a series of points."))
    end
 =#
    if agd.geomdim(src_geom) != 0
        throw(DomainError(source_file, "The shapefile must contain a point."))
    end

    refsys = agd.getspatialref(src_layer)

    dem = agd.readraster(dem_file)

    if agd.importWKT(agd.getproj(dem)) != refsys
        throw(DomainError("The reference systems are not uniform. Aborting analysis."))
    end

 #= *
    nfeature = 0
    features = agd.getgeom.( agd.getlayer(source, 0) )
    for feature in features
        x_source = agd.getx(feature, 0)
        y_source = agd.gety(feature, 0)
        nfeature += 1
        # The assumption is that the lightsource can only project in a downward emisphere (like a streetlight) and so all terrain that is higher than the source will
        #   block the light for the terrain beyond.
        vis_mat = Viewshed.viewshed( dem, x_source, y_source, source_height )
    end
 =#
    x_source = agd.getx(src_geom, 0)
    y_source = agd.gety(src_geom, 0)
    # The assumption is that the lightsource can only project in a downward emisphere (like a streetlight) and so all terrain that is higher than the source will
    #   block the light for the terrain beyond.
    vis_mat = Viewshed.viewshed( dem, x_source, y_source, source_height )
    rows, cols = size(vis_mat)
    noDataValue = -9999.f0
    data = fill(noDataValue, rows, cols)
    @inbounds for r in 1:rows, c in 1:cols
        if dem[r, c] != noDataValue
            if !vis_mat[r, c]
                data[r, c] = 0.f0
            else
                # I₂ = ( d₁^2 / d₂^2 ) * I₁
                #  d₁ = 1.0 -> I₂ = I₁ * ( 1 / d₂^2 ) = I₁ / d₂^2
                data[r, c] = intensity / Viewshed.distance( toCoords(dem, r, c), (x_source, y_source) )^2
            end
        end
    end
    Functions.writeRaster(data, agd.getdriver("GTiff"), agd.getgeotransform(dem), resolution, refsys, noDataValue, output_path)
end

end # module

# ------------------------------------------------------------------------------- TESTING ---------------------------------------------------------------------------------------------

#= CON "GeoArrays.jl"
    using GeoArrays
    const ga = GeoArrays

    rotate_x( xp::In64, yp::In64, xc::In64, yc::In64, θ::In64 ) = round( Int64, (xp - xc)cos(deg2rad(θ)) - (yp - yc)sin(deg2rad(θ)) + xc )
    rotate_y( xp::In64, yp::In64, xc::In64, yc::In64, θ::In64 ) = round( Int64, (xp - xc)sin(deg2rad(θ)) + (yp - yc)cos(deg2rad(θ)) + yc )
    rotate_point( xp::In64, yp::In64, xc::In64, yc::In64, θ::In64 )::Tuple{Int64, Int64} = θ == 0 ? (xp, yp) : ( rotate_x( xp, yp, xc, yc, θ ), rotate_y( xp, yp, xc, yc, θ ) )
    distance( x0::AbstractFloat, y0::AbstractFloat, x1::AbstractFloat, y1::AbstractFloat ) = √( ( x1 - x0 )^2 + ( y1 - y0 )^2 )
    distance( p0::Tuple{T, T}, p1::Tuple{T, T} ) where {T <: AbstractFloat} = √( ( p1[1] - p0[1] )^2 + ( p1[2] - p0[2] )^2 )

    function toCoords(dem, r, c)
        return Tuple{Float64, Float64}(ga.coords(dem, [r, c]))
    end

    function toIndexes(dem, x, y)
        return Tuple{Int64, Int64}(ga.indices(dem, [x, y]))
    end

    # VERSIONE CON "GeoArrays"
    function viewshed( dtm::GeoArray{Float32}, x0::Float64, y0::Float64, h0::Float32 )
        # Source cell
        r0, c0 = toIndexes(dtm, x0, y0)
        # Final point of the right horizontal profile
        rl, cl = size(dtm)
        # Total height of the source accounting for terrain
        th0 = dtm[r0, c0][1] + h0
        # Visibility matrix
        vizmat = falses(size(dtm)...)
        # The source is always visible
        vizmat[r0, c0] = true
        # Check the visibility along 360 lines radially expanding from the source at an angle of 1° between two adjacent.
        for α in 1:89, β in [0, 90, 180, 270]
            # Final point of the profile of the line with agle `α + β`
            rn, cn = rotate_point( rl, c0, r0, c0, α + β )
            Δr = rn - r0
            Δc = cn - c0
            # Values to add to row and column of a cell on the line to pass to another cell of the line
            r_inc, c_inc = (Δr, Δc) ./ max( abs(Δr), abs(Δc) )
            # Indexes of the first cell after the source cell
            r = r0 + r_inc
            c = c0 + c_inc
            r1 = round(Int64, r)
            c1 = round(Int64, c)
            # The first cell after the source is always visible
            vizmat[r1, c1] = true
            # Slope between the source and the first cell
            slope = ( dtm[r1, c1] - dtm[r0, c0] ) / distance( toCoords(dtm, r1, c1), (0.0, 0.0) )
            # Iterate over each cell of the profile on the line
            while r <= rl && r <= rn && c <= cl && c <= cl
                # Calculate the precise indexes of the cell
                rint, cint = round.(Int64, [r, c])
                # Skip cell already known to be visible
                !vizmat[rint, cint] && continue
                # Compute the slope of the new cell
                new_slope = (dtm[rint, cint] - dtm[r0, c0]) / distance( toCoords(dtm, r1, c1), (0.0, 0.0) )
                # If the new slope is greater than the original one the cell is visible
                if new_slope >= slope
                    vizmat[rint, cint] = true
                    slope = new_slope
                end
                # If the cell is higher than the source al cell beyond the current ne will be hidden
                dtm[rint, cint][1] > th0 && break
                r += r_inc
                c += c_inc
            end
        end
        return vizmat
    end

    function run_light( dem, source::Tuple{Float64, Float64}, resolution::Integer, intensity::Real, source_height::Real=0.0, observer_height::Real=1.75, rarefraction::Real=0.14286;
                        output_path::AbstractString=".\\light_intensity.tiff"  )
        if intensity < 0
            throw(DomainError(intenisty, "`intenisty` must be greater than 0."))
        end

        x_source, y_source = source
        # The assumption is that the lightsource can only project in a downward emisphere (like a streetlight) and so all terrain that is higher than the source will
        #   block the light for the terrain beyond.
        vis_mat = viewshed( dem, x_source, y_source, source_height )
        rows, cols = size(vis_mat)
        noDataValue = -9999.f0
        data = fill(noDataValue, rows, cols)
        @inbounds for r in 1:rows, c in 1:cols
            if dem[r, c] != noDataValue
                if !vis_mat[r, c]
                    data[r, c] = 0.f0
                else
                    # I₂ = ( d₁^2 / d₂^2 ) * I₁
                    #  d₁ = 1.0 -> I₂ = I₁ * ( 1 / d₂^2 ) = I₁ / d₂^2
                    data[r, c] = intensity / distance( toCoords(dem, r, c), (x_source, y_source) )^2
                end
            end
        end
        return data
        #   Functions.writeRaster( data, agd.getdriver("GTiff"), agd.getgeotransform(dem), resolution, refsys, noDataValue, output_path, false )
    end

    dtmga = ga.read(split( @__DIR__ , "\\Porting\\")[1] * "\\Mappe\\DTM_32.tiff")
    #   dtm = ga.read(split( @__DIR__ , "\\Porting\\")[1] * "\\Mappe\\DTM_wgs84.tiff")
    src = (726454.9302346368, 5.025993899219433e6)
    #   src = (11.930065824163105,45.425861311724816)
    run_light( dtm, src, 25, 15.0, 1.0f0, output_path="C:\\Users\\DAVIDE-FAVARO\\Desktop\\test_light.tiff" )
=#

# CON "Rasters"

#=   QUESTE FORMULE NON TENGONO CONTO CHE LA "Y" E' DECRESCENTE, MA FORSE NON CAMBIA PERCHE' L'ANGOLO E' LO STESSO MA SPOSTATO/INVERTITO ?
#       forse ho sistemato
rotate_x( xp::Float64, yp::Float64, xc::Float64, yc::Float64, θ::Int64 ) = xc + ( (xp - xc)cos(deg2rad(θ)) - (yp - yc)sin(deg2rad(θ)) )
rotate_y( xp::Float64, yp::Float64, xc::Float64, yc::Float64, θ::Int64 ) = yc - ( (xp - xc)sin(deg2rad(θ)) + (yp - yc)cos(deg2rad(θ)) )
rotate_point( xp::Float64, yp::Float64, xc::Float64, yc::Float64, θ::Int64 ) = θ == 0 ? (xp, yp) : ( rotate_x( xp, yp, xc, yc, θ ), rotate_y( xp, yp, xc, yc, θ ) )
distance( x0::Float64, y0::Float64, x1::Float64, y1::Float64 ) = √( ( x1 - x0 )^2 + ( y1 - y0 )^2 )
distance( p0::Tuple{T, T}, p1::Tuple{T, T} ) where {T <: Float64} = √( ( p1[1] - p0[1] )^2 + ( p1[2] - p0[2] )^2 )
Base.getindex( collection::Raster, x::Float64, y::Float64 ) = collection[X(Near(x)), Y(Near(y))][1]
Base.setindex!( collection::Raster, v, x::Float64, y::Float64, ) = collection[X(Near(x)), Y(Near(y))] .= v
condition( x::Float64, xn::Float64, y::Float64, yn::Float64, θ::Int64 ) = θ == 0 ? (x <= xn && y >= yn) :
                                                                              θ == 90 ? (x >= xn && y >= yn) :
                                                                                θ == 180 ? (x >= xn && y <= yn) :
                                                                                θ == 270 ? (x <= xn && y <= yn) : throw(DomainError(β, "`β` must either be 0, 90, 180 or 270."))

findNearest( dtm::Raster{Float32, 3}, x0::Float64, y0::Float64 ) = findmin( x -> abs(x - x0), dtm.dims[1].val )[2], findmin( y -> abs(y - y0), dtm.dims[2].val )[2]

function viewshed( dtm::Raster{Float32, 3}, x0::Float64, y0::Float64, h0::Float32 )
    # Dimensions of the raster.
        # N.B. Y axis is in decreasing order.
    xmin = dtm.dims[1][1]
    xmax = dtm.dims[1][end]
    ymax = dtm.dims[2][1]
    ymin = dtm.dims[2][end]
    # Displacement from a cell to the next alogn the two axis.
    xstep = dtm.dims[1][2] - xmin
    ystep = dtm.dims[2][2] - ymax
    # Total height of the source accounting for terrain.
    th0 = dtm[x0, y0] + h0
    # Visibility matrix.
    vizmat = Raster(falses(size(dtm)...), dtm.dims)
    # The source is always visible.
    vizmat[x0, y0] = true
    # Check the visibility along 360 lines radially expanding from the source at an angle of 1° between two adjacent ones.
        # Iterating over the four quadrants and then on the angles in each of them makes easier to check the coordinates.
    for α in [0, 90, 180, 270], β in 0:89
        # Final point of the profile of the line with agle `α + β`.
        xn, yn = rotate_point(xmax, y0, x0, y0, α + β)
 #=
        # Diplacement between the center and the final point along the two axis.
        Δx = xn - x0
        Δy = yn - y0
 =#
        # Values to add to row and column of a cell on the line to pass to another cell of the line.
        x_inc, y_inc = ( (xn - x0, yn - y0) ./ max(abs(xn - x0), abs(yn - y0)) ) .* (xstep, ystep)
        # Coordinates of the first cell after the source cell.
        x = x0 + x_inc
        y = y0 + y_inc
        # The first cell after the source is always visible.
        vizmat[x, y] = true
        # Slope between the source and the first cell.
        slope = (dtm[x, y] - dtm[x0, y0]) / distance((x, y), (0.0, 0.0))
        # Go to the following cell in the "profile".
        x += x_inc
        y += y_inc
        # Iterate over each cell of the profile on the line.
        while xmin <= x <= xmax && ymin <= y <= ymax && condition(x, xn, y, yn, α)
            # Skip cells already known to be visible.
            vizmat[x, y] && continue
            # Compute the slope of the new cell.
            new_slope = (dtm[x, y] - dtm[x0, y0]) / distance((x, y), (0.0, 0.0))
            # If the new slope is greater than the original one the cell is visible.
            if new_slope >= slope
                vizmat[x, y] = true
                slope = new_slope
            end
            # If the cell is higher than the source, all cells beyond it will be hidden.
            dtm[x, y] > th0 && break
            # Go to the following cell.
            x += x_inc
            y += y_inc
        end
    end
    return vizmat
end

function run_light( dem::Raster{Float32, 3}, source::Tuple{Float64, Float64}, resolution::Integer, intensity::Real, source_height::Real=0.0, observer_height::Real=1.75, rarefraction::Real=0.14286;
                    output_path::AbstractString=".\\light_intensity.tiff"  )
    if intensity < 0
        throw(DomainError(intenisty, "`intenisty` must be greater than 0."))
    end

    x_source, y_source = source
    # The assumption is that the lightsource can only project in a downward emisphere (like a streetlight) and so all terrain that is higher than the source will
    #   block the light for the terrain beyond.
    vis_mat = viewshed(dem, x_source, y_source, source_height)
    rows, cols = size(vis_mat)
    noDataValue = -9999.f0
    data = fill(noDataValue, rows, cols)
    @inbounds for r in 1:rows, c in 1:cols
        if dem[r, c] != noDataValue
            if !vis_mat[r, c]
                data[r, c] = 0.0f0
            else
                # I₂ = ( d₁^2 / d₂^2 ) * I₁
                #  d₁ = 1.0 -> I₂ = I₁ * ( 1 / d₂^2 ) = I₁ / d₂^2
                data[r, c] = intensity / distance( toCoords(dem, r, c), (x_source, y_source) )^2
            end
        end
    end
    return data
    #   Functions.writeRaster( data, agd.getdriver("GTiff"), agd.getgeotransform(dem), resolution, refsys, noDataValue, output_path, false )
end

x0, y0 = src
xmin = dtmr.dims[1][1]
xmax = dtmr.dims[1][end]
ymax = dtmr.dims[2][1]
ymin = dtmr.dims[2][end]
xstep = dtmr.dims[1][2] - xmin
ystep = dtmr.dims[2][2] - ymax
th0 = dtmr[x0, y0] + 1.0f0

# Punti profilo a 184°
xn, yn = rotate_point(xmax, y0, x0, y0, 184)
Δx = xn - x0
Δy = yn - y0
x_inc, y_inc = ( (Δx, Δy) ./ max(abs(Δx / xstep), abs(Δy / ystep)) )
x = x0 + x_inc
y = y0 + y_inc
arr1 = [ (x, y) for (x, y) in zip(x0 + 2x_inc:x_inc:xn, y0 + 2y_inc:y_inc:yn) ]
plot(arr1)

# Punti profilo a 96°
xn, yn = rotate_point(xmax, y0, x0, y0, 96)
Δx = xn - x0
Δy = yn - y0
x_inc, y_inc = ( (Δx, Δy) ./ max(abs(Δx / xstep), abs(Δy / ystep)) )
x = x0 + x_inc
y = y0 + y_inc
arr2 = [ (x, y) for (x, y) in zip(x0 + 2x_inc:x_inc:xn, y0 + 2y_inc:y_inc:yn) ]
plot!(arr2)

# Punti profilo a 48°
xn, yn = rotate_point(xmax, y0, x0, y0, 48)
Δx = xn - x0
Δy = yn - y0
x_inc, y_inc = ( (Δx, Δy) ./ max(abs(Δx / xstep), abs(Δy / ystep)) )
x = x0 + x_inc
y = y0 + y_inc
arr3 = [ (x, y) for (x, y) in zip(x0 + 2x_inc:x_inc:xn, y0 + 2y_inc:y_inc:yn) ]
plot(arr3)

# Punti profilo a 290°
xn, yn = rotate_point(xmax, y0, x0, y0, 290)
Δx = xn - x0
Δy = yn - y0
x_inc, y_inc = ( (Δx, Δy) ./ max(abs(Δx / xstep), abs(Δy / ystep)) )
x = x0 + x_inc
y = y0 + y_inc
arr4 = [ (x, y) for (x, y) in zip(x0 + 2x_inc:x_inc:xn, y0 + 2y_inc:y_inc:yn) ]
plot!(arr4)

# Forma dei profili
plot([dtmr[xy...] for xy in arr1])
# Slope (La moltiplicazione per 10000 serve a rendere comprensibile la curva)
plot!([ (dtmr[arr1[i]...] - dtmr[x0, y0]) / distance(arr1[i], (x0, y0)) * 10000 for i in 2:length(arr1) ])

plot([dtmr[xy...] for xy in arr2])
plot([dtmr[xy...] for xy in arr3])
plot([dtmr[xy...] for xy in arr4])

plot(dtmr[4000:end, 6000:end])

function viewshed( dtm::Raster{Float32, 3}, x0::Float64, y0::Float64, h0::Float32 )
    # Dimensions of the raster.
        # N.B. Y axis is in decreasing order.
    xmin = dtm.dims[1][1]
    xmax = dtm.dims[1][end]
    ymax = dtm.dims[2][1]
    ymin = dtm.dims[2][end]
    # Displacement from a cell to the next alogn the two axis.
    xstep = dtm.dims[1][2] - xmin
    ystep = dtm.dims[2][2] - ymax
    # Total height of the source accounting for terrain.
    th0 = dtm[x0, y0] + h0
    # Visibility matrix.
    vizmat = Raster(zeros(Int64, size(dtm)...), dtm.dims, missingval=0)
    # The source is always visible.
    vizmat[x0, y0] = 1
    # Check the visibility along 360 lines radially expanding from the source at an angle of 1° between two adjacent ones.
        # Iterating over the four quadrants and then on the angles in each of them makes easier to check the coordinates.
    for α in  1:89, β in [0, 90, 180, 270]
        println("Θ = $(α + β)")
        # Final point of the profile of the line with agle `α + β`.
        xn, yn = rotate_point(xmax, y0, x0, y0, α + β)
        # Diplacement between the center and the final point along the two axis.
        Δx = xn - x0
        Δy = yn - y0
        # Values to add to row and column of a cell on the line to pass to another cell of the line.
        x_inc, y_inc = ( (Δx, Δy) ./ max(abs(Δx / xstep), abs(Δy / ystep)) )
        # Coordinates of the first cell after the source cell.
        x = x0 + x_inc
        y = y0 + y_inc
        # The first cell after the source is always visible.
        vizmat[x, y] = 1
        # Slope between the source and the first cell.
        slope = (dtm[x, y] - dtm[x0, y0]) / distance((x0, y0), (x, y))
        # Iterate over each cell of the profile on the line.
        for (y, x) in zip( (y0 + 2y_inc):y_inc:yn, (x0 + 2x_inc):x_inc:xn )
            # If the coordinates are outside of the raster break.
            ( x < xmin || x > xmax || y < ymin || y > ymax ) && break
            # If the cell has already been visited and set visible and it doesn't hold a missing value continue to the next cell
            ( vizmat[x, y] != 1 && dtm[x, y] == dtm.missingval ) && continue
            # Compute the slope of the new cell.
            new_slope = (dtm[x, y] - dtm[x0, y0]) / Float32(distance((x0, y0), (x, y)))
            # If the new slope is greater than the original one the cell is visible.
            if new_slope >= slope
                vizmat[x, y] = 1
                slope = new_slope
            end
            # If the cell is higher than the source, all cells beyond it will be hidden.
            dtm[x, y] > th0 && break
        end
 #=
        while xmin <= x <= xmax && ymin <= y <= ymax && condition(x, xn, y, yn, α)
            # Skip cells already known to be visible.
            vizmat[x, y] == 1 && continue
            # Compute the slope of the new cell.
            new_slope = (dtm[x, y] - dtm[x0, y0]) / distance((x0, y0), (x, y))
            # If the new slope is greater than the original one the cell is visible.
            if new_slope >= slope
                vizmat[x, y] = 1
                slope = new_slope
            else
                vizmat[x, y] = -1
            end
            # If the cell is higher than the source, all cells beyond it will be hidden.
            dtm[x, y] > th0 && break
            # Go to the following cell.
            x += x_inc
            y += y_inc
        end
 =#
    end
    return vizmat
end

# FUNZIONA ED E' RAGIONEVOLMENTE VELOCE MA NON SONO CERTO CHE IL RISULTATO SIA CORRETTO
function run_light( dem::Raster{Float32, 3}, intensity::Float32, x0::Float64, y0::Float64, h0::Float32 )
    # Dimensions of the raster.
        # N.B. Y axis is in decreasing order.
    xmin = dem.dims[1][1]
    xmax = dem.dims[1][end]
    ymax = dem.dims[2][1]
    ymin = dem.dims[2][end]
    # Displacement from a cell to the next alogn the two axis.
    xstep = dem.dims[1][2] - xmin
    ystep = dem.dims[2][2] - ymax
    # Total height of the source accounting for terrain.
    th0 = dem[x0, y0] + h0
    # Light intensity matrix.
    data = Raster(fill(dem.missingval, size(dem)...), dem.dims, missingval=dem.missingval )
    # The intensity of the light at the source is given
    data[x0, y0] = intensity
    # Check the visibility along 360 lines radially expanding from the source at an angle of 1° between two adjacent ones.
        # Iterating over the four quadrants and then on the angles in each of them makes easier to check the coordinates.
    for α in 1:89, β in [0, 90, 180, 270]
        # Final point of the profile of the line with agle `α + β`.
        xn, yn = rotate_point(xmax, y0, x0, y0, α + β)
        # Diplacement between the center and the final point along the two axis.
        Δx = xn - x0
        Δy = yn - y0
        # Values to add to row and column of a cell on the line to pass to another cell of the line.
        x_inc, y_inc = ( (Δx, Δy) ./ max(abs(Δx / xstep), abs(Δy / ystep)) )
        # Coordinates of the first cell after the source cell.
        x = x0 + x_inc
        y = y0 + y_inc
        # Distance between source and first cell of the profile
        dist = Float32(distance((x0, y0), (x, y)))
        # Slope between the source and the first cell.
        slope = (dem[x, y] - dem[x0, y0]) / dist
        # Value of light intensity at the cell of coordinates (`x`, `y`).
            # I₂ = I₁ * ( d₁^2 / d₂^2 )  ->  d₁ = 1.0  ->  I₂ = I₁ * ( 1 / d₂^2 ) = I₁ / d₂^2.
        data[x, y] = intensity / dist^2
        # Iterate over each cell of the profile on the line.
        for (y, x) in zip( (y0 + 2y_inc):y_inc:yn, (x0 + 2x_inc):x_inc:xn )
            # If the coordinates are outside of the raster break.
            ( x < xmin || x > xmax || y < ymin || y > ymax ) && break
            # If the cell has already been visited and set visible and it doesn't hold a missing value continue to the next cell
            ( ( data[x, y] != 0.0f0 && data[x, y] != dem.missingval ) || dem[x, y] == dem.missingval ) && continue
            # Distance of the cell from the source.
            dist = Float32(distance((x0, y0), (x, y)))
            # Compute the slope of the new cell.
            new_slope = (dem[x, y] - dem[x0, y0]) / dist
            # If the new slope is greater than the original one the cell is visible.
            if new_slope >= slope
                data[x, y] = intensity / dist^2
                slope = new_slope
            else
                data[x, y] = 0.0f0
            end
            # If the cell is higher than the source, all cells beyond it will be hidden.
            dem[x, y] > th0 && break
        end
 #=
        while xmin <= x <= xmax && ymin <= y <= ymax && condition(x, xn, y, yn, α)
            # Skip cells already known to be visible.
            ( ( data[x, y] != 0.0f0 && data[x, y] != dem.missingval ) || dem[x, y] == dem.missingval ) && continue
            # Distance of the cell from the source.
            dist = distance((x, y), (0.0, 0.0))
            # Compute the slope of the new cell.
            new_slope = (dem[x, y] - dem[x0, y0]) / dist
            # If the new slope is greater than the original one the cell is visible.
            if new_slope >= slope
                data[x, y] = intensity / dist^2
                slope = new_slope
            else
                data[x, y] = 0.0f0
            end
            # If the cell is higher than the source, all cells beyond it will be hidden.
            dem[x, y] > th0 && break
            # Go to the following cell.
            x += x_inc
            y += y_inc
        end
 =#
    end
    return data
end

function run_light( dem::Raster{Float32, 3}, intensity::Float32, x0::Float64, y0::Float64, h0::Float32 )
    rmax, cmax = size(dem)
    r0, c0 = findNearest(dem, x0, y0)
    # Total height of the source accounting for terrain.
    th0 = dem[r0, c0] + h0
    # Light intensity matrix.
    data = Raster(fill(dem.missingval, rmax, cmax), dem.dims)
    # The intensity of the light at the source is given
    data[r0, c0] = intensity
    # Check the visibility along 360 lines radially expanding from the source at an angle of 1° between two adjacent ones.
        # Iterating over the four quadrants and then on the angles in each of them makes easier to check the coordinates.
    for α in 0:89, β in [0, 90, 180, 270]
        # Final point of the profile of the line with agle `α + β`.
        rn, cn = rotate_point(rmax, c0, r0, c0, α + β)
        # Diplacement between the center and the final point along the two axis.
        Δr = rn - r0
        Δc = cn - c0
        # Values to add to row and column of a cell on the line to pass to another cell of the line.
        r_inc, c_inc = ( (Δr, Δc) ./ max(abs(Δr), abs(Δc)) )
        # Coordinates of the first cell after the source cell.
        r = r0 + r_inc
        c = c0 + c_inc
        # Distance between source and first cell of the profile
        dist = Float32(distance((r, c), (0.0, 0.0)))
        # Slope between the source and the first cell.
        slope = (dem[r, c] - dem[r0, c0]) / dist
        # Value of light intensity at the cell of coordinates (`x`, `y`).
            # I₂ = I₁ * ( d₁^2 / d₂^2 )  ->  d₁ = 1.0  ->  I₂ = I₁ * ( 1 / d₂^2 ) = I₁ / d₂^2.
        data[r, c] = intensity / dist^2
        # Iterate over each cell of the profile on the line.
        for (r, c) in zip( (c0 + 2c_inc):c_inc:cn, (r0 + 2r_inc):r_inc:rn )
            # If the coordinates are outside of the raster break.
            ( r < rmin || r > rmax || c < cmin || c > cmax ) && break
            # If the cell has already been visited and set visible and it doesn't hold a missing value continue to the nert cell
            ( ( data[r, c] != 0.0f0 && data[r, c] != dem.missingval ) || dem[r, c] == dem.missingval ) && continue
            # Distance of the cell from the source.
            dist = Float32(distance((r, c), (0.0, 0.0)))
            # Compute the slope of the new cell.
            new_slope = (dem[r, c] - dem[r0, c0]) / dist
            # If the new slope is greater than the original one the cell is visible.
            if new_slope >= slope
                data[r, c] = intensity / dist^2
                slope = new_slope
            else
                data[r, c] = 0.0f0
            end
            # If the cell is higher than the source, all cells beyond it will be hidden.
            dem[r, c] > th0 && break
        end
 #=
        while xmin <= x <= xmax && ymin <= y <= ymax && condition(x, xn, y, yn, α)
            # Skip cells already known to be visible.
            ( ( data[x, y] != 0.0f0 && data[x, y] != dem.missingval ) || dem[x, y] == dem.missingval ) && continue
            # Distance of the cell from the source.
            dist = distance((x, y), (0.0, 0.0))
            # Compute the slope of the new cell.
            new_slope = (dem[x, y] - dem[x0, y0]) / dist
            # If the new slope is greater than the original one the cell is visible.
            if new_slope >= slope
                data[x, y] = intensity / dist^2
                slope = new_slope
            else
                data[x, y] = 0.0f0
            end
            # If the cell is higher than the source, all cells beyond it will be hidden.
            dem[x, y] > th0 && break
            # Go to the following cell.
            x += x_inc
            y += y_inc
        end
 =#
    end
    return data
end

using Plots
using ProfileView
using Rasters
using Revise
dtmr = Raster(split( @__DIR__ , "\\Porting\\")[1] * "\\Mappe\\DTM_32.tiff")
src = (726454.9302346368, 5.025993899219433e6)

@code_warntype viewshed(dtmr, src[1], src[2], 1.0f0)
@time vizmat = viewshed(dtmr, src[1], src[2], 1.0f0)

@code_warntype run_light(dtmr, 15.0f0, src[1], src[2], 0.0f0)
@time intmat = run_light(dtmr, 15.0f0, src[1], src[2], 0.0f0)

@code_warntype run_light(dtmr, src, 25, 15.0, 1.0f0, output_path="C:\\Users\\DAVIDE-FAVARO\\Desktop\\test_light.tiff")
intmat = run_light(dtmr, src, 25, 15.0, 1.0f0, output_path="C:\\Users\\DAVIDE-FAVARO\\Desktop\\test_light.tiff")

ProfileView.@profview viewshed(dtmr, src[1], src[2], 1.0f0)
ProfileView.@profview run_light(dtmr, 15.0f0, src[1], src[2], 0.0f0)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

#=
    ENV["PYTHON"] = raw"C:\Users\DAVIDE-FAVARO\.julia\conda\3\python.exe"
    using Pkg
    Pkg.build("PyCall")
    using PyCall
    import PyCall.Conda as cnd

    pushfirst!(PyVector(pyimport("sys")."path"), "C:\\Users\\DAVIDE-FAVARO\\.julia\\conda\\3\\Library")
    pushfirst!(PyVector(pyimport("sys")."path"), @__DIR__)
=#
using PyCall
# Serve per indicare a Python dove trovare i moduli aggiuntivi scaricati con Conda (in questo caso "qgis")
#   pushfirst!(PyVector(pyimport("sys")."path"), "C:\\Users\\DAVIDE-FAVARO\\.julia\\conda\\3\\Library\\python")
# pushfirst!(PyVector(pyimport("sys")."path"), "C:\\Users\\Lenovo\\.julia\\conda\\3\\Library\\python")
pushfirst!(PyVector(pyimport("sys")."path"), "C:\\Users\\DAVIDE-FAVARO\\.julia\\conda\\3\\Library\\python")
dtm = split( @__DIR__ , "\\Porting\\")[1] * "\\Mappe\\DTM_wgs84.tiff"
# output_path = "D:\\Z_Tirocinio_Dati\\test_viewshed.tiff"
output_path = "C:\\Users\\DAVIDE-FAVARO\\Desktop\\DatiRaster\\test_viewshed.tiff"
x_source, y_source = (11.930065824163105,45.425861311724816)
coords = string(x_source)*","*string(y_source)*"[4326]"
source_height = 1.0
observer_height = 1.75
rarefraction = 0.14286

#= ESEMPI DI USO DI Qgis DA Python
https://gis.stackexchange.com/questions/129513/accessing-processing-with-python
https://gis.stackexchange.com/questions/279874/using-qgis-3-processing-algorithms-from-pyqgis-standalone-scripts-outside-of-gu
=#
py"""
import sys
import qgis
from qgis.core import *

# QGIS installation path
# QgsApplication.setPrefixPath("C:\\Program Files\\QGIS 3.16.16\\bin\\qgis-ltr-bin.exe", True)
QgsApplication.setPrefixPath("C:\\Program Files\\QGIS 3.22.0\\bin\\qgis-bin.exe", True)
# Reference the application without using a GUI
qgs = QgsApplication([], False)
# Load providers
qgs.initQgis()

sys.path.append("C:\\Program Files\\QGIS 3.22.0\\apps\\qgis\\plugins")

import processing
from porcessing.core.Processing import Processing
Processing.initialize()

dem = QgsRasterLayer(dtm, "dem")

viewshed = processing.run(
    "grass7:r.viewshed",
    {
        "-c" : True, # Account for earth curvature
        "-r" : False, # Consider atmospheric rarefraction
        "-b" : True,  # Return boolean map ( invisible = 0/false, visible = 1/true )
        "-e" : False, # Assign (non boolean) values to the result ( invisible = NULL, visible = current elevation - viewpoint elevation )
        "input" : dem,
        "output" : output_path,
 # AL POSOTO DI "4326" BISOGNEREBBE METTERE UN COMANDO PER OTTENERE IL CODICE DEL SISTEMA DI RIFERIMENTO UTILIZZATO MA NON TROVO UN COMANDO COSI' IN ArchGDAL
        "coordinates" : coords,
        "observer_elevation" : source_height,
        "target_elevation" : observer_height,
        "max_distance" : -1, # -1 reppresents infinity as maximum possible distance
        "refraction_coeff" : rarefraction,
        #   "memory" : memory,
        "GRASS_RASTER_FORMAT_META" : "",
        "GRASS_RASTER_FORMAT_OPT" : "",
        "GRASS_REGION_CELLSIZE_PARAMETER" : 0,
        #   "GRASS_REGION_PARAMETER" :grass_area,
    }
)

viewshed

# Remove providers and free memory
qgs.exitQgis()
"""

function run_light( dem::AbstractString, source, resolution::Integer, intenisty::Real, source_height::Real=0.0, observer_height::Real=1.75, rarefraction::Real=0.14286,
                    output_path::AbstractString=".\\light_intensity.tiff"  )

    if any( i -> i < 0, intensity )
        throw(DomainError(intenisty, "`intenisty` must be greater than 0."))
    end

    src_layer = agd.getlayer(source, 0)
    src_geoms = agd.getgeom.(collect(src_layer))

    if any( geom -> agd.geomdim(geom) != 0, src_geoms )
        throw(DomainError(source, "`source` must be a point."))
    end

    refsys = agd.getspatialref(src_layer)

    if agd.importWKT(agd.getproj(dem)) != refsys
        throw(DomainError("The reference systems are not uniform. Aborting analysis."))
    end

    rows, cols = size(dem)

    features = agd.getgeom.( agd.getlayer(source, 0) )
    data = fill(0.0f0, rows, cols)
    for feature in features
        x_source = agd.getx(feature, 0)
        y_source = agd.gety(feature, 0)

 # L'IDEALE SAREBBE AVERE GLI IMPORT E I VARI SETUP (CHE AL MOMENTO MANCANO) FUORI DAL CICLO (MANTENEDOLO UN CICLO DI JULIA) E FARE SOLO I CALCOLI ALL'INTERNO
        py"""
        import qgis
        from qgis.core import QgsApplication, QgsProcessingFeedback, QgsVectorLayer

        QgsApplication.setPrefixPath('C:\\Program Files\\QGIS 3.16.16\\bin', True)
        qgs = QgsApplication([], False)
        qgs.initQgis()

        import processing
        from porcessing.core.Processing import Processing
        Processing.initialize()

        viewshed = processing.run(
            "grass7:r.viewshed",
            {
                "-c" : false, # Account for earth curvature
                "-r" : false, # Consider atmospheric rarefraction
                "-b" : true,  # Return boolean map ( invisible = 0/false, visible = 1/true )
                "-e" : false, # Assign (non boolean) values to the result ( invisible = NULL, visible = current elevation - viewpoint elevation )
                "input" : dem,
                "output" : output_path,
         # AL POSOTO DI "4326" BISOGNEREBBE METTERE UN COMANDO PER OTTENERE IL CODICE DEL SISTEMA DI RIFERIMENTO UTILIZZATO MA NON TROVO UN COMANDO COSI' IN ArchGDAL
                "coordinates" : string(x_source)*","*string(y_source)*"[4326]",
                "observer_elevation" : source_height,
                "target_elevation" : observer_height,
                "max_distance" : -1, # -1 reppresents infinity as maximum possible distance
                "refraction_coeff" : rarefraction,
                #   "memory" : memory,
                "GRASS_RASTER_FORMAT_META" : "",
                "GRASS_RASTER_FORMAT_OPT" : "",
                "GRASS_REGION_CELLSIZE_PARAMETER" : 0,
                #   "GRASS_REGION_PARAMETER" :grass_area,
            }
        )
        qgs.exitQgis()
        """

 # RIGHE E COLONNE DELL MATRICE  viewshed DOVREBBERO ESSERE LE STESSE DEL DEM, QUINDI ANCHE COORDINATE MINIME E MASSIME.

        # Matrix of results
        for row in 1:rows, col in 1:cols
            # Coordinates of the point
            x = (col * resolution) + x_min + (resolution / 2)
            y = (row * resolution) + y_min + (resolution / 2)
            # If a point is visible from the source compute the intensity of light at that point.
            if viewshed[row, col] == true
                # "data" contains all zeroes after intialization, so instead of making a case to set its values if it's the first
                #   feature we simply add the maximum between the intensity and 0.
                data[row, col] += max( # Maximum between the new intensity and 0
                    intensity / ( # Diatance^2
                            √( # Distance (x, y) of form (x_source, y_source)
                                (y - y_source)^2 +
                                (x - x_source)^2
                            )
                        )^2,
                    0
                )
            end
        end
    end
    Functions.writeRaster(data, agd.getdriver("GTiff"), geotransform, resolution, refsys, noData, output_path, false)

 #=
     # QUESTO O QUALCHE ALTRO MODO PER OTTENERE I LIMITI DEL RASTER VIEWSHED
        # Define the area of the `dtm` that coincides with the area of the `viewshed`
        minX, maxY, maxX, maxY = getSidesDistances()
        row_begin, col_begin = toIndexes(dtm, minX, minY)
        row_end, col_end = toIndexes(dtm, maxX, maxY)

        valNoData = -9999
        gtiff_driver = agd.getdriver("GTiff")
        target_ds = agd.create( output_path, gtiff_driver, row_end-row_begin, col_end-col_begin, 1, agd.GDAL.GDT_Float32 )
        agd.setgeotransform!(target_ds, [ minX, resolution, 0.0, maxY, 0.0, -resolution ])
        agd.setproj!( target_ds, refsys )
     """ NON SO QUALE SIA IL COMANDO PER SETTARE I METADATI CON `ArchGDAL`
        target_ds.SetMetadata(
            Dict(
                "credits" => "Envifate - Francesco Geri, Oscar Cainelli, Paolo Zatelli, Gianluca Salogni, Marco Ciolli - DICAM Università degli Studi di Trento - Regione Veneto",
                "modulo" => "Analisi inquinamento luminoso",
                "descrizione" => "Simulazione di inquinamento luminoso da sorgente puntuale singola o multipla",
                "srs" => refsys,
                "data" => today()
            )
        )
     """
        band1 = agd.getband(target_ds, 1)
        agd.setnodatavalue!( band1, Float64(valNoData) )
        agd.fillraster!(band1, valNoData)
        band = agd.read(band1)

        # Compute the light intensity in said area of the dtm
        for row in row_begin:row_end
            for col in col_begin:col_end
                if viewshed[row, col] == 1
                    x, y = toCoords(dtm, row, col)
                    dist = √( (y - y_source)^2 + (x - x_source)^2 )
                    new_intensity = intensity / dist

                    res = new_intensity > 0 ? new_intensity : 0
                    if nfeature == 1
                        band[row, col] = res
                    else
                        band[row, col] += res
                    end
                else
                    band[row, col] = 0
                end
            end
        end
    end

    # MANCA SCRIVERE IL RASTER
 =#

end
=#

max_stars_repo_name: DavideFavaro/Tirocinio

max_issues_repo_name: DavideFavaro/Tirocinio

max_forks_repo_name: DavideFavaro/Tirocinio


text:
"""
The syntax for a while is

    while *condition*
        *loop body*
    end
"""

n = 1
while n <= 10
    n += n
    println(n)
end

# My only friends are numbers 😢
friends = [1, 2, 3, 4, 5, 6]

i = 1 # This language better have an iterator/higher-order functions STG
while i <= length(friends)
    friend = friends[i]
    println("Hi $friend, it's nice to see that you're still hanging around these parts.")
    i += 1
end

# Praise
for friend in friends
    println(friend)
end

# Julia has "ranges"
example_range = 1:100 # UnitRange{Int64}
another_range = 'a':'z' # StepRange{Char, Int64}

for letter in 'a':'z'
    print(letter)
end
println()

# Julia also has unicode overloading
for letter = 'a':'z'
    print(letter)
end
println()

# That's a Unicode U+2208 "Element of" symbol
for letter ∈ 'a':'z'
    print(letter)
end
println()

# Make a 5x5 addition table
m, n = 5, 5
table = zeros(m, n)
println(typeof(table)) # Matrix{Float64}, NOT a 2D array
for i in 1:m
    for j in 1:n
        table[i, j] = i + j
    end
end
table

# Same thing but condensed
table = zeros(m, n)
for i in 1:m, j in 1:n
    table[i, j] = i+j
end
table

# More exploration of condensing loops.
combinations_2 = []
for f1 in friends, f2 in friends
    if f1 != f2
        push!(combinations_2, (f1, f2))
    end
end

combinations_3 = []
#for f1 in friends, f2 = friends, f3 ∈ friends
for f1 = friends, f2 = friends, f3 = friends
    if f1 != f2 && f2 != f3 && f1 != f3
        push!(combinations_3, (f1, f2, f3))
    end
end
combinations_3

# COMPREHENSIONS!!!
table = [i+j for i in 1:m, j in 1:n]

for n in 1:10
    powtable = [i^j for i = 1:n, j in 1:n]
    display(powtable)
end

ext: jl

max_stars_repo_name: AustinTSchaffer/DailyProgrammer

max_issues_repo_name: AustinTSchaffer/DailyProgrammer

max_forks_repo_name: AustinTSchaffer/DailyProgrammer


text:
# assert file to test polynomial implementation
using Base.Test
using Polynomial

pNULL = Poly(Float32[])
p0 = Poly([0])
p1 = Poly([0,0,0,0,0,0,0,0,0,0,0,0,0,1])
p2 = Poly([0,0,1,1])
p3 = Poly([0,0,0,0,1,2,1])
p4 = Poly([0,1,3,3,1])
p5 = Poly([0,0,0,0,0,0,0,0,0,0,0,0,0,1,4,6,4,1])
pN = Poly([0,24,15,87,3,276])
pR = Poly([1//1, -2//1, 3//4])
p1000 = Poly(randn(1000))

@test length(pNULL) == 0
@test length(p1000) == 1000
sprint(show, p1000)
sprint(show, pNULL)

@test p3 == Poly([1,2,1])
@test pN*10 == Poly([240, 150, 870, 30, 2760])
@test pN/10 == Poly([2.4, 1.5, 8.7, 0.3, 27.6])
@test 10*pNULL + pN == pN
@test 10*p0 + pN == pN
@test p5 + 2*p1 == Poly([1,4,6,4,3])
@test 10*pNULL - pN == -pN
@test p0 - pN == -pN
@test p5 - 2*p1 == Poly([1,4,6,4,-1])
@test p2*p2*p2 == p4
@test p2^4 == p5
@test pNULL^3 == pNULL
@test pNULL*pNULL == pNULL

@test polyval(pN, -.125) == 276.9609375
@test polyval(pNULL, 10) == 0
@test polyval(p0, -10) == 0
@test polyval(poly([1//2, 3//2]), 1//2) == 0//1
@test polyder(polyint(pN)) == pN
@test polyder(pR) == Poly([2//1, -2//1])
@test polyint(pNULL,1) == p1
@test polyint(Poly(Rational[3, 2, 1])) == Poly(Rational[1, 1, 1, 0])
@test polyder(p3) == Poly([2,2])
@test polyder(p1) == polyder(p0) == polyder(pNULL) == pNULL

@test poly([-1,-1]) == p3
@test roots(p0)==roots(p1)==roots(pNULL)==[]
@test roots(p2) == [-1]
a_roots = copy(pN.a)
@test all(abs(sort(roots(poly(a_roots))) - sort(a_roots)) .< 1e6)
@test length(roots(p5)) == 4
@test roots(pNULL) == []
@test roots(pR) == [1//2, 3//2]

@test pNULL + 2 == p0 + 2 == 2 + p0 == Poly([2])
@test p2 - 2 == -2 + p2 == Poly([1,-1])
@test 2 - p2 == Poly([-1,1])

p0 = Poly([0])
p1 = Poly([1])
p2 = Poly([4, 2, -3, 6, 5])
p3 = Poly([6, 2, -3, 7])
p4 = p2 * p3
@test divrem(p4, p2) == (p3, zero(p3))
@test p3%p2 == p3
@test all((abs((p2/p3 - Poly([2/3,1/9])).a)) .< eps())
@test divrem(p0,p1) == (p0,p0)
@test divrem(p1,p1) == (p1,p0)
@test divrem(p2,p2) == (p1,p0)
@test divrem(pR, pR) == (one(pR), zero(pR))
@test_throws p1/p0
@test_throws divrem(p0,p0)

@test_throws gcd(p1, p2)
p1 = Poly([0,0,1.0,1.0])
p2 = Poly([0,0,0,1.0,2.0,1.0])
@test gcd(p1, p2) == p1
@test gcd(p2, p1) == p1

#Tests for multivariable support
pX = Poly([1, 2, 3, 4, 5])
pS1 = Poly([1, 2, 3, 4, 5], "s")
pS2 = Poly([1, 2, 3, 4, 5], 's')
pS3 = Poly([1, 2, 3, 4, 5], :s)
@test pX != pS1
@test pS1 == pS2
@test pS1 == pS3
@test_throws pS1 + pX
@test_throws pS1 - pX
@test_throws pS1 * pX
@test_throws pS1 / pX
@test_throws pS1 % pX


text:
using FileIO

function fileio_load(f::File{format"SIXEL"}; kwargs...)
    open(f, "r") do io
        sixel_decode(io)
    end
end

function fileio_load(io::Stream{format"SIXEL"}; kwargs...)
    sixel_decode(RGB{N0f8}, io; kwargs...)
end

sixel_decode(::Type{T}, io::Stream{format"SIXEL"}, dec::SDC=default_decoder(); kwargs...) where T =
    sixel_decode(T, io.io, dec; kwargs...)

function fileio_save(f::File{format"SIXEL"}, img::AbstractArray; kwargs...)
    open(f, "w") do io
        sixel_encode(io, img; kwargs...)
    end
end

function fileio_save(io::Stream{format"SIXEL"}, img::AbstractArray; kwargs...)
    sixel_encode(io, img; kwargs...)
end

sixel_encode(io::Stream{format"SIXEL"}, img::AbstractArray, enc::SEC=default_encoder(img); kwargs...) =
    sixel_encode(io.io, img, enc; kwargs...)


text:
module model1d1d
using Base: AbstractFloat
using Dates
using QuadGK
using LinearAlgebra
using ProgressMeter
using Statistics
using LaTeXStrings
# suppressses the GKS QtTerm Windows from appreaing when using plots for savefig
ENV["GKSwstype"] = "nul"
using Plots
using Plots.PlotMeasures
gr()
using Printf
using Colors
using StaticArrays

const Tmax = 2000.0
const Twarn = 1900.0
const Trange = LinRange(0.0,Tmax+1, 10_000)
const sigma = 5.670374419E-8
const dt_max = 0.1

# use for logging inside module files
# initialize using new_log to custom logfile location (name)
m1d1d_logfile = "m1d1dlogfile.log"

include("print_log.jl")
include("front_radiation_types.jl")
include("matprop_types.jl")
include("matprop_functions.jl")
include("heat_transfer_types.jl")
include("heat_transfer_functions.jl")
include("model1d1d_types.jl")
include("mesh_types.jl")
include("varcalc.jl")
include("faster_vec_ops.jl")
# solver function using self written explicit euler methods
# may be instable (should be used for testing only)
# include("solver_1d1d_ee.jl")
# Main solver function using DifferentialEquations.jl
include("solver_1d1d_diffeq.jl")
include("blocking_model.jl")
include("post_process.jl")

export print_log, new_log, test_global_logfilename
export solve_1d1d_diffeq
export SimSettings

end


text:
# Short

"f_0/0"
f_0() = ()

"f_1/1"
f_1(x) = ()

"""f_2/1"""
f_2(x) = ()

# Long

"f_3/1"
function f_3(x)
    ()
end

"""f_4/1"""
function f_4(x)
    ()
end

# Complex type signatures.

"f_5"
f_5{T <: Integer}(a::T, b::Vector{T}, c::Complex{T}) = ()

"f_6"
f_6(::Type{String}, a::Array{Int8,1}) = ()

"f_7"
f_7(n::@compat(Union{Int64,Int64})) = ()

"f_8"
f_8{T<:Integer}(n::T, k::T) = ()

"f_9"
f_9{T<:Real}(::Complex{T}, x::Real) = ()

"f_10"
f_10{N}(dest::Array, src::Array, I::NTuple{N,@compat(Union{Int,AbstractVector})}...) = ()

"f_11"
f_11(V, indexes::NTuple, dims::Dims, linindex::UnitRange{Int}) = ()

"f_12"
f_12{T}(::Type{typejoin()}, ::Type{T}) = ()

"f_13"
f_13{T<:@compat(Union{Float32,Float64})}(::Type{T},i::Integer) = ()

"f_14"
f_14{T<:AbstractFloat}(f, a::T,b::T,c::T...) = ()

"f_15"
f_15(b, z, m = length(b) + 1) = ()

"f_16"
f_16{T<:Integer}(::Type{T}, x::AbstractFloat) = ()

"f_17"
f_17(io::IO, itr::AbstractArray, op, delim, cl, delim_one, compact = false) = ()

"f_18"
f_18(io::IO, items, sep, indent::Int, prec::Int = 0, enclose_operators::Bool = false) = ()

"f_19"
f_19{T,N,P,IV}(V::SubArray{T,N,P,IV}, m::Real) = ()

"f_20"
f_20{T,N,P,IV}(V::SubArray{T,N,P,IV}, I::AbstractArray{Bool,N}) = ()

"f_21"
f_21{T,N,P,IV}(V::SubArray{T,N,P,IV}, v, I::@compat(Union{Real,AbstractVector})...) = ()

"f_22"
f_22(::Type{Ptr{Int8}}, s::String) = ()

"f_23"
f_23(T, matvecA::Function, matvecB::Function, solveSI::Function, n::Integer,
     sym::Bool, cmplx::Bool, bmat::String,
     nev::Integer, ncv::Integer, which::String,
     tol::Real, maxiter::Integer, mode::Integer, v0::Vector) = ()

"f_24"
f_24{T}(::Matrix{T}, A::Bidiagonal{T}) = ()

"f_25"
f_25(f::ANY, t::Array, i, lim::Integer, matching::Array{Any,1}) = ()

"f_26"
f_26(f::Function, types::tup(Vararg{Type})) = ()

"f_27"
f_27(f::Base.Callable, ::tup(), ts::Tuple...) = ()

"f_28"
f_28(x::tup(Any, Vararg{Any})) = ()

"f_29"
f_29(f::Base.Callable, t::tup(Any,Any), s::tup(Any,Any)) = ()

"f_30"
f_30{T, S}(::Type{Matrix{T}}, ::Type{Vector{S}}) = ()

"f_31"
f_31{T <: Integer, S <: Real, N}(::Type{Array{T, N}}, ::S, other = nothing) = ()

"f_32"
f_32{T, S}(::Type{Matrix{T}}, ::Type{Vector{S}}) = ()

"f_33"
f_33{T <: Integer, S <: Real}(::Type{Matrix{T}}, ::Type{Vector{S}}) = ()

# Generic functions.

"g_f_1"
f_1

"""g_f_2"""
f_2

# Generic pre-docs.

"f_34";
:f_34

f_34() = ()

"""f_35""";
:f_35

f_35() = ()

# Grouped docstrings.

"f_36";
(:f_36, Vararg)

f_36() = ()
f_36(x) = ()
f_36(x, y) = ()
f_36(x, y, z) = ()

# [] group syntax.

"f_37";
[:f_37, :f_37!]

f_37() = ()
f_37!() = ()

"f_38/f_39";
[(:f_38, Any), (:f_39, Any)]

f_38(x) = ()
f_39(x) = ()

f_40(x) = ()
f_40!(x) = ()

"f_40";
[f_40, f_40!]

f_41()  = ()
f_41(x) = ()
f_42()  = ()
f_42(x) = ()

"f_41/f_42";
[(f_41,),
 (f_41, Any),
 (f_42,),
 (f_42, Any)]

# Multiline Signatures.

"f_43"
function f_43(
        a,
        b,
        c
        )
    ()
end

f_44() = 44
f_45() = 45

"f_44, f_45";
f_44, f_45


text:
using   DataFrames,
        StatsBase,
        DataFramesMeta,
        Plots,
        StatsPlots

"""
    Nucleotide_Essentials.PlotQuality

Returns a plot of the quality profile of a .fastq or .fastq.gz file

This function plots a visual summary of the distribution of quality scores (automatically detects Phred+33 or Phred+64 encoding) as a function of sequence position for the input fastq file(s).

The plotted lines show summary statistics at each sequence position:

* green is the mean
* dashed red lines are the 25th and 75th quantiles

Supported keyword arguments include:

* ```Input::FastqRecord```: The name of a FastqRecord for plotting
* ```verbose::Bool``` (optional): Whether or not to show some intermediary feedback on the progress of the function (default = false)
* ```outputfigure::Bool``` (optional): Whether or not to output a .png file with the created QualityPlot (default = false)
* ```figurepath::String``` (optional): If outputting a .png figure to file, specify the path to a directory where the file should be written to (default = ```pwd()```)

# Example:
```julia
# A quality profile can be created by supply the path to a .fastq or .fastq.gz file
PlotQuality("path/to/my/file.fastq")
```
"""
function PlotQuality(Input::String, verbose::Bool = false, outputfigure::Bool = false, figurepath::String = pwd())
    if verbose == true
        println(string("Reading file: ", Input))
    end

    if endswith(Input, ".gz") === true
        if verbose == true
            println(string("""|
                        |----Unzipping """, Input))
        end
        run(`gunzip $Input`); #this works
        unzipped = true
    else
        unzipped = false
    end

    if endswith(Input, ".gz") === true
        R1_fixed = replace(Input, ".fastq.gz"=>".fastq")
        R1 = readFastq(R1_fixed)
    else
        R1 = readFastq(Input)
    end
    quals = R1.quality

    if occursin("^", join(quals)) == true && occursin("a", join(quals)) == true && occursin("]", join(quals)) == true && occursin("f", join(quals)) == true
        encoding = 64;
        if verbose == true
            println("""|
            |----Detected Phred+64 quality encoding""")
        end
    end
    if occursin("^", join(quals)) == false && occursin("a", join(quals)) == false && occursin("]", join(quals)) == false && occursin("f", join(quals)) == false
        if verbose == true
            println("""|
            |----Detected Phred+33 quality encoding""")
        end
        encoding = 33;
    end

    if encoding == 33
        scores_33 = Dict(
            ('!') => "0",
            ('\"') => "1",
            ('#') => "2",
            ('\$') => "3",
            ('%') => "4",
            ('&') => "5",
            ('\'') => "6",
            ('(') => "7",
            (')') => "8",
            ('*') => "9",
            ('+') => "10",
            (',') => "11",
            ('-') => "12",
            ('.') => "13",
            ('/') => "14",
            ('0') => "15",
            ('1') => "16",
            ('2') => "17",
            ('3') => "18",
            ('4') => "19",
            ('5') => "20",
            ('6') => "21",
            ('7') => "22",
            ('8') => "23",
            ('9') => "24",
            (':') => "25",
            (';') => "26",
            ('<') => "27",
            ('=') => "28",
            ('>') => "29",
            ('?') => "30",
            ('@') => "31",
            ('A') => "32",
            ('B') => "33",
            ('C') => "34",
            ('D') => "35",
            ('E') => "36",
            ('F') => "37",
            ('G') => "38",
            ('H') => "39",
            ('I') => "40",
            ('J') => "41",
            ('K') => "42"
        )
    end
    if encoding == 64
        scores_64 = Dict(
            ('@') => "0",
            ('A') => "1",
            ('B') => "2",
            ('C') => "3",
            ('D') => "4",
            ('E') => "5",
            ('F') => "6",
            ('G') => "7",
            ('H') => "8",
            ('I') => "9",
            ('J') => "10",
            ('K') => "11",
            ('L') => "12",
            ('M') => "13",
            ('N') => "14",
            ('O') => "15",
            ('P') => "16",
            ('Q') => "17",
            ('R') => "18",
            ('S') => "19",
            ('T') => "20",
            ('U') => "21",
            ('V') => "22",
            ('W') => "23",
            ('X') => "24",
            ('Y') => "25",
            ('Z') => "26",
            ('[') => "27",
            ('\\') => "28",
            (']') => "29",
            ('^') => "30",
            ('_') => "31",
            ('`') => "32",
            ('a') => "33",
            ('b') => "34",
            ('c') => "35",
            ('d') => "36",
            ('e') => "37",
            ('f') => "38",
            ('g') => "39",
            ('h') => "40",
            ('i') => "41",
            ('j') => "42"
        )
    end

    function convert_phred(c)
        if encoding == 33
            scores_33[c]
        else
            scores_64[c]
        end
    end

    max_length = maximum(length.(quals))
    converted_quals = DataFrame([Float64[] for i in 1:max_length], :auto)
    for i in 1:length(quals)
        temp = parse.(Float64, (map(convert_phred, collect(quals[i]))))
        if length(temp) < max_length
            difference = max_length - length(temp)
            for r in 1:difference
                append!(temp, [0.00])
            end
        end
        converted_quals = push!(converted_quals, temp)
    end

    summarized = describe(converted_quals)
    summarized[!, :Cycle] = (1:nrow(summarized))
    summarized.Cycle = 1:nrow(summarized)

    output_long = stack(converted_quals)
    output_long.variable .= replace.(output_long.variable, "x" => "")
    cycles = Vector(output_long.variable)
    cycles = Array(parse.(Float64, cycles))
    output_long.Cycle = cycles
    res = @transform groupby(output_long, :variable) begin
       :quantile_1 = quantile(:value, 0.25)
       :quantile_2 = quantile(:value, 0.5)
       :quantile_3 = quantile(:value, 0.75)
    end
    res.Cycle = cycles
    annotation = nrow(converted_quals)
    # make our actual plot
    if verbose == true
        println("""|
        |----Creating Quality Plot""")
    end
    gr()
    qualityplot = begin
        @df summarized Plots.plot(
            :Cycle,
            :mean,
            legend = false,
            title = R1.filename,
            linecolor = :green,
            titlefontsize = 8,
            xlabel = "Cycle",
            ylabel = "Quality Score")
        @df res Plots.plot!(:Cycle, :quantile_1, linealpha = 0.7, linestyle = :dot, linecolor = :red)
        @df res Plots.plot!(:Cycle, :quantile_3, linealpha = 0.7, linestyle = :dot, linecolor = :red)
    ylims!((0,42))
    #vline!([median(length.(quals))]) # show the median length of the reads
    annotate!(50, 2, Plots.text(string("Reads: ", annotation), 12))
    end

    if outputfigure == true
        fig_out = string(figurepath, "/", R1.filename, "_qualityplot.png")
        if verbose == true
            println(string("""|
            |----Saving output file to """, fig_out))
        end
        savefig(qualityplot, fig_out)
    end

    if unzipped == true
        if verbose == true
            println(string("""|
            |----Cleaning things up"""))
        end
        run(`gzip $R1_fixed`);
    end
    return qualityplot
end


text:
mutable struct ImageSegmentation{N} <: DLPipelines.LearningMethod{ImageSegmentationTask}
    sz::NTuple{N}
    classes::AbstractVector
    downscale::Int
    projectivetransforms::ProjectiveTransforms
    imagepreprocessing::ImagePreprocessing
end

function ImageSegmentation(
        classes::AbstractVector,
        sz=(224, 224);
        augmentations=Identity(),
        downscale=1,
        means=IMAGENET_MEANS,
        stds=IMAGENET_STDS,
        C=RGB{N0f8},
        T=Float32)

    projectivetransforms = ProjectiveTransforms(sz, augmentations=augmentations)
    imagepreprocessing = ImagePreprocessing(means, stds; C=C, T=T)
    return ImageSegmentation(
        sz, classes, downscale,
        projectivetransforms, imagepreprocessing)
end

# ## Core interface

DLPipelines.encode(method::ImageSegmentation, context, sample::NamedTuple) =
    DLPipelines.encode(method, context, Tuple(sample))

function DLPipelines.encode(
        method::ImageSegmentation,
        context,
        sample::Tuple)
    image, mask = sample
    imagec, maskc = run(
        method.projectivetransforms,
        context,
        (DataAugmentation.Image(image), MaskMulti(mask, 1:length(method.classes))))

    x = run(method.imagepreprocessing, context, imagec)

    f = method.downscale
    if f != 1
        newsz = ntuple(i -> round(Int, size(image, i) * 1 / f), ndims(image))
        ytfm = ScaleFixed(newsz) |> DataAugmentation.Crop(newsz, DataAugmentation.FromOrigin()) |> OneHot()
    else
        ytfm = OneHot()
    end
    y = apply(ytfm, MaskMulti(maskc, 1:length(method.classes))) |> itemdata
    return (x, y)
end

function DLPipelines.decodeŷ(method::ImageSegmentation, context, ŷ)
    return onecoldmask(ŷ)
end

function DLPipelines.mocksample(method::ImageSegmentation)
    inputsz = rand.(UnitRange.(method.sz, method.sz .* 2))
    return (
        input = rand(RGB{N0f8}, inputsz),
        target = rand(1:length(method.classes), inputsz)
    )
end

function DLPipelines.mockmodel(method::ImageSegmentation)
    return function segmodel(xs)
        outsz = (
            round.(Int, size(xs)[1:end-1] ./ method.downscale)...,
            length(method.classes),
            size(xs)[end])
        return rand(Float32, outsz)
    end
end

onecoldmask(mask) = reshape(map(I -> I.I[end], argmax(mask; dims=ndims(mask))), size(mask)[1:end-1])

# ## Plotting interface

function plotsample!(f, method::ImageSegmentation, sample)
    image, mask = sample
    f[1, 1] = ax1 = imageaxis(f)
    f[1, 2] = ax2 = imageaxis(f)
    plotimage!(ax1, image)
    plotmask!(ax2, mask, method.classes, )
end

function plotxy!(f, method::ImageSegmentation, (x, y))
    image = invert(method.imagepreprocessing, x)
    mask = decodeŷ(method, Inference(), y)
    f[1, 1] = ax1 = imageaxis(f)
    f[1, 2] = ax2 = imageaxis(f)
    plotimage!(ax1, image)
    plotmask!(ax2, mask, method.classes, )
end

function DLPipelines.methodmodel(method, backbone)
    h, w, ch, b = Flux.outdims(backbone, (method.sz..., 3, 1))
    headupscale = (method.sz[1] / method.downscale) / h
    n_upscale = ceil(Int, log2(headupscale))
    return Chain(
        backbone,
        Models.pixelshufflehead(ch, length(method.classes), n_upscale = n_upscale)
    )

end


text:
module Interfaces

import LinearAlgebra: rank, mul!, ⋅

export id, value, rank, dimension
export ⊕, ⊗, ⋅, add!, sub!, mul!, div!
export expand, expand!, decompose, decompose!, permute

"""
ID.
"""
function id end

"""
Value.
"""
function value end

"Rank."
function rank end

"Dimension."
function dimension end

"Direct sum."
function ⊕ end

"Direct product."
function ⊗ end

"Dot product."
function ⋅ end

"In-place addition."
function add! end

"In-place subtraction."
function sub! end

"In-place multiplication."
function mul! end

"In-place division."
function div! end

"Get the expansion."
function expand end

"In place expansion."
function expand! end

"Decompose."
function decompose end

"In place decomposition."
function decompose! end

"Get the permutation."
function permute end

end # module


text:
using Plasmo
using KaHyPar
using Ipopt
using SchwarzSolver

n_nodes = 100  #number of intervals
N = 3000
k = round(Int64,N/n_nodes)
d = sin.(1:N)

disturbance_matrix = reshape(d, (k, div(length(d), k)))

graph = ModelGraph()
@node(graph,nodes[1:n_nodes])
for (i,node) in enumerate(nodes)
    node_disturbance = disturbance_matrix[:,i]
    M = length(node_disturbance)
    @variable(node, x[1:M+1])
    @variable(node, u[1:M])
    @constraint(node, dynamics[i in 1:M], x[i+1] == x[i] + u[i] + d[i])
    @objective(node, Min, 0.001*sum(x[i]^2 - 2*x[i]*d[i] for i in 1:M) + sum(u[i]^2 for i in 1:M))
end

n1 = getnode(graph,1)
@constraint(n1,n1[:x][1] == 0)                      #First node in planning horizon has initial condition

links = []
for i = 1:n_nodes - 1
    ni = getnode(graph,i)
    nj = getnode(graph,i+1)
    lref = @linkconstraint(graph, ni[:x][end] == nj[:x][1])  #last state in partition i is first state in partition
    push!(links,lref)
end

#Partition the problem
hypergraph,hyper_map = gethypergraph(graph) #create hypergraph object based on graph
partition_vector = KaHyPar.partition(hypergraph,8,configuration = :connectivity,imbalance = imbalance)
partition = Partition(hypergraph,partition_vector,hyper_map)
make_subgraphs!(graph,partition)

overlap = 1
# schwarz_solve(graph,overlap;sub_optimizer = with_optimizer(Ipopt.Optimizer,tol = 1e-12,print_level = 0),max_iterations = 100,tolerance = 1e-10)

schwarz_solve(graph,overlap;sub_optimizer = with_optimizer(Ipopt.Optimizer,tol = 1e-12,print_level = 0),max_iterations = 100,tolerance = 1e-10,
dual_links = [],primal_links = [])

# schwarz_solve(graph,overlap;sub_optimizer = with_optimizer(Gurobi.Optimizer,OutputFlag=0,BarConvTol=1e-10,FeasibilityTol=1e-9,OptimalityTol= 1e-9),max_iterations = 100,tolerance = 1e-10)


text:
using Revise
using BNPRegressionGGA2021
using CSV
using DataStructures
using DataFrames
using Distributions
using Gadfly
using Random
using Statistics
const BNP = BNPRegressionGGA2021

# Define the DGP
function simulate_sample(N0, N1)
    X0 = [ones(N0) rand([0, 1], N0, 5)]
    # X1 = kron([ones(2, 5) [1, 0]], ones(N1, 1))

    X1 = rand(2, 5), rand(2, 4)
    α = 4
    δ = α / √(1 + α^2)
    ω = 1
    ξ = - ω * δ * √(2 / π)
    ϵ = √(1 - δ^2) .* randn(N0) + δ .* abs.(randn(N0)) .+ ξ
    y0 = 0.8 * X0[:, end] + ϵ
    y1 = LinRange(-3, 3, N1) |> collect |> x -> repeat(x, size(X1, 1) ÷ N1)
    return y0, X0, y1, X1
end

N0, N1, Niter = 100, 2, 1;
df = zeros(Bool, Niter, 6);
Random.seed!(1);
y0, X0, y1, X1 = simulate_sample(N0, N1);
plot(x = y0, Geom.density())
m = BNP.DGSBPNormal(; y0, X0, y1, X1);
_, _, chaing = BNP.sample!(m; mcmcsize = 10000, burnin = 5000);
mean(chaing)

# Get the MAP of γ for 100 simulated samples
begin
    N0, N1, Niter = 500, 2, 100;
    df = zeros(Bool, Niter, 6);
    Random.seed!(1);
    for iter in 1:Niter
        y0, X0, y1, X1 = simulate_sample(N0, N1);
        m = BNP.DGSBPNormal(; y0, X0, y1, X1);
        _, _, chaing = BNP.sample!(m; mcmcsize = 10000, burnin = 5000);
        df[iter, :] .= begin
            chaing |>
            DataStructures.counter |>
            collect |>
            x -> sort(x, by = x -> x[2], rev = true) |>
            x -> getindex(x, 1) |>
            x -> getindex(x, 1)
        end
        println(iter)
        println(df[iter, :])
    end
    CSV.write("data/simulation-normal-skewnormal-gamma-bnp.csv", DataFrame(df, :auto))
end

# Save the 100 simulated samples
Random.seed!(1);
df = map(1:100) do iter
    y0, X0, _, _ = simulate_sample(N0, N1);
    df = DataFrame(X0, :auto)
    df[!, :y0] = y0
    df
end;
df = reduce(vcat, df);
CSV.write("data/simulation-normal-skewnormal-data.csv", df);


text:
include("./../para.jl")

struct Spins
	N::Int64
	M::Int64
	A::Array
end

energy(S::Spins) = sum(
+(S.A[2:para.N+2,2:para.M+2] * S.A[1:para.N+1,2:para.M+2]
, S.A[2:para.N+2,2:para.M+2] * S.A[3:para.N+3,2:para.M+2]
, S.A[2:para.N+2,2:para.M+2] * S.A[2:para.N+2,1:para.M+1]
, S.A[2:para.N+2,2:para.M+2] * S.A[2:para.N+2,3:para.M+3]))

max_stars_repo_name: anirbrhm/curious-eavesdropper

max_issues_repo_name: anirbrhm/curious-eavesdropper

max_forks_repo_name: anirbrhm/curious-eavesdropper


text:
using AWS
using DynamoDB
using Base.Test

type Foo
    a
    b
end

# tests on component functions
include("dynamo_json.jl")
include("dynamo_dsl.jl")
include("dynamo.jl")
include("dynamo_row_ops.jl")

# live tests
include("integration_tests.jl")


text:
# 1-D interpolation from irregularly sampled data.

import PyPlot
import Random
import FIRRemez

fig_num = 1
PyPlot.close("all")

Random.seed!(25)

# generate sampling locations.
N_samples = 25
X = convert( Vector{BigFloat}, randn(N_samples).*π )

# generate samples.
Y = sin.(X)

# query.
N_queries = 500
xq_BF = convert(Vector{BigFloat}, collect( LinRange(minimum(X), maximum(X), N_queries) ) )
yq_BF = collect( FIRRemez.Barycentric2nditp(xq_BF[i],X,Y) for i = 1:length(xq_BF) )

xq = convert(Vector{Float64}, xq_BF)
yq = convert(Vector{Float64}, yq_BF)

X_display = convert(Vector{Float64}, X)
Y_display = convert(Vector{Float64}, Y)

PyPlot.figure(fig_num)
fig_num += 1
PyPlot.plot(xq, yq, label = "Barycentric")
PyPlot.plot(X_display, Y_display, label = "samples", "x")
PyPlot.title("Barycentric interpolation")
PyPlot.legend()

max_stars_repo_name: RoyCCWang/FIRRemez

max_issues_repo_name: RoyCCWang/FIRRemez

max_forks_repo_name: RoyCCWang/FIRRemez


text:
export GaussianFilter, GaborWavelet, MorletWavelet

abstract type AbstractFilter end

"""
    gabor1d(ω, ξ, σ)

Calculate the Fourier transform of a gabor wavelet.

# Arguments
 - `ω::Vector{Float64}`: Frequencies over which to calculate the filter.
 - `σ::Float64`: Bandwidth of the filter
 - `ξ::Float64`: Central frequency
"""
function gabor1d(ω::Vector{Float64}, σ::Float64, ξ::Float64)
    return exp.(-(ω.-ξ).^2 ./ (2*σ^2))
end
"""
    gabor1d(N, ξ, σ)

Convenience method for calculating the Fourier transform of a gabor wavelet.

# Arguments
 - `ω::Vector{Float64}`: Frequencies over which to calculate the filter.
 - `σ::Float64`: Bandwidth of the filter
 - `ξ::Float64`: Central frequency
"""
function gabor1d(N::Int64, σ::Float64, ξ::Float64)
    ω = FFTW.fftfreq(N)
    return gabor1d(ω[:], σ, ξ)
end

"""
    gauss1d(ω, σ)

Return a Fourier transformed Gaussian lowpass filter.

# Arguments
 - `ω::Vector{Float64}`: Frequencies over which to calculate the filter.
 - `σ::Float64`: Bandwidth of lowpass filter.
"""
function gauss1d(ω::Vector{Float64}, σ::Float64)
    return exp.(-(ω.^2)./(2*σ^2))
end
"""
    gauss1d(N, σ)

Convenience method for calculate the Fourier transform of a Gaussian lowpass filter.

# Argument
 - `N::Int64`: Length of filter in the time domain.
 - `σ::Float64`: Bandwidth of the lowpass filter.
"""
function gauss1d(N::Int64, σ::Float64)
    ω = collect(0.0:N-1)./N
    return gauss1d(ω[:], σ)
end

"""
    morlet1d(gab, lowpass)

Calculate the Fourier transform of a Morlet wavelet.

# Arguments
 - `gab::Vector{Float64}`: Fourier transformed Gabor wavelet
 - `lowpass::Vector{Float64}`: Fourier transformed Gaussian lowpass filter
"""
function morlet1d(gab::Vector{Float64}, lowpass::Vector{Float64})
    # summation factor to ensure first value is 0
    κ = gab[1]/lowpass[1]
    return gab .- κ.*lowpass
end
"""
    morlet1d(N, ξ, σ)

Calculate the Fourier transform of a Morlet wavelet.

# Arguments
 - `N::Int64`: Length of filter in the time domain.
 - `σ::Float64`: Bandwidth of the lowpass filter.
 - `ξ::Float64`: Central frequency.
"""
function morlet1d(N::Int64, σ::Float64, ξ::Float64)
    gab = gabor1d(N, σ, ξ)
    lowpass = gauss1d(N, σ)
    return morlet1d(gab, lowpass)
end

"""
    GaussianFilter

A structure for Fourier transformed Gaussian lowpass Filters

## Fields: GaussianFilter
 | **Field** | **Description** |
 |:----------|:----------------|
 | :σ        | Filter Bandwidth |
 | :ϕ        | Fourier-transformed Gaussian lowpass filter |
"""
struct GaussianFilter <: AbstractFilter
    σ::Float64 # filter bandwidth

    ϕ::Vector{Float64} # Fourier-transformed lowpass Gaussian filter

    function GaussianFilter(N::Int64, σ::Float64)
        ϕ = gauss1d(N, σ)
        return new(σ, ϕ)
    end
end

"""
    GaborWavelet

A structure for Fourier transformed Gabor Wavelets

## Fields: GaborWavelet
 | **Field** | **Description** |
 |:----------|:----------------|
 | :σ        | Filter Bandwidth |
 | :ξ        | Central Frequency |
 | :j        | Maximal subsampling |
 | :γ        | Fourier-transformed Gabor wavelet |
"""
struct GaborWavelet <: AbstractFilter
    σ::Float64 # filter bandwidth
    ξ::Float64 # central frequency
    j::Int64 # 2^j is the maximal subsampling possible without aliasing

    γ::Vector{Float64} # Fourier-transformed Gabor wavelet

    function GaborWavelet(N::Int64, σ::Float64, ξ::Float64, j::Int64)
        γ = gabor1d(N, σ, ξ)
        return new(σ, ξ, j, γ)
    end
end

"""
    MorletWavelet

A structure for Fourier transformed Morlet Wavelets

## Fields: MorletWavelet
 | **Field** | **Description** |
 |:----------|:----------------|
 | :σ        | Filter Bandwidth |
 | :ξ        | Central Frequency |
 | :j        | Maximal subsampling |
 | :ψ        | Fourier-transformed Morlet wavelet |
"""
struct MorletWavelet <: AbstractFilter
    σ::Float64 # filter bandwidth
    ξ::Float64 # central frequency
    j::Int64 # 2^j is the maximal subsampling possible without aliasing

    ψ::Vector{Float64} # Fourier-transformed Morlet wavelet

    function MorletWavelet(N::Int64, σ::Float64, ξ::Float64, j::Int64)
        ψ = morlet1d(N, σ, ξ)
        return new(σ, ξ, j, ψ)
    end
    function MorletWavelet(γ::GaborWavelet, ϕ::GaussianFilter)
        ψ = morlet1d(γ.γ, ϕ.ϕ)
        return new(γ.σ, γ.ξ, γ.j, ψ)
    end
end


text:
module TestPI

using Test, PIMD

@testset "PI" begin
    β, m, V = 1., 1., harmonicPotential
    system = OneDimensionalQuantumSystem(β, m, V)
    @test system.potential(10.) ≈ 50.
    @test system.force(10.) ≈ -10.
    @test system.hessian(10.) ≈ 1.
end

end


text:
using HorizonSideRobots
include("BackPath.jl")

function left(side::HorizonSide)
    return(HorizonSide((Int(side) + 1) % 4))
end

function try_move!(r::Robot,side::HorizonSide)::Bool
    ort_side = left(side)
    n=0
    while isborder(r,side)
        if !isborder(r,ort_side)
            move!(r,ort_side)
            n += 1
        else
            break
        end
    end
    if isborder(r,side)
        movements!(r,inverse(ort_side),n)
        return false
    end
    move!(r,side)
    while (isborder(r,inverse(ort_side)) && n != 0)
        move!(r,side)
    end
    movements!(r,inverse(ort_side),n)
    return true
end

max_stars_repo_name: Mananhina/mirea-progs

max_issues_repo_name: Mananhina/mirea-progs

max_forks_repo_name: Mananhina/mirea-progs


text:
using DataFrameMacros
using DataFrames
using Test

@testset "transform & select" begin
    df = DataFrame(a = 1:100, b = rand('a':'e', 100), c = randn(100))

    df2 = @transform(df, :a + 1)
    @test df2 == transform(df, :a => ByRow(x -> x + 1))
    df3 = @transform(df, :d = :a + 1)
    @test df3 == transform(df, :a => ByRow(x -> x + 1) => :d)
    df4 = @transform(df, :d = :a + 1, :e = :c * 2)
    @test df4 == transform(df,
        :a => ByRow(x -> x + 1) => :d,
        :c => ByRow(x -> x * 2) => :e
    )
    df5 = @transform(df, :d = :a + 1, :e = :c * 2, :f = :a * :c - (:c / :a))
    @test df5 == transform(df,
        :a => ByRow(x -> x + 1) => :d,
        :c => ByRow(x -> x * 2) => :e,
        [:a, :c] => ByRow((x, y) -> x * y - (y / x)) => :f
    )

    col1 = :a
    col2 = :c
    df6 = @transform(df, :d = $col1 + $col2 * :a)
    @test df6 == transform(df, [col1, col2, :a] => ByRow((x, y, z) -> x + y * z) => :d)

    s = "d"
    df7 = @transform(df, Symbol(s) = :a + 1)
    @test df7 == transform(df, :a => ByRow(x -> x + 1) => :d)

    df8 = @transform(df, AsTable = (x = "$(:a) $(:b)", y = :c * 2))

    @test df8 == transform(df,
        [:a, :b, :c] =>
            ByRow((a, b, c) -> (x = "$a $b", y = c * 2)) => AsTable
    )

    df9 = @transform(df, :a)
    @test df9 == transform(df, :a)

    df10 = @select(df, :a * 2)
    @test df10 == select(df, :a => ByRow(x -> x * 2))

    df11 = @select(df, :x = :a * 2)
    @test df11 == select(df, :a => ByRow(x -> x * 2) => :x)

    myfunc(x) = x + 2
    df12 = @transform(df, myfunc(:a))
    @test "a_myfunc" in names(df12)
    @test df12 == transform(df, :a => ByRow(myfunc))

    df13 = @transform(df, :x = @c sum(:a) .+ :c)
    @test df13 == transform(df, [:a, :c] => ((x, y) -> sum(x) .+ y) => :x)

    df14 = DataFrame(a = ['a', 'b', missing])
    df15 = @transform(df14, @m uppercase(:a))
    # @test df15 == transform(df14, :a => ByRow(passmissing(uppercase)))
end

@testset "combine" begin
    df = DataFrame(a = 1:100, b = rand('a':'e', 100), c = randn(100))

    df1 = @combine(df, sum(:a))
    @test df1 == combine(df, :a => sum)

    df2 = @combine(df, sum(:a) + 3)
    @test df2 == combine(df, :a => (x -> sum(x) + 3))

    df3 = @combine(df, :x = sum(:a) + 3)
    @test df3 == combine(df, :a => (x -> sum(x) + 3) => :x)

    df4 = @combine(df, :y = @r 3 * :a)
    @test df4 == combine(df, :a => ByRow(x -> 3 * x) => :y)
end

@testset "subset" begin
    df = DataFrame(a = [-1, 1, -2, 2, -3, 3], b = randn(6))
    df2 = @subset(df, abs(:a) < 2)
    @test df2 == subset(df, :a => ByRow(x -> abs(x) < 2))
end

@testset "groupby" begin
    df = DataFrame(a = [-1, 1, -2, 2, -3, 3], b = randn(6))
    gdf = @groupby(df, abs(:a))

    gdf2 = @groupby(df, :a = abs(:a))

    gdf3 = @groupby(df, :x = abs(:a))
end

@testset "sort" begin
    df = DataFrame(a = rand(1:10, 100), b = rand('a':'z', 100))
    df2 = @sort(df, :a)
    @test df2 == sort(df, :a)

    df3 = @sort(df, -:a)
    @test df3 == sort(df, :a, rev = true)

    df4 = @sort(df, -:a, :b)
    @test df4 == sort(df, [order(:a, rev = true), :b])

    df5 = @sort(df, :a; rev = true)
    @test df5 == sort(df, :a, rev = true)
end

@testset "escaping symbols" begin
    df = DataFrame(color = [:red, :green, :blue])
    df2 = @transform(df, :color == $:red)
    @test df2[:, 2] == transform(df, :color => ByRow(==(:red)))[:, 2]
end

@testset "mutating" begin
    df = DataFrame(a = 1:100, b = rand('a':'e', 100), c = randn(100))
    @transform!(df, :a = :a * 2)
    @sort!(df, :b)
    @subset!(df, :a > 20)
    @select!(df, :a, :c)
end

@testset "string and int column specification" begin
    df = DataFrame([Symbol("a column") => [1, 2, 3], :b => [4, 5, 6]])
    df2 = @transform(df, :c = $"a column" * $2)
    #"
    @test df2 == transform(df, ["a column", names(df)[2]] => ByRow(*) => :c)
end

module Mod1
    module Mod2
        export func2
        func2(x) = 3x
    end
    using .Mod2
    export Mod2

    export func
    func(x) = 2x
end

@testset "modules" begin

    using .Mod1

    df = DataFrame(a = [1, 2, 3])
    df2 = @transform(df, :b = Mod1.func(:a), :c = Mod1.Mod2.func2(:a))
    @test df2 == transform(df,
        :a => ByRow(x -> Mod1.func(x)) => :b,
        :a => ByRow(x -> Mod1.Mod2.func2(x)) => :c,
    )
end

@testset "simple formula" begin
    df = DataFrame(a = 1:10)
    df2 = @transform(df, :b = :a + :a)
    @test df2 == transform(df, :a => ByRow(x -> x + x) => :b)

    using .Mod1

    # check that column names are correct for simple function calls
    df3 = @transform(df, Mod1.func(:a), Mod1.Mod2.func2(:a))
    @test df3 == transform(df,
        :a => ByRow(Mod1.func),
        :a => ByRow(Mod1.Mod2.func2),
    )
end

@testset "getproperty dot syntax" begin
    df = DataFrame(a = [(;b = 1), (;b = 3)])
    df2 = @transform(df, :c = :a.b)
    @test df2 == transform(df, :a => ByRow(x -> x.b) => :c)
end

@testset "block syntax" begin
    df = DataFrame(a = 1:3, b = 'a':'c')
    df2 = @transform df begin
        :c = :a + 1
        :d = string(:b) ^ 2
    end
    @test df2 == @transform(df, :c = :a + 1, :d = string(:b) ^ 2)
end

@testset "table shortcut @t" begin
    df = DataFrame(a = 1:3, b = 4:6)
    df2 = @transform(df, @t begin
        x = :a + :b
        :y = x * 2
        :z = x + 4
    end)
    @test df2 == transform(df, [:a, :b] => ByRow((a, b) -> begin
        x = a + b
        vary = x * 2
        varz = x + 4
        (y = vary, z = varz)
    end) => AsTable)

    df3 = @transform(df, @t begin
        if iseven(:a)
            :z = :a
        else
            :z = :b
        end
    end)

    # tuple destructuring
    df4 = @transform(df, @t begin
        x = :a + :b
        :y, :z = x * 2, x + 4
    end)
    @test df4 == df2

    # tuple destructuring with non symbols mixed in
    df5 = @transform(df, @t begin
        x = :a + :b
        :y, qqq, :z = x * 2, "hello", x + 4
    end)
    @test df5 == df2 == df4
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 ImageDataDisk <: SwaggerModel
    lun::Any # spec type: Union{ Nothing, Int32 } # spec name: lun
    snapshot::Any # spec type: Union{ Nothing, SubResource } # spec name: snapshot
    managedDisk::Any # spec type: Union{ Nothing, SubResource } # spec name: managedDisk
    blobUri::Any # spec type: Union{ Nothing, String } # spec name: blobUri
    caching::Any # spec type: Union{ Nothing, String } # spec name: caching
    diskSizeGB::Any # spec type: Union{ Nothing, Int32 } # spec name: diskSizeGB
    storageAccountType::Any # spec type: Union{ Nothing, StorageAccountType } # spec name: storageAccountType

    function ImageDataDisk(;lun=nothing, snapshot=nothing, managedDisk=nothing, blobUri=nothing, caching=nothing, diskSizeGB=nothing, storageAccountType=nothing)
        o = new()
        validate_property(ImageDataDisk, Symbol("lun"), lun)
        setfield!(o, Symbol("lun"), lun)
        validate_property(ImageDataDisk, Symbol("snapshot"), snapshot)
        setfield!(o, Symbol("snapshot"), snapshot)
        validate_property(ImageDataDisk, Symbol("managedDisk"), managedDisk)
        setfield!(o, Symbol("managedDisk"), managedDisk)
        validate_property(ImageDataDisk, Symbol("blobUri"), blobUri)
        setfield!(o, Symbol("blobUri"), blobUri)
        validate_property(ImageDataDisk, Symbol("caching"), caching)
        setfield!(o, Symbol("caching"), caching)
        validate_property(ImageDataDisk, Symbol("diskSizeGB"), diskSizeGB)
        setfield!(o, Symbol("diskSizeGB"), diskSizeGB)
        validate_property(ImageDataDisk, Symbol("storageAccountType"), storageAccountType)
        setfield!(o, Symbol("storageAccountType"), storageAccountType)
        o
    end
end # type ImageDataDisk

const _property_map_ImageDataDisk = Dict{Symbol,Symbol}(Symbol("lun")=>Symbol("lun"), Symbol("snapshot")=>Symbol("snapshot"), Symbol("managedDisk")=>Symbol("managedDisk"), Symbol("blobUri")=>Symbol("blobUri"), Symbol("caching")=>Symbol("caching"), Symbol("diskSizeGB")=>Symbol("diskSizeGB"), Symbol("storageAccountType")=>Symbol("storageAccountType"))
const _property_types_ImageDataDisk = Dict{Symbol,String}(Symbol("lun")=>"Int32", Symbol("snapshot")=>"SubResource", Symbol("managedDisk")=>"SubResource", Symbol("blobUri")=>"String", Symbol("caching")=>"String", Symbol("diskSizeGB")=>"Int32", Symbol("storageAccountType")=>"StorageAccountType")
Base.propertynames(::Type{ ImageDataDisk }) = collect(keys(_property_map_ImageDataDisk))
Swagger.property_type(::Type{ ImageDataDisk }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ImageDataDisk[name]))}
Swagger.field_name(::Type{ ImageDataDisk }, property_name::Symbol) =  _property_map_ImageDataDisk[property_name]

const _allowed_ImageDataDisk_caching = ["None", "ReadOnly", "ReadWrite"]

function check_required(o::ImageDataDisk)
    (getproperty(o, Symbol("lun")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ImageDataDisk }, name::Symbol, val)
    if name === Symbol("caching")
        Swagger.validate_param(name, "ImageDataDisk", :enum, val, _allowed_ImageDataDisk_caching)
    end
end


text:
# Incidence matrix for element connections in the grid
function create_incidence_matrix(g::Grid)
    cell_containing_node = Dict{Int, Set{Int}}()
    for (cellid, cell) in enumerate(g.cells)
        for v in cell.nodes
            if !haskey(cell_containing_node, v)
                cell_containing_node[v] = Set{Int}()
            end
            push!(cell_containing_node[v], cellid)
        end
    end

    I, J, V = Int[], Int[], Bool[]
    for (_, cells) in cell_containing_node
        for cell1 in cells # All these cells have a neighboring node
            for cell2 in cells
                # if true # cell1 != cell2
                if cell1 != cell2
                    push!(I, cell1)
                    push!(J, cell2)
                    push!(V, true)
                end
            end
        end
    end

    incidence_matrix = sparse(I, J, V)
    return incidence_matrix
end

# Greedy algorithm for coloring a grid such that no two cells with the same node
# have the same color
function greedy_coloring(incidence_matrix, cells=1:size(incidence_matrix, 1))
    cell_colors = Dict{Int,Int}(i => 0 for i in cells) # Zero represents no color set yet
    occupied_colors = Set{Int}()
    final_colors = Vector{Int}[]
    total_colors = 0
    for cellid in cells
        empty!(occupied_colors)
        # loop over neighbors
        for r in nzrange(incidence_matrix, cellid)
            cell_neighbour = incidence_matrix.rowval[r]
            cell_neighbour in cells || continue # Only care about the subset given in cells
            color = cell_colors[cell_neighbour]
            if color != 0
                push!(occupied_colors, color)
            end
        end

        # occupied colors now contains all the colors we are not allowed to use
        free_color = 0
        for attempt_color in 1:total_colors
            if attempt_color ∉ occupied_colors
                free_color = attempt_color
                break
            end
        end
        if free_color == 0 # no free color found, need to bump max colors
            total_colors += 1
            free_color = total_colors
            push!(final_colors, Int[])
        end
        @assert free_color != 0
        cell_colors[cellid] = free_color
        push!(final_colors[free_color], cellid)
    end
    return final_colors
end

# See Appendix A in https://www.math.colostate.edu/%7Ebangerth/publications/2013-pattern.pdf
function workstream_coloring(incidence_matrix)
    ###################
    # 1. Partitioning #
    ###################
    zones = Set{Int}[]
    ## Zone 1: Just the first element
    push!(zones, Set{Int}(1))
    Z = 2
    Z0 = Set{Int}() # Dummy zone
    n_visited = 1
    ## Zone N: All elements with connection to elements in Zone N-1
    while true
        s = Set{Int}()
        # Loop over all elements in previous zone and add their neigbouring elements
        # unless they are in any of the previous 2 zones.
        for c in get(zones, Z-1, Z0)
            for r in nzrange(incidence_matrix, c)
                cell_neighbour = incidence_matrix.rowval[r]
                if !(cell_neighbour in get(zones, Z-2, Z0) || cell_neighbour in get(zones, Z-1, Z0))
                    push!(s, cell_neighbour)
                end
            end
        end
        push!(zones, s)
        n_visited += length(s)
        if n_visited >= size(incidence_matrix, 1)
            break
        end
        Z += 1
    end

    ###############
    # 2. Coloring #
    ###############
    # TODO: The reference uses DSATUR algorithm instead of greedy
    # TODO: Zones can be colorized in parallel
    zone_colors = Vector{Vector{Int}}[greedy_coloring(incidence_matrix, z) for z in zones]

    ################
    # 3. Gathering #
    ################
    Nodd,  Zodd  = findmax(x -> isodd(x)  ? length(zone_colors[x]) : typemin(Int), 1:length(zone_colors))
    Neven, Zeven = findmax(x -> iseven(x) ? length(zone_colors[x]) : typemin(Int), 1:length(zone_colors))
    N = Nodd + Neven
    final_colors = append!(zone_colors[Zodd], zone_colors[Zeven]) # Reuse these for output
    color_sizes = map(length, final_colors)
    used_for_zone = Set{Int}()
    for Z in 1:length(zone_colors)
        (Z == Zodd || Z == Zeven) && continue
        zone_color_vectors = zone_colors[Z]
        odd = isodd(Z)

        empty!(used_for_zone)

        for local_color in sortperm(zone_color_vectors; by=length, rev=true)
            cond = odd ? (x -> x <= Nodd) : (x -> x > Nodd)
            _, global_color = findmin(x -> (cond(x) && x ∉ used_for_zone) ? color_sizes[x] : typemax(Int), 1:N)
            push!(used_for_zone, global_color)
            append!(final_colors[global_color], zone_color_vectors[local_color])
            map!(length, color_sizes, final_colors)
        end
    end

    # Maybe nice to sort?
    foreach(sort!, final_colors)

    return final_colors
end

@enum ColoringAlgorithm GREEDY WORKSTREAM

"""
    create_coloring(g::Grid; alg::ColoringAlgorithm)

Create a coloring of the cells in grid `g` such that no neighboring cells
have the same color.

Returns a vector of vectors with cell indexes, e.g.:

```julia
ret = [
   [1, 3, 5, 10, ...], # cells for color 1
   [2, 4, 6, 12, ...], # cells for color 2
]
```

Two different algorithms are available, specified with the `alg` keyword argument:
 - `alg = Ferrite.WORKSTREAM` (default): Three step algorithm from
   [*WorkStream*](https://www.math.colostate.edu/%7Ebangerth/publications/2013-pattern.pdf)
   , albeit with a greedy coloring in the second step.
 - `alg = Ferrite.GREEDY`: greedy algorithm that works well for structured grid such as
   e.g. grids from `generate_grid`.

The resulting colors can be visualized using [`vtk_cell_data_colors`](@ref).

!!! note "Cell to color mapping"
    In a previous version of Ferrite this function returned a dictionary mapping
    cell ID to color numbers as the first argument. If you need this mapping you
    can create it using the following construct:
    ```julia
    colors = create_coloring(...)
    cell_colormap = Dict{Int,Int}(
        cellid => color for (color, cellids) in enumerate(final_colors) for cellid in cellids
    )
    ```
"""
function create_coloring(g::Grid; alg::ColoringAlgorithm=WORKSTREAM)
    incidence_matrix = create_incidence_matrix(g)
    if alg === WORKSTREAM
        return workstream_coloring(incidence_matrix)
    elseif alg === GREEDY
        return greedy_coloring(incidence_matrix)
    else
        error("impossible")
    end
end

"""
    vtk_cell_data_colors(vtkfile, cell_colors, name="coloring")

Write cell colors (see [`create_coloring`](@ref)) to a VTK file for visualization.
"""
function vtk_cell_data_colors(vtkfile, cell_colors::AbstractVector{<:AbstractVector{<:Integer}}, name="coloring")
    color_vector = zeros(sum(length, cell_colors))
    for (i, cells_color) in enumerate(cell_colors)
        for cell in cells_color
            color_vector[cell] = i
        end
    end
    vtk_cell_data(vtkfile, color_vector, name)
end


text:
# Clean up tools!

using GPUBenchmarks
db = GPUBenchmarks.get_database()
using GPUBenchmarks: name, timestamp

timestamp = last(sort(unique(timestamp.(db))))
single_ts = GPUBenchmarks.timestamp.(db, timestamp)
GPUBenchmarks.update_database!(single_ts)


text:
if length(ARGS) < 1 error("L argument missing") end
L=parse(Int, ARGS[1])
# L = 4
N=2*L^2
nbonds = Int(3*N/2)

# Build honeycomb lattice
println("Building honeycomb lattice...")
basis = [[0.,0.],[0.,1.]] # bond length = 1
a1 = [cos(deg2rad(30)), sin(deg2rad(30)) + 1]
a2 = a1.*[-1,1]

# site positions
positions = Vector{Vector{Float64}}(undef, N)
pi = 1
for ix in 0:L-1
    for iy in 0:L-1
        global positions, pi
        positions[pi] = basis[1]+(iy)*a1+ix*(a1-a2)
        positions[pi+1] = basis[2]+(iy)*a1+ix*(a1-a2)
        pi += 2
    end
end

# unit cells neighbortable
ucs = reshape(1:L^2, (L,L)) # unitcells form square lattice
downright = circshift(ucs, (+1,-1)) # down right bond
up = circshift(ucs, (-1,0)) # up "right" bond
ucneighs = vcat(downright[:]', up[:]') # column = unitcell linidx

getsite(uc::Int, site::Int) = 2 * (uc-1) + site

# bonds
bonds = Vector{Vector{Int}}(undef, nbonds) # from, to, type
bondvecs = Vector{Vector{Float64}}(undef, nbonds) # to_pos - from_pos
ibond = 1
for uc in 1:L^2
    global bonds, bondvecs, ibond
    # inner bond
    s1 = getsite(uc, 1)
    s2 = getsite(uc, 2)
    bonds[ibond] = [s1, s2, 0]
    bondvecs[ibond] = [0.0, 1.0]
    ibond += 1

    # down-right bond
    trg = getsite(ucneighs[1,uc],2)
    bonds[ibond] = [s1, trg, 1]
    bondvecs[ibond] = [0.8660254037844387, 0.5]
    ibond += 1

    # up-"right" bond
    trg = getsite(ucneighs[2,uc],1)
    bonds[ibond] = [s2, trg, 2]
    bondvecs[ibond] = [0.8660254037844387, -0.5]
    ibond += 1
end

# Create ALPS like xml file
println("Creating ALPS like XML file...")
using LightXML

# create an empty XML document
xdoc = XMLDocument()

# create & attach a graph (root) node
graph = create_root(xdoc, "GRAPH")
nedges = Int(3*N/2)
set_attributes(graph, dimension=2, vertices=N, edges=nedges)

function addvertex!(node, id, typ, coordinates)
    v = new_child(node, "VERTEX")
    set_attributes(v, Dict("id"=>id, "type"=>typ))

    vc = new_child(v, "COORDINATE")
    add_text(vc, "$(join(coordinates, " "))")
end

function addedge!(node, src, trg, id, typ, vector)
    edge = new_child(node, "EDGE")
    set_attributes(edge, id=id)
    set_attributes(edge, source=src, target=trg)
    set_attributes(edge, Dict("type"=>typ, "vector"=>join(vector, " ")))
end

# add vertices to graph
for (vi, pos) in enumerate(positions)
    addvertex!(graph, vi, Int(iseven(vi)), pos)
end

# add edges to graph
for (ibond, b) in enumerate(bonds)
    addedge!(graph, b[1], b[2], ibond, b[3], bondvecs[ibond])
end

# s = string(xdoc);
# print(xdoc)

save_file(xdoc, "honeycomb_L_$(L)_W_$(L).xml");

# Create plot
println("Creating a plot of the lattice...")
using PyPlot

# bondcolors=["C7", "C8", "C9"]
bondcolors=fill("black", 3)
sitecolors=["white", "gray"]

# plot bonds
for (ibond, b) in enumerate(bonds)
    src = positions[b[1]]
    trg = positions[b[2]]

    # don't show pbc bonds
    b[2] > b[1] || continue
    b[2] - b[1] < 2*L || continue

    plot([src[1], trg[1]], [src[2], trg[2]], color=bondcolors[b[3]+1], linewidth=2.0)
    gca()[:set_axis_off]()
end

# plot sites
for (s, pos) in enumerate(positions)
    stype = Int(iseven(s))+1
#     plot(pos..., marker="o", markersize=13, color=sitecolors[stype])
    plot(pos..., marker="o", markersize=13, color=sitecolors[stype], mec="black", mew=1.5)
end
# scatter(getindex.(positions, 1), getindex.(positions, 2), c=1:length(positions), cmap="Blues")

axis("equal")
savefig("honeycomb_L_$(L)_W_$(L).pdf");

println("Done.")

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
using OrdinaryDiffEq, Test
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear

@testset "saveat Tests" begin
  prob = prob_ode_linear

  sol =solve(prob,DP5(),dt=1//2^(2),save_everystep=false)
  sol2=solve(prob,DP5(),dt=1//2^(2),save_everystep=false,saveat=[1/2])

  @test sort!(symdiff(sol.t,sol2.t)) == [0.0,1/2,1.0]

  sol2=solve(prob,DP5(),dt=1//2^(2),save_everystep=false,saveat=[0.0,1/2,1.0])

  @test symdiff(sol.t,sol2.t) == [1/2]

  sol2=solve(prob,DP5(),dt=1//2^(2),save_everystep=false,saveat=1/2)

  @test sort!(symdiff(sol.t,sol2.t)) == [1/2]

  sol3=solve(prob,DP5(),dt=1//2^(2),save_everystep=false,saveat=[1/2],tstops=[1/2])

  @test sol3.t == [1/2]

  sol3=solve(prob,DP5(),dt=1//2^(2),saveat=[0.0,1/2,1.0],tstops=[1/2])

  @test sol3.t == [0.0,1/2,1.0]

  sol3=solve(prob,DP5(),dt=1//2^(2),saveat=1/10,tstops=[1/2])

  @test sol3.t == collect(0.0:0.1:1.00)

  sol =solve(prob,RK4(),dt=1/2^(2),save_everystep=true,adaptive=false)
  sol2=solve(prob,RK4(),dt=1/2^(2),save_everystep=true,adaptive=false,saveat=[.125,.6,.61,.8])

  @test symdiff(sol.t,sol2.t) == [.125,.6,.61,.8]

  sol =solve(prob,Rosenbrock32(),dt=1/2^(2),save_everystep=true)
  sol2=solve(prob,Rosenbrock32(),dt=1/2^(2),save_everystep=true,saveat=[.125,.6,.61,.8])

  @test symdiff(sol.t,sol2.t) == [.125,.6,.61,.8]

  sol =solve(prob,GenericTrapezoid(),dt=1/2^(2),save_everystep=true)
  sol2=solve(prob,GenericTrapezoid(),dt=1/2^(2),save_everystep=true,saveat=[.125,.6,.61,.8])

  @test symdiff(sol.t,sol2.t) == [.125,.6,.61,.8]

  prob = prob_ode_2Dlinear

  sol =solve(prob,DP5(),dt=1//2^(2),save_everystep=true)
  sol2=solve(prob,DP5(),dt=1//2^(2),save_everystep=true,saveat=[0.0,1/2,1.0])

  @test symdiff(sol.t,sol2.t) == [1/2]

  sol =solve(prob,RK4(),dt=1/2^(2),save_everystep=true,adaptive=false)
  sol2=solve(prob,RK4(),dt=1/2^(2),save_everystep=true,adaptive=false,saveat=[0.0,.125,.6,.61,.8,1.0])

  @test symdiff(sol.t,sol2.t) == [.125,.6,.61,.8]

  sol =solve(prob,Rosenbrock32(),dt=1/2^(2),save_everystep=true)
  sol2=solve(prob,Rosenbrock32(),dt=1/2^(2),save_everystep=true,saveat=[.125,.6,.61,.8])

  @test symdiff(sol.t,sol2.t) == [.125,.6,.61,.8]

  sol =solve(prob,GenericTrapezoid(),dt=1/2^(2),save_everystep=false)
  sol2=solve(prob,GenericTrapezoid(),dt=1/2^(2),saveat=[.125,.6,.61,.8])

  @test sort!(symdiff(sol.t,sol2.t)) == [0.0,.125,.6,.61,.8,1.0]

  sol=solve(prob,GenericTrapezoid(),dt=1/2^(2),save_everystep=true,dense=false,saveat=[0,.125,.6,.61,.8])

  @test !(sol.t[2] ≈ 0)

  # Test Iterators

  sol2=solve(prob,DP5(),dt=1//2^(2),save_everystep=false,dense=false,saveat=0:1//100:1)

  @test sol2.t ≈ collect(0:1//100:1)

  sol2=solve(prob,DP5(),dt=1//2^(2),save_everystep=false,dense=false,saveat=range(0,stop=1,length=100))

  @test sol2.t ≈ range(0,stop=1,length=100)

  f = (du,u,p,t) -> prob.f(du,u,p,t)
  prob2 = ODEProblem(f,vec(prob.u0),prob.tspan,1.01)

  sol2=solve(prob2,DP5(),dt=1//2^(2),saveat=.1,save_idxs=1:2:5)

  for u in sol2.u
    @test length(u) == 3
  end

  sol2=solve(prob2,DP5(),dt=1//2^(2),saveat=.1,save_idxs=1:2:5,save_everystep=true)

  sol=solve(prob2,DP5(),dt=1//2^(2),save_start=false)

  @test sol.t[1] == 1//2^(2)
end


text:
using Test
using LogicCircuits
using ProbabilisticCircuits

@testset "Entropy and KLD" begin

    pc1, vtree = load_struct_prob_circuit(
                    zoo_psdd_file("simple2.1.psdd"), zoo_vtree_file("simple2.vtree"))
    pc2, vtree = load_struct_prob_circuit(
                    zoo_psdd_file("simple2.2.psdd"), zoo_vtree_file("simple2.vtree"))
    pc3, vtree = load_struct_prob_circuit(
                    zoo_psdd_file("simple2.3.psdd"), zoo_vtree_file("simple2.vtree"))

    @test entropy(pc1) ≈ 1.2899219826090118
    @test entropy(pc2) ≈ 0.9359472745536583

    @test kl_divergence(pc1, pc2) ≈ 0.5672800167911778
    @test kl_divergence(pc2, pc3) ≈ 0.38966506

end