text:
"""
The `Messages` module defines structs that represent [`HTTP.Request`](@ref)
and [`HTTP.Response`](@ref) Messages.

The `Response` struct has a `request` field that points to the corresponding
`Request`; and the `Request` struct has a `response` field.
The `Request` struct also has a `parent` field that points to a `Response`
in the case of HTTP Redirect.

The Messages module defines `IO` `read` and `write` methods for Messages
but it does not deal with URIs, creating connections, or executing requests.

The `read` methods throw `EOFError` exceptions if input data is incomplete.
and call parser functions that may throw `HTTP.ParsingError` exceptions.
The `read` and `write` methods may also result in low level `IO` exceptions.

### Sending Messages

Messages are formatted and written to an `IO` stream by
[`Base.write(::IO,::HTTP.Messages.Message)`](@ref) and or
[`HTTP.Messages.writeheaders`](@ref).

### Receiving Messages

Messages are parsed from `IO` stream data by
[`HTTP.Messages.readheaders`](@ref).
This function calls [`HTTP.Messages.appendheader`](@ref) and
[`HTTP.Messages.readstartline!`](@ref).

The `read` methods rely on [`HTTP.IOExtras.unread!`](@ref) to push excess
data back to the input stream.

### Headers

Headers are represented by `Vector{Pair{String,String}}`. As compared to
`Dict{String,String}` this allows [repeated header fields and preservation of
order](https://tools.ietf.org/html/rfc7230#section-3.2.2).

Header values can be accessed by name using
[`HTTP.Messages.header`](@ref) and
[`HTTP.Messages.setheader`](@ref) (case-insensitive).

The [`HTTP.Messages.appendheader`](@ref) function handles combining
multi-line values, repeated header fields and special handling of
multiple `Set-Cookie` headers.

### Bodies

The `HTTP.Message` structs represent the Message Body as `Vector{UInt8}`.

Streaming of request and response bodies is handled by the
[`HTTP.StreamLayer`](@ref) and the [`HTTP.Stream`](@ref) `<: IO` stream.
"""

module Messages

export Message, Request, Response, HeaderSizeError,
       reset!,
       iserror, isredirect, ischunked, issafe, isidempotent,
       header, hasheader, setheader, defaultheader, appendheader,
       mkheaders, readheaders, headerscomplete, readtrailers, writeheaders,
       readstartline!, writestartline,
       bodylength, unknown_length,
       load, payload

import ..HTTP

using ..Pairs
using ..@warn
using ..IOExtras
using ..Parsers
import ..Parsers: headerscomplete, reset!

const unknown_length = typemax(Int)

abstract type Message end

"""
    Response <: Message

Represents a HTTP Response Message.

- `version::VersionNumber`
   [RFC7230 2.6](https://tools.ietf.org/html/rfc7230#section-2.6)

- `status::Int16`
   [RFC7230 3.1.2](https://tools.ietf.org/html/rfc7230#section-3.1.2)
   [RFC7231 6](https://tools.ietf.org/html/rfc7231#section-6)

- `headers::Vector{Pair{String,String}}`
   [RFC7230 3.2](https://tools.ietf.org/html/rfc7230#section-3.2)

- `body::Vector{UInt8}`
   [RFC7230 3.3](https://tools.ietf.org/html/rfc7230#section-3.3)

- `request`, the `Request` that yielded this `Response`.
"""

mutable struct Response <: Message
    version::VersionNumber
    status::Int16
    headers::Headers
    body::Vector{UInt8}
    request::Message

    function Response(status::Int=0, headers=[]; body=UInt8[], request=nothing)
        r = new()
        r.version = v"1.1"
        r.status = status
        r.headers = mkheaders(headers)
        r.body = body
        if request != nothing
            r.request = request
        end
        return r
    end
end

Response(s::Int, body::AbstractVector{UInt8}) = Response(s; body=body)
Response(s::Int, body::AbstractString) = Response(s, Vector{UInt8}(body))

Response(bytes) = parse(Response, bytes)

function reset!(r::Response)
    r.version = v"1.1"
    r.status = 0
    if !isempty(r.headers)
        empty!(r.headers)
    end
    if !isempty(r.body)
        empty!(r.body)
    end
end

"""
    Request <: Message

Represents a HTTP Request Message.

- `method::String`
   [RFC7230 3.1.1](https://tools.ietf.org/html/rfc7230#section-3.1.1)

- `target::String`
   [RFC7230 5.3](https://tools.ietf.org/html/rfc7230#section-5.3)

- `version::VersionNumber`
   [RFC7230 2.6](https://tools.ietf.org/html/rfc7230#section-2.6)

- `headers::Vector{Pair{String,String}}`
   [RFC7230 3.2](https://tools.ietf.org/html/rfc7230#section-3.2)

- `body::Vector{UInt8}`
   [RFC7230 3.3](https://tools.ietf.org/html/rfc7230#section-3.3)

- `response`, the `Response` to this `Request`

- `parent`, the `Response` (if any) that led to this request
  (e.g. in the case of a redirect).
   [RFC7230 6.4](https://tools.ietf.org/html/rfc7231#section-6.4)
"""

mutable struct Request <: Message
    method::String
    target::String
    version::VersionNumber
    headers::Headers
    body::Vector{UInt8}
    response::Response
    parent
end

Request() = Request("", "")

function Request(method::String, target, headers=[], body=UInt8[];
                 version=v"1.1", parent=nothing)
    r = Request(method,
                target == "" ? "/" : target,
                version,
                mkheaders(headers),
                body,
                Response(),
                parent)
    r.response.request = r
    return r
end

Request(bytes) = parse(Request, bytes)

mkheaders(h::Headers) = h
mkheaders(h)::Headers = Header[string(k) => string(v) for (k,v) in h]

"""
    issafe(::Request)

https://tools.ietf.org/html/rfc7231#section-4.2.1
"""

issafe(r::Request) = r.method in ["GET", "HEAD", "OPTIONS", "TRACE"]

"""
    isidempotent(::Request)

https://tools.ietf.org/html/rfc7231#section-4.2.2
"""

isidempotent(r::Request) = issafe(r) || r.method in ["PUT", "DELETE"]

"""
    iserror(::Response)

Does this `Response` have an error status?
"""

iserror(r::Response) = r.status != 0 && r.status != 100 && r.status != 101 &&
                       (r.status < 200 || r.status >= 300) && !isredirect(r)

"""
    isredirect(::Response)

Does this `Response` have a redirect status?
"""
isredirect(r::Response) = r.status in (301, 302, 307, 308)

"""
    statustext(::Response) -> String

`String` representation of a HTTP status code. e.g. `200 => "OK"`.
"""

statustext(r::Response) = Base.get(STATUS_MESSAGES, r.status, "Unknown Code")

"""
    header(::Message, key [, default=""]) -> String

Get header value for `key` (case-insensitive).
"""
header(m::Message, k, d="") = header(m.headers, k, d)
header(h::Headers, k::String, d::String="") = getbyfirst(h, k, k => d, lceq)[2]
lceq(a,b) = lowercase(a) == lowercase(b)

"""
    hasheader(::Message, key) -> Bool

Does header value for `key` exist (case-insensitive)?
"""
hasheader(m, k::String) = header(m, k) != ""

"""
    hasheader(::Message, key, value) -> Bool

Does header for `key` match `value` (both case-insensitive)?
"""
hasheader(m, k::String, v::String) = lowercase(header(m, k)) == lowercase(v)

"""
    setheader(::Message, key => value)

Set header `value` for `key` (case-insensitive).
"""
setheader(m::Message, v) = setheader(m.headers, v)
setheader(h::Headers, v::Pair) = setbyfirst(h, Pair{String,String}(v), lceq)

"""
    defaultheader(::Message, key => value)

Set header `value` for `key` if it is not already set.
"""

function defaultheader(m, v::Pair)
    if header(m, first(v)) == ""
        setheader(m, v)
    end
    return
end

"""
    ischunked(::Message)

Does the `Message` have a "Transfer-Encoding: chunked" header?
"""

ischunked(m) = any(h->(lowercase(h[1]) == "transfer-encoding" &&
                       endswith(lowercase(h[2]), "chunked")),
                   m.headers)

"""
    appendheader(::Message, key => value)

Append a header value to `message.headers`.

If `key` is `""` the `value` is appended to the value of the previous header.

If `key` is the same as the previous header, the `value` is [appended to the
value of the previous header with a comma
delimiter](https://stackoverflow.com/a/24502264)

`Set-Cookie` headers are not comma-combined because [cookies often contain
internal commas](https://tools.ietf.org/html/rfc6265#section-3).
"""

function appendheader(m::Message, header::Pair{String,String})
    c = m.headers
    k,v = header
    if k == ""
        c[end] = c[end][1] => string(c[end][2], v)
    elseif k != "Set-Cookie" && length(c) > 0 && k == c[end][1]
        c[end] = c[end][1] => string(c[end][2], ", ", v)
    else
        push!(m.headers, header)
    end
    return
end

"""
    httpversion(::Message)

e.g. `"HTTP/1.1"`
"""

httpversion(m::Message) = "HTTP/$(m.version.major).$(m.version.minor)"

"""
    writestartline(::IO, ::Message)

e.g. `"GET /path HTTP/1.1\\r\\n"` or `"HTTP/1.1 200 OK\\r\\n"`
"""

function writestartline(io::IO, r::Request)
    write(io, "$(r.method) $(r.target) $(httpversion(r))\r\n")
    return
end

function writestartline(io::IO, r::Response)
    write(io, "$(httpversion(r)) $(r.status) $(statustext(r))\r\n")
    return
end

"""
    writeheaders(::IO, ::Message)

Write `Message` start line and
a line for each "name: value" pair and a trailing blank line.
"""

function writeheaders(io::IO, m::Message)
    writestartline(io, m)                 # FIXME To avoid fragmentation, maybe
    for (name, value) in m.headers        # buffer header before sending to `io`
        write(io, "$name: $value\r\n")
    end
    write(io, "\r\n")
    return
end

"""
    write(::IO, ::Message)

Write start line, headers and body of HTTP Message.
"""

function Base.write(io::IO, m::Message)
    writeheaders(io, m)
    write(io, m.body)
    return
end

function Base.String(m::Message)
    io = IOBuffer()
    write(io, m)
    String(take!(io))
end

#Like https://github.com/JuliaIO/FileIO.jl/blob/v0.6.1/src/FileIO.jl#L19 ?
load(m::Message) = payload(m, String)

function payload(m::Message)::Vector{UInt8}
    enc = lowercase(first(split(header(m, "Transfer-Encoding"), ", ")))
    return enc in ["", "identity", "chunked"] ? m.body : decode(m, enc)
end

payload(m::Message, ::Type{String}) =
    hasheader(m, "Content-Type", "ISO-8859-1") ? iso8859_1_to_utf8(payload(m)) :
                                                 String(payload(m))

function decode(m::Message, encoding::String)::Vector{UInt8}
    if enc == "gzip"
        # Use https://github.com/bicycle1885/TranscodingStreams.jl ?
    end
    @warn "Decoding of HTTP Transfer-Encoding is not implemented yet!"
    return m.body
end

"""
    readstartline!(::Parsers.Message, ::Message)

Read the start-line metadata from Parser into a `::Message` struct.
"""

function readstartline!(m::Parsers.Message, r::Response)
    r.version = VersionNumber(m.major, m.minor)
    r.status = m.status
    return
end

function readstartline!(m::Parsers.Message, r::Request)
    r.version = VersionNumber(m.major, m.minor)
    r.method = m.method
    r.target = m.target
    return
end

"""
Arbitrary limit to protect against denial of service attacks.
"""
const header_size_limit = 0x10000

struct HeaderSizeError <: Exception end

"""
    readheaders(::IO, ::Parser, ::Message)

Read headers (and startline) from an `IO` stream into a `Message` struct.
Throw `EOFError` if input is incomplete.
"""

function readheaders(io::IO, parser::Parser, message::Message)

    n = 0
    while !headerscomplete(parser) && !eof(io)
        bytes = readavailable(io)
        n += length(bytes)
        excess = parseheaders(parser, bytes) do h
            appendheader(message, h)
        end
        unread!(io, excess)
        n -= length(excess)
        if n > header_size_limit
            throw(HeaderSizeError())
        end
    end
    if !headerscomplete(parser)
        throw(EOFError())
    end
    readstartline!(parser.message, message)
    return message
end

"""
    headerscomplete(::Message)

Have the headers been read into this `Message`?
"""

headerscomplete(r::Response) = r.status != 0 && r.status != 100
headerscomplete(r::Request) = r.method != ""

"""
    readtrailers(::IO, ::Parser, ::Message)

Read trailers from an `IO` stream into a `Message` struct.
"""

function readtrailers(io::IO, parser::Parser, message::Message)
    if messagehastrailing(parser)
        readheaders(io, parser, message)
    end
    return message
end

"""
"The presence of a message body in a response depends on both the
 request method to which it is responding and the response status code.
 Responses to the HEAD request method never include a message body [].
 2xx (Successful) responses to a CONNECT request method (Section 4.3.6 of
 [RFC7231]) switch to tunnel mode instead of having a message body.
 All 1xx (Informational), 204 (No Content), and 304 (Not Modified)
 responses do not include a message body.  All other responses do
 include a message body, although the body might be of zero length."
[RFC7230 3.3](https://tools.ietf.org/html/rfc7230#section-3.3)
"""

bodylength(r::Response)::Int =
                 r.request.method == "HEAD" ? 0 :
                     r.status in [204, 304] ? 0 :
    (l = header(r, "Content-Length")) != "" ? parse(Int, l) :
                                              unknown_length

"""
"The presence of a message body in a request is signaled by a
 Content-Length or Transfer-Encoding header field.  Request message
 framing is independent of method semantics, even if the method does
 not define any use for a message body."
[RFC7230 3.3](https://tools.ietf.org/html/rfc7230#section-3.3)
"""

bodylength(r::Request)::Int =
    ischunked(r) ? unknown_length :
                   parse(Int, header(r, "Content-Length", "0"))

"""
    readbody(::IO, ::Parser) -> Vector{UInt8}

Read message body from an `IO` stream.
"""

function readbody(io::IO, parser::Parser, m::Message)
    if ischunked(m)
        body = IOBuffer()
        while !bodycomplete(parser) && !eof(io)
            data, excess = parsebody(parser, readavailable(io))
            write(body, data)
            unread!(io, excess)
        end
        m.body = take!(body)
    else
        l = bodylength(m)
        m.body = read(io, l)
        if l != unknown_length && length(m.body) < l
            throw(EOFError())
        end
    end
end

function Base.parse(::Type{T}, str::AbstractString) where T <: Message
    bytes = IOBuffer(str)
    p = Parser()
    r = Request()
    m::T = T == Request ? r : r.response
    readheaders(bytes, p, m)
    readbody(bytes, p, m)
    readtrailers(bytes, p, m)
    if ischunked(m) && !messagecomplete(p)
        throw(EOFError())
    end
    return m
end

"""
    set_show_max(x)

Set the maximum number of body bytes to be displayed by `show(::IO, ::Message)`
"""

set_show_max(x) = global body_show_max = x
body_show_max = 1000

"""
    bodysummary(bytes)

The first chunk of the Message Body (for display purposes).
"""
bodysummary(bytes) = view(bytes, 1:min(length(bytes), body_show_max))

function compactstartline(m::Message)
    b = IOBuffer()
    writestartline(b, m)
    strip(String(take!(b)))
end

function Base.show(io::IO, m::Message)
    if get(io, :compact, false)
        print(io, compactstartline(m))
        if m isa Response
            print(io, " <= (", compactstartline(m.request::Request), ")")
        end
        return
    end
    println(io, typeof(m), ":")
    println(io, "\"\"\"")
    writeheaders(io, m)
    summary = bodysummary(m.body)
    write(io, summary)
    if length(m.body) > length(summary)
        println(io, "\n⋮\n$(length(m.body))-byte body")
    end
    print(io, "\"\"\"")
    return
end

const STATUS_MESSAGES = (()->begin
    @assert ccall(:jl_generating_output, Cint, ()) == 1
    v = fill("Unknown Code", 530)
    v[100] = "Continue"
    v[101] = "Switching Protocols"
    v[102] = "Processing"                            # RFC 2518 => obsoleted by RFC 4918
    v[200] = "OK"
    v[201] = "Created"
    v[202] = "Accepted"
    v[203] = "Non-Authoritative Information"
    v[204] = "No Content"
    v[205] = "Reset Content"
    v[206] = "Partial Content"
    v[207] = "Multi-Status"                          # RFC 4918
    v[300] = "Multiple Choices"
    v[301] = "Moved Permanently"
    v[302] = "Moved Temporarily"
    v[303] = "See Other"
    v[304] = "Not Modified"
    v[305] = "Use Proxy"
    v[307] = "Temporary Redirect"
    v[400] = "Bad Request"
    v[401] = "Unauthorized"
    v[402] = "Payment Required"
    v[403] = "Forbidden"
    v[404] = "Not Found"
    v[405] = "Method Not Allowed"
    v[406] = "Not Acceptable"
    v[407] = "Proxy Authentication Required"
    v[408] = "Request Time-out"
    v[409] = "Conflict"
    v[410] = "Gone"
    v[411] = "Length Required"
    v[412] = "Precondition Failed"
    v[413] = "Request Entity Too Large"
    v[414] = "Request-URI Too Large"
    v[415] = "Unsupported Media Type"
    v[416] = "Requested Range Not Satisfiable"
    v[417] = "Expectation Failed"
    v[418] = "I'm a teapot"                        # RFC 2324
    v[422] = "Unprocessable Entity"                # RFC 4918
    v[423] = "Locked"                              # RFC 4918
    v[424] = "Failed Dependency"                   # RFC 4918
    v[425] = "Unordered Collection"                # RFC 4918
    v[426] = "Upgrade Required"                    # RFC 2817
    v[428] = "Precondition Required"               # RFC 6585
    v[429] = "Too Many Requests"                   # RFC 6585
    v[431] = "Request Header Fields Too Large"     # RFC 6585
    v[440] = "Login Timeout"
    v[444] = "nginx error: No Response"
    v[495] = "nginx error: SSL Certificate Error"
    v[496] = "nginx error: SSL Certificate Required"
    v[497] = "nginx error: HTTP -> HTTPS"
    v[499] = "nginx error or Antivirus intercepted request or ArcGIS error"
    v[500] = "Internal Server Error"
    v[501] = "Not Implemented"
    v[502] = "Bad Gateway"
    v[503] = "Service Unavailable"
    v[504] = "Gateway Time-out"
    v[505] = "HTTP Version Not Supported"
    v[506] = "Variant Also Negotiates"             # RFC 2295
    v[507] = "Insufficient Storage"                # RFC 4918
    v[509] = "Bandwidth Limit Exceeded"
    v[510] = "Not Extended"                        # RFC 2774
    v[511] = "Network Authentication Required"     # RFC 6585
    v[520] = "CloudFlare Server Error: Unknown"
    v[521] = "CloudFlare Server Error: Connection Refused"
    v[522] = "CloudFlare Server Error: Connection Timeout"
    v[523] = "CloudFlare Server Error: Origin Server Unreachable"
    v[524] = "CloudFlare Server Error: Connection Timeout"
    v[525] = "CloudFlare Server Error: Connection Failed"
    v[526] = "CloudFlare Server Error: Invalid SSL Ceritificate"
    v[527] = "CloudFlare Server Error: Railgun Error"
    v[530] = "Site Frozen"
    return v
end)()

end # module Messages


text:
"""
label(node)

Get the label a node.
"""
label(p::AbstractNode) = p.label

"""
label!(node, v)

Set the string `v` as label of `node`.
"""
function label!(p::AbstractNode, v::String)
	p.label = v

	return nothing
end

"""
	getspecies(node)

Get the number of the species assigned to `node`.
"""
getspecies(p) = p.species

"""
	setspecies!(node, i)

Assign the species number `i` to `node`.
"""
function setspecies!(p::AbstractNode, i::Int)
	p.species = i

	return nothing
end

"""
	annotation(node, k)

Get the value of the annotation with key `k` in `node`.
"""
annotation(p::AbstractNode, k::AbstractString) = p.annotations[k]

"""
	annotate!(node, key, value)

Set the `value` for an annotation `key` in `node`.

Annotation keys can only be `AbstractString`s.
"""
function annotate!(p::AbstractNode, k::AbstractString, v::Any)
	p.annotations[k] = v

	return nothing
end

"""
	deannotate!(node, key)

Remove an annotation `key` from a `node`.
"""
function deannotate!(p::AbstractNode, k::AbstractString)
	delete!(p.annotations, k)

	return nothing
end

ext: jl


text:
"""
Used to cache various things the Gibbs sampler needs
"""
mutable struct GibbsSamplerState

    bn::BayesNet
    key_constructor_name_order::Array{Symbol,1}
    max_cache_size::Union{Int, Nothing}
    markov_blanket_cpds_cache::Dict{Symbol, Array{CPD}}
    markov_blanket_cache::Dict{Symbol, Vector{Symbol}}
    finite_distrbution_is_cacheable::Dict{Symbol, Bool}
    finite_distribution_cache::Dict{String, Array{Float64, 1}}

    function GibbsSamplerState(
        bn::BayesNet,
        max_cache_size::Union{Int, Nothing}=nothing
        )

        a = rand(bn)
        markov_blankets = Dict{Symbol, Vector{Symbol}}()
        is_cacheable = Dict{Symbol, Bool}()
        markov_blanket_cpds_cache = Dict{Symbol, Array{CPD}}()
        for name in names(bn)
            markov_blankets[name] = Symbol[ele for ele in markov_blanket(bn, name)]
            is_cacheable[name] = reduce(&, [hasfinitesupport(get(bn, mb_name)(a)) for mb_name in markov_blankets[name]])
            markov_blanket_cpds_cache[name] = get(bn, markov_blanket(bn, name))
        end

        new(bn, names(bn), max_cache_size,
            markov_blanket_cpds_cache,
            markov_blankets,
            is_cacheable,
            Dict{String, Array{Float64, 1}}(),
            )
    end
end

"""
Helper to sample_posterior_finite

Modifies a and gss
"""
function get_finite_distribution!(gss::GibbsSamplerState, varname::NodeName, a::Assignment, support::AbstractArray)

    is_cacheable = gss.finite_distrbution_is_cacheable[varname]
    key = ""
    if is_cacheable

        key = join([string(a[name]) for name in gss.markov_blanket_cache[varname]], ",")
        key = join([string(varname), key], ",")

        if haskey(gss.finite_distribution_cache, key)
            return gss.finite_distribution_cache[key]
        end

    end

    markov_blanket_cpds = gss.markov_blanket_cpds_cache[varname]
    posterior_distribution = zeros(length(support))
    for (index, domain_element) in enumerate(support)
        a[varname] = domain_element
        # Sum logs for numerical stability
        posterior_distribution[index] = exp(sum([logpdf(cpd, a) for cpd in markov_blanket_cpds]) +
                                                 logpdf(get(gss.bn, varname), a)) # because A is not in its own Markov blanket
    end
    posterior_distribution = posterior_distribution / sum(posterior_distribution)
    if is_cacheable && ( gss.max_cache_size == nothing || length(gss.finite_distribution_cache) < gss.max_cache_size )
        gss.finite_distribution_cache[key] = posterior_distribution
    end
    return posterior_distribution
end

"""
Helper to sample_posterior
Should only be called if the variable associated with varname is discrete

set a[varname] ~ P(varname | not varname)

Modifies both a and gss
"""
function sample_posterior_finite!(gss::GibbsSamplerState, name::NodeName, a::Assignment, support::AbstractArray)

   posterior_distribution = get_finite_distribution!(gss, name, a, support)

   # Adapted from Distributions.jl, credit to its authors
   p = posterior_distribution
   n = length(p)
   i = 1
   c = p[1]
   u = rand()
   while c < u && i < n
       c += p[i += 1]
   end

   a[name] = support[i]
end

"""
Implements Metropolis-Hastings with a normal distribution proposal with mean equal to the previous value
of the variable "varname" and stddev equal to 10 times the standard deviation of the distribution of the target
variable given its parents ( var_distribution should be get(bn, varname)(a) )

MH will go through nsamples iterations.  If no proposal is accepted, the original value will remain

This function expects that a[varname] is within the support of the distribution, it will not check to make sure this is true

Helper to sample_posterior
Should only be used to sampling continuous distributions

set a[varname] ~ P(varname | not varname)

Modifies a and caches in gss
"""
function sample_posterior_continuous!(
    gss::GibbsSamplerState,
    varname::NodeName,
    a::Assignment,
    var_distribution::ContinuousUnivariateDistribution;
    MH_iterations::Integer=10,
    )

    # Random Walk Metropolis Hastings
    markov_blanket_cpds = gss.markov_blanket_cpds_cache[varname]
    stddev = std(var_distribution) * 10.0
    previous_sample_scaled_true_prob = exp(sum([logpdf(cpd, a) for cpd in markov_blanket_cpds]))
    proposal_distribution = Normal(a[varname], stddev) # TODO why does calling this constructor take so long?

    for sample_iter = 1:MH_iterations

        # compute proposed jump
        current_value = a[varname]
        proposed_jump = rand(proposal_distribution)
        if ~ insupport(var_distribution, proposed_jump)
            continue # reject immediately, zero probability
        end

        # Compute acceptance probability
        a[varname] = proposed_jump
        proposed_jump_scaled_true_prob = exp(sum([logpdf(cpd, a) for cpd in markov_blanket_cpds]))
        # Our proposal is symmetric, so q(X_new, X_old) / q(X_old, X_new) = 1
        # accept_prob = min(1, proposed_jump_scaled_true_prob/previous_sample_scaled_true_prob)
        accept_prob = proposed_jump_scaled_true_prob/previous_sample_scaled_true_prob # min operation is unnecessary

        # Accept or reject and clean up
        if rand() <= accept_prob
            # a[varname] = proposed_jump
            previous_sample_scaled_true_prob = proposed_jump_scaled_true_prob
            proposal_distribution = Normal(proposed_jump, stddev)
        else
            a[varname] = current_value
        end
    end

    # a[varname] is set in the above for loop
end

"""
set a[varname] ~ P(varname | not varname)

Modifies a and caches in gss
"""
function sample_posterior!(gss::GibbsSamplerState, varname::Symbol, a::Assignment)
    bn = gss.bn
    cpd = get(bn, varname)
    distribution = cpd(a)
    if hasfinitesupport(distribution)
        sample_posterior_finite!(gss, varname, a, support(distribution))
    elseif typeof(distribution) <: DiscreteUnivariateDistribution
        error("Infinite Discrete distributions are currently not supported in the Gibbs sampler")
    else
        # sample_posterior_continuous_adaptive!(gss, varname, a, distribution)
        sample_posterior_continuous!(gss, varname, a, distribution)
    end
end

"""
The main loop associated with Gibbs sampling
Returns a data frame with nsamples samples

Supports the various parameters supported by gibbs_sample
Refer to gibbs_sample for parameter meanings
"""
function gibbs_sample_main_loop(
    gss::GibbsSamplerState,
    nsamples::Integer,
    thinning::Integer,
    start_sample::Assignment,
    consistent_with::Assignment,
    variable_order::Union{Vector{Symbol}, Nothing},
    time_limit::Union{Int, Nothing},
    )

    start_time = now()

    bn = gss.bn
    a = start_sample
    if variable_order == nothing
         v_order = names(bn)
    else
         v_order = coalesce(variable_order)
    end

    # v_order = [varname for varname in v_order if ~haskey(consistent_with, varname)]
    v_order_query_only = Symbol[]
    for varname in v_order
        if ~haskey(consistent_with, varname)
            push!(v_order_query_only, varname)
        end
    end
    v_order = v_order_query_only

    t = Dict{Symbol, Vector{Any}}()
    for name in v_order
        t[name] = Any[]
    end

    for sample_iter in 1:nsamples
        if (time_limit!=nothing) && ((now() - start_time).value > time_limit)
            break
        end

        if variable_order == nothing
            v_order = shuffle!(v_order)
        end

        # skip over thinning samples
        for skip_iter in 1:thinning
            for varname in v_order
                 sample_posterior!(gss, varname, a)
            end

            if variable_order == nothing
                v_order = shuffle!(v_order)
            end
        end

        for varname in v_order
            sample_posterior!(gss, varname, a)
            push!(t[varname], a[varname])
        end

    end

    return convert(DataFrame, t), (now() - start_time).value
end

"""
Implements Gibbs sampling. (https://en.wikipedia.org/wiki/Gibbs_sampling)
For finite variables, the posterior distribution is sampled by building the exact distribution.
For continuous variables, the posterior distribution is sampled using Metropolis Hastings MCMC.
Discrete variables with infinite support are currently not supported.
The Gibbs Sampler only supports CPDs that return Univariate Distributions. (CPD{D<:UnivariateDistribution})

bn:: A Bayesian Network to sample from.  bn should only contain CPDs that return UnivariateDistributions.

nsamples: The number of samples to return.

burn_in:  The first burn_in samples will be discarded.  They will not be returned.
The thinning parameter does not affect the burn in period.
This is used to ensure that the Gibbs sampler converges to the target stationary distribution before actual samples are drawn.

thinning: For every thinning + 1 number of samples drawn, only the last is kept.
Thinning is used to reduce autocorrelation between samples.
Thinning is not used during the burn in period.
e.g. If thinning is 1, samples will be drawn in groups of two and only the second sample will be in the output.

time_limit: The number of milliseconds to run the algorithm.
The algorithm will return the samples it has collected when either nsamples samples have been collected or time_limit
milliseconds have passed.  If time_limit is null then the algorithm will run until nsamples have been collected.
This means it is possible that zero samples are returned.

error_if_time_out: If error_if_time_out is true and the time_limit expires, an error will be raised.
If error_if_time_out is false and the time limit expires, the samples that have been collected so far will be returned.
	This means it is possible that zero samples are returned.  Burn in samples will not be returned.
If time_limit is null, this parameter does nothing.

consistent_with: the assignment that all samples must be consistent with (ie, Assignment(:A=>1) means all samples must have :A=1).
Use to sample conditional distributions.

max_cache_size:  If null, cache as much as possible, otherwise cache at most "max_cache_size"  distributions

variable_order: variable_order determines the order of variables changed when generating a new sample.
If null use a random order for every sample (this is different from updating the variables at random).
Otherwise should be a list containing all the variables in the order they should be updated.

initial_sample:  The inital assignment to variables to use.  If null, the initial sample is chosen by
briefly running rand(bn, get_weighted_dataframe).
"""
function gibbs_sample(bn::BayesNet, nsamples::Integer, burn_in::Integer;
        thinning::Integer=0,
        consistent_with::Assignment=Assignment(),
        variable_order::Union{Vector{Symbol}, Nothing}=nothing,
        time_limit::Union{Int, Nothing}=nothing,
        error_if_time_out::Bool=true,
        initial_sample::Union{Assignment, Nothing}=nothing,
        max_cache_size::Union{Int, Nothing}=nothing
        )
    # check parameters for correctness
    nsamples > 0 || throw(ArgumentError("nsamples parameter less than 1"))
    burn_in >= 0 || throw(ArgumentError("Negative burn_in parameter"))
    thinning >= 0 || throw(ArgumentError("Negative thinning parameter"))
    if variable_order != nothing
        v_order = coalesce(variable_order)
        bn_names = names(bn)
        for name in bn_names
            name in v_order || throw(ArgumentError("Gibbs sample variable_order must contain all variables in the Bayes Net"))
        end
        for name in v_order
            name in bn_names || throw(ArgumentError("Gibbs sample variable_order contains a variable not in the Bayes Net"))
        end
    end
    if time_limit != nothing
        time_limit > 0 || throw(ArgumentError(join(["Invalid time_limit specified (", time_limit, ")"])))
    end
    if initial_sample != nothing
        init_sample = coalesce(initial_sample)
        for name in names(bn)
            haskey(init_sample, name) || throw(ArgumentError("Gibbs sample initial_sample must be an assignment with all variables in the Bayes Net"))
        end
        for name in keys(consistent_with)
            init_sample[name] == consistent_with[name] || throw(ArgumentError("Gibbs sample initial_sample was inconsistent with consistent_with"))
        end
        pdf(bn, init_sample) > 0 || throw(ArgumentError("Gibbs sample initial_sample has a pdf value of zero"))
    end

    gss = GibbsSamplerState(bn, max_cache_size)

    # Burn in
    # for burn_in_initial_sample use get_weighted_dataframe, should be consistent with the varibale consistent_with
    if initial_sample == nothing
        rand_samples = get_weighted_dataframe(bn, 50, consistent_with)
    	if reduce(|, isnan.(convert(Array{AbstractFloat}, rand_samples[!,:p])))
    		error("Gibbs Sampler was unable to find an inital sample with non-zero probability, please supply an inital sample")
    	end
        burn_in_initial_sample = sample_weighted_dataframe(rand_samples)
    else
        burn_in_initial_sample = initial_sample
    end
    burn_in_samples, burn_in_time = gibbs_sample_main_loop(gss, burn_in, 0, burn_in_initial_sample,
                                         consistent_with, variable_order, time_limit)

    # Check that more time is available
    remaining_time = nothing
    if time_limit != nothing
        remaining_time = time_limit - burn_in_time
        if error_if_time_out
            remaining_time > 0 || error("Time expired during Gibbs sampling")
        end
    end

    # Real samples
    main_samples_initial_sample = burn_in_initial_sample
    if burn_in != 0 && size(burn_in_samples)[1] > 0
        main_samples_initial_sample = Assignment()
        for varname in names(bn)
            main_samples_initial_sample[varname] = haskey(consistent_with, varname) ?
                                                       consistent_with[varname] : burn_in_samples[end, varname]
        end
    end
    samples, samples_time = gibbs_sample_main_loop(gss, nsamples, thinning,
                               main_samples_initial_sample, consistent_with, variable_order, remaining_time)
    combined_time = burn_in_time + samples_time
    if error_if_time_out && time_limit!=nothing
        combined_time < time_limit || error("Time expired during Gibbs sampling")
    end

    # Add in columns for variables that were conditioned on
    evidence = Dict{Symbol, Array{Float64,1}}()
    for varname in keys(consistent_with)
        evidence[varname] = ones(size(samples)[1]) * consistent_with[varname]
    end
    evidence = DataFrame(evidence)
    return hcat(samples, evidence)
end

"""
The GibbsSampler type houses the parameters of the Gibbs sampling algorithm.  The parameters are defined below:

burn_in:  The first burn_in samples will be discarded.  They will not be returned.
The thinning parameter does not affect the burn in period.
This is used to ensure that the Gibbs sampler converges to the target stationary distribution before actual samples are drawn.

thinning: For every thinning + 1 number of samples drawn, only the last is kept.
Thinning is used to reduce autocorrelation between samples.
Thinning is not used during the burn in period.
e.g. If thinning is 1, samples will be drawn in groups of two and only the second sample will be in the output.

time_limit: The number of milliseconds to run the algorithm.
The algorithm will return the samples it has collected when either nsamples samples have been collected or time_limit
milliseconds have passed.  If time_limit is null then the algorithm will run until nsamples have been collected.
This means it is possible that zero samples are returned.

error_if_time_out: If error_if_time_out is true and the time_limit expires, an error will be raised.
If error_if_time_out is false and the time limit expires, the samples that have been collected so far will be returned.
        This means it is possible that zero samples are returned.  Burn in samples will not be returned.
If time_limit is null, this parameter does nothing.

consistent_with: the assignment that all samples must be consistent with (ie, Assignment(:A=>1) means all samples must have :A=1).
Use to sample conditional distributions.

max_cache_size:  If null, cache as much as possible, otherwise cache at most "max_cache_size"  distributions

variable_order: variable_order determines the order of variables changed when generating a new sample.
If null use a random order for every sample (this is different from updating the variables at random).
Otherwise should be a list containing all the variables in the order they should be updated.

initial_sample:  The inital assignment to variables to use.  If null, the initial sample is chosen by
briefly using a LikelihoodWeightedSampler.
"""
mutable struct GibbsSampler <: BayesNetSampler

    evidence::Assignment
    burn_in::Int
    thinning::Int
    variable_order::Union{Vector{Symbol}, Nothing}
    time_limit::Union{Int, Nothing}
    error_if_time_out::Bool
    initial_sample::Union{Assignment, Nothing}
    max_cache_size::Union{Int, Nothing}

    function GibbsSampler(evidence::Assignment=Assignment();
        burn_in::Int=100,
        thinning::Int=0,
        variable_order::Union{Vector{Symbol}, Nothing}=nothing,
        time_limit::Union{Int, Nothing}=nothing,
        error_if_time_out::Bool=true,
        initial_sample::Union{Assignment, Nothing}=nothing,
        max_cache_size::Union{Int, Nothing}=nothing
        )

        new(evidence, burn_in, thinning, variable_order, time_limit, error_if_time_out, initial_sample, max_cache_size)
    end
end

"""
Implements Gibbs sampling. (https://en.wikipedia.org/wiki/Gibbs_sampling)
For finite variables, the posterior distribution is sampled by building the exact distribution.
For continuous variables, the posterior distribution is sampled using Metropolis Hastings MCMC.
Discrete variables with infinite support are currently not supported.
The Gibbs Sampler only supports CPDs that return Univariate Distributions. (CPD{D<:UnivariateDistribution})

Sampling requires a GibbsSampler object which contains the parameters for Gibbs sampling.
See the GibbsSampler documentation for parameter details.
"""
function Base.rand(bn::BayesNet, sampler::GibbsSampler, nsamples::Integer)

    return gibbs_sample(bn, nsamples, sampler.burn_in; thinning=sampler.thinning,
		consistent_with=sampler.evidence, variable_order=sampler.variable_order,
		time_limit=sampler.time_limit, error_if_time_out=sampler.error_if_time_out,
		initial_sample=sampler.initial_sample, max_cache_size=sampler.max_cache_size)
end

"""
    NOTE: this is inefficient. Use rand(bn, GibbsSampler, nsamples) whenever you can
"""
function Random.rand!(a::Assignment, bn::BayesNet, sampler::GibbsSampler)
    df = rand(bn, sampler, 1)
    for name in names(bn)
        a[name] = df[1, name]
    end
    a
end


text:
@testset "Misc/MPolyFactor" begin

  function check_factorization(num::Integer, a)
    f = factor(a)
    g = f.unit
    for i in f.fac
      num -= i[2]
      g *= i[1]^i[2]
    end
    @test num <= 0
    @test g == a
  end

  R, (x, y, z) = PolynomialRing(Hecke.QQ, ["x", "y", "z"])
  check_factorization(6, x^24-y^24*z^12)

  R, (x, y, z) = PolynomialRing(Hecke.QQ, ["x", "y", "z"])
  check_factorization(2, ((x^2+y^2+z^2+2)*(x+1)*(y+2)*(z+3)+x*y*z)*(x^2+y^2+z^2+3))

  R, (x, y, z, t) = PolynomialRing(Hecke.QQ, ["x", "y", "z", "t"])
  for n in 1:8
    count = 0
    for i in 1:2:n
        if n%i == 0
            count += 1
        end
    end
    check_factorization(count, ((1+x+y+z+t)^n+1)*((1+x+y+z+t)^n+2))
  end

  R, (x, y, z, u, v) = PolynomialRing(Hecke.QQ, ["x", "y", "z", "u", "v"])
  check_factorization(2, (y^2*v*(z+v)*x+v+y)*(v^2*(y+z)*x+z^2))

  QA, A = PolynomialRing(Hecke.QQ, "A")
  K, a = number_field(A^3 - 2, "a")
  R, (x, y, z) = PolynomialRing(K, ["x", "y", "z"])
  check_factorization(0, R(0))
  check_factorization(0, R(1))
  check_factorization(0, R(2))
  check_factorization(6, 2*x^2*y*(1+y)*(x^3-2))
  check_factorization(7, 2*x^2*y*(1+y)*(2*x^3-1)*(a*x+a^2*y+1)*(a^2*x+a*y+1)^2)
  check_factorization(2, x^6-2*y^3)
  check_factorization(2, (y*z*x^2+z-1)*(x+a*z^3*y))
  check_factorization(3, (x+y+a)*(x+a*y+1)*(x+y+2*a))
  check_factorization(3, (x^2+y^2+a)*(x^3+a*y^3+1)*(x^2+y+2*a))
  check_factorization(8, (1+x+a*y+x*y)^3*(1+x+a^2*y+x*y^2)^5)
  check_factorization(4, (x-(a+1)^4*y+z)*(a*x+y+z)*(x+a*y+z)*(x+y+a*z))

end


text:
struct ModulationFilterbank{F<:FilterCoefficients,T<:Real,U<:Real} <: Filterbank
    filters::Vector{F}
    center_frequencies::Vector{T}
    fs::U
end

function make_modulation_filter(w0, Q)
    W0 = tan(w0/2)
    B0 = W0/Q
    b = [B0, 0, -B0]
    a = [(1 + B0 + W0^2), (2*W0^2 - 2), (1 - B0 + W0^2)]
    BiquadFilter(b[1], b[2], b[3], a[1], a[2], a[3])
end

function modulation_filterbank(mf, fs, q)
    ModulationFilterbank([make_modulation_filter(w0, q) for w0 in 2*pi*mf/fs], mf, fs)
end

function compute_modulation_cfs(min_cf, max_cf, n)
    spacing_factor = (max_cf/min_cf)^(1/(n-1))
    cfs = zeros(n)
    cfs[1] = min_cf
    for k=2:n
        cfs[k] = cfs[k-1]*spacing_factor
    end
    return cfs
end


text:
using ..CryptoGroups: static

struct ECPoint{P<:AbstractPoint, S} <: AbstractPoint # The same contract is satisfied thus a subtype
    p::P

    function ECPoint{P}(order::Integer, cofactor::Integer; name=nothing) where P <: AbstractPoint
        svars = static(; order, cofactor, name)
        return ECPoint{P, svars}
    end

    function ECPoint{P, S}(x::P) where {P <: AbstractPoint, S}

        @assert oncurve(x) "Point is not in curve"
        # A test with cofactor also here
        new{P, S}(x)
    end

    function ECPoint(p::P, order::Integer, cofactor::Integer; name=nothing) where P <: AbstractPoint
        EP = ECPoint{P}(order, cofactor; name)
        return EP(p)
    end
end

order(::Type{ECPoint{P, S}}) where {P <: AbstractPoint, S} = convert(Integer, S.order)
cofactor(::Type{ECPoint{P, S}}) where {P <: AbstractPoint, S} = convert(Integer, S.cofactor)

#modulus(::Type{ECPoint{P}}) where P <: AbstractPoint = modulus(P) # Only makes sense for point in prime fields!

name(::Type{ECPoint}) = nothing
name(::Type{ECPoint{P}}) where P <: AbstractPoint = nothing
name(::Type{ECPoint{P, S}}) where {P <: AbstractPoint, S} = convert(Symbol, S.name)

eq(::Type{ECPoint{P, S}}) where {P <: AbstractPoint, S} = eq(P)
field(::Type{ECPoint{P, S}}) where {P <: AbstractPoint, S} = field(P)

### May need to do epoint seperatelly
function Base.show(io::IO, ::Type{P}) where P <: ECPoint
    if @isdefined P
        if name(P) == nothing
            Base._show_type(io, P)
        else
            print(io, name(P))
        end
    else
        print(io, "ECPoint")
    end
end

function Base.show(io::IO, p::P) where P <: ECPoint
    show(io, P)
    print(io, " <| (")
    show(io, gx(p))
    print(io, ", ")
    show(io, gy(p))
    print(io, ")")
end

function Base.display(::Type{P}) where P <: ECPoint
    show(P)
    ### I could be more precise on the Point
    ### Like ECPoint{AffinePoint{<:Weierstrass, <:F2GNB}, S}
    if name(P) != nothing
        print(" (alias for ECPoint{<:AffinePoint, ::static_ECPoint})")
    end
end

function Base.:+(x::P, y::P) where P <: ECPoint
    @assert x.p != y.p "AffinePoint's can't be equal. Multiply by 2 instead."
    # Also assertion that x.p != -x.p
    return P(x.p + y.p)
end

Base.:*(x::P, n::Integer) where P <: ECPoint = P(x.p * n)
Base.:*(n::Integer, x::ECPoint) = x * n

<|(::Type{ECPoint{P, S}}, x) where {P <: AbstractPoint, S} = ECPoint{P, S}(P <| x)
validate(p::P) where P <: ECPoint = order(P) * p.p == zero(p.p)

oncurve(p::ECPoint) = oncurve(p.p)

Base.:(==)(x::ECPoint, y::ECPoint) = x.p == y.p

gx(p::ECPoint) = gx(p.p)
gy(p::ECPoint) = gy(p.p)


text:
first = 1
second = 2.1
string = "I am a string"
boolvar = true

println(first)
println(second)
println(string)
println(boolvar)

const MYNAME = "Richard"
println(MYNAME)

function f1()
    x::Int64 = 21
    println(x)
end
f1()

function f2()
    y = Int64(1234)
    z = Int(2)
    println(typeof(y))
    println(typeof(z))
end
f2()

function f3()
    println(typemax(Int64))
    println(typemin(Int64))
end
f3()

println(Sys.WORD_SIZE)
println(zero(Int64))

function f4()
    x1 = 'x'
    println(Int(x1))
    println(Char(120))
end
f4()

myStr = """
Hello from
the
other
side of the room
"""
println(myStr)
println(length(myStr))
println(myStr[1:10])

for c in myStr
    print(Int64(c))
end
println()

x2 = Int64(12)
x3 = Int64(119)
println("The result of $x2 + $x3 is $(x2 + x3)")

max_stars_repo_name: rpillar/JuliaExamples

max_issues_repo_name: rpillar/JuliaExamples

max_forks_repo_name: rpillar/JuliaExamples


text:
"""Module with functions for the computation of profiles and viewshed."""
module Viewshed

"""
    rotate_x( xp::Number, yp::Number, xc::Number, yc::Number, θ::Number )::Int64

Find the cartesian `x` coordinate of a point in `xp` after a rotation of `θ` degrees around a point in `xc`.

Used as part of the `rotate_point` function.
"""
rotate_x( xp::Number, yp::Number, xc::Number, yc::Number, θ::Number )::Int64 = round( Int64, (xp - xc)cos(deg2rad(θ)) - (yp - yc)sin(deg2rad(θ)) + xc )

"""
    rotate_x( xp::Number, yp::Number, xc::Number, yc::Number, θ::Number )::Int64

Find the cartesian `y` coordinate of a point in `yp` after a rotation of `θ` degrees around a point in `yc`.

Used as part of the `rotate_point` function.
"""
rotate_y( xp::Number, yp::Number, xc::Number, yc::Number, θ::Number ) = round( Int64, (xp - xc)sin(deg2rad(θ)) + (yp - yc)cos(deg2rad(θ)) + yc )

"""
    rotate_point( xp::Number, yp::Number, xc::Number, yc::Number, θ::Number )::Tuple{Int64, Int64}

Rotate a point of cartesian coordinates (`xp`, `yp`) around a point (`xc`, `yc`) by `θ` degrees.
"""
rotate_point( xp::Number, yp::Number, xc::Number, yc::Number, θ::Number ) = θ == 0 ? (xp, yp) : ( rotate_x( xp, yp, xc, yc, θ ), rotate_y( xp, yp, xc, yc, θ ) )

# SAREBBE PIU' CORRETTO METTERE dtm::Matrix{T} E DENTRO LA TUPLE DI RITORNO COME PRIMO ELEMENTO Vector{T}
"""
    getprofile( dtm::AbstractArray, r0::Integer, c0::Integer, rn::Integer, cn::Integer )::Tuple{ Vector{T}, Vector{ Tuple{Real, Real} } } where {T <: Number}

Return a `Tuple` of `Vectors` containing the values and the coordinates of all cells on the line between the point at indexes (`r0`, `c0`) and the one at (`rn`, `cn`) of matrix `dtm`.

The function is an implementation of the DDA (Digital Differential Analyser) algorithm for line rasterization.
"""
function getprofile( dtm::AbstractArray, r0::Integer, c0::Integer, rn::Integer, cn::Integer )
    Δr = rn - r0
    Δc = cn - c0
    steps = max( abs(Δr), abs(Δc) )
    r_inc = Δr / steps
    c_inc = Δc / steps
    r_end = r0 + r_inc * steps
    c_end = c0 + c_inc * steps
    heigths_profile = []
    coords_profile = []
    for (r, c) in zip( r0:r_inc:r_end, c0:c_inc:c_end )
        rint, cint = round.(Int64, [r, c])
        push!( heigths_profile, dtm[rint, cint][1] )
        push!( coords_profile, Tuple(GeoArrays.coords( dtm, [rint, cint])) )
    end
    return heigths_profile, coords_profile
end

"""
    function getprofile( dtm::AbstractArray, r0::Integer, c0::Integer, rn::Integer, cn::Integer, h0::Real )

Return a `Tuple` of `Vectors` containing the values and the coordinates of all cells on the line between the point at indexes (`r0`, `c0`) and the one at (`rn`, `cn`) of matrix `dtm`.
The algorithm will stop at the first cell with value greater than the sum of `h0` and the value of cell (`r0`, `c0`).

The function is an implementation of the DDA (Digital Differential Analyser) algorithm for line rasterization.
"""
function getprofile( dtm::AbstractArray, r0::Integer, c0::Integer, rn::Integer, cn::Integer, h0::Real )
    Δr = rn - r0
    Δc = cn - c0
    # Height of the light source
    th0 = dtm[r0, c0] + h0
    # Number of steps to reach the final point
    steps = max( abs(Δr), abs(Δc) )
    r_inc = Δr / steps
    c_inc = Δc / steps
    #   r_end = r0 + r_inc * steps
    #   c_end = c0 + c_inc * steps
    heigths_profile = []
    coords_profile = []
    for (r, c) in zip( r0:r_inc:rn, c0:c_inc:cn )
        rint, cint = round.(Int64, [r, c])
        push!( heigths_profile, dtm[rint, cint][1] )
        push!( coords_profile, Tuple(GeoArrays.coords( dtm, [rint, cint])) )
        if dtm[rint, cint] >= th0
            break
        end
    end
    return heigths_profile, coords_profile
end

#= CREDO SI POSSANO TOGLIERE
    """
        generate_profile( dtm::AbstractArray, x0::Real, y0::Real, xn::Real, yn::Real )::Tuple{ Vector{T}, Vector{ Tuple{Real, Real} } } where {T <: Number}

    Return the profile, formed by a `Vector` of heights and one of coordinates, that spans from the point of coordinates (`x0`, `y0`) to the one at (`xn`, `yn`).
    """
    function generate_profile( dtm::AbstractArray, x0::Real, y0::Real, xn::Real, yn::Real )
        r0, c0 = toIndexes(dtm, x0, y0)
        rn, cn = toIndexes(dtm, xn, yn)
        return getprofile(dtm, r0, c0, rn, cn)
    end

    function generate_profile( dtm::AbstractArray, x0::Real, y0::Real, h0::Real )
        r0, c0 = toIndexes(dtm, x0, y0)
        return getprofile(dtm, r0, c0, h0)
    end
=#

# Given the number of points in a profile, the indexes of the source and the raster, generate 360 profiles
"""
    generate_profiles( dtm, x0::Real, y0::Real, points_number::Integer )::Tuple{ Vector{ Vector{T} }, Vector{ Vector{ Tuple{Integer, Integer} } } } where {T <: Number}

Given a raster `dtm`, a source point of coordinates (`x0`, `y0`) and `points_number` number of desired points for each profile,
return two `Vector`s containing, respectively, the `Vector`s of the heights and the corresponding indexes of the 360 profiles of `dtm`,
obtained as radii with the point as center.
"""
function generate_profiles( dtm::AbstractArray, x0::Real, y0::Real, points_number::Integer )
    values_profiles = []
    indexes_profiles = []
    r0, c0 = toIndexes(dtm, x0, y0)
    rn = r0 + points_number
    for α in [0, 90, 180, 270], β in 0:89
        rn, cn = rotate_point(rn, c0, r0, c0, α+β)
        values, indexes = getprofile(dtm, r0, c0, rn, cn)
        push!(values_profiles, values)
        push!(indexes_profiles, indexes)
    end
    return values_profiles, indexes_profiles
end

"""
Generate the profiles from the source point to the first point higher than the source (the one that will hide all the following)
"""
function generate_profiles( dtm::AbstractArray, x0::Real, y0::Real, h0::Real )
    values_profiles = []
    indexes_profiles = []
    r0, c0 = toIndexes(dtm, x0, y0)
    rn = size(dtm, 1)
    for α in [0, 90, 180, 270], β in 0:89
        rn, cn = rotate_point(rn, c0, r0, c0, α+β)
        values, indexes = getprofile(dtm, r0, c0, rn, cn, h0)
        push!(values_profiles, values)
        push!(indexes_profiles, indexes)
    end
    return values_profiles, indexes_profiles
end

# Given a profile, return only the visible points from the source, if `result == :visible`, all the invisible ones otherwise
"""
    veiwshed( profile::AbstractVector, result::Symbol=:visible )::AbstractVector

Return a `Vector` of all the visible points of `profile`, if `result` is set to `:visible`, or of all the hidden points if set to `:hidden`.

The functions assumes the origin point, from which the viewshed is computed, to be in the first cell of the array.
"""
function veiwshed( profile::AbstractVector, result::Symbol=:visible )::AbstractVector
    if result ∉ [:visible, :hidden]
        throw(DomainError(result, "result must either be `:visible` or `:hidden`."))
    end
    if result == :visible
        visible = [ profile[1], profile[2] ]
        slope = ( profile[2][2] - profile[1][2] ) / ( profile[2][1] - profile[1][1] )
        for i in 3:length(profile)
            new_slope = ( profile[i][2] - profile[1][2] ) / ( profile[i][1] - profile[1][1] )
            if new_slope >= slope
                push!( visible, profile[i] )
                slope = new_slope
            end
        end
        return visible
    else
        hidden = []
        for i in 3:length(profile)
            new_slope = ( profile[i][2] - profile[1][2] ) / ( profile[i][1] - profile[1][1] )
            if new_slope < slope
                push!( nonvisible, profile[i] )
                slope = new_slope
            end
        end
        return hidden
    end
end

distance( x0, y0, x1, y1 ) = √( ( x1 - x0 )^2 + ( y1 - y0 )^2 )
distance( p0, p1 ) = √( ( p1[1] - p0[1] )^2 + ( p1[2] - p0[2] )^2 )

"""
    viewshed( dtm, x0::Real, y0::Real, h0::Real )

Return the visibility matrix containing boolean values indicating wether a cell is reached by the light generated by a source positioned at (`x0`, `y0`) and with relative height of `h0`.
The matrix will have the same size of `dtm`.
"""
function viewshed( dtm, x0::Real, y0::Real, h0::Real )
    # Source cell
    r0, c0 = toIndexes(dtm, x0, y0)
    # Final point of the right horizontal profile
    rm = size(dtm, 1)
    # Total height of the source accounting for terrain
    th0 = dtm[r0, c0] + 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 0:89, β in [0, 90, 180, 270]
        # Final point of the profile of the line with agle `α + β`
        rn, cn = rotate_point( rm, c0, r0, c0, α + β )
        Δr, Δc = (rn - r0, 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 = (rn - r0, cn - c0) ./ max( abs(Δr), abs(Δc) )
        # Indexes of the first cell after the source cell
        r1, c1 = round.(Int64, (r0, c0) .+ c_inc)
        # 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.(Ref(dtm), r1, c1), (0, 0) )
        # Iterate over each cell of the profile on the line
        for (r, c) in zip(r1:r_inc:rn, c1:c_inc:cn)
            # 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.(Ref(dtm), r1, c1), (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] > th0 && break
        end
    end
    return vizmat
end

#= References
    function getprofile( dtm::AbstractArray, r0::Integer, c0::Integer, rn::Integer, cn::Integer, h0::Real )
        Δr = rn - r0
        Δc = cn - c0
        th0 = dtm[r0, c0] + h0
        steps = max( abs(Δr), abs(Δc) )
        r_inc = Δr / steps
        c_inc = Δc / steps
        heigths_profile = []
        coords_profile = []
        for (r, c) in zip( r0:r_inc:rn, c0:c_inc:cn )
            rint, cint = round.(Int64, [r, c])
            push!( heigths_profile, dtm[rint, cint][1] )
            push!( coords_profile, Tuple(GeoArrays.coords( dtm, [rint, cint])) )
            if dtm[rint, cint] >= th0
                break
            end
        end
        return heigths_profile, coords_profile
    end

    function generate_profiles( dtm::AbstractArray, x0::Real, y0::Real, h0::Real )
        values_profiles = []
        indexes_profiles = []
        r0, c0 = toIndexes(dtm, x0, y0)
        rm = size(dtm, 1)
        for α in [0, 90, 180, 270], β in 0:89
            rn, cn = rotate_point(rm, c0, r0, c0, α+β)
            values, indexes = getprofile(dtm, r0, c0, rn, cn, h0)
            push!(values_profiles, values)
            push!(indexes_profiles, indexes)
        end
        return values_profiles, indexes_profiles
    end

    function veiwshed( profile::AbstractVector )::AbstractVector
        visible = [ profile[1], profile[2] ]
        slope = ( profile[2][2] - profile[1][2] ) / ( profile[2][1] - profile[1][1] )
        for i in 3:length(profile)
            new_slope = ( profile[i][2] - profile[1][2] ) / ( profile[i][1] - profile[1][1] )
            if new_slope >= slope
                push!( visible, profile[i] )
                slope = new_slope
            end
        end
        return visible
    end
=#

# === OLD VIEWSHED TEST ===========================================================================================================================

profile1 = [
    (0,12),
    (25,12),
    (50,8),
    (75,10),
    (100,13),
    (125,11),
    (150,15),
    (175,14),
    (200,23)
]
profile2 = [
    (0,500),
    (25,12),
    (50,8),
    (75,10),
    (100,9),
    (125,11),
    (150,15),
    (175,14),
    (200,23)
]
profile3 = [
    (0,12),
    (25,12),
    (50,8),
    (75,10),
    (100,16),
    (125,11),
    (150,500),
    (175,14),
    (200,23)
]
profile = profile3

visible = [ profile[1], profile[2] ]
vali = abs( ( profile[2][2] - profile[1][2] ) / ( profile[2][1] - profile[1][1] ) )
for i in 3:length(profile)
    println("vali: $vali")
    val = ( profile[i][2] - profile[i-1][2] ) / ( profile[i][1] - profile[i-1][1] )
    print("$(profile[i]): $val")
    if val >= abs(vali)
        print("    PUSH")
        push!( visible, profile[i] )
    end
    println("\n")
    vali = vali < 0 ? vali+val : vali-val
end
visible

# ============================================================================================================================================

max_stars_repo_name: DavideFavaro/Tirocinio

max_issues_repo_name: DavideFavaro/Tirocinio

max_forks_repo_name: DavideFavaro/Tirocinio


text:
include("../gcc_common.jl")
build_and_upload_gcc(v"11.0.0-iains")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
@safetestset "isemail" begin
    using Validator

    @test Validator.isemail("foo@bar.com") === true
    @test Validator.isemail("x@x.au") === true
    @test Validator.isemail("foo@bar.com.au") === true
    @test Validator.isemail("foo+bar@bar.com") === true
    @test Validator.isemail("hans.m端ller@test.com") === true
    @test Validator.isemail("hans@m端ller.com") === true
    @test Validator.isemail("test|123@m端ller.com") === true
    @test Validator.isemail("test123+ext@gmail.com") === true
    @test Validator.isemail("some.name.midd.leNa.me.and.locality+extension@GoogleMail.com") === true
    @test Validator.isemail("test.1@gmail.com") === true
    @test Validator.isemail("test@1337.com") === true

    @test Validator.isemail("invalidemail@") === false
    @test Validator.isemail("invalid.com") === false
    @test Validator.isemail("@invalid.com") === false
    @test Validator.isemail("foo@bar.com.") === false
    @test Validator.isemail("somename@ｇｍａｉｌ.com") === false
    @test Validator.isemail("foo@bar.co.uk.") === false
    @test Validator.isemail("z@co.c") === false
    @test Validator.isemail("ｇｍａｉｌｇｍａｉｌｇｍａｉｌｇｍａｉｌｇｍａｉｌ@gmail.com") === false
    @test Validator.isemail("test1@invalid.co m") === false
    @test Validator.isemail("multiple..dots@stillinvalid.com") === false
    @test Validator.isemail("test123+invalid! sub_address@gmail.com") === false

    # Domain Blacklist
    @test Validator.isemail("email@foo.gmail.com", false, false, ["gmail.com", "foo.bar.com"]) === true

    @test Validator.isemail("foo+bar@gmail.com", false, false, ["gmail.com", "foo.bar.com"]) === false
    @test Validator.isemail("email@foo.bar.com", false, false, ["gmail.com", "foo.bar.com"]) === false

    # Ignore max length
    @test Validator.isemail("very-long-email-that-exceed-the-max-charactersvery-long-email-that-exceed-the-max-characters@bar.com") === false

    @test Validator.isemail("very-long-email-that-exceed-the-max-charactersvery-long-email-that-exceed-the-max-characters@bar.com", true) === true

    # Allow Ip Domain
    @test Validator.isemail("email@[123.123.123.123]", false, true) === true
    @test Validator.isemail("email@255.255.255.255", false, true) === true

    @test Validator.isemail("email@0.0.0.256", false, true) === false
    @test Validator.isemail("email@26.0.0.256", false, true) === false
    @test Validator.isemail("email@[266.266.266.266]", false, true) === false
end


text:
using Deprecations
using Test

@testset "Deprecations.jl" begin
    # Write your own tests here.
end


text:
"""
    MovingWindowSolver <: AbstractMovingWindowSolver

A composite type for the moving-window algorithm.

# Constructors
```julia
MovingWindowSolver(parallelsolver, budget, adaptive, τ, Δτ)
MovingWindowSolver(parallelsolver; budget, adaptive, τ, Δτ)
MoWA(args...; kwargs...)
```

# Arguments
- `parallelsolver :: AbstractTimeParallelSolver`
- `adaptive :: Union{AbstractAdaptiveParameters, Nothing}`
- `budget :: Union{AbstractBudget, Nothing}` : budget available.
- `τ :: Real` : time-window length.
- `Δτ :: Real` : time-window shifting.
"""
mutable struct MovingWindowSolver{parallelsolver_T<:AbstractTimeParallelSolver, adaptive_T<:Union{AbstractAdaptiveParameters, Nothing}, budget_T<:Union{AbstractBudget, Nothing}, τ_T<:Real, Δτ_T<:Real} <: AbstractMovingWindowSolver
    parallelsolver::parallelsolver_T
    adaptive::adaptive_T
    budget::budget_T
    τ::τ_T
    Δτ::Δτ_T
end
function MovingWindowSolver(parallelsolver::AbstractTimeParallelSolver; adaptive::Union{AbstractAdaptiveParameters, Nothing}=nothing, budget::Union{AbstractBudget, Nothing}=nothing, τ::Real, Δτ::Real=τ)
    @↓ P = parallelsolver
    if Δτ < τ/P || Δτ > τ
        error("Select τ/P ≤ Δτ ≤ τ.")
    end
    return MovingWindowSolver(parallelsolver, adaptive, budget, τ, Δτ)
end
@doc (@doc MovingWindowSolver) MoWA(args...; kwargs...) = MovingWindowSolver(args...; kwargs...)


text:
#!/usr/bin/env julia

export pe16

pe16(n::Integer=big(2)^1000)=digits(n)|>sum

if !haskey(ENV,"PROJECT_EULER_WITHOUT_RESULT")
    pe16()|>print
end

max_stars_repo_path: 16. Power digit sum/16.jl

max_stars_repo_name: miRoox/ProjectEuler

max_issues_repo_path: 16. Power digit sum/16.jl

max_issues_repo_name: miRoox/ProjectEuler

max_forks_repo_path: 16. Power digit sum/16.jl

max_forks_repo_name: miRoox/ProjectEuler


text:
using CSV, DataFrames, Dates, Plots
include("src/Hymod.jl")

@time begin
    obsPath = "../example_data/LA011201_obs.csv"
    obs = CSV.read(obsPath)
    # obs.isodate = Hymod.Utils.parseDates(obs, format="m/d/Y",dateCol=:date)

    forcingPath = "../example_data/LA011201_forcings.csv"
    forcings = CSV.read(forcingPath)

    forcings.pet = Hymod.hargreaves(forcings,tminCol=:tmin,tmaxCol=:tmax,dtCol=:isodate)

    pars = Hymod.get_random_params()
    q = Hymod.simulate(forcings,precipCol=:precip, petCol=:pet; pars...)
end

max_stars_repo_name: jsta/hymod

max_issues_repo_name: jsta/hymod

max_forks_repo_name: jsta/hymod


text:
export detx

"""
`detx(A::AbstractMatrix{T})` is an exact determinant of the matrix.
Here `T` can be any kind of integer, rational, `Mod`, or `GF2`.

I hope to expand this to `Polynomial`s.
"""
function detx(A::AbstractMatrix{T}) where {T<:IntegerX}
    # @info "Using IntegerX detx{$T}"
    try
        return T(det(A // 1))
    catch
        A = big.(A)
        return detx(A)
    end
end

function detx(A::AbstractMatrix{T}) where {T<:RationalX}
    # @info "Using RationalX detx{$T}"
    try
        return det(A)
    catch
        return det(big.(A))
    end
end

function detx(A::AbstractMatrix{T}) where {T}
    # @info "Using detx! on matrix of type $T"
    r, c = size(A)
    B = Matrix{Any}(undef, r, c)
    for i = 1:r
        for j = 1:c
            B[i, j] = A[i, j]
        end
    end
    # B = copy(A)
    try
        return detx!(B)
    catch
        @warn "Using cofactor expansion to calculate determinant; may be very slow."
        return cofactor_det(B)  # if all else fails!
    end
end

# detx! is the work behind det when we can't use Julia's det.
# it can modify the matrix so this is not exposed to the general public.

function detx!(A::AbstractMatrix{T}) where {T}
    r, c = size(A)
    @assert r == c "Matrix must be square"

    if r == 0
        return 1
    end

    if r == 1
        return A[1, 1]
    end

    if r == 2
        return A[1, 1] * A[2, 2] - A[1, 2] * A[2, 1]
    end

    col = A[:, 1]  # first column
    if all(col .== 0)
        return 0
    end

    k = findfirst(col .!= 0)  # we pivot here

    sign_factor = 1
    if k > 1
        sign_factor = -sign_factor
        row_swap!(A, 1, k)
    end

    factor = A[1, 1] * sign_factor  # multiply by this at the end

    row_scale!(A, 1, 1 // A[1, 1])

    for i = 2:r
        row_add_mult!(A, 1, -A[i, 1], i)
    end

    return detx!(A[2:end, 2:end]) * factor

end


text:
module Issue34

using Interpolations, Base.Test

A = rand(1:20, 100, 100)

# In #34, this incantation throws
itp = interpolate(A, BSpline(Quadratic(Flat())), OnCell())
# Sanity check that not only don't throw, but actually interpolate
for i in 1:size(A,1), j in 1:size(A,2)
    @test_approx_eq itp[i,j] A[i,j]
end

end


text:
using Test
using HDF5
using DSGE

path = dirname(@__FILE__)

m = AnSchorfheide()
Γ0, Γ1, C, Ψ, Π = eqcond(m)
stake = 1 + 1e-6
#G1, C, impact, fmat, fwt, ywt, gev, eu, loose = gensys(Γ0, Γ1, C, Ψ, Π, stake)
G1, C, impact, eu = gensys(Γ0, Γ1, C, Ψ, Π, stake)

file = "$path/../reference/gensys.h5"
G1_exp = h5read(file, "G1_gensys")
C_exp = h5read(file,"C_gensys")
impact_exp = h5read(file, "impact")
eu_exp = h5read(file, "eu")

@testset "Check gensys outputs match reference" begin
    @test @test_matrix_approx_eq G1_exp G1
    @test @test_matrix_approx_eq C_exp C
    @test @test_matrix_approx_eq impact_exp impact

    @test isequal(eu_exp, eu)
end


text:
# Construction of independent buckets of containment boxes

# Given as input a list `randomLineArray` of pairs of 2D points, the function
# `containment2DBoxes` returns, in the same order, the list of _containment boxes_ of the
# input lines. A _containment box_ of a geometric object of dimension $d$ is defined as the
# minimal $d$-cuboid, equioriented with the reference frame, that contains the object. For a
# 2D line it is given by the tuple $(x1,y1,x2,y2)$, where $(x1,y1)$ is the point of minimal
# coordinates, and $(x2,y2)$ is the point of maximal  coordinates.

using LAR

larview(randomLines(1000,0.1)...)

randomLineArray = randomLines(200,0.2);
lines = lar2hpc(randomLineArray...);
yellow = p.COLOR(p.YELLOW)(lines);
boxes = containment2DBoxes(randomLineArray...);
rects = [box2rect(box) for box in boxes];
polylines = p.AA(polyline)(rects);
cyan = p.COLOR(p.CYAN)(p.STRUCT(polylines));
p.VIEW(p.STRUCT([yellow,cyan]))


text:
using IterTools
using Combinatorics

## Tools: functions for multiple iterators
function comb2col(comb, n)
  #given a possible combination of 1:n, returns a matrix with 1 at the choosen values
  temp = zeros(Int, n, 1)
  temp[comb] = 1
  return temp
end

function matcheck(mat, m)
  #returns false if any row of the matrix is zero or any row or column has a sum greater than m
  a = sum(mat, 2) #sum of each row
  if minimum(a) == 0 #zero row
    return false
  end
  if maximum(a) > m #row with sum over m
    return false
  end
  if maximum(sum(mat,1)) > m #col with a sum over m
    return false
  end
  return true
end

## Angle matrix 1st row
function acol1(n)
  #returns an iterator for the first (straight edge) column of the angle matrix
  iter = [zeros(Int, n, 1)] #first element
  for i = 1:n-3 #loop through number of angles containing edge
    comb = combinations(1:n-3, i) #combinations 1:n, choose i
    temp = imap(x->comb2col(x .+ 3, n), comb) #convert to 1s and 0s. the .+3 shifts the combination to ensure that the first 3 elements are zero
    iter = chain(iter, temp) #build the chain
  end
  return iter
end

## Angle matrix top half
function rowpadzeros(comb, k, numzeros)
  #returns the combination with 1 at the choosen values, then padded with zeros
  rowbase = comb2col(comb, k)' #non-zero part of the rowbase. transposed to make a row
  return [zeros(Int, 1, numzeros) rowbase]
end

function aupperrow(n, k)
  #returns the iterator for a angle matrix row with total lenth n.
  #the first n-k elements are always zero
  iter1 = [zeros(Int, 1, n)] #first element
  comb = combinations(1:k, 1) #choose one of the last n-k rows
  iter = imap(x->rowpadzeros(x, k, n-k), comb)
  return chain(iter1, iter)
end

function aupper(n)
  #returns an iterator for the upper triangular half of the angle matrix
  allrows = (aupperrow(n, n-i+1) for i in 1:n) #each row
  a = product(allrows...)
  b = imap(x -> vcat(x...), a)
  return b
end

## Angle matrix bottom half
immutable AngRowGen
  #generates the possible row with sum less than m
  n::Int #total row length
  m::Int #row maximum sum
  l::Int #non-zero row length
end

function Base.start(a::AngRowGen)
  row = zeros(Int, 1, a.n)
  row[1] = -1 #so the first iteration gives all zeros
  return row
end

function Base.next(a::AngRowGen, row)
  for i in 1:a.l #iterate through element, starting with first non-zero element
    if sum(row) < a.m #still room in the raw
      row[i] += 1 #increment and return
      return row, row
    else #no room
      row[i] = 0
    end
  end
  return row, row
end

function Base.done(a::AngRowGen, row)
  return row[a.l] == a.m
end

function alower(n, m)
  #returns an iterator for the upper triangular half of the angle matrix
  allrows = (AngRowGen(n, m, i) for i in 1:n-1)
  a = product(allrows...)
  b = imap(x->vcat(zeros(Int, 1, n), x...), a) #build matrix with zero row on top
  return chain([zeros(Int, n, n)], b) #because product is removing the all zero option
end

## Build anlge matrix
function anglemat(n, m)
  #returns an iterator of the possible angle matrices
  a = product(aupper(n), alower(n, m)) #all upper and lower half combinations
  b = imap(x->x[1] + x[2], a) #add upper and lower half
  c = Iterators.filter(x->matcheck(x, m), b) #filter out zero rows and rows and columns summing over m
  d = product(acol1(n), c) #combine first column and rest of matrix
  f = imap(x->hcat(x...), d)
  return f
end

## Side matrix row
immutable SideRowGen
  #generates the possible row with sum less than m
  n::Int #total row length
  m::Int #row maximum sum
  l::Int #sum of elements from l to end is less than floor(Int, sqrt(m))
end

function Base.start(d::SideRowGen)
  row = zeros(Int, 1, d.n)
  return row
end

function Base.next(d::SideRowGen, row)
  flsqm = floor(Int, sqrt(d.m))
  for i in 1:d.n #iterate through each element
    #(still room in the row) and (element is in beginning of or there is room at the end)
    if (sum(row) < d.m) && ((i < d.l) || sum(row[d.l:end]) < flsqm)
      row[i] += 1
      return row, row
    else #no room
      row[i] = 0
    end
  end
  return row, row
end

function Base.done(d::SideRowGen, row)
  flsqm = floor(Int, sqrt(d.m))
  if row[end] == flsqm #last element at max
    if d.l == 1 #avoid looking for zero index
      return row[end] == flsqm #last unconstrained value + last value
    else
      return row[d.l-1] + row[end] == d.m #last unconstrained value + last value
    end
  end
  return false
end

## Build side matrix
function sidemat(n, m)
  #returns an iterator of the possible side matricies

  ###### correct version. unexpected behavior in product is skipping the first value
  #allrows = (SideRowGen(n, m, i) for i in 1:n) #build all row iterators
  #a = product(allrows...) #combine all row iterators
  ###### end correct version

  ###### begin workaround
  allrows = (SideRowGen(n, m, i) for i in 2:n) #build all row iterators
  firstrow = ([zeros(Int, 1, i) 1 zeros(Int, 1, n-i-1)] for i = 0:n-1)
  a = product(firstrow, allrows...) #combine all row iterators
  ###### end workaround

  b = imap(x->vcat(x...),a)
  c = Iterators.filter(x->matcheck(x,m), b) ##filter out zero rows and rows and columns summing over m
  #(zero rows are impossible, this could be sped up)
  return c
end

max_stars_repo_name: ericbumbalough/Reptiles

max_issues_repo_name: ericbumbalough/Reptiles

max_forks_repo_name: ericbumbalough/Reptiles


text:
abstract type AbstractIndex; end

"""
    accelerate(a, IndexType)

Return an `AcceleratedArray` wrapping `a` using the acceleration index of type `T`.

This operation will not modify `a` but the acceleration index will be invalidated (and
become unsafe to use) if `a` is modified directly after the index is constructed. (See also
`accelerate!`).
"""
function accelerate(a::AbstractArray, ::Type{T}) where {T <: AbstractIndex}
	AcceleratedArray(a, T(a))
end

"""
    accelerate!(a, IndexType)

Return an `AcceleratedArray` wrapping `a` using the acceleration index of type `T`.

Depending on the index type, this operation may also modify `a`. For example a `SortIndex`
will `sort!` the array `a` to maximize cache efficiency. The acceleration index will be
invalidated (and become unsafe to use) if `a` is modified directly after the index is
constructed. (See also `accelerate`).
"""
function accelerate!(a::AbstractArray, ::Type{T}) where {T <: AbstractIndex}
	accelerate(a, T)
end


text:
immutable GaussianCopula
    cov::Array{Float64, 2}
    _chol::Array{Float64, 2}
    # GaussianCopula(cov::Array{Float64,2}) = new(cov, chol(cov, Val{:L}))
    # two argument constructor is allowed for easy copy-paste for debugging
end
GaussianCopula(cov::Array{Float64,2}) = GaussianCopula(cov, ctranspose(chol(cov)))
GaussianCopula(ndim::Int, rho::Float64) = GaussianCopula(rho*ones(ndim, ndim) + (1.0-rho)*diagm(ones(ndim)))

ndim(c::GaussianCopula) = size(c.cov,1)

normal_pdf(x) = exp(-(x.^2)/2.0)/sqrt(2*pi) # assumes zero mean, unit variance
normal_cdf(x) = 0.5*(1.0 + erf(x/sqrt(2.0))) # assumes zero mean, unit variance

function sample_mvn(rng::AbstractRNG, cov_chol::Array{Float64})
    return cov_chol*randn(rng, size(cov_chol,1))
end

function rand(rng::AbstractRNG, c::GaussianCopula, sample::Vector{Float64}=Array(Float64, ndim(c)))
    return normal_cdf(sample_mvn(rng, c._chol))
end
rand(c::GaussianCopula) = rand(Base.GLOBAL_RNG, c)

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
module JuliaDBQuery

using JuliaDB
using OnlineStats
using JSON3
using Dates

const table_cache = Dict()

# yes, this is bad, I'll be adding stuff to our query parser rather than this hack
Base.:(<)(date::Date,s::String) = date < Date(s)
Base.:(>)(date::Date,s::String) = date > Date(s)
Base.:(<=)(date::Date,s::String) = date <= Date(s)
Base.:(>=)(date::Date,s::String) = date >= Date(s)
Base.:(==)(date::Date,s::String) = date == Date(s)

@inline function _all(fs, x)
    result = true
    @inbounds for i in fs
        result = result && i(x)
    end
    return result
end

"loads a table from disk, if it hasn't already loaded, returns the table either way"
function load_table(table_name)
    return get!(table_cache, table_name) do
        load(table_name)
    end
end

function run_basic_query(table, q::String)
    run_basic_query(table,JSON3.read(q,Dict))
end

function run_pipeline_query(table, q::String)
    run_pipeline_query(table, JSON3.read(q,Dict))
end

function run_basic_query(table, q::Dict{String,Any})
    _filter(table, q)
end

function run_pipeline_query(table, q::Array{Dict{String, Any}})
    # we do this first, so that we can work on any parsing errors before we start processing
    # yeah, it isn't the julia way, but it is better this way ;).

    foldl(q; init=table) do step
        # each step is a Dict, with a key describing the operation
        local ks = keys(step)
        local vs = values(step)

        if length(ks) != 1
            throw(ArgumentError("each step in a pipline must have only one instruction \n" * JSON2.write(step)))
        end
        local k = ks[1]
        local v = vs[1]
        if k == "$filter"
            filter(table,v)
        elseif k == "$group"
            group(table,v)
        elseif k == "$fields"
            fields(table,v)
        else k == "$distinct"
            distinct(table,v)
        end
    end
    reduce(pipeline_steps, table)
end

# due to name collisions, I've called anything which is a method for handling
# query fragments, with an _ on the front eg, _filter, is the method which handles
# $filter fragments.

"""
this takes a spec like
{
    home:{
        \$eq:"gore"
    },
    age:{
        \$gt:25,
        \$lt:45
    },
    \$or:[
        {
            job:{
                \$eq:"programmer"
            }
        }, {
            retired:{
                \$eq:true
            }
        }
    ]
}
and returns a function which you can apply to a row, giving a true or false
"""
function predicate_for(q::Dict)
    local symbols = Set{Symbol}()
    functions = map(keys(q) |> collect) do k
        local v = q[k]
        # this if, returns the function for the current branch in the json
        if k == "\$and"
            local fs_with_symbols = map(predicate_for,v) |> collect
            local fs = getindex.(fs_with_symbols,1)
            local symbol_tuples = getindex.(fs_with_symbols,2)
            symbols = symbols ∪ reduce(∪, symbol_tuples)
            r -> _all(fs,r)
        elseif k == "\$or"
            local fs_with_symbols = map(predicate_for,v) |> collect
            local fs = getindex.(fs_with_symbols,1)
            local symbol_tuples = getindex.(fs_with_symbols,2)
            symbols = reduce(∪, symbol_tuples; init=symbols)

            function(r) return any(map(f -> f(r),fs)) end
        else
            # we need to know the fields which we need to select on
            # so add it to the list
            local field = Symbol(k)
            push!(symbols, field)

            # make a function for each predicate
            local fs = map(keys(v) |> collect) do vk
                local vv = v[vk]
                lookup_predicate(vk,field,vv)
            end
            if (length(fs) == 1)
                return fs[1]
            end
            # if the row handles all of the predicates, then the row should return true
            r -> _all(fs,r)
        end
    end
    # if the row handles all of the predicates, then the row should return true
    local symbols_tuple = ((symbols |> collect)...,)
    if (length(functions) == 1)
        return functions[1], symbols_tuple
    end
    return r -> _all(functions,r), symbols_tuple
end

function lookup_predicate(name, field, q)
    # if name == "$and"
    #     _and(q)
    # elseif name == "$or"
    #     _or(q)
    return if name == "\$not"
        _not(field, q)
    elseif name == "\$gt"
        _gt(field, q)
    elseif name == "\$gte"
        _gte(field, q)
    elseif name == "\$lt"
        _lt(field, q)
    elseif name == "\$lte"
        _lte(field, q)
    elseif name == "\$eq"
        _eq(field, q)
    elseif name == "\$in"
        _in(field, q)
    end
end

#
# function _or(q)
#     return x -> x[field] ∈ q
# end

@inline function _eq(field, q)
    return x -> x[field] == q
end

@inline function _not(field, q)
    return x -> !x[field]
end

@inline function _gt(field, q)
    return x -> x[field] > q
end

@inline function _gte(field, q)
    return x -> x[field] >= q
end

@inline function _lt(field, q)
    return x -> x[field] < q
end

@inline function _lte(field, q)
    return x -> x[field] <= q
end

@inline function _in(field, q::Array)
    return x -> x[field] ∈ q
end

@inline function _filter(table, q::Dict)
    @time local f, symbols = predicate_for(q)
    println("filtering for ", symbols)
    @time filter(f, table ; select=symbols)
end

function _group(table, q::Dict)
end

function _fields(table, q::Dict)
end

function _min(table, q::Dict)
end

function _max(table, q::Dict)
end

function _avg(table, q::Dict)
end

function _sum(table, q::Dict)
end

function _distinct(table, q::Dict)
end

function _js(table,q::Dict)
end

function _julia(table,q::Dict)
end

function _r(table,q::Dict)
end

end


text:
# localtestmemaloccrc1.jl
using MicrochipRTDM
using FTD2XX

projectdirectory = "c:/Users/Chris/MPLABXProjects/lbcmcp/PWM_RTDM.x"
mapdict = microchip_parsemap(projectdirectory,"production")
uartconfig1 = UARTConfiguration(921600,8,1,"n",readtimeout = 1000,writetimeout = 1000)
io = open(FT_DeviceIndex(0),uartconfig1)

const x = zeros(UInt64,10)
const buffer = Array(UInt8,1)
function f()
  i = UInt16(0)
  crc = MicrochipRTDM.rtdm_crc(0xffff,i)
  for i::UInt16 in 1:255
    crc = MicrochipRTDM.rtdm_crc(crc,b"asdfg")
  end
end

g() = for i in 1:1000; f(); end
g()

Profile.clear_malloc_data()
g()

close(io)
quit()


text:
using Statistics
using FFTW
using Plots
using BenchmarkTools
using Profile
using LinearAlgebra
using Measures
using ProgressMeter
using Plots
theme(:dark)
using Revise
push!(LOAD_PATH, pwd()*"/main")
using DHC_2DUtils
using DHC_tests

function fink_filter_bank_p(c, L; nx=256, wd=2, t=1, shift=false, Omega=false, safety_on=true, wd_cutoff=1, p=2)
    #c     - sets the scale sampling rate (1 is dyadic, 2 is half dyadic)
    #L     - number of angular bins (usually 8*t or 16*t)
    #wd    - width of the wavelets (default 1, wd=2 for a double covering)
    #t    - plane coverage (default 1, full 2pi 2)
    #shift - shift in θ by 1/2 of the θ spacing
    #Omega - true= append Omega filter (all power beyond Nyquist) so the sum of filters is 1.0

    # -------- assertion errors to make sure arguments are reasonable
    #@test wd <= L/2

    # -------- set parameters
    dθ   = t*π/L
    θ_sh = shift ? dθ/2 : 0.0
    dx   = nx/2-1

    im_scale = convert(Int8,log2(nx))
    # -------- number of bins in radial direction (size scales)
    J = (im_scale-3)*c + 1
    normj = 1/sqrt(c)

    # -------- allocate output array of zeros
    filt      = zeros(nx, nx, J*L+(Omega ? 2 : 1))
    psi_index = zeros(Int32, J, L)
    psi_ind_in= zeros(Int32, J*L+(Omega ? 2 : 1), 2)
    psi_ind_L = zeros(Int32, J*L+(Omega ? 2 : 1))
    theta     = zeros(Float64, L)
    j_value   = zeros(Float64, J)
    info=Dict{String,Any}()

    # -------- compute the required wd
    j_rad_exp = zeros(J)
    for j_ind = 1:J
        j = (j_ind-1)/c
        jrad  = im_scale-j-2
        j_rad_exp[j_ind] = 2^(jrad)
    end

    wd_j = max.(ceil.(wd_cutoff.*L./(t.*π.*j_rad_exp)),wd)

    if !safety_on
        wd_j.=wd
    end

    # loop over wd from small to large
    ## there is some uneeded redundancy added by doing this esp in l loop
    for wd in sort(unique(wd_j))
        # -------- allocate theta and logr arrays
        θ    = zeros(nx, nx)
        logr = zeros(nx, nx)

        wdθ  = wd*dθ
        norm = 1.0/((wd).^(1/p))
        # -------- loop over l
        for l = 0:L-1
            θ_l        = dθ*l+θ_sh
            theta[l+1] = θ_l

        # -------- allocate anggood BitArray
            anggood = falses(nx, nx)

        # -------- loop over pixels
            for x = 1:nx
                sx = mod(x+dx,nx)-dx -1    # define sx,sy so that no fftshift() needed
                for y = 1:nx
                    sy = mod(y+dx,nx)-dx -1
                    θ_pix  = mod(atan(sy, sx)+π -θ_l, 2*π)
                    θ_good = abs(θ_pix-π) <= wdθ

                    # If this is a pixel we might use, calculate log2(r)
                    if θ_good
                        anggood[y, x] = θ_good
                        θ[y, x]       = θ_pix
                        r2            = sx^2 + sy^2
                        logr[y, x]    = 0.5*log2(max(1,r2))
                    end
                end
            end
            angmask = findall(anggood)
        # -------- compute the wavelet in the Fourier domain
        # -------- the angular factor is the same for all j
            F_angular = norm .* cos.((θ[angmask].-π).*(L/(2*wd*t))).^(2/p)

        # -------- loop over j for the radial part
        #    for (j_ind, j) in enumerate(1/c:1/c:im_scale-2)
            j_ind_w_wd = findall(wd_j.==wd)
            for j_ind in j_ind_w_wd
                j = (j_ind-1)/c
                j_value[j_ind] = 1+j  # store for later
                jrad  = im_scale-j-2
                Δj    = abs.(logr[angmask].-jrad)
                rmask = (Δj .<= 1) #deprecating the 1/c to 1, constant width

        # -------- radial part
                F_radial = normj .* cos.(Δj[rmask] .* (π/2)).^(2/p) #deprecating c*π/2 to π/2
                ind      = angmask[rmask]
        #      Let's have these be (J,L) if you reshape...
        #        f_ind    = (j_ind-1)*L+l+1
                f_ind    = j_ind + l*J
                filt[ind, f_ind] = F_radial .* F_angular[rmask]
                psi_index[j_ind,l+1] = f_ind
                psi_ind_in[f_ind,:] = [j_ind-1,l]
                psi_ind_L[f_ind] = 1
            end
        end
    end

    # -------- phi contains power near k=0 not yet accounted for
    filter_power = (sum(filt.^p, dims=3))[:,:,1]

    # -------- for plane half-covered (t=1), add other half-plane
    if t == 1
        filter_power .+= circshift(filter_power[end:-1:1,end:-1:1],(1,1))
    end

    # -------- compute power required to sum to 1.0
    i0 = round(Int16,nx/2-2)
    i1 = round(Int16,nx/2+4)
    center_power = 1.0 .- fftshift(filter_power)[i0:i1,i0:i1]
    zind = findall(center_power .< 1E-15)
    center_power[zind] .= 0.0  # set small numbers to zero
    phi_cen = zeros(nx, nx)
    phi_cen[i0:i1,i0:i1] = (center_power).^(1/p)

    # -------- before adding ϕ to filter bank, renormalize ψ if t=1
    if t==1 filt .*= 2.0^(1/p) end  # double power for half coverage

    # -------- add result to filter array
    phi_index  = J*L+1
    filt[:,:,phi_index] .= fftshift(phi_cen)
    psi_ind_in[phi_index,:] = [J,0]
    psi_ind_L[phi_index] = 0

    if Omega     # append a filter containing the rest (outside Nyquist)
        filter_power += filt[:,:,phi_index].^p
        edge_power    = 1.0 .- filter_power
        zind          = findall(edge_power .< 1E-15)
        edge_power[zind]     .= 0.0  # set small numbers to zero
        Omega_index           = J*L+2
        info["Omega_index"]   = Omega_index
        filt[:,:,Omega_index] = (edge_power).^(1/p)
        psi_ind_in[Omega_index,:] = [J,1]
        psi_ind_L[Omega_index] = 0
    end

    # -------- metadata dictionary
    info["npix"]         = nx
    info["j_value"]      = j_value
    info["theta_value"]  = theta
    info["psi_index"]    = psi_index
    info["phi_index"]    = phi_index
    info["J_L"]          = psi_ind_in
    info["t"]            = t
    info["wd"]           = wd_j
    info["wd_cutoff"]    = wd_cutoff
    info["fs_center_r"]  = j_rad_exp
    info["psi_ind_L"]    = psi_ind_L
    info["p"]            = p

    return filt, info
end

testtest_p = fink_filter_bank_p(1,8)
test_bank = fink_filter_bank(1,8)
test_bank .== test_p

using PyCall
using PyPlot
plt.matplotlib.style.use("dark_background")

plot_filter_bank_QA(test_p[1], test_p[2]; fname="/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/scratch_AKS/images/filter_bank_QA_p2.png")

test_p1 = fink_filter_bank_p(1,8,p=1)

disc = plot_filter_bank_QA(test_p1[1], test_p1[2]; fname="/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/scratch_AKS/images/filter_bank_QA_p1.png",p=1)

heatmap(disc)

test_pp5 = fink_filter_bank_p(1,8,p=0.5)

disc = plot_filter_bank_QA(test_pp5[1], test_pp5[2]; fname="/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/scratch_AKS/images/filter_bank_QA_pp5.png",p=0.5)

heatmap(disc)

test_p3 = fink_filter_bank_p(1,8,p=3)

disc = plot_filter_bank_QA(test_p3[1], test_p3[2]; fname="/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/scratch_AKS/images/filter_bank_QA_p3.png",p=3)

heatmap(disc)

function DHC_compute_p(image::Array{Float64,2}, filter_hash::Dict, filter_hash2::Dict=filter_hash;
    doS2::Bool=true, doS20::Bool=false, norm=true, iso=false, FFTthreads=2, normS1::Bool=false, normS1iso::Bool=false,
    p=2)
    # image        - input for WST
    # filter_hash  - filter hash from fink_filter_hash
    # filter_hash2 - filters for second order.  Default to same as first order.
    # doS2         - compute S2 coeffs
    # doS20        - compute S2 coeffs
    # norm         - scale to mean zero, unit variance
    # iso          - sum over angles to obtain isotropic coeffs

    # Use 2 threads for FFT
    FFTW.set_num_threads(FFTthreads)

    # array sizes
    (Nx, Ny)  = size(image)
    if Nx != Ny error("Input image must be square") end
    (Nf, )    = size(filter_hash["filt_index"])
    if Nf == 0  error("filter hash corrupted") end
    @assert Nx==filter_hash["npix"] "Filter size should match npix"
    @assert Nx==filter_hash2["npix"] "Filter2 size should match npix"
    @assert (normS1 && normS1iso) != 1 "normS1 and normS1iso are mutually exclusive"

    # allocate coeff arrays
    out_coeff = []
    S0  = zeros(Float64, 2)
    S1  = zeros(Float64, Nf)
    if doS2  S2  = zeros(Float64, Nf, Nf) end  # traditional 2nd order
    if doS20 S20 = zeros(Float64, Nf, Nf) end  # real space correlation
    anyM2 = doS2 | doS20

    # allocate image arrays for internal use
    im_rd_0_1  = Array{Float64, 3}(undef, Nx, Ny, Nf)

    ## 0th Order
    S0[1]   = mean(image)
    norm_im = image.-S0[1]
    S0[2]   = sum(abs.(norm_im).^p)/(Nx*Ny)
    if norm
        norm_im ./= (Nx*Ny*S0[2]).^(1/p)
    else
        norm_im = copy(image)
    end

    append!(out_coeff,S0[:])

    ## 1st Order
    im_fd_0 = fft(norm_im)  # total power=1.0

    # unpack filter_hash
    f_ind   = filter_hash["filt_index"]  # (J, L) array of filters represented as index value pairs
    f_val   = filter_hash["filt_value"]

    zarr = zeros(ComplexF64, Nx, Ny)  # temporary array to fill with zvals

    # make a FFTW "plan" for an array of the given size and type
    P = plan_ifft(im_fd_0)  # P is an operator, P*im is ifft(im)

    ## Main 1st Order and Precompute 2nd Order
    for f = 1:Nf
        f_i = f_ind[f]  # CartesianIndex list for filter
        f_v = f_val[f]  # Values for f_i
        if length(f_i) > 0
            zarr[f_i] = f_v.*im_fd_0[f_i]
            im_rd_0_1[:,:,f] .= abs.(P*zarr)
            S1[f] = sum(im_rd_0_1[:,:,f].^(p))  # image power
            #zero out the intermediate arrays
            zarr[f_i] .= 0
        end
    end

    append!(out_coeff, iso ? filter_hash["S1_iso_mat"]*S1 : S1)

    if normS1iso
        S1iso = vec(reshape(filter_hash["S1_iso_mat"]*S1,(1,:))*filter_hash["S1_iso_mat"])
    end

    Mat2 = filter_hash["S2_iso_mat"]
    if doS2
        f_ind2   = filter_hash2["filt_index"]  # (J, L) array of filters represented as index value pairs
        f_val2   = filter_hash2["filt_value"]

        ## Traditional second order
        for f1 = 1:Nf
            thisim = fft(im_rd_0_1[:,:,f1])  # Could try rfft here
            # Loop over f2 and do second-order convolution
            if normS1
                normS1pwr = S1[f1]
            elseif normS1iso
                normS1pwr = S1iso[f1]
            else
                normS1pwr = 1
            end
            for f2 = 1:Nf
                f_i = f_ind2[f2]  # CartesianIndex list for filter
                f_v = f_val2[f2]  # Values for f_i
                zarr[f_i] = f_v .* thisim[f_i]
                S2[f1,f2] = sum(abs.(P*zarr).^(p))/normS1pwr
                #zero out the intermediate arrays
                zarr[f_i] .= 0
            end
        end
        append!(out_coeff, iso ? Mat2*S2[:] : S2[:])
    end

    # Real domain 2nd order
    if doS20
        Amat = reshape(im_rd_0_1, Nx*Ny, Nf)
        S20  = Amat' * Amat
        append!(out_coeff, iso ? Mat2*S20[:] : S20[:])
    end

    return out_coeff
end

test = zeros(16,16,6)
for i=1:6
    a = rand(16,16)
    test[:,:,i] .= a
end

test

filter_hash = fink_filter_hash(1,8)
f_ind   = filter_hash["filt_index"]  # (J, L) array of filters represented as index value pairs
f_val   = filter_hash["filt_value"]

temp_mat = zeros(Float64, 256, 256)
f_i = f_ind[10]  # CartesianIndex list for filter
f_v = f_val[10]  # Values for f_i

temp_mat[f_i] = f_v
temp_mat
@benchmark temp_mat[f_i] = f_v

@benchmark for i = 1:length(f_i)
    ind       = f_i[i]
    temp_mat[ind] = f_v[i]        # filter*image in Fourier domain
end

a = rand(256,256)
temp_mat[f_i] = f_v
temp_mat
@benchmark temp_mat[f_i] = f_v.*a[f_i]

@benchmark for i = 1:length(f_i)
    ind       = f_i[i]
    temp_mat[ind] = f_v[i]*a[ind]       # filter*image in Fourier domain
end

a = rand(256,256)

original = DHC_compute(a,filter_hash)

new = DHC_compute_p(a,filter_hash)

@benchmark DHC_compute_p(a,filter_hash)

@benchmark DHC_compute(a,filter_hash)

6/.178

# so we are slower by slightly less than a factor of Nf... not sure we can help that

function fink_filter_hash_p(c, L; nx=256, wd=2, t=1, shift=false, Omega=false, safety_on=true, wd_cutoff=1, p=2)
    # -------- compute the filter bank
    filt, hash = fink_filter_bank_p(c, L; nx=nx, wd=wd, t=t, shift=shift, Omega=Omega, safety_on=safety_on, wd_cutoff=wd_cutoff, p=p)

    # -------- list of non-zero pixels
    flist = fink_filter_list(filt)

    # -------- pack everything you need into the info structure
    hash["filt_index"] = flist[1]
    hash["filt_value"] = flist[2]

    # -------- compute matrix that projects iso coeffs, add to hash
    S1_iso_mat = S1_iso_matrix(hash)
    hash["S1_iso_mat"] = S1_iso_mat
    S2_iso_mat = S2_iso_matrix(hash)
    hash["S2_iso_mat"] = S2_iso_mat
    hash["num_iso_coeff"] = size(S1_iso_mat)[1] + size(S2_iso_mat)[1] + 2
    hash["num_coeff"] = size(S1_iso_mat)[2] + size(S2_iso_mat)[2] + 2

    return hash
end

filter_hash = fink_filter_hash(1,8)
filter_hashp2 = fink_filter_hash_p(1,8,p=2)

filter_hashp2 == filter_hash

filter_hashp1 = fink_filter_hash_p(1,8,p=1)

@benchmark DHC_compute_p(a,filter_hashp1,p=1)

filter_hashp3 = fink_filter_hash_p(1,8,p=3)

@benchmark DHC_compute_p(a,filter_hashp3,p=3)

filter_hashpp5 = fink_filter_hash_p(1,8,p=0.5)

@benchmark DHC_compute_p(a,filter_hashpp5,p=0.5)

b = [-1, 2, 3, 4]

norm(b)

norm(b,1)

norm(b,3)

norm(b,1/3)

sum(abs.(b).^1)

sqrt(sum(abs.(b).^2))

@benchmark sum(abs.(b).^1)

@benchmark norm(b,1)

@benchmark sum(abs.(b).^2)

@benchmark norm(b,2)^2

Profile.clear()
@profile DHC_compute_p(a,filter_hashp3,p=3)

Juno.profiler()

b = rand(ComplexF64,256,256)

sum(abs.(b).^3)

@benchmark sum(abs.(b).^3)

@benchmark norm(b,3)^3

@benchmark sum(abs2.(b).^(3/2))

sum(abs2.(b).^(3/2))

ref = fink_filter_bank(1,8)
test = fink_filter_bank_p(1,8)

ref == test

a = rand(256,256)

original = DHC_compute(a,filter_hash)

new = DHC_compute_p(a,filter_hash)

original[2:5]./new[2:5]

256*256

original == new

original ≈ new

nearlysame(x, y) = x ≈ y || (isnan(x) & isnan(y))
nearlysame(A::AbstractArray, B::AbstractArray) = all(map(nearlysame, A, B))

nearlysame(original[1:6*8+1],new[1:6*8+1])

maximum(abs.(new.-original))

original[100:105]./new[100:105]

new[100:105]

sum(original[3:end])

sum(original[3:6*8+1])

function rod_image(xcen, ycen, length, pa, fwhm; nx=256)
    # returns image of a rod with some (x,y) position, length,
    #   position angle, and FWHM in an (nx,nx) image.
    rodimage = zeros(nx,nx)

    x=0
    y=0
    sig = fwhm/2.355
    dtor = π/180
    # -------- define a unit vector in direction of rod at position angle pa
    ux = sin(pa*dtor)   #  0 deg is up
    uy = cos(pa*dtor)   # 90 deg to the right

    for i=1:nx
        for j=1:nx
            x=i-nx/2+xcen
            y=j-nx/2+ycen

            # -------- distance parallel and perpendicular to
            dpara =  ux*x + uy*y
            dperp = -uy*x + ux*y

            if abs(dpara)-length <0
                dpara= 0
            end
            dpara = abs(dpara)
            dpara = min(abs(dpara-length),dpara)

            rodimage[i,j] = exp(-(dperp^2+dpara^2)/(2*sig^2))
        end
    end
    rodimage ./= sqrt(sum(rodimage.^2))
    return rodimage
end

rod_test = rod_image(1,1,30,35,8)

heatmap(rod_test)

original = DHC_compute(rod_test,filter_hash)

new = DHC_compute_p(rod_test,filter_hash)

original ≈ new

sum(original[3:end])

sum(new[3:end])

sum(new[3:end])-2 < 1e-6
sum(original[3:end])-2 < 1e-6

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
export Butterfly

module Butterfly

using CSV
using ..Empirikos: MultinomialSummary, TruncatedPoissonSample

const DATA = joinpath(@__DIR__, "butterfly.txt")

function load_table()
    CSV.File(DATA)
end

function ebayes_samples()
    tbl = load_table()
    MultinomialSummary(TruncatedPoissonSample.(tbl.x), tbl.y)
end

end


text:
@testset "Georeferencing" begin
  table = Table(x=rand(3), y=[1,2,3], z=["a","b","c"])
  tuple = (x=rand(3), y=[1,2,3], z=["a","b","c"])

  # explicit domain types
  d = georef(table, PointSet(rand(2,3)))
  @test domain(d) isa PointSet
  d = georef(tuple, PointSet(rand(2,3)))
  @test domain(d) isa PointSet
  d = georef(table, CartesianGrid(3))
  @test domain(d) isa CartesianGrid
  d = georef(tuple, CartesianGrid(3))
  @test domain(d) isa CartesianGrid

  # coordinates of point set
  d = georef(table, rand(2,3))
  @test domain(d) isa PointSet
  d = georef(tuple, rand(2,3))
  @test domain(d) isa PointSet

  # coordinates names in table
  d = georef(table, (:x,:y))
  @test domain(d) isa PointSet

  # coordinates names in named tuple
  d = georef((a=rand(10), b=rand(10), c=rand(10)), (:b,:c))
  @test domain(d) isa PointSet
  @test Tables.columnnames(values(d)) == (:a,)

  # regular grid data
  d = georef(tuple)
  @test domain(d) == CartesianGrid(3)
  tuple2D = (x=rand(10,10), y=rand(Int,10,10))
  d = georef(tuple2D)
  @test domain(d) == CartesianGrid(10,10)
  d = georef(tuple2D, (1.,2.), (3.,4.))
  @test domain(d) == CartesianGrid((10,10), (1.,2.), (3.,4.))
  d = georef(tuple2D, spacing=(3.,4.))
  @test domain(d) == CartesianGrid((10,10), (0.,0.), (3.,4.))
  tuple3D = (x=rand(10,10,10), y=rand(10,10,10))
  d = georef(tuple3D)
  @test domain(d) == CartesianGrid(10,10,10)
end