text:
export read_non_aligned_sequence, write_mapped_contacts, score_BM, PPV

##############################################################
"""
    function read_non_aligned_sequence(seq::String)

    (seq[might be a path]) --> sequence

    USE:
    Given the path of a non aligned sequence (with '.' and insertions)
    returns the sequence. Takes as input .fasta format, or sequence itself,
    as well as the sequence without the path.

 	INPUT:
    `seq`: amino acid sequence or path to it

    OUTPUT:
    String

"""

function read_non_aligned_sequence(seq::String)
    isfasta(file_path::String) = readdlm(file_path)[1][1] == '>'
    seq =   if ispath(seq)
                if isfasta(seq)
                    readdlm(seq)[2]
                else
                    readdlm(seq)
                end

            else
                seq
            end
    return seq
end

##############################################################
"""
    function afa2model(seq::String)

    (seq) --> dictionary

    USE:
    Given a .afa non aligned sequence coming from a HMM
    returns a dictionary that maps its positions to the matched states (model).

    Example: yAC returns the following dictionary

    position       -->  position
    in the model   -->  in the original sequence

    2              -->  1
    3              -->  2

 	INPUT:
    - `seq`: amino acid sequence with insertions

    OUTPUT:
    Dictionary

"""

function afa2model(seq::String)
    seq = read_non_aligned_sequence(seq)
    p_new =  Dict{Int64,Int64}()
    k_lower_and_upper = 0
    k_gap_and_upper = 0
    for lettr in seq
        if islowercase(lettr)
            k_lower_and_upper += 1
        elseif lettr == '-'
            k_gap_and_upper +=1
        elseif isuppercase(lettr)
            k_gap_and_upper +=1
            k_lower_and_upper += 1
            p_new[k_lower_and_upper] = k_gap_and_upper
        end
    end

    if length(p_new) != length([lettr for lettr in seq if isuppercase(lettr)])
        error("Error: there is some inconsistency in the dictionary, verify input.")
    end
    return p_new
end

##############################################################
"""
    function write_mapped_contacts(afa::String, old_distances_path:String; path_out::String = "")
    function write_mapped_contacts(afa::String, old_distances::Array{Float64, 2}; path_out::String = "")

    (afa, old_distances_path, path_out) --> Nothing

    USE:
    Given a .afa sequence with insertions aligned to the distances file
    writes a new distance file with mapped position numbers.

 	INPUT:
    - `afa`: amino acid sequence with insertions

    - `old_distance_path`: path of amino acid distances from crystal data
    - `old_distances`: array with pairs of positions and residue distances

    - `path_out`: output path, if not declared will be the same of "old_distances_path"

    OUTPUT:
    Nothing

"""

function write_mapped_contacts(afa::String, old_distances_path::String; path_out::String = "")
    println("MAKE SURE THAT THE .afa FILE AND THE DISTANCE FILE REFER TO THE SAME SEQUENCE OR ARE ALIGNED.\n")

    afa_sequence = read_non_aligned_sequence(afa)
    l_afa = length(afa_sequence)
    l_gaps = length([lettr for lettr in afa_sequence if lettr == '-'])
    l_upper = length([lettr for lettr in afa_sequence if isuppercase(lettr)])
    full_length = l_afa - l_gaps

    println(
"The length of the full sequence is          $(full_length).
The length of the model-aligned sequence is $(l_upper).
The gaps in the model-aligned sequence are  $(l_gaps).")

    pdb2wt = afa2model(afa)
    old_distances = readdlm(old_distances_path)

    # count how many pairs with distances there will be in the new file a cool julia oneliner
    n_distances = length([i for i in 1:size(old_distances, 1) if haskey(pdb2wt, old_distances[i, 1]) && haskey(pdb2wt, old_distances[i, 2]) ])

    # new file for amino acid distances
    new_distances = zeros(Float64, n_distances, 3)

    # fill new matrix with positions mapped to our models
    l = 0
    for i in 1:size(old_distances, 1)
        if haskey(pdb2wt, old_distances[i, 1]) && haskey(pdb2wt, old_distances[i, 2])
            l +=1
            new_distances[l, :] = [pdb2wt[old_distances[i, 1]], pdb2wt[old_distances[i, 2]] , old_distances[i, 3]]
        end
    end

    path_out =  path_out == "" ? join(split(old_distances_file, '.')[1:end-1], '.')*"_mapped_model.txt" : path_out
    writedlm(path_out, new_distances)
end

function write_mapped_contacts(afa::String, old_distances::Array{Float64, 2}; path_out::String = "")
    println("MAKE SURE THAT THE .afa FILE AND THE DISTANCE FILE REFER TO THE SAME SEQUENCE OR ARE ALIGNED.\n")

    afa_sequence = read_non_aligned_sequence(afa)
    l_afa = length(afa_sequence)
    l_gaps = length([lettr for lettr in afa_sequence if lettr == '-'])
    l_upper = length([lettr for lettr in afa_sequence if isuppercase(lettr)])
    full_length = l_afa - l_gaps

    println(
"The length of the full sequence is          $(full_length).
The length of the model-aligned sequence is $(l_upper).
The gaps in the model-aligned sequence are  $(l_gaps).")

    pdb2wt = afa2model(afa)

    # count how many pairs with distances there will be in the new file a cool julia oneliner
    n_distances = length([i for i in 1:size(old_distances, 1) if haskey(pdb2wt, old_distances[i, 1]) && haskey(pdb2wt, old_distances[i, 2]) ])

    # new file for amino acid distances
    new_distances = zeros(Float64, n_distances, 3)

    # fill new matrix with positions mapped to our models
    l = 0
    for i in 1:size(old_distances, 1)
        if haskey(pdb2wt, old_distances[i, 1]) && haskey(pdb2wt, old_distances[i, 2])
            l +=1
            new_distances[l, :] = [pdb2wt[old_distances[i, 1]], pdb2wt[old_distances[i, 2]] , old_distances[i, 3]]
        end
    end

    path_out =  path_out == "" ? join(split(old_distances_file, '.')[1:end-1], '.')*"_mapped_model.txt" : path_out
    writedlm(path_out, new_distances)
end
;

##############################################################
"""
    function PPV(scores::Array{Float64,2}, distances::Array{Float64,2}; amino_dist = 4, atom_dist = 8)

    (scores, distances) --> PPV

    USE:
    Given a list of residue pairs with a score (arbitrarly defined) and a list of experimental residue distances
    returns the fraction of pairs - out of the highest scoring ones - that are considered "contacts".
    The definition of contacts depends on linear amino acid distance and physical distance between residues.

    INPUT:
    - `scores`: list of protein positions with a score
    - `distances`: list of protein positions with physical distance
    - `amino_dist`: threshold in terms of linear amino acid distance
    - `atom_dist`: threshold in terms of physical distance

    OUTPUT:
    PPV vector

"""

function PPV(scores::Array{Float64,2}, distances::Array{Float64,2}; amino_dist = 4, atom_dist = 8)
    if size(scores[1,:],1)!=3 || size(distances[1,:],1)!=3
        error("contactprediction.jl - PPV: `scores` and `distances` should be in the format `i j value`.")
    end

    scores = sortslices(scores, dims=1, rev=true, by=x->x[3])
    d = Dict{Tuple{Float64,Float64},Float64}()
    for l in 1:size(distances,1)
        d[(distances[l,1], distances[l,2])] = distances[l,3]
    end
    tp = 0
    np = 0
    TP = []
    cc=0
    for r in 1:size(scores,1)
        if abs(scores[r,1] - scores[r,2]) > amino_dist && get(d, (scores[r,1], scores[r,2]), -1.)>=0.
            np +=1
            if d[(scores[r,1], scores[r,2])]<atom_dist
                tp+=1
            else
                # cc<10?println("$(scores[r,1]), $(scores[r,2])"):Void
                cc+=1
            end
            push!(TP,tp/np)
        end
    end

    return TP
end

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

"""
    score_BM(J::Array{Float64,2}; APC::Bool = true, gap::Bool = false)
    score_BM(path_param::AbstractString; APC::Bool = true, gap::Bool = false)

    (J) --> scores
    (path_param) --> scores

    USE:
    Compute Frobenius norm of `q x q` blocks in matrix `J` and adds APC correction.
    This is to use the DCA model as a contact prediction predictor.

    INPUT:
    - `J`: coupling matrix in `q x q x L x L` format
    - `path_param`: path of fields and couplings parameters
    - `APC`: apply the famous APC correction. Default to `true`.
    - `gap`: Remove the state `21` from the Frobenius norm. Default to `false`.

    OUTPUT:
    scores

"""

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

function score_BM(path_param::AbstractString; APC::Bool = true, gap::Bool = false)
    h, J = extract_params(path_param)
    q = size(J, 1)
    L = size(J, 4)

    println("Size of input matrix is ",size(J))

    S = zeros(Float64, L,L)
    for i in 1:L-1
        for j in i+1 : L
            S[i,j] = gap ? sqrt(sum(J[:, : ,i, j].^2)) : sqrt(sum(J[1:20, 1:20, i, j].^2))
            S[j,i] = S[i,j]
        end
    end

    F = zeros(Float64, convert(Int64, L*(L-1)/2),3)
    pos = 1
    for i in 1:L
        for j in i+1:L
            F[pos,:] = [i,j,S[i,j] - ( sum(S[:,j])/L * sum( S[:,i])/L ) / (sum(S)/(L^2)) * APC]
            pos += 1
        end
    end

    return F
end

function score_BM(J::Array{Float64,4}; APC::Bool = true, gap::Bool = false)
    q = size(J, 1)
    L = size(J, 4)

    println("Size of input matrix is ",size(J))

    S = zeros(Float64, L,L)
    for i in 1:L-1
        for j in i+1 : L
            S[i,j] = gap ? sqrt(sum(J[:, : ,i, j].^2)) : sqrt(sum(J[1:20, 1:20, i, j].^2))
            S[j,i] = S[i,j]
        end
    end

    F = zeros(Float64, convert(Int64, L*(L-1)/2),3)
    pos = 1
    for i in 1:L
        for j in i+1:L
            F[pos,:] = [i,j,S[i,j] - ( sum(S[:,j])/L * sum( S[:,i])/L ) / (sum(S)/(L^2)) * APC]
            pos += 1
        end
    end

    return F
end

max_stars_repo_name: matteobisardi/KitMSA

max_issues_repo_name: matteobisardi/KitMSA

max_forks_repo_name: matteobisardi/KitMSA


text:
testfiles = [
    "moiety.jl"
    "ullmann.jl"
    "findreplace.jl"
    "examples.jl"
]

@assert VERSION.major == 1
@assert VERSION.minor ≥ 6

using Test, Documenter, PoreMatMod
PoreMatMod.banner()

for testfile ∈ testfiles
    @info "Running test/$testfile"
    @time include(testfile)
end

if VERSION.minor ≥ 7
    @time doctest(PoreMatMod)
end

@info "Done."

ext: jl


text:
module TestSetMacro

module Clone
    import Accessors: setmacro, opticmacro, modifymacro
    macro optic(ex)
        opticmacro(identity, ex)
    end
    macro set(ex)
        setmacro(identity, ex)
    end
    macro modify(f, ex)
        modifymacro(identity, f, ex)
    end
end#module Clone

using Accessors: Accessors
using Test
using .Clone: Clone

using StaticArrays: @SMatrix
using StaticNumbers

@testset "setmacro, opticmacro isolation" begin

    # test that no symbols like `IndexLens` are needed:
    Clone.@optic(_                                         )
    Clone.@optic(_.a                                       )
    Clone.@optic(_[1]                                      )
    Clone.@optic(first(_)                                  )
    Clone.@optic(_[end]                                    )
    Clone.@optic(_[static(1)]                              )
    Clone.@optic(_.a[1][end, end-2].b[static(1), static(1)])

    @test Accessors.@optic(_.a) === Clone.@optic(_.a)
    @test Accessors.@optic(_.a.b) === Clone.@optic(_.a.b)
    @test Accessors.@optic(_.a.b[1,2]) === Clone.@optic(_.a.b[1,2])

    o = (a=1, b=2)
    @test Clone.@set(o.a = 2) === Accessors.@set(o.a = 2)
    @test Clone.@set(o.a += 2) === Accessors.@set(o.a += 2)

    @test Clone.@modify(x -> x+1, o.a) === Accessors.@modify(x -> x+1, o.a)

    m = @SMatrix [0 0; 0 0]
    m2 = Clone.@set m[end-1, end] = 1
    @test m2 === @SMatrix [0 1; 0 0]
    m3 = Clone.@set(first(m) = 1)
    @test m3 === @SMatrix[1 0; 0 0]
end

end#module


text:
function Base.rand(::Type{SurfaceRangeEnvelope}, presences::T) where {T <: SimpleSDMLayer}
    @assert SimpleSDMLayers._inner_type(presences) <: Bool
    n_occ = sum(presences)
    iszero(n_occ) && throw(ArgumentError("The presences layer is empty"))
    _msk = mask(presences, presences)
    acceptable_cells = setdiff(keys(presences), keys(_msk))
    rlat = extrema([p[2] for p in keys(_msk)])
    rlon = extrema([p[1] for p in keys(_msk)])
    filter!(p -> rlat[1] <= p[2] <= rlat[2], acceptable_cells)
    filter!(p -> rlon[1] <= p[1] <= rlon[2], acceptable_cells)
    if length(acceptable_cells) < n_occ
        @warn "There are fewer acceptable pseudo-absences than occurrences"
    end
    pa_places = sample(acceptable_cells, min(length(acceptable_cells), n_occ), replace=false)
    pa = similar(presences, Bool)
    pa[pa_places] = fill(true, length(pa_places))
    return pa
end


text:
#

function get_singlepage(documentpath; pagesize=300, pagetoken::String="")
    @show "Firestore.get_singlepage", documentpath, pagesize, pagetoken
    base_url = FIRESTORE_URL
    headers = Dict(auth_header())
    project_id = projectid()
    firestore_path = "/projects/$project_id/databases/(default)/documents"
    url = "$base_url$firestore_path/$documentpath"
    query = Dict{String,Any}("pageSize" => pagesize)
    if !isempty(pagetoken)
        query["pageToken"] = pagetoken
    end
    http_response = HTTP.get(url, headers; query=query)
    response_string = String(http_response.body)
    JSON.parse(response_string)
end

function get_collection(documentpath; pagesize=300)
    @show "Firestore.get_collection", documentpath, pagesize
    docs = []
    pagetoken = ""
    while true
        response = get_singlepage(
            documentpath,
            pagesize=pagesize,
            pagetoken=pagetoken,
        )
        isempty(response) && return docs
        push!(docs, response["documents"]...)
        !haskey(response, "nextPageToken") && break
        pagetoken = response["nextPageToken"]
    end
    docs
end

get_document(documentpath) = get_singlepage(documentpath)

function get(path...)
    @show "Firestore.get", path
    full_path = join(path, '/')
    length(path) % 2 == 0 && return convert_map(get_document(full_path))
    map(convert_map, get_collection(full_path))
end

"""
get_request(url::String)

# Get request which is general in nature and gets any document or collection

# Example
```julia
init("[FIREBASE_ADMIN_SDK].json")
res = get_request("/firebase_test/firebase_get") # document get
res = get_request("/firebase_test/firebase_get/firebase_get_collection") # collection fetch
```
"""
function get_request(url::String)
    header = Dict(auth_header())
    project_id = projectid()
    url = "$(FIRESTORE_URL)/projects/$project_id/databases/(default)/documents$url"
    println(url)
    http_response = HTTP.get(url, header)
    response_string = String(http_response.body)
    JSON.parse(response_string)
end

body = Dict(
    "newTransaction" => Dict(),
    # "transaction" => "",
    # "readTime" => "",
    "documents" => [
    "projects/fir-jl-457eb/databases/(default)/documents/firebase_test/firebase_get",
    "projects/fir-jl-457eb/databases/(default)/documents/firebase_test/firebase_get?collectionID=firebase_get_collection",
    "projects/fir-jl-457eb/databases/(default)/documents/firebase_test/firebase_get/firebase_get_collection/test",
    ],
    "mask" => Dict(
        "fieldPaths" => [
            "first1",
        ]),

)

# mask gets particular fields looks like
# weird field in mask is just ignored and nothing happens
"""
firestore_batchget(batch=body)

Gets multiple documents
Documents returned by this method are not guaranteed
to be returned in the same order that they were requested.

# Example
```julia
init("[FIREBASE_ADMIN_SDK].json")
res = firestore_batchget()
```
"""
function firestore_batchget(batch=body)
    header = Dict(auth_header())
    project_id = projectid()
    println(batch)
    body = JSON.json(batch)
    println("Batch:", body)
    final_url = "$(FIRESTORE_URL)/projects/$project_id/databases/(default)/documents:batchGet"
    http_response =  HTTP.request(
        "POST",
        final_url,
        header,
        body,
    )
    response_string = String(http_response.body)
    JSON.parse(response_string)
end

body = Dict(
    "writes" => [
        Dict(
            "update" => Dict(
                "name" => "projects/fir-jl-457eb/databases/(default)/documents/firebase_test/firebase_get",
                "fields" => Dict(
                    "first" => Dict(
                        "stringValue" => "ashwani",
                        ),
                        "last" => Dict(
                            "stringValue" => "rathee",
                            )
                )
            ),
            "updateMask" => Dict(
                "fieldPaths" => [
                    "first",
                    "last",
                ]
            )
        )
    ],
    "labels" => Dict(
        string => "ABC"
    )
)

# study types of writes from here https://firebase.google.com/docs/firestore/reference/rest/v1/Write
# study the effect of labelss, and update masks
# Labels associated with this batch write.
# mask gets particular fields looks like
# weird field in mask is just ignored and nothing happens
"""
firestore_batchwrite(batch=body)

Applies a batch of write operations.
The documents.batchWrite method does not apply the write operations atomically
and can apply them out of order. Method does not allow more than one write per document.
Each write succeeds or fails independently. See the BatchWriteResponse for the success
status of each write.

If you require an atomically applied set of writes, use documents.commit instead.

# Example
```julia
init("[FIREBASE_ADMIN_SDK].json")
res = firestore_batchwrite()
```
"""
function firestore_batchwrite(batch=body)
    header = Dict(auth_header())
    project_id = projectid()
    println(batch)
    body = JSON.json(batch)
    println("Batch:", body)
    final_url = "$(FIRESTORE_URL)/projects/$project_id/databases/(default)/documents:batchWrite"
    http_response =  HTTP.request(
        "POST",
        final_url,
        header,
        body,
    )
    response_string = String(http_response.body)
    JSON.parse(response_string)
end

body = Dict(
    "options" => Dict(

    )
)

# study transaction options and what is transactions used for here
"""
firestore_beginTransaction()

Starts a new transaction.
"options"
The options for the transaction. Defaults to a read-write transaction.

# Example
```julia
init("[FIREBASE_ADMIN_SDK].json")
firestore_beginTransaction()
```
"""
function firestore_beginTransaction()
    header = Dict(auth_header())
    project_id = projectid()
    # println(batch)
    # body = JSON.json(batch)
    # println("Batch:", body)
    final_url = "$(FIRESTORE_URL)/projects/$project_id/databases/(default)/documents:beginTransaction"
    http_response =  HTTP.request(
        "POST",
        final_url,
        header,
        # body,
    )
    response_string = String(http_response.body)
    JSON.parse(response_string)
end

body = Dict(
    "writes" => [
        Dict(
            "update" => Dict(
                "name" => "projects/fir-jl-457eb/databases/(default)/documents/firebase_test/firebase_get",
                "fields" => Dict(
                    "first" => Dict(
                        "stringValue" => "ashwani",
                        ),
                        "last" => Dict(
                            "stringValue" => "rathee",
                            )
                )
            ),
            "updateMask" => Dict(
                "fieldPaths" => [
                    "first",
                    "last",
                ]
            )
        )
    ],

    "transaction" => "EZTktxl/jJ9GIlkA3ModskgtczlV2HUqlqGQbh6U8YClsJvU07MTX5/8q0Fhd3hAvi0aSvm4NLF2H8FNiEzspxJLVM4NOZVDsClQgpuLqwkA9uFrQM81VXxzYzGLLvU0hssZZg=="
)

# commits a particular transaction that was started or is happening
# this one is particularly interesting
"""
firestore_commit(transactionid="",body =body)

Commits a transaction, while optionally updating documents.

# Example
```julia
init("[FIREBASE_ADMIN_SDK].json")
data = firestore_beginTransaction()
firestore_commit(data["transaction"])
```
"""
function firestore_commit(transactionid="", body=body)
    header = Dict(auth_header())
    project_id = projectid()
    body["transaction"] = transactionid
    body = JSON.json(body)
    println("Body:", body)
    final_url = "$(FIRESTORE_URL)/projects/$project_id/databases/(default)/documents:commit"
    http_response =  HTTP.request(
        "POST",
        final_url,
        header,
        body,
    )
    response_string = String(http_response.body)
    JSON.parse(response_string)
end

    # Create Document
# body = """{
#         "fields": {
#           "test": {
#             "stringValue": "1"
#           }
#         }
#       }
#       """
body = Dict(
    "fields" => Dict(
        "test" => Dict(
            "stringValue" => "1"
        )
    )
)
# above :The request body contains an instance of Document.

# how two add field to the
# dpcument id is supposed to be a query parameter as of now it's hardcoded
"""
firestore_createdoc(url, body)

# Example

```julia
init("[FIREBASE_ADMIN_SDK].json")
firestore_createdoc("/firebase_test/firebase_post/two")
firestore_createdoc("/firebase_test/firebase_post/three";query ="hello") # with document id name in query
```
"""
function firestore_createdoc(url, body=body;query=query)
    header = Dict(auth_header())
    project_id = projectid()
    body = JSON.json(body)
    println("Body:", body)
    query = Dict{String,Any}("documentId" => query)
    # mask part doesn't work yet
        # query = JSON.json(query)
        # println(query)
    final_url = "$(FIRESTORE_URL)/projects/$project_id/databases/(default)/documents$url"
    http_response = HTTP.post(
        final_url,
        header,
        body;query=query
    )
    response_string = String(http_response.body)
    JSON.parse(response_string)
end

# implement query properly for documentid and mask in here

"""
firestore_deletedoc(url)

Deletes a document

# Example
```julia
init("[FIREBASE_ADMIN_SDK].json")
firestore_createdoc("/firebase_test/firebase_post/three";query ="hello") # creates the document and then we delete it
firestore_deletedoc("/firebase_test/firebase_post/three/hello")
```
"""
function firestore_deletedoc(url)
    header = Dict(auth_header())
    project_id = projectid()
    final_url = "$(FIRESTORE_URL)/projects/$project_id/databases/(default)/documents$url?currentDocument.exists=true"
    http_response = HTTP.delete(
        final_url,
        header,
    )
    response_string = String(http_response.body)
    JSON.parse(response_string)
end

"""
firestore_getdoclist(url)

Lists documents

# Example
```julia
init("[FIREBASE_ADMIN_SDK].json")
firestore_getdoclist("/firebase_test")
```
"""
function firestore_getdoclist(url)
    header = Dict(auth_header())
    project_id = projectid()
    final_url = "$(FIRESTORE_URL)/projects/$project_id/databases/(default)/documents$url"
    http_response =  HTTP.get(
        final_url,
        header,
    )
    response_string = String(http_response.body)
    JSON.parse(response_string)
end

"""
firestore_getcollectionids(url;query="")

Lists all the collection IDs underneath a document.

# Example
```julia
init("[FIREBASE_ADMIN_SDK].json")
firestore_getcollectionids("/firebase_test/firebase_post")
```
"""
function firestore_getcollectionids(url;query="")
    header = Dict(auth_header())
    project_id = projectid()
    final_url = "$(FIRESTORE_URL)/projects/$project_id/databases/(default)/documents$url:listCollectionIds"
    pagesize = 300
    pagetoken = ""
    if !isempty(pagetoken)
        query["pageToken"] = pagetoken
    end
    query = Dict{String,Any}("pageSize" => pagesize, "pageToken" => pagetoken)

    # this query content is supposed to be request body
    http_response = HTTP.post(
        final_url,
        header;
        query=query,
    )
    response_string = String(http_response.body)
    JSON.parse(response_string)
end

# Doesn't work really
body = Dict(
    "partitionCount" => "1",
    "pageToken" => "",
    "pageSize" => 300,
    "structuredQuery" => Dict(
        "where" => Dict(
            "fieldFilter" => Dict(
                "field" => Dict("fieldPath" => "test"),
                "op" => "EQUAL",
                "value" => Dict("integerValue" => 1)
            )
        ),
        "from" => [Dict("collectionId" => "new1")]
    )
)

"""
firestore_partitionquery(url)

Partitions a query by returning partition cursors that can be used
to run the query in parallel. The returned partition cursors are
split points that can be used by documents.runQuery as starting/end points for the query results.
"""
function firestore_partitionquery(url, body=body)
    header = Dict(auth_header())
    project_id = projectid()
    final_url = "$(FIRESTORE_URL)/projects/$project_id/databases/(default)/documents$url:partitionQuery"
    body = JSON.json(body)
    println("Body:", body)
    http_response = HTTP.post(
        final_url,
        header,
        body
        ;
    )
    response_string = String(http_response.body)
    JSON.parse(response_string)
end

body = Dict(
    "fields" => Dict(
        "test1" => Dict(
            "stringValue" => "1"
        )
    )
)

# need to implement query parameters like updatemask ,mask currentDocument

"""
firestore_patch()

Updates or inserts a document.

# Example

```julia
init("[FIREBASE_ADMIN_SDK].json")
firestore_patch("/firebase_test/firebase_patch")
```
"""
function firestore_patch(url, body=body)
    header = Dict(auth_header())
    project_id = projectid()
    final_url = "$(FIRESTORE_URL)/projects/$project_id/databases/(default)/documents$url"
    body = JSON.json(body)
    println("Body:", body)
    # query = Dict{String,Any}("documentId" => query)
    # mask part doesn't work yet
        # query = JSON.json(query)
        # println(query)
    http_response = HTTP.patch(
        final_url,
        header,
        body
        ;
    )
    response_string = String(http_response.body)
    JSON.parse(response_string)
end

"""
firestore_rollback(transactionid,body=body)

Rolls back a transaction.

# Example

```julia
init("[FIREBASE_ADMIN_SDK].json")
data = firestore_beginTransaction()
# firestore_commit(data["transaction"])
firestore_rollback(data["transaction"])
```

# Output
If successful, the response body will be empty.

"""

body = Dict("transaction" => "")
function firestore_rollback(transactionid, body=body)
    header = Dict(auth_header())
    project_id = projectid()
    final_url = "$(FIRESTORE_URL)/projects/$project_id/databases/(default)/documents:rollback"
    body["transaction"] = transactionid
    body = JSON.json(body)
    println("Body:", body)

    http_response = HTTP.post(
        final_url,
        header,
        body
        ;
    )
    response_string = String(http_response.body)
    JSON.parse(response_string)

end

body = Dict(
    "newTransaction" => Dict(),
    "structuredQuery" => Dict(
        "select" => Dict(
         "fields" => [
          Dict(
           "fieldPath" => "test"
          ),
         ]
        )
    )
)

body = Dict(
    "newTransaction" => Dict(),
    "structuredQuery" => Dict(
        "where" => Dict(
            "fieldFilter" => Dict(
                "field" => Dict("fieldPath" => "test"),
                "op" => "EQUAL",
                "value" => Dict("stringValue" => "1")
            )
        ),
        "from" => [Dict("collectionId" => "new")]
    )
)
# very specific as shown above is successfully run and it works by choosing 1/3 or 2/3 as the variables are changed

"""
firestore_runquery()

Runs a query.

# Example
```julia
init("[FIREBASE_ADMIN_SDK].json")
firestore_runquery("/firebase_test/firebase_post")
```
"""
function firestore_runquery(url, body=body)
    header = Dict(auth_header())
    project_id = projectid()
    final_url = "$(FIRESTORE_URL)/projects/$project_id/databases/(default)/documents$url:runQuery"
    body = JSON.json(body)
    println("Body:", body)
    println("FINAL_URL:", final_url)
    http_response = HTTP.post(
        final_url,
        header,
        body
    )
    response_string = String(http_response.body)
    JSON.parse(response_string)
end


text:
immutable Net
    layers::Vector{Layer}

    function Net(layers::Vector{Layer})
        for i in 1:length(layers) - 1
            @assert inputSize(layers[i+1]) == outputSize(layers[i])
        end
        new(layers)
    end
end

gradientStates(net::Net) = map(gradientStates, net.layers)

function forward(net::Net, inputActivation::Vol)
    act = forward(net.layers[1], inputActivation)
    for i in 2:length(net.layers)
        act = forward(net.layers[i], act)
    end
    invariants(net)
    return act
end

function backward(net::Net, class::Int)
    invariants(net)
    numLayers = length(net.layers)
    loss = backward(net.layers[numLayers], class)
    for i = numLayers-1:-1:1
        backward(net.layers[i])
    end
    invariants(net)
    return loss
end

function invariants(net::Net)
    for i in 1:length(net.layers) - 1
        @assert is(net.layers[i].outputActivation, net.layers[i+1].inputActivation)
        @assert is(net.layers[i].outputActivation.w, net.layers[i+1].inputActivation.w)
        @assert is(net.layers[i].outputActivation.dw, net.layers[i+1].inputActivation.dw)
        if i + 2 <= length(net.layers)
            @assert !is(net.layers[i].outputActivation, net.layers[i+2].inputActivation)
            @assert !is(net.layers[i].outputActivation.w, net.layers[i+2].inputActivation.w)
            @assert !is(net.layers[i].outputActivation.dw, net.layers[i+2].inputActivation.dw)
        end
    end
end

function train(net::Net, inputActivation::Vol, class::Int)
    forward(net, inputActivation)
    result = backward(net, class)
    invariants(net)
    return result
end

function predict(net::Net, inputActivation::Vol)
    forward(net, inputActivation)
    _, index = findmax(last(net.layers).outputActivation.w)
    invariants(net)
    return index
end


text:
#!/usr/bin/julia

# Trizen
# 17 September 2016
# https://github.com/trizen

# A decently efficient algorithm for computing the results of the Kempner/Smarandache function.

# See also: https://projecteuler.net/problem=549
#           https://en.wikipedia.org/wiki/Kempner_function
#           http://mathworld.wolfram.com/SmarandacheFunction.html

# ∑S(i) for 2 ≤ i ≤ 10^2 == 2012
# ∑S(i) for 2 ≤ i ≤ 10^6 == 64938007616
# ∑S(i) for 2 ≤ i ≤ 10^8 == 476001479068717

using Primes

function smarandache(n::Int64, cache)

    isprime(n) && return n

    f = factor(n)

    count = 0
    distinct = true
    for v in values(f)
        count += v
        if (distinct && v != 1)
            distinct = false
        end
    end

    distinct && return maximum(keys(f))

    if (length(f) == 1)

        k = collect(keys(f))[1]

        (count <= k) && return k*count

        if haskey(cache, n)
            return cache[n]
        end

        max = k*count
        ff  = factorial(BigInt(max - k))

        while (ff % n == 0)
            max -= k
            ff /= max
        end

        cache[n] = max
        return max
    end

    arr = Int64[]

    for (k,v) in f
        push!(arr, v == 1 ? k : smarandache(k^v, cache))
    end

    maximum(arr)
end

#
## Tests
#

function test()
    cache = Dict{Int64, Int64}()

    limit = 10^2
    sumS = 0

    for k in 2:limit
        sumS += smarandache(k, cache)
    end

    println("∑S(i) for 2 ≤ i ≤ $limit == $sumS")

    if (limit == 100 && sumS != 2012)
        warn("However that is incorrect! (expected: 2012)")
    end
end

test()

max_stars_repo_name: trizen/julia-scripts

max_issues_repo_name: trizen/julia-scripts

max_forks_repo_name: trizen/julia-scripts


text:
using EstimationStatistics
using Test

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


text:
module StoppingInterface

using LinearAlgebra, Stopping

using SolverCore
include("stopping_to_jso_stats.jl")
export status_stopping_to_stats, stopping_to_stats

using JSOSolvers, NLPModels
include("stopping_to_jsosolvers.jl")

using NLPModelsIpopt
include("stopping_to_nlpmodels_ipopt.jl")

using Requires
is_knitro_installed = false
include("stopping_to_nlpmodels_knitro.jl")

end # module


text:
function derivative(f, x::Union{Number,AbstractArray{<:Number}},
                    integrator)
    local d
    alg = unwrap_alg(integrator, true)
    if get_current_alg_autodiff(integrator.alg, integrator.cache)
      d = ForwardDiff.derivative(f, x)
    else
      d = DiffEqDiffTools.finite_difference_gradient(f, x, alg.diff_type, eltype(x), Val{false})
    end
    d
end

function derivative!(df::AbstractArray{<:Number}, f, x::Union{Number,AbstractArray{<:Number}}, fx::AbstractArray{<:Number}, integrator::DEIntegrator)
    if alg_autodiff(integrator.alg)
        ForwardDiff.derivative!(df, f, fx, x)
    else
        RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
        DiffEqDiffTools.finite_difference!(df, f, x, integrator.alg.diff_type, RealOrComplex, fx)
    end
    nothing
end

function jacobian(f, x,
                  integrator::DiffEqBase.DEIntegrator)
    local J
    alg = unwrap_alg(integrator, true)
    isarray = typeof(x) <: AbstractArray
    if get_current_alg_autodiff(integrator.alg, integrator.cache)
      if isarray
        J = ForwardDiff.jacobian(f,x)
      else
        J = ForwardDiff.derivative(f,x)
      end
    else
      if isarray
        J = DiffEqDiffTools.finite_difference_jacobian(f, x, alg.diff_type, eltype(x), Val{false})
      else
        J = DiffEqDiffTools.finite_difference_derivative(f, x, alg.diff_type, eltype(x))
      end
    end
    J
end

function jacobian!(J::AbstractMatrix{<:Number}, f, x::AbstractArray{<:Number}, fx::AbstractArray{<:Number}, integrator::DEIntegrator, jac_config)
    if alg_autodiff(integrator.alg)
      ForwardDiff.jacobian!(J, f, fx, x, jac_config)
    else
      DiffEqDiffTools.finite_difference_jacobian!(J, f, x, jac_config)
    end
    nothing
end

function build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  if !has_jac(f)
    if alg_autodiff(alg)
      jac_config = ForwardDiff.JacobianConfig(uf,du1,uprev,ForwardDiff.Chunk{determine_chunksize(u,alg)}())
    else
      if alg.diff_type != Val{:complex}
        jac_config = DiffEqDiffTools.JacobianCache(tmp,du1,du2,alg.diff_type)
      else
        jac_config = DiffEqDiffTools.JacobianCache(Complex{eltype(tmp)}.(tmp),Complex{eltype(du1)}.(du1),nothing,alg.diff_type)
      end
    end
  else
    jac_config = nothing
  end
  jac_config
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

# based on code from https://www.ilikebigbits.com/2015_03_04_plane_from_points.html
"""
    plane_from_points(points::Vector{SVector{3, Float64}}) ->  centroid, normal

Estimate the best fitting plane for a set of points in 3D.

See utilities.jl for a simpler and possibly more numerically stable version.
"""
function plane_from_points2(points::Vector{SVector{3, Float64}})
    if length(points) < 3
        return nothing      # At least three points required
    end

    sum = Base.sum(points)
    centroid = Vec3(sum * (1.0 / (length(points))))

    # Calc full 3x3 covariance matrix, excluding symmetries:
    xx = 0.0
    xy = 0.0
    xz = 0.0
    yy = 0.0
    yz = 0.0
    zz = 0.0

    for p in points
        r = p - centroid;
        xx += r[1] * r[1]
        xy += r[1] * r[2]
        xz += r[1] * r[3]
        yy += r[2] * r[2]
        yz += r[2] * r[3]
        zz += r[3] * r[3]
    end

    det_x = yy*zz - yz*yz
    det_y = xx*zz - xz*xz
    det_z = xx*yy - xy*xy

    det_max = max(det_x, det_y, det_z)
    if det_max <= 0.0
        return nothing      # The points don't span a plane
    end

    # Pick path with best conditioning:
    if      (det_max == det_x)
        dir = Vec3(det_x,xz*yz - xy*zz,xy*yz - xz*yy)
    elseif  (det_max == det_y)
        dir = Vec3(xz*yz - xy*zz, det_y, xy*xz - yz*xx)
    else
        dir = Vec3(xy*yz - xz*yy,xy*xz - yz*xx,det_z)
    end

    return centroid, normalize(dir)
end

function csg_sphere(;radius = 10.0)
    sph = Sphere(radius)
    csg = leaf(sph, Transform(Vec3(0.0, 0.0, radius)))
    csg = csg(identitytransform())
    return csg
end

function csg_cylinder(;radius = 10.0, height = 50.0, added_rotation = rotationY(π/2.0))
    cyl = Cylinder(radius, height)
    csg = leaf(cyl, Transform(Vec3(0.0, 0.0, radius)) * added_rotation)
    csg = csg(identitytransform())
    return csg
end

function csg_plane()
    plane = Plane(unitZ3() * -1, zero(Vec3))
    csg = leaf(plane, identitytransform())
    csg = csg(identitytransform())
    return csg
end


text:
using Random

## Binary string fitnesses

function onemax(instance:: Dict{String, }, point::BitArray{1})
    sum(point)
end

function zeromax(instance:: Dict{String, }, point::BitArray{1})
    length(point) - sum(point)
end

function fractionmax(instance:: Dict{String, }, fraction::Float64)
    return function _F(instance, point)
        ## length(candidate)*fraction
        lastZero = Int(instance["size"]*fraction)
        ZeroMaxFitness(instance, point[1:lastZero]) + OneMaxFitness(instance, point[lastZero+1:end])
    end
end

max_stars_repo_name: ahrajabi/DiscO

max_issues_repo_name: ahrajabi/DiscO

max_forks_repo_name: ahrajabi/DiscO


text:
@testset "interpolations" begin

for interpolation in (Lagrange{1, RefCube, 1}(),
                      Lagrange{1, RefCube, 2}(),
                      Lagrange{2, RefCube, 1}(),
                      Lagrange{2, RefCube, 2}(),
                      Lagrange{2, RefTetrahedron, 1}(),
                      Lagrange{2, RefTetrahedron, 2}(),
                      Lagrange{3, RefCube, 1}(),
                      Serendipity{2, RefCube, 2}(),
                      Serendipity{3, RefCube, 2}(),
                      Lagrange{3, RefTetrahedron, 1}(),
                      Lagrange{3, RefTetrahedron, 2}(),)

    # Test of utility functions
    ndim = Ferrite.getdim(interpolation)
    r_shape = Ferrite.getrefshape(interpolation)
    func_order = Ferrite.getorder(interpolation)
    @test typeof(interpolation) <: Interpolation{ndim,r_shape,func_order}
    @test typeof(Ferrite.getlowerdim(interpolation)) <: Interpolation{ndim-1}
    @test typeof(Ferrite.getlowerorder(interpolation)) <: Interpolation{ndim,r_shape,func_order-1}

    n_basefuncs = getnbasefunctions(interpolation)
    x = rand(Tensor{1, ndim})
    f = (x) -> Ferrite.value(interpolation, Tensor{1, ndim}(x))
    @test vec(ForwardDiff.jacobian(f, Array(x))') ≈
           reinterpret(Float64, Ferrite.derivative(interpolation, x))
    @test sum(Ferrite.value(interpolation, x)) ≈ 1.0

    coords = Ferrite.reference_coordinates(interpolation)
    for node in 1:n_basefuncs
        N_node = Ferrite.value(interpolation, coords[node])
        for k in 1:node
            if k == node
                @test N_node[k] ≈ 1.0
            else
                @test N_node[k] ≈ 0.0
            end
        end
    end
end

end


text:
"""
    GaussLik <: Likelihood

Gaussian, a.k.a. Normal, likelihood
```math
p(y | f, σ) = 1 / √(2πσ²) \\exp(-(x - f)²/(2σ²)),
```
where standard deviation ``σ`` is a non-fixed hyperparameter and ``f`` is the latent
Gaussian process.
"""
mutable struct GaussLik <: Likelihood
    "Standard deviation"
    σ::Float64
    "Priors for likelihood parameters"
    priors::Array

    """
        GaussLik(lσ::Float64)

    Create `GaussLik` with standard deviation `exp(lσ)`.
    """
    GaussLik(lσ::Float64) = new(exp(lσ), [])
end

# log of probability density
function log_dens(gauss::GaussLik, f::VecF64, y::VecF64)
    return (-0.5 * log(2 * pi) - log(gauss.σ)) .- 0.5 * ((y - f) / gauss.σ).^2
end

# derivative of log pdf wrt latent function
function dlog_dens_df(gauss::GaussLik, f::VecF64, y::VecF64)
    return [(yi-fi)/gauss.σ^2 for (fi,yi) in zip(f,y)]
end

# derivative of log pdf wrt to parameters
function dlog_dens_dθ(gauss::GaussLik, f::VecF64, y::VecF64)
    return gauss.σ*[-1/gauss.σ + 1/gauss.σ^3*(yi-fi).^2 for (fi,yi) in zip(f,y)]
end

#mean and variance under likelihood
mean_lik(gauss::GaussLik, f::VecF64) = f
var_lik(gauss::GaussLik, f::VecF64) = ones(length(f))*gauss.σ^2

function set_params!(gauss::GaussLik, hyp::VecF64)
    length(hyp) == 1 || throw(ArgumentError("Gaussian/Normal likelihood has only one free parameter"))
    gauss.σ = exp(hyp[])
end

get_params(gauss::GaussLik) = Float64[log(gauss.σ)]
num_params(gauss::GaussLik) = 1

#Computes the predictive mean and variance
function predict_obs(gauss::GaussLik, fmean::VecF64, fvar::VecF64)
    return fmean, fvar + gauss.σ^2
end


text:
if abspath(PROGRAM_FILE) == @__FILE__
    (a, b, c) = [1, 2, 3]
    @assert(a == 1)
    @assert(b == 2)
    @assert(c == 3)
    pairs = [(([2, 3, 4], 5, 6), ([7, 8, 9], 10, 11))]
    for (((x1, y1, z1), v1, m1), ((x2, y2, z2), v2, m2)) in pairs
        @assert(x1 == 2)
        @assert(y1 == 3)
        @assert(z1 == 4)
        @assert(v1 == 5)
        @assert(m1 == 6)
        @assert(x2 == 7)
        @assert(y2 == 8)
        @assert(z2 == 9)
        @assert(v2 == 10)
        @assert(m2 == 11)
    end
end


text:
export Unit,
       Choice,
       𝕀,
       ℝ,
       BinaryChoice
export PrimitiveParam, Param

# # Primitives
# These are primitive parameters.

abstract type PrimitiveParam end
(pp::PrimitiveParam)(id, ω) = ω[Member(id, p)]

"Nondeterministic choice of true or false"
struct BinaryChoice{T <: Integer} <: PrimitiveParam end
Base.eltype(::Type{BinaryChoice{T}}) where T = T

"A Real-valued variable of type `T` in the unit interval: [0, 1]"
struct Unit{T} <: PrimitiveParam end
Base.eltype(::Type{Unit{T}}) where T = T
const 𝕀 = Unit

"The Real numbers"
struct Reals{T} <: PrimitiveParam end
Base.eltype(::Type{Reals{T}}) where T = T
const ℝ = Reals

# #

# # Families of parameters
# # Non primtiives
"Parameter family"
abstract type Param end

"Parametric choice of values in collection"
struct Choice{T} <: Param
  of::T
end
function (c::Choice)(id, ω)
  # isempty(c.of) && throw(ArgumentError("range must be non-empty"))
  c.of[(Real{Int}()(id, ω) % length(c.of)) + 1 ]
end

"One of many other parametric choices"
struct Disjunction{XS} <: Param
  xs::XS
end

function (c::Disjunction)(id, ω)
  i = Choice(1:length(x.xs))(id, ω)
  c.xs[i](id, ω)
end

"Interval `[a, b]`"
struct Interval{T} <: Param
  a::T
  b::T
end

struct Permutation
  n
end

struct Shuffle
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
### A Pluto.jl notebook ###
# v0.14.4

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 42b4e6fe-a897-11eb-2676-f7fd96f35a22
begin
	using NeidSolarScripts
	using RvSpectMLBase
	using EchelleInstruments
	using EchelleCCFs
	using RvSpectML
	using JLD2, FileIO
	using PlutoUI
	using Plots, ColorSchemes
	using Statistics, NaNMath
	using Polynomials
end

# ╔═╡ 5313b674-a897-11eb-3820-218f26a14d4d
begin
	ccf_path = ".."
	ccf_fn = "daily_ccfs.jld2"
	@assert isfile(joinpath(ccf_path,ccf_fn))
	@assert filesize(joinpath(ccf_path,ccf_fn)) > 0
end

# ╔═╡ 37a6c623-e901-4f8b-8989-835958471549
begin
	data = load(joinpath(ccf_path,ccf_fn))
	num_obs = size(data["order_ccfs"],3)
	num_order_idx = size(data["order_ccfs"],2)
	num_v = size(data["order_ccfs"],1)
	num_pix = size(data["mean_clean_flux"],1)
	num_pix_continuum = size(data["mean_clean_flux_continuum_normalized"],1)
	date_str = match(r"neidL1_(\d+)T\d+\.fits",data["manifest"].Filename[1]).captures[1]
	size(data["order_ccfs"])
end

# ╔═╡ d1fe0408-cccd-45c9-a5ef-c13a42386e44
begin
	(vmin_default, vmax_default) = extrema(data["v_grid"])
	order_ccfs_norm = data["order_ccfs"]./maximum(data["order_ccfs"],dims=1)
	order_ccfs_mean = mean(order_ccfs_norm,dims=3)
	ccf_order_depths = 1.0 .- minimum(order_ccfs_norm,dims=1)./maximum(order_ccfs_norm,dims=1)
	ccf_mean_order_depth = vec(mean(ccf_order_depths,dims=2))
end;

# ╔═╡ 609e3d45-e1a5-4060-97a8-5b668f8137c1
md"""
### Pixel ranges for figures plotting mean flux.
Colum  (min, max) $(@bind xmin NumberField(1:num_pix; default=1))  $(@bind xmax NumberField(1:num_pix; default=num_pix))
Show uncertainties?
Raw $(@bind plot_raw_errs CheckBox(default=false))
Continuum $(@bind plot_continum_errs CheckBox(default=false))
"""

# ╔═╡ 28aced5a-3100-4523-9f1a-f7b14016791b
md"""
#### Order index & Observation for plotting single CCF.

Order index  $(@bind order_idx NumberField(1:num_order_idx; default=floor(Int64,median(1:num_order_idx))))
Obs index  $(@bind obs_idx NumberField(1:num_obs; default=1))
"""

# ╔═╡ abbc4880-fedd-4bdd-a1db-b94e8a4d5654
md" vmin $(@bind vmin_kms NumberField(-100:100; default=-30))
vmax $(@bind vmax_kms NumberField(-100:100; default=30))
Ready? $(@bind ready_to_plot CheckBox(default=false))
"

# ╔═╡ 6c469209-c7e0-41f2-97e3-db587b894e56
begin
	vmin = vmin_kms *1000;
	vmax = vmax_kms *1000;
end

# ╔═╡ ab0fc462-a897-11eb-1727-e71788522c02
if !all(isnan.(data["order_ccfs"][:,order_idx,obs_idx]))
  plt1 = plot(data["v_grid"],data["order_ccfs"][:,order_idx,obs_idx], ms=1.2,legend=:none)
  xlims!(plt1,vmin,vmax)
	title!("CCF Order Idx " * string(data["orders_to_use"][order_idx]) * " Obs " * string(obs_idx))
  #ylims!(plt1,extrema(f_filtered[λmin.<=lambda.<=λmax]))
  plt1
end

# ╔═╡ c403d4ff-3035-4bfc-9b9e-3f33ed0a6b7a
md"""
### Ranges for figures plotting multiple CCFs.
Order index (min, max, step) $(@bind order_idx_lo NumberField(1:num_order_idx; default=1))  $(@bind order_idx_hi NumberField(1:num_order_idx; default=num_order_idx))
$(@bind order_idx_step NumberField(1:num_order_idx; default=1))
"""

# ╔═╡ d84af1b1-1023-4493-8b0f-90c9030be4bd
md"Scale CCF depths? $(@bind scale_order_ccf_depths CheckBox(default=false))  "

# ╔═╡ 7dc1bf77-174b-4614-8aea-15cad25c0335
begin
  plt2 = plot()
  pal2 = palette(:berlin, length(order_idx_lo:order_idx_step:order_idx_hi))
  for (i,ord_idx) in enumerate(order_idx_lo:order_idx_step:order_idx_hi)
		if all(isnan.(order_ccfs_norm[:,ord_idx,obs_idx])) continue end
  		y = scale_order_ccf_depths ? 1. .- (1.0 .-order_ccfs_norm[:,ord_idx,obs_idx])./ccf_order_depths[1,ord_idx,obs_idx] :
		order_ccfs_norm[:,ord_idx,obs_idx]
		#plot!(plt2,data["v_grid"],order_ccfs_norm[:,ord_idx,obs_idx], ms=1.2,color=pal2[i],legend=:none)
		plot!(plt2,data["v_grid"],y, ms=1.2,color=pal2[i],legend=:none)
  end
  xlims!(plt2,vmin,vmax)
  xlabel!(plt2,"v (m/s)")
  ylabel!(plt2,scale_order_ccf_depths ? "Scaled CCFs" : "Normalized CCFs")
  title!(plt2,"CCF Order Idxs " * string(data["orders_to_use"][order_idx_lo]) * "-" * string(data["orders_to_use"][order_idx_hi]) * " Obs " * string(obs_idx))
  #ylims!(plt2,extrema(f_filtered[λmin.<=lambda.<=λmax]))
  plt2
end

# ╔═╡ ceae0ff5-f650-4dbf-81e5-26864eb1adcc
begin
	local z = scale_order_ccf_depths ? 1. .- (1.0 .-order_ccfs_norm[:,order_idx_lo:order_idx_hi,obs_idx])./reshape(ccf_order_depths[1,order_idx_lo:order_idx_hi,obs_idx],1,length(order_idx_lo:order_idx_hi)) :
		order_ccfs_norm[:,order_idx_lo:order_idx_hi,obs_idx]
	plt5 = heatmap(data["v_grid"],data["orders_to_use"][order_idx_lo:order_idx_hi],z')
	xlims!(plt5,vmin,vmax)
end

# ╔═╡ df7571c9-8b3f-4503-8bbe-42a77f0f9c0e
md"""Obs index range (min, max, step) $(@bind obs_idx_lo NumberField(1:num_obs; default=1))
$(@bind obs_idx_hi NumberField(1:num_obs; default=num_obs))
$(@bind obs_idx_step NumberField(1:num_obs; default=1))
"""

# ╔═╡ 468fe2c8-7c4c-45f2-9577-a2cecba64d2d
begin
  pal3 = palette(:berlin, length(obs_idx_lo:obs_idx_step:obs_idx_hi))
  plt3 = plot()
  for (i,obs_id) in enumerate(obs_idx_lo:obs_idx_step:obs_idx_hi)
		if all(isnan.(order_ccfs_norm[:,order_idx,i])) continue end
		plot!(plt3,data["v_grid"],order_ccfs_norm[:,order_idx,i].-order_ccfs_mean[:,order_idx], ms=1.2,c=pal3[i],legend=:none)
	end
  xlims!(plt3,vmin,vmax)
  #ylims!(plt2,extrema(f_filtered[λmin.<=lambda.<=λmax]))
  xlabel!(plt3,"v (m/s)")
  ylabel!(plt3,"Normalized CCF - Time ave")
  title!(plt3,"CCF Order Idx " * string(order_idx) * "   Obs " * string(obs_idx_lo) * "-" * string(obs_idx_hi) )
  plt3
end

# ╔═╡ 0e4e2de1-d706-4019-bece-e24367b38d33
begin
	plt4 = heatmap(data["v_grid"],obs_idx_lo:obs_idx_hi,(order_ccfs_norm[:,order_idx,obs_idx_lo:obs_idx_step:obs_idx_hi].-order_ccfs_mean[:,order_idx])')
	xlabel!(plt4,"v (m/s)")
	ylabel!(plt4,"Observation number")
end

# ╔═╡ ed5e02d9-cfba-4ee6-85f2-3678cb1899f9

begin
	if eltype(data["normalization_anchors"][1]) <:AbstractVector
	num_orders_continuum = length(first(data["normalization_anchors"]))
	anchor_density = zeros(num_pix,num_orders_continuum)
	sigma_pix = 3
	for i in 1:num_obs
		for j in 1:num_orders_continuum
			anchors_here = data["normalization_anchors"][i][j]
			n_anchors_here = length(anchors_here)
			for k in 1:n_anchors_here
				anchor_density[:,j] .+= map(l->exp(-0.5*(anchors_here[k]-l)^2/sigma_pix^2),1:num_pix)./(sigma_pix*sqrt(2π))
			end
		end
	end
	end
end;

# ╔═╡ d251eb17-a164-4a36-90f3-d1d5aa87405c
begin
	plt = plot(legend=:bottomleft)
	#xmin = 1800;  xmax = 9500;
	pixlo = searchsortedfirst(1:size(data["mean_clean_flux_continuum_normalized"],1),xmin)
	pixhi = searchsortedlast(1:size(data["mean_clean_flux_continuum_normalized"],1),xmax)
	yscale_raw = NaNMath.maximum(data["mean_clean_flux"][pixlo:pixhi,data["orders_to_use"][order_idx]])
	yscale_sed = NaNMath.maximum(data["mean_clean_flux_sed_normalized"][pixlo:pixhi,data["orders_to_use"][order_idx]])
	yscale_continuum = NaNMath.maximum(data["mean_clean_flux_continuum_normalized"][pixlo:pixhi,data["orders_to_use"][order_idx]])

	plot!(plt,pixlo:pixhi,
		data["mean_clean_flux_sed_normalized"][pixlo:pixhi,data["orders_to_use"][order_idx]]./yscale_sed, color=:green, label="SED Normalized")
	plot!(plt,pixlo:pixhi,
			data["mean_clean_flux_continuum_normalized"][pixlo:pixhi,data["orders_to_use"][order_idx]]./yscale_continuum, color=:red, label="Continuum Normalized")
	if plot_continum_errs
		plot!(plt,pixlo:pixhi,
			data["mean_clean_flux_continuum_normalized"][pixlo:pixhi,data["orders_to_use"][order_idx]]./yscale_continuum.+sqrt.(data["mean_clean_var_continuum_normalized"][pixlo:pixhi,data["orders_to_use"][order_idx]])./yscale_continuum,color=:red, ls=:dot,label=:none)
		plot!(plt,pixlo:pixhi,
			data["mean_clean_flux_continuum_normalized"][pixlo:pixhi,data["orders_to_use"][order_idx]]./yscale_continuum.-sqrt.(data["mean_clean_var_continuum_normalized"][pixlo:pixhi,data["orders_to_use"][order_idx]])./yscale_continuum,color=:red, ls=:dot,label=:none)
	else
	end
	plot!(plt,pixlo:pixhi,data["mean_clean_flux"][pixlo:pixhi,data["orders_to_use"][order_idx]]./yscale_raw, color=:blue, label="Raw")
	if plot_raw_errs
		plot!(plt,pixlo:pixhi,data["mean_clean_flux"][pixlo:pixhi,data["orders_to_use"][order_idx]]./yscale_raw.+sqrt.(data["mean_clean_var"][pixlo:pixhi,data["orders_to_use"][order_idx]])./yscale_raw
			, color=:blue, ls=:dot, label=:none)
		plot!(plt,pixlo:pixhi,data["mean_clean_flux"][pixlo:pixhi,data["orders_to_use"][order_idx]]./yscale_raw.-sqrt.(data["mean_clean_var"][pixlo:pixhi,data["orders_to_use"][order_idx]])./yscale_raw
			, color=:blue, ls=:dot, label=:none)
	end
	if eltype(data["normalization_anchors"][1]) <:AbstractVector
		plot!(plt,anchor_density[:,order_idx]./maximum(anchor_density[:,order_idx]),color=:orange,label="Anchor density")
	else
		anchor_idx = data["normalization_anchors"][order_idx]
		if length(anchor_idx) >0
			scatter!(plt,anchor_idx,ones(length(anchor_idx)),color=:magenta,label="Anchors")
			scatter!(plt,anchor_idx,data["mean_clean_flux"][anchor_idx,data["orders_to_use"][order_idx]]./yscale_raw,color=:blue,ms=2,label=:none)
		end
	end
	xlims!(pixlo,pixhi)
	#ylims!(0.0,1.1)
	plt
end

# ╔═╡ eff49841-4793-486a-8ed0-1ebe6deb4759
begin
	anchors_merged = []
	if eltype(data["normalization_anchors"][1]) <:AbstractVector
	threshold = maximum(anchor_density[:,order_idx])/2
	for i in 2:(num_pix-1)
		if anchor_density[i,order_idx] > threshold
			if (anchor_density[i,order_idx]>anchor_density[i+1,order_idx]) && (anchor_density[i,order_idx]>anchor_density[i-1,order_idx])
				push!(anchors_merged, i)
			end
		end
	end
	end
	anchors_merged
end

# ╔═╡ 3e6f6f21-878d-4618-9865-c1e656bc3104
begin
	if eltype(data["normalization_anchors"][1]) <:AbstractVector
	plt12 = scatter(anchor_density[:,order_idx])
	xlabel!(plt12,"Pixel")
	ylabel!(plt12,"Anchor density")
	#xlims!(plt12,900,1100)
	plt12
	end
end

# ╔═╡ c7a9d874-d823-45ff-8214-da400d01f383
collect(keys(data))

# ╔═╡ e3a2d62e-2d0f-4306-bb26-474445ce894e
data["orders_to_use"]

# ╔═╡ 1d3b78eb-983a-4f74-b064-e88a387dc5a1
data["mean_clean_flux"]

# ╔═╡ e0d97999-78cc-4e89-9534-66ca44dffee5
data["mean_clean_flux_sed_normalized"]

# ╔═╡ 4e718a48-6b04-441d-a0f6-17af907780ed
data["mean_clean_flux_continuum_normalized"]

# ╔═╡ c15b050a-f0de-4079-8fcb-912630722788
all(isnan.(data["mean_clean_flux_continuum_normalized"][1,:]))

# ╔═╡ a0103542-f745-4913-a8cf-a3bddd1631e4
length(data["mean_clean_flux_continuum_normalized"][:,1])

# ╔═╡ d9782c42-a2f7-41ae-a324-0cc86a56f393
data["order_ccfs"]

# ╔═╡ 9bafebc2-daae-42ea-9f68-4020a590a030
if eltype(data["normalization_anchors"][1]) <:AbstractVector
	map(i->data["normalization_anchors"][i][21],1:length(data["normalization_anchors"]))
end

# ╔═╡ 2735c043-d57d-4782-8ad1-a1254e373fe0
keys(data)

# ╔═╡ 82f4ceb0-617a-434e-97de-c33d500958f4
begin
	mrv = map(ord_idx->RVFromCCF.MeasureRvFromCCFTemplate(v_grid=data["v_grid"],template=vec(mean(data["order_ccfs"],dims=3)[:,ord_idx,:]), frac_of_width_to_fit = 1.0, measure_width_at_frac_depth=0.5),1:num_order_idx)
	begin
		rv_matrix = zeros(num_order_idx, num_obs)
		σ_rv_matrix = zeros(num_order_idx, num_obs)
		for obs_idx in 1:num_obs
			for ord_idx in 1:num_order_idx
				(rv, σ_rv) = mrv[ord_idx](data["v_grid"],data["order_ccfs"][:,ord_idx,obs_idx],data["order_ccf_vars"][:,ord_idx,obs_idx])
				rv_matrix[ord_idx,obs_idx] = rv
				σ_rv_matrix[ord_idx,obs_idx] = σ_rv
			end
		end
	end
end

# ╔═╡ c7ec68a3-6a2b-4cad-b1ff-2cb9b23cc499
begin
	scatter(σ_rv_matrix[:,1])
	scatter!(map(i->std(rv_matrix[i,:]),1:21))
end

# ╔═╡ e7541e93-8b7a-4c26-901e-d887931b8103
order_idx_for_rvs = vcat(1:5,6,8,9,10,17,21)
#order_idx_for_rvs = vcat(1:11,19:21)
#order_idx_for_rvs = vcat(1:5,8:10,19:21)

# ╔═╡ 1a8e2e8d-97a9-4634-9121-0f33266227f8
scatter(data["orders_to_use"][1:21],map(i->std(rv_matrix[i,:]),1:21))

# ╔═╡ e7906830-5dfd-42ea-ad8e-948b7a3138db
begin
	rvs = map(i->NaNMath.sum(rv_matrix[order_idx_for_rvs,i]./σ_rv_matrix[order_idx_for_rvs,i].^2)./NaNMath.sum(1.0 ./σ_rv_matrix[order_idx_for_rvs,i].^2),1:num_obs)
	rvs .+= data["manifest"][!,"Δv_diff_ext"]
	σ_rvs = sqrt.(map(i->1.0 ./NaNMath.sum(1.0 ./σ_rv_matrix[order_idx_for_rvs,i].^2),1:num_obs))
end;

# ╔═╡ 9d38c13c-fbba-4fa2-b528-feae03551f54
(t_binned, rvs_binned) = bin_times_and_rvs_max_Δt(times=data["manifest"].bjd,rvs=rvs,Δt_threshold= 5/(24*60))

# ╔═╡ aac89b26-bd5d-402e-9b1c-4ab6a70dfe51
(NaNMath.std(rvs), NaNMath.mean(σ_rvs), std(rvs_binned))

# ╔═╡ 493b5c40-214a-41cd-aaaf-02609704ad49
begin
	daily_fit = Polynomials.fit(data["manifest"].bjd.-mean(data["manifest"].bjd),rvs,1)
	daily_slope = Polynomials.fit(t_binned,rvs_binned,1)[1]
end

# ╔═╡ 2921145f-df6a-4b14-9faf-46a27e8f9057
#std(rvs.-daily_fit.(data["manifest"].bjd.-mean(data["manifest"].bjd))),
binned_rms_minus_linear = std(rvs_binned.-daily_fit.(t_binned.-mean(data["manifest"].bjd)))

# ╔═╡ facd440e-0e06-4824-9da0-cc7dd549db91
begin
	plt11 = scatter(data["manifest"].bjd.-mean(data["manifest"].bjd),rvs,yerr=σ_rvs, ms=2.0,label=:none)
	scatter!(plt11,t_binned.-mean(data["manifest"].bjd),rvs_binned, ms=4,label="Binned")
	plot!(plt11,data["manifest"].bjd.-mean(data["manifest"].bjd),daily_fit.(data["manifest"].bjd.-mean(data["manifest"].bjd)),label="Daily slope = " * string(round(Polynomials.fit(t_binned,rvs_binned,1)[1]/24,digits=3)) * "m/s/hr")
	xlabel!(plt11,"Time (d)")
	ylabel!(plt11,"RV (m/s)")
	title!(plt11,date_str * "  RMS_5min-line = " * string(round(binned_rms_minus_linear,digits=3)) * "m/s")
	plt11
end

# ╔═╡ 8c9d4dbf-c36e-4eee-b49f-ed0be110b55b
scatter(t_binned,rvs_binned)

# ╔═╡ 01d28d16-6acc-4209-893f-cd8cb730922d
begin
	scatter(rv_matrix[1,:])
	scatter!(rv_matrix[2,:])
	scatter!(rv_matrix[3,:])
	scatter!(rv_matrix[4,:])
	scatter!(rv_matrix[5,:])
	scatter!(rv_matrix[6,:])
	scatter!(rv_matrix[8,:])
	scatter!(rv_matrix[9,:])
	scatter!(rv_matrix[10,:])
end

# ╔═╡ 859a9d38-f3d5-42ea-ba91-1673ded404ad
data["manifest"][!,"Δv_diff_ext"]

# ╔═╡ Cell order:
# ╠═42b4e6fe-a897-11eb-2676-f7fd96f35a22
# ╠═5313b674-a897-11eb-3820-218f26a14d4d
# ╠═37a6c623-e901-4f8b-8989-835958471549
# ╠═d1fe0408-cccd-45c9-a5ef-c13a42386e44
# ╟─609e3d45-e1a5-4060-97a8-5b668f8137c1
# ╟─d251eb17-a164-4a36-90f3-d1d5aa87405c
# ╟─28aced5a-3100-4523-9f1a-f7b14016791b
# ╟─abbc4880-fedd-4bdd-a1db-b94e8a4d5654
# ╟─6c469209-c7e0-41f2-97e3-db587b894e56
# ╟─ab0fc462-a897-11eb-1727-e71788522c02
# ╟─c403d4ff-3035-4bfc-9b9e-3f33ed0a6b7a
# ╟─d84af1b1-1023-4493-8b0f-90c9030be4bd
# ╟─7dc1bf77-174b-4614-8aea-15cad25c0335
# ╟─ceae0ff5-f650-4dbf-81e5-26864eb1adcc
# ╟─df7571c9-8b3f-4503-8bbe-42a77f0f9c0e
# ╟─468fe2c8-7c4c-45f2-9577-a2cecba64d2d
# ╟─0e4e2de1-d706-4019-bece-e24367b38d33
# ╟─ed5e02d9-cfba-4ee6-85f2-3678cb1899f9
# ╟─eff49841-4793-486a-8ed0-1ebe6deb4759
# ╟─3e6f6f21-878d-4618-9865-c1e656bc3104
# ╠═c7a9d874-d823-45ff-8214-da400d01f383
# ╠═e3a2d62e-2d0f-4306-bb26-474445ce894e
# ╠═1d3b78eb-983a-4f74-b064-e88a387dc5a1
# ╠═e0d97999-78cc-4e89-9534-66ca44dffee5
# ╠═4e718a48-6b04-441d-a0f6-17af907780ed
# ╠═c15b050a-f0de-4079-8fcb-912630722788
# ╠═a0103542-f745-4913-a8cf-a3bddd1631e4
# ╠═d9782c42-a2f7-41ae-a324-0cc86a56f393
# ╠═9bafebc2-daae-42ea-9f68-4020a590a030
# ╠═2735c043-d57d-4782-8ad1-a1254e373fe0
# ╠═82f4ceb0-617a-434e-97de-c33d500958f4
# ╠═c7ec68a3-6a2b-4cad-b1ff-2cb9b23cc499
# ╠═e7541e93-8b7a-4c26-901e-d887931b8103
# ╠═1a8e2e8d-97a9-4634-9121-0f33266227f8
# ╠═e7906830-5dfd-42ea-ad8e-948b7a3138db
# ╠═9d38c13c-fbba-4fa2-b528-feae03551f54
# ╠═aac89b26-bd5d-402e-9b1c-4ab6a70dfe51
# ╠═493b5c40-214a-41cd-aaaf-02609704ad49
# ╠═2921145f-df6a-4b14-9faf-46a27e8f9057
# ╠═facd440e-0e06-4824-9da0-cc7dd549db91
# ╠═8c9d4dbf-c36e-4eee-b49f-ed0be110b55b
# ╠═01d28d16-6acc-4209-893f-cd8cb730922d
# ╠═859a9d38-f3d5-42ea-ba91-1673ded404ad


text:
struct WRectangle <: Widget
    geometry::Quadrangle
    z::Int
    cb::WidgetCallbacks
end

AbstractGUI.zindex(w::WRectangle) = w.z
AbstractGUI.vertex_data(w::WRectangle) = w.geometry
Base.in(p::Point, w::WRectangle) = p in w.geometry
AbstractGUI.callbacks(w::WRectangle) = w.cb


text:
primitive type F2 <: Signed 8 end

F2(x::Int8) = reinterpret(F2, x)
F2(x::Number) = F2(Int8(x))
Core.Int8(x::F2) = reinterpret(Int8, x)
Core.Int32(x::F2) = Int32(Int8(x))
Core.Int64(x::F2) = Int64(Int8(x))
Core.Bool(x::F2) = ~iszero(x)
Base.show(io::IO, x::F2) = print(io, Int8(x))
#Base.:+(a::F2, b::F2) = F2(mod(Int8(a) + Int8(b), Int8(2)))
Base.:+(a::F2, b::F2) = F2(Int8(a) ⊻ Int8(b))
Base.:-(a::F2, b::F2) = a + b
Base.:*(a::F2, b::F2) = F2(Int8(a) & Int8(b))

function standard_basis(T, dim::Integer)
    return Matrix{T}(I, (dim, dim))
end

function standard_basis(example::Array)
    T = eltype(example)
    dim, _ = size(example)
    return standard_basis(T, dim)
end

function swaprow!(x, i, j)
    for k in axes(x, 2)
        x[i, k], x[j, k] = x[j, k], x[i, k]
    end
    return x
end

function swapcol!(x, i, j)
    for k in axes(x, 1)
        x[k, i], x[k, j] = x[k, j], x[k, i]
    end
    return x
end

nonzero = (!) ∘ iszero

function smith_normal_form(M::Matrix{F2})
    @inbounds begin
    M = copy(M)
    nrows, ncols = size(M)
    U = standard_basis(F2, nrows)
    V = standard_basis(F2, ncols)
    for c in 1:min(nrows, ncols)
        found = false
        k = c
        l = c
        while !found
            if nonzero(M[k, l])
                swaprow!(M, c, k)
                swaprow!(U, c, k)
                swapcol!(M, c, l)
                swapcol!(V, c, l)
                found = true
            elseif l == ncols && k == nrows
                # no nonzeros in row or column
                found = true
            elseif l == ncols
                l = c
                k += 1
            else
                l += 1
            end
        end
        for i in findall(nonzero, view(M, c+1:nrows, c)) .+ c
            M[i, :] += M[c, :]
            U[i, :] += U[c, :]
        end
        for j in findall(nonzero, view(M, c, c+1:ncols)) .+ c
            M[:, j] += M[:, c]
            V[:, j] += V[:, c]
        end
    end
    end
    return U, M, V
end

function smith_normal_form(M::Union{BitMatrix, Matrix{Bool}})
    @inbounds begin
    M = copy(M)
    nrows, ncols = size(M)
    U = standard_basis(Bool, nrows)
    V = standard_basis(Bool, ncols)
    for c in 1:min(nrows, ncols)
        found = false
        k = c
        l = c
        while !found
            if nonzero(M[k, l])
                swaprow!(M, c, k)
                swaprow!(U, c, k)
                swapcol!(M, c, l)
                swapcol!(V, c, l)
                found = true
            elseif l == ncols && k == nrows
                # no nonzeros in row or column
                found = true
            elseif l == ncols
                l = c
                k += 1
            else
                l += 1
            end
        end
        for i in findall(identity, view(M, c+1:nrows, c)) .+ c
            for j in c:ncols
                M[i, j] = M[i, j] ⊻ M[c, j]
            end
            U[i, :] = U[i, :] .⊻ U[c, :]
        end
        for j in findall(identity, view(M, c, c+1:ncols)) .+ c
            for i in c:nrows
                M[i, j] = M[i, j] ⊻ M[i, c]
            end
            V[:, j] = V[:, j] .⊻ V[:, c]
        end
    end
    end
    return U, M, V
end

function rowreduction(M::Matrix{F2})
    @inbounds begin
    M = copy(M)
    nrows, ncols = size(M)
    U = standard_basis(F2, nrows)
    pivotcols = []
    pivotcol = 1 :: Int64
    for r in 1:nrows
        if pivotcol > ncols # check
            return U, M, pivotcols
        end
        i = r :: Int64
        while iszero(M[i, pivotcol])
            i += 1
            if i > nrows
                i = r
                pivotcol += 1
                if pivotcol > ncols #check
                    return U, M, pivotcols
                end
            end
        end
        swaprow!(M, i, r)
        swaprow!(U, i, r)
        # findall(nonzero, view(M, :, pivotcol))
        for j in 1:nrows
            if nonzero(M[j, pivotcol]) && j != r
                for k in pivotcol:ncols
                    M[j, k] = M[j, k] + M[r, k] #check
                end
                U[j, :] += U[r, :]
            end
        end
        push!(pivotcols, pivotcol)
        pivotcol += 1
    end
    end
    return U, M, pivotcols
end

function rowreduction(M::Union{BitMatrix, Matrix{Bool}})
    @inbounds begin
    M = copy(M)
    nrows, ncols = size(M)
    U = standard_basis(Bool, nrows)
    pivotcols = []
    pivotcol = 1 :: Int64
    for r in 1:nrows
        if pivotcol > ncols # check
            return U, M, pivotcols
        end
        i = r :: Int64
        while !M[i, pivotcol]
            i += 1
            if i > nrows
                i = r
                pivotcol += 1
                if pivotcol > ncols #check
                    return U, M, pivotcols
                end
            end
        end
        swaprow!(M, i, r)
        swaprow!(U, i, r)
        for j in findall(view(M, :, pivotcol))
            if j != r
                for k in pivotcol:ncols
                    M[j, k] = M[j, k] ⊻ M[r, k] #check
                end
                U[j, :] = U[j, :] .⊻ U[r, :]
            end
        end
        push!(pivotcols, pivotcol)
        pivotcol += 1
    end
    end
    return U, M, pivotcols
end

function issquare(M)
    nrows, ncols = size(M)
    square = nrows == ncols
    return square
end

function isfullrank(M)
    _, echelonform, _ = rowreduction(M)
    fullrank = isfullrank(M, echelonform)
    return fullrank
end

function isfullrank(M, echelonform)
    fullrank = all(any(map(nonzero, echelonform), dims=2))
    return fullrank
end

function inverse(M)
    inverse, echelonform, _ = rowreduction(M)
    square = issquare(M)
    fullrank = isfullrank(M, echelonform)
    if ~square || ~fullrank
        throw(ErrorException("Matrix is not invertible"))
    else
        return inverse
    end
end

function kernel_image(M)
    U, M, V = smith_normal_form(M)
    kernel = V[:, [i[2] for i in findall(!, any(nonzero.(M), dims=1))]]
    image = inverse(U)[:, [i[1] for i in findall(any(nonzero.(M), dims=2))]]
    return kernel, image
end

function extend_basis_relative(U, V)#::FilteredBasis
    dim, nelements = size(U)
    all_vectors = hcat(U, V)
    transform, echelonform, pivotcols = rowreduction(all_vectors)
    extended_basis = all_vectors[:, pivotcols]
    return [U, extended_basis[:, (nelements + 1):end]]
end

function extend_basis(U)#::FilteredBasis
    return extend_basis_relative(U, standard_basis(U))
end

function extend_basis(Us...)#::FilteredBasis
    Us = [Us...]
    ambient_basis = standard_basis(first(Us))
    push!(Us, ambient_basis)
    descending = reverse(Us)
    subspace = pop!(descending)
    filtered_bases = [subspace]
    while ~isempty(descending)
        superspace = pop!(descending)
        _, extension = extend_basis_relative(subspace, superspace)
        push!(filtered_bases, extension)
        subspace = hcat(subspace, extension)
    end
    return filtered_bases
end

max_stars_repo_name: osumray/ChainComplexes

max_issues_repo_name: osumray/ChainComplexes

max_forks_repo_name: osumray/ChainComplexes


text:
function initstate(interpreter::AbstractInterpreter,
                   domain::AbstractedDomain, problem::GenericProblem)
    concrete_state = initstate(ConcreteInterpreter(), domain.domain, problem)
    absdom = AbstractedDomain(domain, interpreter)
    return abstractstate(absdom, concrete_state)
end

function initstate(interpreter::AbstractInterpreter,
                   domain::AbstractedDomain, objtypes::AbstractDict)
    concrete_state = initstate(ConcreteInterpreter(), domain.domain, objtypes)
    absdom = AbstractedDomain(domain, interpreter)
    return abstractstate(absdom, concrete_state)
end

function initstate(interpreter::AbstractInterpreter,
                   domain::AbstractedDomain, objtypes::AbstractDict, fluents)
    concrete_state = initstate(ConcreteInterpreter(),
                               domain.domain, objtypes, fluents)
    absdom = AbstractedDomain(domain, interpreter)
    return abstractstate(absdom, concrete_state)
end