text:
module TGridReader
#=
This module can read fluent .msh files in ascii formating
=#
include("TGridTypes.jl")

export TgridMeshImport

function TgridMeshImport(filename::String)
    io = Base.open(filename, "r")
    sections = loopGridSections(io)

    myTGrid = evalSections(sections)
    return myTGrid
end

function loopGridSections(io)
    sections = gridSectionStruct[]
    iobuff = IOBuffer(UInt8[], read=true, write=true, append=true)
    openbrackets = 0
    closedbrackets= 0

    for line in eachline(io)

            for s in line
                if s == '('
                    if openbrackets > 0
                        push!(sections, gridSectionStruct(openbrackets, read(iobuff,String)))
                        flush(iobuff)
                    end
                    openbrackets = openbrackets + 1
                elseif s == ')'
                    closedbrackets = closedbrackets + 1
                else
                    write(iobuff, s);
                end

                if openbrackets == closedbrackets
                    write(iobuff, "\n")
                    push!(sections, gridSectionStruct(openbrackets, read(iobuff,String)))
                    flush(iobuff)
                    openbrackets = 0
                    closedbrackets = 0
                end
            end

            if openbrackets != closedbrackets
                write(iobuff, "\n")
            end
    end

    if openbrackets != closedbrackets
        @error "Unbalances brackets in file"
    end

    return sections
end

function evalSections(sections::Array{gridSectionStruct, 1})
    mainGridID = -1
    dimensions = -1
    header = ""
    nodesSections = NodeSection[]
    cellSections = CellSection[]
    faceSections = FaceSection[]
    namedZones = NamedZone[]
    zoneID = -1
    firstIdx = 0
    lastIdx = 0
    type = -1
    cellType = -1
    faceType = -1
    i = 0

    for (k, sec) in enumerate(sections)

        if sec.gridIDx == 1
            s = split(sec.content, " ")
            contentArr = s[s .!= ""]

            mainGridID = parse(Int64, contentArr[1])

            if  length(contentArr) > 1

                if mainGridID == xfgridsections["xf-dimension"]
                    dimensions = parse(Int64, contentArr[2])
                end

                if mainGridID == xfgridsections["xf-header"]
                    header = contentArr[2:end]
                end

            end
            i = 1

        else
            if mainGridID == xfgridsections["xf-comment"]
                #empty

            elseif mainGridID == xfgridsections["xf-node"]
                s = split(sec.content, " ")
                contentArr = s[s .!= ""]

                if i == 1
                    s = split(sec.content, " ")
                    contentArr = s[s .!= ""]
                    zoneID = parse(Int64, contentArr[1], base=16)

                    if zoneID == 0
                        firstIdx = 1
                        lastIdx = parse(Int64, contentArr[3], base=16)
                        push!(nodesSections, NodeSection(zoneID, firstIdx, lastIdx, nothing))
                    else
                        firstIdx = parse(Int64, contentArr[2], base=16)
                        lastIdx = parse(Int64, contentArr[3], base=16)
                    end
                elseif i==2
                    s = split(sec.content, r" |\n")
                    contentArr = s[s .!= ""]

                    if dimensions == 3
                        x = parse.(Float64, contentArr[1:3:end])
                        y = parse.(Float64, contentArr[2:3:end])
                        z = parse.(Float64, contentArr[3:3:end])
                    elseif dimensions == 2
                        x = parse.(Float64, contentArr[1:2:end])
                        y = parse.(Float64, contentArr[2:2:end])
                    else
                        @error "wrong size of dimensions in node section"
                    end

                    if dimensions == 3
                        push!(nodesSections, NodeSection(zoneID, firstIdx, lastIdx, [x y z]))
                    elseif dimensions == 2
                        push!(nodesSections, NodeSection(zoneID, firstIdx, lastIdx, [x y]))
                    end
                else
                    @warn "Wrong size of subsections in node section"
                end

            elseif mainGridID == xfgridsections["xf-periodic-face"]
                # not implemented jet
                @warn "xf-periodic-face no implemented jet"

            elseif mainGridID == xfgridsections["xf-cell"]
                s = split(sec.content, " ")
                contentArr = s[s .!= ""]

                if i == 1
                    s = split(sec.content, " ")
                    contentArr = s[s .!= ""]
                    zoneID = parse(Int64, contentArr[1], base=16)

                    if zoneID == 0
                        firstIdx = 1
                        lastIdx = parse(Int64, contentArr[3], base=16)
                        type = 0
                        cellType = -1
                    else
                        firstIdx = parse(Int64, contentArr[2], base=16)
                        lastIdx = parse(Int64, contentArr[3], base=16)
                        type = parse(Int64, contentArr[4])
                        cellType = parse(Int64, contentArr[5])
                    end

                    if cellType != 0
                        push!(cellSections, CellSection(zoneID, firstIdx,
                        lastIdx, type, cellType, nothing))
                    end

                elseif i==2
                    if cellType == 0
                        s = split(sec.content, r" |\n")
                        contentArr = s[s .!= ""]

                        mixedCellTypes = parse.(Int, contentArr)
                        push!(cellSections, CellSection(zoneID, firstIdx,
                            lastIdx, type, cellType, mixedCellTypes))
                    else
                        @warn "This should not occur possible parsing error of cell section..."
                    end

                else
                    @warn "Wrong size of subsections in cell section"
                end
            elseif mainGridID == xfgridsections["xf-face"]

                s = split(sec.content, " ")
                contentArr = s[s .!= ""]

                if i == 1
                    s = split(sec.content, " ")
                    contentArr = s[s .!= ""]
                    #println(contentArr[1])
                    zoneID = parse(Int64, contentArr[1], base=16)

                    if zoneID == 0
                        firstIdx = 1
                        lastIdx = parse(Int64, contentArr[3], base=16)
                        type = 0
                        faceType = -1
                        push!(faceSections, FaceSection(zoneID, firstIdx,
                        lastIdx, type, faceType, nothing))
                    else
                        firstIdx = parse(Int64, contentArr[2], base=16)
                        lastIdx = parse(Int64, contentArr[3], base=16)
                        type = parse(Int64, contentArr[4])
                        faceType = parse(Int64, contentArr[5])
                    end

                elseif i==2
                    localFaceType = 0
                    s = split(sec.content, "\n")
                    contentArr = s[s .!= ""]
                    faces = Array{Face,1}(undef, length(contentArr))

                    if faceType == 0
                        for (ii, ss) in enumerate(contentArr)
                            sss = split(ss, " ")
                            contentArr1 = sss[sss .!= ""]
                            localFaceType = parse(Int64, contentArr1[1])

                            faceNodes = facetype_nodes_per_face[localFaceType]

                            nodes = parse(Int64, contentArr1[2:(1+faceNodes)], base=16)
                            cr = parse(Int64, contentArr1[(2+faceNodes)], base=16)
                            cl = parse(Int64, contentArr1[(3+faceNodes)], base=16)

                            faces[ii] =  Face(localFaceType, nodes, cr, cl)
                        end
                    else
                        localFaceType = faceType
                        faceNodes = facetype_nodes_per_face[localFaceType]

                        for (ii, ss) in enumerate(contentArr)
                            sss = split(ss, " ")
                            contentArr1 = sss[sss .!= ""]

                            nodes = parse.(Int64, contentArr1[1:faceNodes], base=16)
                            cr = parse(Int64, contentArr1[(1+faceNodes)], base=16)
                            cl = parse(Int64, contentArr1[(2+faceNodes)], base=16)

                            faces[ii] =  Face(localFaceType, nodes, cr, cl)
                        end

                    end

                    push!(faceSections, FaceSection(zoneID, firstIdx,
                    lastIdx, type, faceType, faces))
                else
                    @warn "Wrong size of subsections in face section"
                end

            elseif mainGridID == xfgridsections["xf-face-tree"]
                # not implemented jet
                @warn "xf-face-tree no implemented jet"

            elseif mainGridID == xfgridsections["xf-cell-tree"]
                # not implemented jet
                @warn "xf-cell-tree no implemented jet"

            elseif mainGridID == xfgridsections["xf-face-parents"]
                # not implemented jet
                @warn "xf-face-parents no implemented jet"

            elseif mainGridID == xfgridsections["xf-zone"]
                s = split(sec.content, " ")
                contentArr = s[s .!= ""]

                if i == 1
                    s = split(sec.content, " ")
                    contentArr = s[s .!= ""]
                    zoneID = parse(Int64, contentArr[1], base=16)
                    zoneType = contentArr[2]
                    zoneName = contentArr[3]

                    push!(namedZones, NamedZone(zoneID, zoneType, zoneName))
                elseif i==2

                else
                    @warn "Wrong size of subsections in face section"
                end

            else
                @warn "no fitting parsing for given main ID"
            end

            i = i + 1;

        end
    end

  myTGrid= TGrid(header, dimensions, cellSections, faceSections, nodesSections, namedZones)

return myTGrid
end

end

max_stars_repo_name: c-schubert/linearMeshIO

max_issues_repo_name: c-schubert/linearMeshIO

max_forks_repo_name: c-schubert/linearMeshIO


text:
## Utility functions for getting type names and modules (from IncrementalInference)
_getmodule(t::T) where T = T.name.module
_getname(t::T) where T = T.name.name

convertPackedType(t::Union{T, Type{T}}) where {T <: FunctorInferenceType} = getfield(_getmodule(t), Symbol("Packed$(_getname(t))"))
function convertStructType(::Type{PT}) where {PT <: PackedInferenceType}
    # see #668 for expanded reasoning.  PT may be ::UnionAll if the type is of template type.
    ptt = PT isa DataType ? PT.name.name : PT
    moduleName = PT isa DataType ? PT.name.module : Main
    symbolName = Symbol(string(ptt)[7:end])
    getfield(moduleName, symbolName)
end

##==============================================================================
## Sorting
##==============================================================================
#Natural Sorting less than

# Adapted from https://rosettacode.org/wiki/Natural_sorting
# split at digit to not digit change
splitbynum(x::AbstractString) = split(x, r"(?<=\D)(?=\d)|(?<=\d)(?=\D)")
#parse to Int
numstringtonum(arr::Vector{<:AbstractString}) = [(n = tryparse(Int, e)) != nothing ? n : e for e in arr]
#natural less than
function natural_lt(x::T, y::T) where T <: AbstractString
    xarr = numstringtonum(splitbynum(x))
    yarr = numstringtonum(splitbynum(y))
    for i in 1:min(length(xarr), length(yarr))
        if typeof(xarr[i]) != typeof(yarr[i])
            return isa(xarr[i], Int)
        elseif xarr[i] == yarr[i]
            continue
        else
            return xarr[i] < yarr[i]
        end
    end
    return length(xarr) < length(yarr)
end

natural_lt(x::Symbol, y::Symbol) = natural_lt(string(x),string(y))

"""
    $SIGNATURES

Convenience wrapper for `Base.sort`.
Sort variable (factor) lists in a meaningful way (by `timestamp`, `label`, etc), for example `[:april;:x1_3;:x1_6;]`
Defaults to sorting by timestamp for variables and factors and using `natural_lt` for Symbols.
See Base.sort for more detail.

Notes
- Not fool proof, but does better than native sort.

Example

`sortDFG(ls(dfg))`
`sortDFG(ls(dfg), by=getLabel, lt=natural_lt)`

Related

ls, lsf
"""
sortDFG(vars::Vector{<:DFGNode}; by=getTimestamp, kwargs...) = sort(vars; by=by, kwargs...)
sortDFG(vars::Vector{Symbol}; lt=natural_lt, kwargs...) = sort(vars; lt=lt, kwargs...)

##==============================================================================
## Validation of session, robot, and user IDs.
##==============================================================================
global _invalidIds = ["USER", "ROBOT", "SESSION", "VARIABLE", "FACTOR", "ENVIRONMENT", "PPE", "DATA_ENTRY", "FACTORGRAPH"]

global _validLabelRegex = r"^[a-zA-Z]\w*$"

"""
$(SIGNATURES)

Returns true if the label is valid for a session, robot, or user ID.
"""
function isValidLabel(id::Union{Symbol, String})::Bool
    if typeof(id) == Symbol
        id = String(id)
    end
    return all(t -> t != uppercase(id), _invalidIds) && match(_validLabelRegex, id) != nothing
end

"""
    $SIGNATURES

Small utility to return `::Int`, e.g. `0` from `getVariableLabelNumber(:x0)`

Examples
--------
```julia
getVariableLabelNumber(:l10)          # 10
getVariableLabelNumber(:x1)           # 1
getVariableLabelNumber(:x1_10, "x1_") # 10
```

DevNotes
- make prefix Regex based for longer -- i.e. `:apriltag578`, `:lm1_4`

"""
getVariableLabelNumber(vs::Symbol, prefix=string(vs)[1]) = parse(Int, string(vs)[(length(prefix)+1):end])

ext: jl


text:
using PMDPs

using StaticArrays
using POMDPs
using MCTS, DiscreteValueIteration
using StatsBase
using Plots
using POMDPSimulators
using D3Trees
using POMDPPolicies
using Random

using Traceur

edges = create_edges(5, 3, [50,60,70,80,90])
products = create_continuous_products(edges)
λ = create_λ(Float64[10,3,3,5,4], products)
mdp = PMDP(edges, products, λ)

# @requirements_info ValueIterationSolver() mdp
@requirements_info MCTSSolver() mdp State{5}(SA[1,1,1,1,1], 89, SA[0,0,1,1,1])

solver = MCTSSolver(n_iterations=100,
                    depth=100,
                    exploration_constant=10.0,
                    enable_tree_vis=true)
planner = solve(solver, mdp)
s = initialstate(mdp, Random.MersenneTwister(4))
a = action(planner, s)
println("Action:")
println(a)
println("State:")
println(s)

s_actions = zeros(500)
for i in 1:10
    a = action(planner, s)
    s_actions[i]=a
end

# println(s_actions)

display(bar(countmap(s_actions)))
println(countmap(s_actions))
println(mean(s_actions))
# display(plot(s_actions))

rand_policy = RandomPolicy(mdp)

# hr = HistoryRecorder(max_steps=100)
# history = simulate(hr, mdp,  )

initial_state = State{5}(SA[5,5,5,5,5], 0, SA[0,0,0,0,0])
rollout_sim = RolloutSimulator(max_steps=10)
r_mcts = simulate(rollout_sim, mdp, planner, initial_state)
r_rand = simulate(rollout_sim, mdp, rand_policy, initial_state)

initial_state = State{5}(SA[5,5,5,5,5], 0, SA[0,0,0,0,0])
hr = HistoryRecorder(max_steps=30)
h_mcts = simulate(hr, mdp, planner, initial_state)
h_rand = simulate(hr, mdp, rand_policy, initial_state)

@show state_hist(h_mcts)
@show collect(action_hist(h_mcts))


text:
"""
    mutable struct Structure_Si  <: Str

This contains information for model constructing of Si and its compounds.

# Supported list:
- structure_si()
- structure_si3n4()
- structure_si3n4_ort()
- structure_sio2()
"""
mutable struct Structure_Si <: Str
    atom_vec::Matrix{Float64}
    cell_vec::Matrix{Float64}
    atom_type
    atom_name
    atom_charge
    para_mass
    num_atoms
    num_atom_types
    vec_type_id
    res_name
end

"""
    structure_si()

Do this will generate a Famliy_Si type which contains all information needed to build a Si moedel.
"""
function structure_si()
    atom_vec = [
        0       0       0
        0       0.5     0.5
        0.5     0       0.5
        0.5     0.5     0
        0.25    0.25    0.25
        0.25    0.75    0.75
        0.75    0.25    0.75
        0.75    0.75    0.25
    ]
    cell_vec = [
        5.4310      0           0
        0           5.4310      0
        0           0           5.4310
    ]
    atom_type = [1 1 1 1 1 1 1]
    atom_charge = 0 .* atom_type
    para_mass = [28.085501] # Si: 28.085501. Unit: g/mol
    atom_name = split("Si")
    num_atoms = length(atom_type)
    num_atom_types = length(atom_name)
    res_name = "SI"
    Structure_Si(atom_vec, cell_vec, atom_type, atom_name, atom_charge, para_mass, num_atoms, num_atom_types, [1, 2], res_name)
end

"""
    structure_si3n4()

Do this will generate a Famliy_Si type which contains all information needed to build a Si3N4 moedel.
"""
function structure_si3n4()
    atom_vec = [
        -0.0977580565637542 0.195516113127508   0
        0.0588368288671674  0.307869462739661   0
        -0.0535176961032620 0.576818144016926   0
        0.327654573699065   0.151274131074550   0
        0.403949967083608   0                   0.500000000000000
        0.438645801493259   0.309541840068124   0
        0.749567567747191   0.272293435934189   0.500000000000000
        0.0180020380857845  0.645233580181295   0.500000000000000
        0.129048914298138   0.803653323477456   0.500000000000000
        0.510352849669919   0.378109310535080   0.500000000000000
        0.397866659130036   0.647057991812345   0.500000000000000
        0.554461544560957   0.759411341424498   0.500000000000000
        -0.291760564557441  0.683850293038517   0
        0.0514815980747357  0.957207969090818   0
    ]
    cell_vec = [
        7.595  0         0
        3.7975 6.577463  0
        0      0         2.902
    ]
    atom_type = [2 1 2 2 1 1 1 1 2 2 1 2 1 1]
    atom_charge = [-0.575925, 0.7679]
    para_mass = [14.0067 28.0855] # N: 14.0067. Si: 28.085501. Unit: g/mol
    atom_name = split("N Si")
    num_atoms = length(atom_type)
    num_atom_types = length(atom_name)
    res_name = "SIN"
    Structure_Si(atom_vec, cell_vec, atom_type, atom_name, atom_charge, para_mass, num_atoms, num_atom_types, [1, 2], res_name)
end

function structure_si3n4_pdb()
    atom_vec = [
        -0.0977580565637542 0.195516113127508   0
        0.0588368288671674  0.307869462739661   0
        -0.0535176961032620 0.576818144016926   0
        0.327654573699065   0.151274131074550   0
        0.403949967083608   0                   0.500000000000000
        0.438645801493259   0.309541840068124   0
        0.749567567747191   0.272293435934189   0.500000000000000
        0.0180020380857845  0.645233580181295   0.500000000000000
        0.129048914298138   0.803653323477456   0.500000000000000
        0.510352849669919   0.378109310535080   0.500000000000000
        0.397866659130036   0.647057991812345   0.500000000000000
        0.554461544560957   0.759411341424498   0.500000000000000
        -0.291760564557441  0.683850293038517   0
        0.0514815980747357  0.957207969090818   0
    ]
    cell_vec = [
        7.595  0         0
        3.7975 6.577463  0
        0      0         2.902
    ]
    atom_type = [2 1 2 2 1 1 1 1 2 2 1 2 1 1]
    _atom_charge = [-0.575925 0.7679]
    _para_mass = [14.0067 28.0855] # N: 14.0067. Si: 28.085501. Unit: g/mol
    _atom_name = split("N SI")
    atom_charge, para_mass, atom_name = [], [], []
    for (atom, type) in enumerate(atom_type)
        atom_charge =  vcat(atom_charge, _atom_charge[type])
        para_mass = vcat(para_mass, _para_mass[type])
        atom_name = vcat(atom_name, join([_atom_name[type], 'A'-1+length(findall(x->x==type,atom_type[1:atom]))]))
    end
    atom_type = Array(1:length(atom_type))
    num_atoms = length(atom_type)
    num_atom_types = length(atom_name)
    res_name = "SIN"
    Structure_Si(atom_vec, cell_vec, atom_type, atom_name, atom_charge, para_mass, num_atoms, num_atom_types, [1, 2], res_name)
end

"""
    structure_si3n4_ort()

Do this will generate a Famliy_Si type which contains all information needed to build a Si3N4 moedel with orthogonal unit cell.
"""
function structure_si3n4_ort()
    atom_vec = [
    0                       0.597031356376194   0
    0.000658327847268042    0.521394509327430   0.500000000000000
    0.0308755760368664      0.401100527527099   0.500000000000000
    0.0967083607636604      0.119153727824721   0
    0.126793943383805       0                   0
    0.127583936800527       0.923220010332222   0.500000000000000
    0.146872942725477       0.363320112468571   0
    0.190125082290981       0.676164982283494   0
    0.296115865701119       0.710448618483487   0.500000000000000
    0.309479921000658       0.175330240356416   0
    0.318104015799868       0.844922570673702   0.500000000000000
    0.331599736668861       0.309804192546631   0
    0.437327188940092       0.344011811814916   0.500000000000000
    0.482422646477946       0.657540834015205   0.500000000000000
    0.500000000000000       0.0970328006978962  0
    0.500658327847268       0.0213959536491320  0.500000000000000
    0.530875576036866       0.901099083205397   0.500000000000000
    0.596708360763660       0.619152283503018   0
    0.626793943383805       0.499998555678298   0
    0.627583936800527       0.423221454653924   0.500000000000000
    0.646872942725477       0.863318668146869   0
    0.690125082290981       0.176166426605196   0
    0.796115865701119       0.210450062805189   0.500000000000000
    0.809479921000659       0.675328796034714   0
    0.818104015799869       0.344924014995404   0.500000000000000
    0.831599736668861       0.809802748224929   0
    0.937327188940092       0.844010367493214   0.500000000000000
    0.982422646477946       0.157542278336908   0.500000000000000
    ]
    cell_vec = [
    7.595   0           0
    0       13.154964   0
    0       0           2.902
    ]
    atom_type = [2 1 2 2 1 2 1 1 2 1 1 2 1 1 2 1 2 2 1 2 1 1 2 1 1 2 1 1]
    atom_charge = 0 .* atom_type
    para_mass = [14.0067 28.0855] # N: 14.0067. Si: 28.085501. Unit: g/mol
    atom_name = split("N Si")
    num_atoms = length(atom_type)
    num_atom_types = length(atom_name)
    res_name = "SIN"
    Structure_Si(atom_vec, cell_vec, atom_type, atom_name, atom_charge, para_mass, num_atoms, num_atom_types, [1, 2], res_name)
end

function structure_si3n4_ort_pdb()
    atom_vec = [
    0                       0.597031356376194   0
    0.000658327847268042    0.521394509327430   0.500000000000000
    0.0308755760368664      0.401100527527099   0.500000000000000
    0.0967083607636604      0.119153727824721   0
    0.126793943383805       0                   0
    0.127583936800527       0.923220010332222   0.500000000000000
    0.146872942725477       0.363320112468571   0
    0.190125082290981       0.676164982283494   0
    0.296115865701119       0.710448618483487   0.500000000000000
    0.309479921000658       0.175330240356416   0
    0.318104015799868       0.844922570673702   0.500000000000000
    0.331599736668861       0.309804192546631   0
    0.437327188940092       0.344011811814916   0.500000000000000
    0.482422646477946       0.657540834015205   0.500000000000000
    0.500000000000000       0.0970328006978962  0
    0.500658327847268       0.0213959536491320  0.500000000000000
    0.530875576036866       0.901099083205397   0.500000000000000
    0.596708360763660       0.619152283503018   0
    0.626793943383805       0.499998555678298   0
    0.627583936800527       0.423221454653924   0.500000000000000
    0.646872942725477       0.863318668146869   0
    0.690125082290981       0.176166426605196   0
    0.796115865701119       0.210450062805189   0.500000000000000
    0.809479921000659       0.675328796034714   0
    0.818104015799869       0.344924014995404   0.500000000000000
    0.831599736668861       0.809802748224929   0
    0.937327188940092       0.844010367493214   0.500000000000000
    0.982422646477946       0.157542278336908   0.500000000000000
    ]
    cell_vec = [
    7.595   0           0
    0       13.154964   0
    0       0           2.902
    ]
    atom_type = [2 1 2 2 1 2 1 1 2 1 1 2 1 1 2 1 2 2 1 2 1 1 2 1 1 2 1 1]
    _atom_charge = [-0.575925 0.7679]
    _para_mass = [14.0067 28.0855] # N: 14.0067. Si: 28.085501. Unit: g/mol
    _atom_name = split("N SI")
    atom_charge, para_mass, atom_name = [], [], []
    for (atom, type) in enumerate(atom_type)
        atom_charge =  vcat(atom_charge, _atom_charge[type])
        para_mass = vcat(para_mass, _para_mass[type])
        atom_name = vcat(atom_name, join([_atom_name[type], 'A'-1+length(findall(x->x==type,atom_type[1:atom]))]))
    end
    atom_type = Array(1:length(atom_type))
    num_atoms = length(atom_type)
    num_atom_types = length(atom_name)
    res_name = "SIN"
    Structure_Si(atom_vec, cell_vec, atom_type, atom_name, atom_charge, para_mass, num_atoms, num_atom_types, [1, 2], res_name)
end

"""
    structure_sio2()

Do this will generate a Famliy_Si type which contains all information needed to build a SiO2 moedel.
"""
function structure_sio2()
    atom_vec = [
        0.196866211329851   0.196866211329851   0
        0.136601044596223   0                   0.178468624064479
        0.596625150662917   0.596625150662917   0.500000000000000
        0.656890317396545   0.793491361992768   0.678468624064479
        0.0966251506629169  0.696866211329852   0.250000000000000
        0.293491361992768   0.636601044596224   0.428468624064479
        0.696866211329852   0.0966251506629169  0.750000000000000
        0.500000000000000   0.156890317396545   0.928468624064479
        0.156890317396545   0.500000000000000   0.0715313759355211
        0.636601044596224   0.293491361992768   0.571531375935521
        0                   0.136601044596223   0.821531375935521
        0.793491361992768   0.656890317396545   0.321531375935521
    ]
    cell_vec = [
        4.9780  0           0
        0       4.9780      0
        0       0           6.9480
    ]
    atom_type = [2 1 2 1 2 1 2 1 1 1 1 1]
    atom_charge = 0 .* atom_type
    para_mass = [15.999400 28.0855] # O: 15.999400. Si: 28.085501. Unit: g/mol
    atom_name = split("O Si")
    num_atoms = length(atom_type)
    num_atom_types = length(atom_name)
    res_name = "SIO"
    Structure_Si(atom_vec, cell_vec, atom_type, atom_name, atom_charge, para_mass, num_atoms, num_atom_types, [1, 2], res_name)
end


text:
@test_throws DimensionMismatch splitobs((X, rand(149)))
@test_throws DimensionMismatch splitobs((X, rand(149)), obsdim=:last)

@testset "typestability" begin
    @testset "Int" begin
        @test_throws ArgumentError splitobs(10, 0.)
        @test_throws ArgumentError splitobs(10, 1.)
        @test_throws ArgumentError splitobs(10, (0.2,0.0))
        @test_throws ArgumentError splitobs(10, (0.2,0.8))
        @test_throws MethodError splitobs(10, 0.5, ObsDim.Undefined())
        @test typeof(@inferred(splitobs(10))) <: NTuple{2}
        @test eltype(@inferred(splitobs(10))) <: UnitRange
        @test typeof(@inferred(splitobs(10, 0.5))) <: NTuple{2}
        @test typeof(@inferred(splitobs(10, (0.5,0.2)))) <: NTuple{3}
        @test eltype(@inferred(splitobs(10, 0.5))) <: UnitRange
        @test eltype(@inferred(splitobs(10, (0.5,0.2)))) <: UnitRange
        @test eltype(@inferred(splitobs(10, at=0.5))) <: UnitRange
    end
    for var in vars
        @test_throws ArgumentError splitobs(var, 0.)
        @test_throws ArgumentError splitobs(var, 1.)
        @test_throws ArgumentError splitobs(var, (0.2,0.0))
        @test_throws ArgumentError splitobs(var, (0.2,0.8))
        @test_throws MethodError splitobs(var, 0.5, ObsDim.Undefined())
        @test typeof(@inferred(splitobs(var))) <: NTuple{2}
        @test eltype(@inferred(splitobs(var))) <: SubArray
        @test typeof(@inferred(splitobs(var, 0.5))) <: NTuple{2}
        @test typeof(@inferred(splitobs(var, (0.5,0.2)))) <: NTuple{3}
        @test eltype(@inferred(splitobs(var, 0.5))) <: SubArray
        @test eltype(@inferred(splitobs(var, (0.5,0.2)))) <: SubArray
        @test typeof(@inferred(splitobs(var, 0.5, ObsDim.Last()))) <: NTuple{2}
        @test typeof(@inferred(splitobs(var, 0.5, ObsDim.First()))) <: NTuple{2}
        @test eltype(@inferred(splitobs(var, 0.5, ObsDim.First()))) <: SubArray
        #@test_throws ErrorException @inferred(splitobs(var, at=0.5))
        #@test_throws ErrorException @inferred(splitobs(var, obsdim=:last))
        @test_throws ErrorException @inferred(splitobs(var, obsdim=1))
    end
    for tup in tuples
        @test_throws ArgumentError splitobs(tup, 0.)
        @test_throws ArgumentError splitobs(tup, 1.)
        @test_throws ArgumentError splitobs(tup, (0.2,0.0))
        @test_throws ArgumentError splitobs(tup, (0.2,0.8))
        @test typeof(@inferred(splitobs(tup, 0.5))) <: NTuple{2}
        @test typeof(@inferred(splitobs(tup, (0.5,0.2)))) <: NTuple{3}
        @test eltype(@inferred(splitobs(tup, 0.5))) <: Tuple
        @test eltype(@inferred(splitobs(tup, (0.5,0.2)))) <: Tuple
        @test typeof(@inferred(splitobs(tup, 0.5, ObsDim.Last()))) <: NTuple{2}
        @test eltype(@inferred(splitobs(tup, 0.5, ObsDim.Last()))) <: Tuple
        @test_throws ErrorException @inferred(splitobs(tup, obsdim=:last))
    end
end

@testset "Int" begin
    @test splitobs(10) == (1:7,8:10)
    @test splitobs(10, 0.5) == (1:5,6:10)
    @test splitobs(10, (0.5,0.3)) == (1:5,6:8,9:10)
    @test splitobs(150) == splitobs(150, 0.7)
    @test splitobs(150, at=0.5) == splitobs(150, 0.5)
    @test splitobs(150, at=(0.5,0.2)) == splitobs(150, (0.5,0.2))
    @test nobs.(splitobs(150)) == (105,45)
    @test nobs.(splitobs(150, at=(.2,.3))) == (30,45,75)
    @test nobs.(splitobs(150, at=(.1,.2,.3))) == (15,30,45,60)
    # tests if all obs are still present and none duplicated
    @test sum(sum.(getobs.(splitobs(150)))) == 11325
    @test sum(sum.(splitobs(150,at=.1))) == 11325
    @test sum(sum.(splitobs(150,at=(.2,.1)))) == 11325
    @test sum(sum.(splitobs(150,at=(.1,.4,.2)))) == 11325
    @test sum.(splitobs(150)) == (5565, 5760)
end

println("<HEARTBEAT>")

@testset "Array, SparseArray, and SubArray" begin
    for var in (Xs, ys, vars...)
        @test splitobs(var) == splitobs(var, 0.7, ObsDim.Last())
        @test splitobs(var, at=0.5) == splitobs(var, 0.5, ObsDim.Last())
        @test splitobs(var, at=(0.5,0.2)) == splitobs(var, (0.5,0.2), ObsDim.Last())
        @test splitobs(var, obsdim=1) == splitobs(var, 0.7, ObsDim.First())
        @test nobs.(splitobs(var)) == (105,45)
        @test nobs.(splitobs(var, at=(.2,.3))) == (30,45,75)
        @test nobs.(splitobs(var, at=(.2,.3), obsdim=:last)) == (30,45,75)
        @test nobs.(splitobs(var, at=(.1,.2,.3))) == (15,30,45,60)
    end
    @test nobs.(splitobs(X', obsdim=1),obsdim=1) == (105,45)
    # tests if all obs are still present and none duplicated
    @test sum(vec.(sum.(getobs.(splitobs(sparse(X1))),dims=2))) == fill(11325,10)
    @test sum(vec.(sum.(splitobs(X1),dims=2))) == fill(11325,10)
    @test sum(vec.(sum.(splitobs(X1,at=.1),dims=2))) == fill(11325,10)
    @test sum(vec.(sum.(splitobs(X1,at=(.2,.1)),dims=2))) == fill(11325,10)
    @test sum(vec.(sum.(splitobs(X1,at=(.1,.4,.2)),dims=2))) == fill(11325,10)
    @test sum(vec.(sum.(getobs.(splitobs(sparse(X1),at=(.2,.1))),dims=2))) == fill(11325,10)
    @test sum(vec.(sum.(splitobs(X1',obsdim=1),dims=1))) == fill(11325,10)
    @test sum.(splitobs(Y1)) == (5565, 5760)
    @test sum.(getobs.(splitobs(sparse(Y1)))) == (5565, 5760)
    @test sum.(splitobs(Y1, obsdim=:first)) == (5565, 5760)
end

println("<HEARTBEAT>")

@testset "Tuple of Array, SparseArray, and SubArray" begin
    for tup in ((Xs,ys), (X,ys), (Xs,y), (Xs,Xs), (XX,X,ys), (X,yv), (Xv,y), tuples...)
        @test_throws MethodError splitobs(tup, 0.5, ObsDim.Undefined())
        @test_throws MethodError splitobs(tup..., 0.5)
        @test_throws MethodError splitobs(tup...)
        @test all(map(x->(typeof(x)<:Tuple), splitobs(tup)))
        @test all(map(x->(typeof(x)<:Tuple), splitobs(tup,at=0.5)))
        @test nobs.(splitobs(tup)) == (105,45)
        @test nobs.(splitobs(tup, at=(.2,.3))) == (30,45,75)
        @test nobs.(splitobs(tup, at=(.2,.3), obsdim=:last)) == (30,45,75)
        @test nobs.(splitobs(tup, at=(.1,.2,.3))) == (15,30,45,60)
    end
    @test nobs.(splitobs((X',y), obsdim=1),obsdim=1) == (105,45)
    # tests if all obs are still present and none duplicated
    # also tests that both paramter are split disjoint
    train,test = splitobs((X1,Y1,X1))
    @test vec(sum(train[1],dims=2)+sum(test[1],dims=2)) == fill(11325,10)
    @test vec(sum(train[3],dims=2)+sum(test[3],dims=2)) == fill(11325,10)
    @test sum(train[2]) + sum(test[2]) == 11325
    @test all(train[1]' .== train[2])
    @test all(train[3]' .== train[2])
    @test all(test[1]' .== test[2])
    @test all(test[3]' .== test[2])
    train,test = splitobs((X1',Y1), obsdim=1)
    @test vec(sum(train[1],dims=1)) == fill(5565,10)
    @test vec(sum(test[1],dims=1)) == fill(5760,10)
    @test sum(train[2]) == 5565
    @test sum(test[2]) == 5760
    @test all(train[1] .== train[2])
    @test all(test[1] .== test[2])
    train,test = splitobs((sparse(X1),Y1),at=0.2)
    @test vec(sum(getobs(train[1]),dims=2)+sum(getobs(test[1]),dims=2)) == fill(11325,10)
    @test sum(train[2]) + sum(test[2]) == 11325
    @test all(getobs(train[1])' .== train[2])
    @test all(getobs(test[1])' .== test[2])
end

@testset "ObsView" begin
    A = ObsView(X)
    b, c = @inferred splitobs(A, .7)
    @test b isa ObsView
    @test c isa ObsView
    @test b == A[1:105]
    @test c == A[106:end]
end

@testset "BatchView" begin
    A = BatchView(X, 5)
    b, c = @inferred splitobs(A, .6)
    @test b isa BatchView
    @test c isa BatchView
    @test b == A[1:18]
    @test c == A[19:end]
end


text:
using DataFrames, DrWatson, LinearAlgebra, Parameters, Printf
using Query, SparseArrays, Statistics, StatsPlots

include(srcdir("load.jl"))

function results(col, fields::Symbol...)
    tempcol = Symbol(string(col) * "-sparse")

    df = collect_results(datadir("info", string(col)); verbose=false)
    disallowmissing!(df)

    select!(df, [:drug, :trial, :gh, :gw, :phase, fields..., col])

    df[!, :gridsize] = (df[:, :gh] .* df[:, :gw]).^2
    df[!, tempcol] = (sparse∘linearize).(df[:, col])
    df[!, :nnz] = nnz.(df[:, tempcol])
    df[!, :eigen] = (real ∘ last ∘ first ∘ eigen ∘ linearize).(df[:, col])

    select!(df, Not(col))
    rename!(df, Dict(tempcol => col))
    df
end

macro infoplot(df, xcol::QuoteNode, ycol::QuoteNode, xlabel::String, ylabel::String,
               legend::Union{QuoteNode,Bool}=:(:topright))
    quote
        kf = $df |> @filter(_.drug == "K-gluconate")
        knf = $df |> @filter(_.drug == "K-gluconate:Nifedipine")

        marker=(:circ, 4, 1.0)
        color=[3 4]
        α=0.5

        kplot = @df kf dotplot($xcol, $ycol, group=:phase, marker=marker, color=color, label="")
        @df kf violin!(kplot, $xcol, $ycol, group=:phase, α=α, color=color)

        knplot = @df knf dotplot($xcol, $ycol, group=:phase, marker=marker, color=color, label="")
        @df knf violin!(knplot, $xcol, $ycol, group=:phase, α=α)

        xmin, xmax = extrema($df[:,$xcol])

        plot(kplot, knplot, layout=2, title=["K-gluconate" "K-gluconate + Nifedipine"], titleloc=:left,
             xlabel=$xlabel, xticks=xmin:xmax, ylabel=$ylabel, size=(1000,500), legend=$legend)
    end
end

macro infoplot(xcol::QuoteNode, ycol::QuoteNode, xlabel::String, ylabel::String, legend::Union{QuoteNode,Bool}=:(:topright))
    :(df -> @infoplot df $xcol $ycol $xlabel $ylabel $legend)
end

function mutualinfo()
    @time mi = results(:mi, :lag)
    groups = groupby(mi, [:gh, :gw])

    for key in keys(groups)
        @unpack gh, gw = NamedTuple(key)

        @info "plotting mutual information for " gh gw

        basepath = projectdir("plots", "info", (@sprintf "%dx%d" gh gw), "mi")
        mkpath(basepath)

        g = groups[key]

        p = g |> @infoplot :lag :eigen "Lag" "Largest Eigenvalue"
        savefig(p, joinpath(basepath, "eigen.svg"))

        p = g |> @infoplot :lag :nnz "Lag" "Number of Nonzero Values"
        savefig(p, joinpath(basepath, "nnz.svg"))

        p = g |> @map({_.drug, _.trial, _.phase, _.lag, μ=mean(_.mi)}) |>
                 DataFrame |>
                 @infoplot :lag :μ "Lag" "Average MI"
        savefig(p, joinpath(basepath, "mean.svg"))

        p = g |> @map({_.drug, _.trial, _.phase, _.lag, μ=mean(nonzeros(_.mi))}) |>
                 @map({_.drug, _.trial, _.phase, _.lag, μ=isnan(_.μ) ? zero(_.μ) : _.μ}) |>
                 DataFrame |>
                 @infoplot :lag :μ "Lag" "Average Nonzero MI" :bottomright
        savefig(p, joinpath(basepath, "mean-nz.svg"))
    end
end

function transferentropy()
    @time te = results(:te, :k)
    groups = groupby(te, [:gh, :gw])

    for key in keys(groups)
        @unpack gh, gw = NamedTuple(key)

        @info "plotting tranfer entropy for " gh gw

        basepath = projectdir("plots", "info", (@sprintf "%dx%d" gh gw), "te")
        mkpath(basepath)

        g = groups[key]

        p = g |> @infoplot :k :eigen "History Length" "Largest Eigenvalue" :bottomright
        savefig(p, joinpath(basepath, "eigen.svg"))

        p = g |> @infoplot :k :nnz "History Length" "Number of Nonzero Values" :bottomright
        savefig(p, joinpath(basepath, "nnz.svg"))

        p = g |> @map({_.drug, _.trial, _.phase, _.k, μ=mean(_.te)}) |>
                 DataFrame |>
                 @infoplot :k :μ "History Length" "Average TE" :bottomright
        savefig(p, joinpath(basepath, "mean.svg"))

        p = g |> @map({_.drug, _.trial, _.phase, _.k, μ=mean(nonzeros(_.te))}) |>
                 @map({_.drug, _.trial, _.phase, _.k, μ=isnan(_.μ) ? zero(_.μ) : _.μ}) |>
                 DataFrame |>
                 @infoplot :k :μ "History Length" "Average Nonzero TE" :bottomright
        savefig(p, joinpath(basepath, "mean-nz.svg"))
    end
end

mutualinfo()
transferentropy()

max_stars_repo_name: ELIFE-ASU/dglmoore

max_issues_repo_name: ELIFE-ASU/dglmoore

max_forks_repo_name: ELIFE-ASU/dglmoore


text:
using PETSc2
using Base.Test

function not_isapprox(args...; kwargs...)
  return !isapprox(args...; kwargs...)
end

import MPI
println("pwd = ", pwd())
println("PETSC_DIR = ", ENV["PETSC_DIR"])
if !MPI.Initialized()
  MPI.Init()
end

comm = MPI.COMM_WORLD
comm_size = MPI.Comm_size(MPI.COMM_WORLD)
comm_rank = MPI.Comm_rank(MPI.COMM_WORLD)

# write your own tests here

# define some arrays to make looping easier
#vec_norms = [PETSC_NORM_1, PETSC_NORM_2, PETSC_NORM_INFINITY]

# skip 1 norm untill Petsc 1 complex 1 norm fix is released
vec_norms = [ NORM_2, NORM_INFINITY]
julia_vec_norms = [ 2, Inf]
# for some reason only the first four vector format works
#vec_formats = [VECSEQ, VECMPI, VECSTANDARD, VECSHARED, VECSEQCUSP, VECMPICUSP, VECCUSP, VECSEQVIENNACL, VECMPIVIENNACL, VECVIENNACL, VECNEST]  # no pthread

if comm_size > 1  # this is a parallel run
  vec_formats = [ VECMPI, VECSTANDARD]
else  # serial run
  vec_formats = [VECSEQ]
end

PetscInitialize(["-ksp_monitor","-malloc","-malloc_debug","-malloc_dump"]);
#PetscInitialize(["-ksp_monitor","-malloc","-malloc_debug"]);

# size of the system owned by this process (ie. local sys size)
sys_size_local = PetscInt(3)
sys_size_global = sys_size_local*comm_size

# indicies of the vector owned by this process
# create these with smallest precision, so they can be promoted
tmp2 = convert(Array{Float32,1}, collect(1.0:3))
tmp = convert(Array{Float32, 2}, [1.0 2.0 3; 4 5 7; 7 8 9])

tmp3 = convert(Array{Complex64, 1}, [1.0 + 0im; 2.0 + 1.0im; 3.0 + 2.0im])
tmp4 = convert( Array{Complex64, 2}, [1.0 + 1im   2 + 2im  3 + 3im; 4 + 4im  5 + 5im 7 + 7im; 7 + 7im 8 + 8im 9 + 9im])

# Have both local and global versions of system for testing
A_julia = zeros(PetscScalar, sys_size_local, sys_size_local)
rhs = zeros(PetscScalar, sys_size_local)
A_julia_global = zeros(PetscScalar, comm_size*sys_size_local, comm_size*sys_size_local)
rhs_global = zeros(PetscScalar, comm_size*sys_size_local)
# convert to arrays of the proper Petsc type
# this facilitates testing with the different Petsc build options

# define tolerances for KSP
if PetscReal == Float64
  rtol = PetscReal(1e-5)
  abstol = PetscReal(1e-14)
else  # single preciiosn
  rtol = PetscReal(1e-5)
  abstol = PetscReal(1e-6)
end
dtol = PetscReal(1e5)
maxits = PetscInt(50)

if PetscScalar <: Real

  for i=1:sys_size_local
    rhs[i] = convert(PetscScalar, tmp2[i])

    for j=0:(comm_size - 1)  # global
      idx = i + j*sys_size_local
      rhs_global[idx] = convert(PetscScalar, tmp2[i])
    end

  end

  for i=1:sys_size_local
    for j=1:sys_size_local
      A_julia[i,j] = convert(PetscScalar, tmp[i,j])

      for k=0:(comm_size - 1)
        idxm = i + k*sys_size_local
        idxn = j + k*sys_size_local
        A_julia_global[idxm, idxn] = convert(PetscScalar, tmp[i,j] )
      end
    end
  end

end  # end if statement

if PetscScalar <: Complex

  for i=1:sys_size_local
    rhs[i] = convert(PetscScalar, tmp3[i])

    for j=0:(comm_size - 1)  # global
      idx = i + j*sys_size_local
      rhs_global[idx] = convert(PetscScalar, tmp3[i])
    end

  end

  for i=1:sys_size_local
    for j=1:sys_size_local
      A_julia[i,j] = convert(PetscScalar, tmp4[i,j])

      for k=0:(comm_size - 1)
        idxm = i + k*sys_size_local
        idxn = j + k*sys_size_local
        A_julia_global[idxm, idxn] = convert(PetscScalar, tmp4[i,j] )
      end

    end
  end
end
#

if comm_rank == 0
  println("rhs_global = ", rhs_global)
  println("A_julia_global = ", A_julia_global)
end

#println("rhs = ", rhs)
#println("A_julia = ", A_julia)
x_julia = A_julia\rhs
#println("x_julia = ", x_julia)

include("test_vec.jl")
include("test_vec_interface.jl")
include("test_mat.jl")
include("test_mat_interface.jl")
include("test_ksp.jl")
include("test_options.jl")
#include("test_prealloc.jl")

PetscFinalize()


text:
module GenTeX

# include("example.jl")
include("latex.jl")
include("cache.jl")

end # module


text:
import RegularExpressions
using Documenter: doctest

doctest(RegularExpressions)


text:
include("../examples/baker.jl")
using LinearAlgebra
using JLD
using SharedArrays
using Distributed
function sens(s,nSteps)
    u = 2*pi*rand(2)
    u_trj = step(u, s, nSteps)
    x, y = view(u_trj,1,:), view(u_trj,2,:)
    du_trj = dstep(u_trj, s)
    ddu_trj = d2step(u_trj,s)
    vs = zeros(2)
    q = rand(2)
    q /= norm(q)

    pp = pert(u_trj, 2)
    q1 = rand(2)
    vs1 = zeros(2)

    dJds_st = 0.
    dJds_ust = 0.
    N = 12
    nSteps = nSteps + 1
    g = zeros(nSteps)
    J = cos.(4*y)

    # Df represents the derivative of f
    # on the
    # unstable manifold.
    # Xu = a q
    # d = 2
    Xu = zeros(2, nSteps)
    a = zeros(nSteps)
    Da = zeros(nSteps)
    Dq = zeros(2)
    Dvs = zeros(2)
    Dvs1 = zeros(2)
    s2 = s[2]
    # nSteps large number.
    for i = 1:nSteps-1
        dui = du_trj[:,:,i] # for large systems, can't store Jacobian.
        ppi = pp[:,i]
        xi, yi = x[i], y[i]
        #dppi = [cos(xi) 0; cos(xi)*sin(yi) sin(xi)*cos(yi)]
        sxi, cxi = sin(xi), cos(xi)
        syi, cyi = sin(2*yi), cos(2*yi)
        du11 = 2 + s2*syi*cxi/2
        du12 = s2*sxi*cyi
        dppi = [syi*cxi/2 cyi*sxi; 0. 0]*
        [1/du11 -2*du12/du11; 0 2.]
        q1 .= dui*q
        z = norm(q1)
        q1 ./= z

        z2 = z*z

        vs1 .= dui*vs + ppi
        a[i+1] = dot(vs1, q1)
        vs1 .-= a[i+1]*q1

        d2q = reshape([dot(ddu_trj[:,j,i],
                    q) for j=1:4],2,2)
        Dq = d2q*q/z2 + dui*Dq/z2
        dzdx = dot(z2*Dq, q1)
        Dq .= Dq .- dot(Dq,q1)*q1

        Dvs1 .= d2q*vs/z + dui*Dvs/z + dppi*q1/z
        Da[i+1] = dot(vs1, Dq) + dot(Dvs1, q1)
        Dvs1 .= Dvs1 - Da[i+1]*q1 - a[i+1]*Dq
        Delta_Da = dot(Dvs1, q1) + dot(vs1, Dq)
        Dvs1 .-= Delta_Da*q
        Da[i+1] += Delta_Da

        g[i+1] = g[i]/z - dzdx/z2
        dJdu = [0., -4*sin(4*yi)]
        dJds_st += dot(dJdu, vs)/nSteps

        vs .= vs1
        q .= q1
        Dvs .= Dvs1

    end
    for n = 1:N
        J_shift = J[n+1:end]
        nJ = length(J_shift)
        g_shift = g[2:nJ+1]
        Da_shift = Da[2:nJ+1]
        a_shift = a[2:nJ+1]
        dJds_ust -= dot(J_shift,Da_shift)/nJ
        dJds_ust -= dot(J_shift,g_shift.*a_shift)/nJ
    end
    @show dJds_st
    return dJds_st + dJds_ust
end
function get_sens(s)
    nSteps = 500000
    # J = cos(4y)
    n_exps = size(s)[2]
    dJds = zeros(n_exps)
    n_rep = 16
    dJds_proc = SharedArray{Float64}(n_rep)
    for k=1:n_exps
        sk = s[:,k]
        dJds_proc .= 0.
        t = @distributed for j=1:n_rep
            dJds_proc[j] = sens(sk,nSteps)/n_rep
        end
        wait(t)
        dJds[k] = sum(dJds_proc)
        @show dJds[k]
    end
    save("../data/sens/dJds_s2.jld", "s2",
         s[2,:], "dJds", dJds)
end

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
using sample
import other
include("ss")

include("another.jl")
include("subdir/yetanother.jl")

max_stars_repo_name: oxinabox/julia-vim-omnicomplete

max_issues_repo_name: oxinabox/julia-vim-omnicomplete

max_forks_repo_name: oxinabox/julia-vim-omnicomplete


text:
"""
Distributed Co-occurence (DisC) structure for
document embedding using word vectors.

# References
  * [Arora et al. ICLR 2018 "A compressed sensing
     view of unsupervised text embeddings, bag-on-n-grams
     and LSTMs"](https://openreview.net/pdf?id=B1e5ef-C-)
"""
struct DisCEmbedder{T,S} <: WordVectorsEmbedder{T,S}
    embeddings::EmbeddingsLibrary{T,S}
    n::Int
    config::NamedTuple
end

function DisCEmbedder(embeddings::EmbeddingsLibrary{T,S},
                      config;
                      n=DEFAULT_DISC_NGRAM,
                      kwargs...
                     ) where {T<:AbstractFloat, S<:AbstractString}
    return DisCEmbedder(embeddings, n, config)
end

# Dimensionality function
function dimensionality(embedder::DisCEmbedder)
    # Embedding corresponding to each n-gram embedding
    # concatenated vertically
    return embedder.n * size(embedder.embeddings)[1]
end

# Sentence embedding function - returns a `embedder.dim`×1 matrix
function sentences2vec(embedder::DisCEmbedder,
                       document_embedding::Vector{Matrix{T}};
                       kwargs...) where {T,S}
    if isempty(document_embedding)
        return zeros(T, dimensionality(embedder), 0)
    else
        return __dist_cooc(document_embedding,
                           embedder.n,
                           dimensionality(embedder))
    end
end

function __dist_cooc(document_embedding::Vector{Matrix{T}},
                     n::Int,
                     dim::Int
                    ) where {T<:AbstractFloat}
    # Pre-allocate output embeddings
    X = zeros(T, dim, length(document_embedding))
    m = Int(dim/n)
    @inbounds @simd for j in 1:length(document_embedding)
        for k in 1:n
            # For each column `j`, `m` values are filled by
            # the inner loop with the `k`-gram embeddings
            X[(k-1)*m+1:k*m, j] = __k_gram_prodsum_embedding(document_embedding[j], k)
        end
    end
    return X
end

function __k_gram_prodsum_embedding(A::Matrix{T}, k::Int) where {T}
    n = size(A, 2)		# number of embedded words
    k = clamp(k, 1, n)  # so it does not fail if too few embeddings
    X = zeros(T, size(A, 1), n-k+1)  # output
    col = 1
    @inbounds @simd for j in 1:n-k+1
        # Multiply `k` 'consecutive' embeddings element-wise
        X[:, col] = prod(A[:,j:j+k-1], dims=2)
        col+= 1
    end
    # Add all embedding products together element-wise and return
    return 1/k * vec(sum(X, dims=2))
end


text:
module FileManagement

export
    get_iters,
    get_time,
    get_final_time,
    get_grid,
    get_fields,
    get_statistics,
    get_wind_stress,
    get_surface_wave_parameters,
    get_parameter,
    get_multiple_parameters,
    set_from_file!,
    calculate_horizontal_average_timeseries,
    collect_horizontal_averages

using JLD2

function get_iters(filename)
    file = jldopen(filename)
    iters = parse.(Int, keys(file["timeseries/t"]))
    close(file)
    return iters
end

function get_time(filename, i)
    file = jldopen(filename)
    t = file["timeseries/t/$i"]
    close(file)
    return t
end

get_final_time(filename) = get_time(filename, get_iters(filename)[end])

function get_grid(filename)
    file = jldopen(filename)

    Nx = file["grid/Nx"]
    Ny = file["grid/Ny"]
    Nz = file["grid/Nz"]

    Lx = file["grid/Lx"]
    Ly = file["grid/Ly"]
    Lz = file["grid/Lz"]

    close(file)

    grid = RegularRectilinearGrid(size=(Nx, Ny, Nz), length=(Lx, Ly, Lz))

    return grid
end

function get_fields(filename, i)
    file = jldopen(filename)

    t = file["timeseries/t/$i"]
    u = file["timeseries/u/$i"]
    v = file["timeseries/v/$i"]
    w = file["timeseries/w/$i"]
    b = file["timeseries/b/$i"]

    close(file)

    return t, u, v, w, b
end

function get_statistics(filename, i)
    t, u, v, w, b = get_fields(filename, i)
    grid = get_grid(filename)
    b′, U, V, B, Bz, E, w², e = calculate_statistics(grid, u, v, w, b)
    return t, b′, U, V, B, Bz, E, w², e
end

function get_averages(filename, i)

    U = mean(u, dims=(1, 2))
    V = mean(v, dims=(1, 2))
    B = mean(b, dims=(1, 2))
    W² = mean(w.^2, dims=(1, 2))

    return t, U, V, W², B
end

function get_wind_stress(filename)
    file = jldopen(filename)
    τ = abs(file["boundary_conditions/Qᵘ₀"])
    close(file)
    return τ
end

function get_surface_wave_parameters(filename)
    file = jldopen(filename)

    aˢʷ = file["surface_waves/aˢʷ"]
    kˢʷ = file["surface_waves/kˢʷ"]

    close(file)

    return aˢʷ, kˢʷ
end

function get_parameter(filename, group, parameter_name)
    parameter = nothing

    jldopen(filename) do file

        if parameter_name ∈ keys(file["$group"])
            parameter = file["$group/$parameter_name"]
        end

    end

    return parameter
end

get_multiple_parameters(filename, group, names...) =
    Tuple(get_parameter(filename, group, name) for name in names)

function set_from_file!(model, filename; i=length(get_iters(filename)))

    file = jldopen(filename)

    # Load initial condition from file
    iter = parse(Int, keys(file["timeseries/t"])[i])

    u₀ = file["timeseries/u/$iter"]
    v₀ = file["timeseries/v/$iter"]
    w₀ = file["timeseries/w/$iter"]
    b₀ = file["timeseries/b/$iter"]

    close(file)

    array_type = typeof(model.velocities.u.data.parent)

    # Set initial condition
    model.velocities.u.data.parent .= array_type(u₀)
    model.velocities.v.data.parent .= array_type(v₀)
    model.velocities.w.data.parent .= array_type(w₀)
    model.tracers.b.data.parent .= array_type(b₀)

    return nothing
end

function calculate_horizontal_average_timeseries(directory; part=0)

    part > 9 && error("Part must be less than 10.")

    filenames = cd(() -> glob("*fields*"), directory)
    sortby(filename) = parse(Int, filename[end-5:end-5])
    sort!(filenames, by=sortby)
    part = part == 0 ? length(filenames) : part

    t, U, V, B, Bz, w², E = [[] for i=1:7]

    filepath = joinpath(directory, filenames[part])

    grid = get_grid(filepath)

    for filename in filenames
        filepath = joinpath(directory, filename)
        iters = get_iters(filepath)

        for iter in iters
            tᵢ, u, v, w, b = get_fields(filepath, iter)

            push!(t, tᵢ)

            b′, Uᵢ, Vᵢ, Bᵢ, Bzᵢ, Eᵢ, w²ᵢ, e = calculate_statistics(grid, u, v, w, b)

            push!(U, Uᵢ)
            push!(V, Vᵢ)
            push!(B, Bᵢ)
            push!(Bz, Bzᵢ)
            push!(w², w²ᵢ)
            push!(E, Eᵢ)
        end
    end

    ii = sortperm(t)

     t = t[ii]
     U = U[ii]
     V = V[ii]
     B = B[ii]
    Bz = Bz[ii]
    w² = w²[ii]
     E = E[ii]

    S = [sqrt.(Ui.^2 .+ Vi.^2) for (Ui, Vi) in zip(U, V)]

    return t, U, V, S, B, Bz, w², E
end

function collect_horizontal_averages(filename)
    iters = get_iters(filename)
    grid = get_grid(filename)

    cell_quantity_names = (:U, :V, :b, :U², :V², :E, :νₑ, :κₑ_b)
    cell_quantities = NamedTuple{cell_quantity_names}(Tuple(zeros(grid.Nz, length(iters))
                                                            for i = 1:length(cell_quantity_names)))

    face_quantity_names = (:W², :W³, :wu, :wv, :wb, :τ₁₃, :τ₂₃, :τ₃₃, :q₃_b)
    face_quantities = NamedTuple{face_quantity_names}(Tuple(zeros(grid.Nz+1, length(iters))
                                                            for i = 1:length(face_quantity_names)))

    t = zeros(length(iters))

    file = jldopen(filename)

    for (j, iter) in enumerate(iters)
        for i = 1:length(cell_quantities)
            c = cell_quantities[i]
            name = propertynames(cell_quantities)[i]
            c[:, j] .= file["timeseries/$name/$iter"][:][2:end-1]
        end

        for i = 1:length(face_quantities)
            f = face_quantities[i]
            name = propertynames(face_quantities)[i]
            f[:, j] .= file["timeseries/$name/$iter"][:][2:end]
        end

        t[j] = file["timeseries/t/$iter"]
    end

    close(file)

    return merge((t=t,), cell_quantities, face_quantities)
end

end # module


text:
todisplay = VERSION <= VersionNumber("1.2") ? true : false

using LinearAlgebra
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using CAGD
if todisplay
    using ViewerGL
    GL = ViewerGL
	include("../views.jl")
end

function twocubegrids(n=1, m=1, p=1, atol = 1e-6)
    V,(VV,EV,FV,CV) = Lar.cuboidGrid([n,m,p],true)
    mybox = (V,CV,FV,EV)

    twocubs = Lar.Struct([mybox, Lar.t(.3,.4,.5), Lar.r(pi/5,0,0), Lar.r(0,0,pi/12), mybox])

    V,CV,FV,EV = Lar.struct2lar(twocubs)
    if todisplay  GL.VIEW([ GL.GLGrid(V,FV, GL.COLORS[1]) ])  end

    cop_EV = convert(Lar.ChainOp, Lar.coboundary_0(EV::Lar.Cells))
    cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);
    model = CAGD.Model(V)
    CAGD.addModelCells!(model, 1, cop_EV)
    CAGD.addModelCells!(model, 2, cop_FE)

    split_model = CAGD.facesplitting(model)
    congr_model = CAGD.mergeModelVertices(split_model, signed_merge=true)
    if todisplay  displayModel(congr_model)  end
    gift_model = deepcopy(congr_model)
    FC, bicon_comps = CAGD.tgw(congr_model, 3, atol = atol)
    CAGD.addModelCells!(gift_model, 3, convert(Lar.ChainOp, FC'))
    if todisplay  viewExplode(gift_model)  end
end


text:
# TODO: User supplied elemtype??
ValueInfoProto(name::String, inshape, elemtype=Float32) =
ValueInfoProto(
    name=name,
    _type=TypeProto(
        tensor_type=TypeProto_Tensor(inshape, elemtype)
    )
)

TypeProto_Tensor(inshape, elemtype) = TypeProto_Tensor(
    elem_type=tp_tensor_elemtype(elemtype),
    shape=TensorShapeProto(inshape)
)
TypeProto_Tensor(::Missing, elemtype) = TypeProto_Tensor(
    elem_type=tp_tensor_elemtype(elemtype)
)

TensorShapeProto(shape) = TensorShapeProto(dim=[tsp_d(s) for s in reverse(shape)])
tsp_d(::Missing) = TensorShapeProto_Dimension()
tsp_d(n::Integer) = TensorShapeProto_Dimension(dim_value=n)
tsp_d(s::String) = TensorShapeProto_Dimension(dim_param=s)
tsp_d(s::Symbol) = tsp_d(string(s))

tp_tensor_elemtype(i::Integer) = i
tp_tensor_elemtype(::Missing) = TensorProto_DataType.UNDEFINED
tp_tensor_elemtype(::Type{Float32}) = TensorProto_DataType.FLOAT
tp_tensor_elemtype(::Type{Float64}) = TensorProto_DataType.DOUBLE

TensorProto(x::Number, name ="") = TensorProto([x], name)

TensorProto(t::AbstractArray{Float64,N}, name ="") where N = TensorProto(
    dims=collect(reverse(size(t))),
    data_type=TensorProto_DataType.DOUBLE,
    double_data = reshape(from_nnlib(t), :),
    name=name)

TensorProto(t::AbstractArray{Float32,N}, name ="") where N = TensorProto(
    dims=collect(reverse(size(t))),
    data_type=TensorProto_DataType.FLOAT,
    float_data = reshape(from_nnlib(t), :),
    name=name)

TensorProto(t::AbstractArray{Float16,N}, name ="") where N = TensorProto(t, TensorProto_DataType.FLOAT16, name)

TensorProto(t::AbstractArray{Int64,N}, name ="") where N = TensorProto(
    dims=collect(reverse(size(t))),
    data_type=TensorProto_DataType.INT64,
    int64_data = reshape(from_nnlib(t), :),
    name=name)

TensorProto(t::AbstractArray{Int32,N}, name ="") where N = TensorProto(
    dims=collect(reverse(size(t))),
    data_type=TensorProto_DataType.INT32,
    int32_data = reshape(from_nnlib(t), :),
    name=name)

TensorProto(t::AbstractArray{Int8,N}, name ="") where N = TensorProto(t, TensorProto_DataType.INT8, name)

TensorProto(t::AbstractArray, data_type::Int32, name) = TensorProto(
    dims=collect(reverse(size(t))),
    data_type=data_type,
    raw_data = reinterpret(UInt8, reshape(from_nnlib(t), :)),
    name=name)

AttributeProto(p::Pair) = AttributeProto(first(p), last(p))
AttributeProto(name::Symbol, v) = AttributeProto(string(name), v)

AttributeProto(name::String, i::Int64) = AttributeProto(
    name=name,
    _type = AttributeProto_AttributeType.INT,
    i = i
)

AttributeProto(name::String, f::Float32) = AttributeProto(
    name=name,
    _type = AttributeProto_AttributeType.FLOAT,
    f = f
)

AttributeProto(name::String, floats::AbstractVector{Float32}) = AttributeProto(
    name=name,
    _type = AttributeProto_AttributeType.FLOATS,
    floats = floats
)

AttributeProto(name::String, f::Float64) = AttributeProto(
    name=name,
    _type = AttributeProto_AttributeType.FLOAT,
    f = Float32(f)
)

AttributeProto(name::String, i::NTuple{N, Int64}) where N = AttributeProto(name, collect(i))

AttributeProto(name::String, i::AbstractVector{Int64}) = AttributeProto(
    name=name,
    _type = AttributeProto_AttributeType.INTS,
    ints = i
)

AttributeProto(name::String, str::AbstractString) = AttributeProto(
    name=name,
    _type = AttributeProto_AttributeType.STRING,
    s = Vector{UInt8}(str)
)

AttributeProto(name::String, strings::AbstractVector{<:AbstractString}) = AttributeProto(
    name=name,
    _type = AttributeProto_AttributeType.STRINGS,
    strings = Vector{UInt8}.(strings)
)

AttributeProto(name::String, tensor::TensorProto) = AttributeProto(
    name=name,
    _type = AttributeProto_AttributeType.TENSOR,
    t = tensor
)


text:
module test_jive_mockup

using Jive # Mock @mockup
using Test

module Goods
struct Foo
end
function f(::Foo)
    10
end
function g(::Foo)
    10
end
end # module Goods

@mockup module Goods
function f(::Foo)
    20
end
end # @mockup module Goods
@test Goods.f(Goods.Foo()) == 10
@test Mock.Goods.f(Mock.Goods.Foo()) == 20
@test Mock.Goods.g(Mock.Goods.Foo()) == 10

Goods3 = @mockup module Goods
function g(::Foo)
    30
end
end # @mockup module Goods
@test Goods.f(Goods.Foo()) == 10
@test Mock.Goods.f(Mock.Goods.Foo()) == 10
@test Mock.Goods.g(Mock.Goods.Foo()) == 30
@test Goods3 isa Module
@test Goods3.g === Mock.Goods.g

end # module test_jive_mockup


text:
using ITensors

function run(; maxdim::Int,
               nsweeps::Int = 5,
               outputlevel::Int = 0,
               conserve_qns::Bool = false,
               use_sparse_mpo::Bool = false,
               N::Int = 100)
  sites = siteinds("S=1", N; conserve_qns = conserve_qns)
  ampo = AutoMPO()
  for j in 1:N-1
    ampo .+= 0.5, "S+", j, "S-", j+1
    ampo .+= 0.5, "S-", j, "S+", j+1
    ampo .+=      "Sz", j, "Sz", j+1
  end
  HD = MPO(ampo,sites)
  if use_sparse_mpo
    println("Splitting blocks of H (= making MPO sparse)")
    H = splitblocks(linkinds, HD)
  else
    println("Using dense MPO")
    H = HD
  end
  psi0 = productMPS(sites, n -> isodd(n) ? "↑" : "↓")
  sweeps = Sweeps(nsweeps)
  maxdims = min.(maxdim, [10, 20, 100, maxdim])
  maxdim!(sweeps, maxdims...)
  cutoff!(sweeps, 0.0)
  t = @elapsed begin
  energy, ψ = dmrg(H, psi0, sweeps;
                   outputlevel = outputlevel)
  end
  println(t)
  return
end


text:
## ----------------------------------------------------------------------------
# marginal distributions
let
    objider = iJR.BIOMASS_IDER
    size = [300, 250]
    method2 = ME_MAX_POL

    # Iders
    model_iders, Hd_iders = [iJR.BIOMASS_IDER], ["D"]
    for Hd_met in CONC_IDERS
        model_met = Hd_mets_map[Hd_met]
        model_exch = Hd_rxns_map[Hd_met]
        push!(model_iders, model_exch)
        push!(Hd_iders, string("u", Hd_met))
    end

    for (model_ider, Hd_ider) in zip(model_iders, Hd_iders)
        ps = Plots.Plot[]
        ps2 = Plots.Plot[]
        for exp in EXPS
            p = plot(title = string(Hd_ider, " exp: ", exp))
            p2 = plot(title = string(Hd_ider, " exp: ", exp))
            margin, m, M = -Inf, Inf, -Inf
            Hd_av = Hd.val(Hd_ider, exp)

            # EP
            for method in ALL_METHODS
                color = method_colors[method]

                datfile = ME_INDEX[method, :DFILE, exp]
                dat = deserialize(datfile)
                model = dat[:model]
                objidx = ChU.rxnindex(model, objider)
                exp_beta = dat[:exp_beta]
                epouts = dat[:epouts]
                epout = epouts[exp_beta]
                ep_av = ChU.av(model, epout, model_ider)
                ep_va = sqrt(ChU.va(model, epout, model_ider))
                # fbaout = dat[:fbaout]

                # ChP.plot_marginal!(p, model, [epout, fbaout], model_exch; legend = false)
                ChP.plot_marginal!(p, model, [epout], model_ider;
                    legend = false, color, alpha = 0.6, lw = 5)

                m = minimum([m, ep_av, Hd_av])
                M = maximum([M, ep_av, Hd_av])
                margin = maximum([margin, 3 * ep_va])

                if method == method2
                    for (beta, epout) in sort(epouts; by = first)
                        ep_av = ChU.av(model, epout, model_ider)
                        ep_va = sqrt(ChU.va(model, epout, model_ider))

                        alpha = 0.15
                        color = method_colors[method]
                        ChP.plot_marginal!(p2, model, epout, model_ider;
                            legend = false, color, alpha, lw = 1)

                        if beta == exp_beta
                            # @info "At" Hd_ider exp
                            ChP.plot_marginal!(p2, model, epout, model_ider;
                                legend = false, color,
                                alpha = 1.0, lw = 3
                            )
                            break
                        end
                    end
                    push!(ps2, p2)
                end

            end
            # Experimental
            vline!(p, [Hd_av]; label = "", lw = 6, color = :black, alpha = 0.3)
            vline!(p2, [Hd_av]; label = "", lw = 6, color = :black, alpha = 0.3)

            plot!(p; xlim = [m - margin, M + margin], size)
            plot!(p2; xlim = [m - margin, M + margin], size)
            push!(ps, p)
        end

        for k in [:xi, :D, :sGLC]
            p = plot(;title = Hd_ider, size)
            xticks =  (EXPS, string.(EXPS))
            p = bar!(p, EXPS, Hd.val(k); title = k, label = "", xticks)
            push!(ps, p)
            push!(ps2, p)
        end

        pname = string(Hd_ider, "_marginals")
        mysavefig(ps, pname)
        pname = string(Hd_ider, "_marginals_vs_beta")
        mysavefig(ps2, pname)
    end

end

## ----------------------------------------------------------------------------
# marginals v2
let
    objider = iJR.BIOMASS_IDER
    size = [300, 250]

    # Iders
    model_iders, Hd_iders = [iJR.BIOMASS_IDER], ["D"]
    for Hd_met in CONC_IDERS
        model_met = Hd_mets_map[Hd_met]
        model_exch = Hd_rxns_map[Hd_met]
        push!(model_iders, model_exch)
        push!(Hd_iders, string("u", Hd_met))
    end

    for (model_ider, Hd_ider) in zip(model_iders, Hd_iders)
        marg_params = (;xlabel = string(Hd_ider), yaxis = nothing, ylabel = "prob", titlefont = 8)

        epps = Plots.Plot[]
        exps = Plots.Plot[]
        for method in ALL_METHODS
            expp = plot(;title = string("Experimental"), marg_params...)
            epp = plot(;title = string(" MaxEnt: \n", method), marg_params...)
            margin, m, M = -Inf, Inf, -Inf

            # EP
            for exp in EXPS
                Hd_av = Hd.val(Hd_ider, exp)
                color = exp_colors[exp]

                datfile = ME_INDEX[method, :DFILE, exp]
                dat = deserialize(datfile)
                model = dat[:model]
                objidx = ChU.rxnindex(model, objider)
                epouts = dat[:epouts]
                exp_beta = maximum(keys(epouts))
                epout = epouts[exp_beta]
                ep_av = ChU.av(model, epout, model_ider)
                ep_va = sqrt(ChU.va(model, epout, model_ider))

                # bounds = ChU.bounds(model, model_ider)
                ChP.plot_marginal!(epp, model, [epout], model_ider;
                    legend = false, color, alpha = 0.8, lw = 3
                )

                m = minimum([m, ep_av, Hd_av])
                M = maximum([M, ep_av, Hd_av])
                margin = maximum([margin, 3 * ep_va])

                # Experimental
                vline!(expp, [Hd_av]; label = "", lw = 3, color, alpha = 0.8)

            end

            map([expp, epp]) do p
                plot!(p; xlim = [m - margin, M + margin], size)
            end

            push!(epps, epp)
            push!(exps, expp)
        end

        extras = Plots.Plot[]
        for k in [:xi, :D, :sGLC]
            p = plot(;title = "Experimental", size,
                xlabel = "rep", ylabel = string(k), titlefont = 8
            )
            xticks =  (EXPS, string.(EXPS))
            vals = [Hd.val(k, exp) for exp in EXPS]
            color = [exp_colors[exp] for exp in EXPS]
            p = bar!(p, EXPS, vals; label = "", xticks, color)
            push!(extras, p)
        end

        ps = Plots.Plot[exps; epps; extras]
        layout = (3, 3)
        pname = string(Hd_ider, "_marginals_v2")
        mysavefig(ps, pname; layout)

    end # for (model_ider, Hd_ider)

end


text:
using TensorTrains

w=2
L=5
h,V = hubbard_2D(w,L;U=4,w_pbc=false,L_pbc=false)
H_tto = hV_to_mpo(h,V,tol=1e-5)
ψ_0 = half_filling(w*L)
E_als,ψ_r = als_eigsolv(H_tto,ψ_0,sweep_schedule=[2,5],rmax_schedule=[16,32])
E_mals,ψ_tt,r = mals_eigsolv(H_tto,ψ_r,tol=1e-6,sweep_schedule=[2,5],rmax_schedule=[32,64])
#Fiedler order
σ_fiedler = fiedler_order(ψ_r)
h_f = h[σ_fiedler,σ_fiedler]
V_f = V[σ_fiedler,σ_fiedler,σ_fiedler,σ_fiedler]
H_tto = hV_to_mpo(h_f,V_f,tol=1e-5)
E_fiedler,ψ_fiedler,r = mals_eigsolv(H_tto,ψ_r,tol=1e-6,sweep_schedule=[3],rmax_schedule=[64])
#BWPO order
σ_bwpo = bwpo_order(ψ_r)
h_b = h[σ_bwpo,σ_bwpo]
V_b = V[σ_bwpo,σ_bwpo,σ_bwpo,σ_bwpo]
H_tto = hV_to_mpo(h_b,V_b,tol=1e-5)
E_bwpo,ψ_bwpo,r = mals_eigsolv(H_tto,ψ_r,tol=1e-6,sweep_schedule=[3],rmax_schedule=[64])

max_stars_repo_name: msdupuy/Tensor-Train-Julia

max_issues_repo_name: msdupuy/Tensor-Train-Julia

max_forks_repo_name: msdupuy/Tensor-Train-Julia


text:
using Statistics, SpecialFunctions

"""
    tnmom1c(c, a, b)

Computes ⟨(x - c)⟩ for the standard normal distribution truncated to the
interval [a, b]. This is more numerically stable than computing the mean
of x and subtracting c. Returns NaN if a > b.
"""
function tnmom1c(c::Real, a::Real, b::Real)
    if !(a ≤ b) || isnan(c)
        return oftype(middle(a, b), NaN)
    elseif a == b || isinf(c)
        return middle(a, b) - c
    elseif abs(a) > abs(b)
        return -tnmom1c(-c, -b, -a)
    elseif c ≤ 0
        #= At this point tnmean(a,b) ≥ 0.
        Therefore the subtraction is fine. =#
        return tnmean(a, b) - c
    elseif isinf(a) && isinf(b)
        return oftype(middle(a, b), -c)
    end

    @assert a < b
    @assert a ≤ 0 ≤ b || 0 ≤ a ≤ b
    @assert c > 0

    exΔm1 = expm1(middle(a, b) * (a - b))

    if a ≤ 0 ≤ b
        exa = √(2/π) * exp(-a^2 / 2)
        exb = √(2/π) * exp(-b^2 / 2)
        exa ≈ exb && return -c
        era = erf(a / √2)
        erb = erf(b / √2)
        era ≈ erb && return middle(a - c, b - c)

        m1 = (exa - exb) / (era - erb)
        c ≈ m1 && return zero(m1)

        fa = exa + c * era
        fb = exb + c * erb
        return (fa - fb) / (erb - era)
    elseif 0 < a < b
        era = erfcx(a / √2)
        erb = erfcx(b / √2)
        if xerfcx_asym_thresh(a / √2, 3)
            fa = (1 - c/a) * √(2/π) + √2 * c/a * xerfcx_pi(a / √2)
        else
            fa = √(2/π) - c * era
        end
        if xerfcx_asym_thresh(b / √2, 3)
            fb = (1 - c/b) * √(2/π) + √2 * c/b * xerfcx_pi(b / √2)
        else
            fb = √(2/π) - c * erb
        end
        fa = (1 - c/a) * √(2/π) + √2 * c/a * xerfcx_pi(a / √2)
        fb = (1 - c/b) * √(2/π) + √2 * c/b * xerfcx_pi(b / √2)
        exΔ = exp(middle(a, b) * (a - b))
        return (fa - fb * exΔ) / (era - erb * exΔ)
    end
end

"""
    tnmom1c(c, a, b, μ, σ)

Computes ⟨(x - c)⟩ for the normal distribution with mean μ and standard
deviation σ, truncated to the interval [a, b]. This is more numerically stable
than computing the mean of x and subtracting c.
"""
function tnmom1c(c, a, b, μ, σ)
    α = (a - μ) / σ
    β = (b - μ) / σ
    γ = (c - μ) / σ
    return σ * tnmom1c(γ, α, β)
end


text:
import Base.convert

export BED

type BED <: Region
    contig :: String
    range :: Range1{Int}
    name :: String
    score :: Float64
    strand :: Char
end

register_file_extension!(BED, ".bed")
register_file_extension!(BED, ".bed.gz")

function Base.read(io :: IO, ::Type{BED})
    while !eof(io)
        try
            return convert(BED, readline(io))
        catch
            continue
        end
    end
    throw(EOFError())
end

function Base.convert(::Type{BED}, s :: String)
    fields = split(strip(s), '\t')
    name = if length(fields)>3 fields[4] else "." end
    score = if length(fields)>4 float(fields[5]) else NaN end
    strand = if length(fields)>5 fields[6][1] else '.' end
    BED(fields[1], int(fields[2]):int(fields[3]),
                   name,
                   score,
                   convert(Strand, strand))
end

function Base.convert(::Type{String}, bed :: BED)
    "$(bed.contig)\t$(bed.range[1])\t$(bed.range[end])\t$(bed.name)\t$(bed.score)\t$(bed.strand)\n"
end


text:
@testset "PosteriorGP" begin
    rng = MersenneTwister(123456)
    N_cond = 3
    N_a = 5
    N_b = 6

    # Specify prior.
    m = CustomMean(sin)
    k = Matern32Kernel()
    f = GP(m, k)

    # Sample from prior.
    x = collect(range(-1.0, 1.0; length=N_cond))
    fx = f(x, 1e-15)
    y = rand(rng, fx)

    # Construct posterior.
    f_post = posterior(fx, y)

    # Verify that posterior collapses around observations.
    @test mean(f_post, x) ≈ y
    @test cov_diag(f_post, x) ≈ zeros(N_cond) rtol=1e-14 atol=1e-14

    # Check interface is implemented fully and consistently.
    a = collect(range(-1.0, 1.0; length=N_a))
    b = randn(rng, N_b)
    abstractgp_interface_tests(f_post, a, b)

    #Check sequential conditioning posterior
    X = rand(5)
    y = rand(5)

    f1 = GP(SqExponentialKernel())
    p_fx1 = posterior(f1(X[1:3], 0.1), y[1:3])
    p_p_fx1 = posterior(p_fx1(X[4:5], 0.1), y[4:5])

    f2 = GP(SqExponentialKernel())
    p_fx2 = posterior(f2(X, 0.1), y)

    @test p_p_fx1.data.C.U ≈ p_fx2.data.C.U atol=1e-5
    @test p_p_fx1.data.α ≈ p_fx2.data.α atol=1e-5
    @test p_p_fx1.data.x ≈ p_fx2.data.x atol=1e-5
    @test p_p_fx1.data.δ ≈ p_fx2.data.δ atol=1e-5

end


text:
using OmegaGrav
using MatterPower
using PyCall
using Dierckx
using OrdinaryDiffEq
using Tables, CSV
# %% Call the python wrapper for CLASS, `classy`, via PyCall
classy = pyimport("classy")
# Create an instance of the CLASS wrapper
cosmo = classy.Class()
# Create a dictionary of the cosmological parameters
deg_ncdm = true # 3 massive neutrinos with degenerate mass?
params = Dict(
   "output" => "mPk",
   "P_k_max_h/Mpc" => 100,
   "z_pk" => 10,
   "non linear" => "halofit",
   "A_s" => 2.097e-9,
   "n_s" => 0.9652,
   "k_pivot" => 0.05,
   "h" => 0.6737,
   "omega_b" => 0.02233,
   "omega_cdm" => 0.1198,
   "N_ncdm" => 1,
)
if deg_ncdm
   push!(params, "m_ncdm" => 0.02, "deg_ncdm" => 3, "N_ur" => 0.00641)
   mν = params["m_ncdm"] * params["deg_ncdm"]
else
   push!(params, "m_ncdm" => 0.06, "N_ur" => 2.0328)
   mν = params["m_ncdm"]
end
h0 = params["h"]
Ωcb = (params["omega_b"] + params["omega_cdm"]) / h0^2
Ων = mν / 93.14 / h0^2 #  Ωνh^2 = ∑mν / (93.14 eV)
Ωm = Ωcb + Ων
# Set the parameters and run the CLASS code
cosmo.set(params)
cosmo.compute()

# %% Compute the analytical linear theory result for K/W
# Reference: Equation (12) of Davis, Miller, White, ApJ, 490, 63 (1997).
zini = 9
a = 1/(1+zini):0.01:1.0 # scale factor
nred = length(a)
sol_δθ = setup_growth(Ωm, 1 - Ωm)
KovW = zeros(nred)
for i = 1:nred
   KovW[i] = -2 / 3 * (sol_δθ(a[i])[2] / sol_δθ(a[i])[1])^2 * a[i] / Ωm
end

# %% Compute `Ωgrav = Ωm * W / 2` at various redshifts
Ωgpk = zeros(nred)
Ωglin = zeros(nred)
for ired = 1:nred
   z = 1 / a[ired] - 1
   # Define functions to return linear and non-linear power spectra
   # Note: The CLASS code takes wavenumbers in units of 1/Mpc (no h) and
   # return power spectra in units of Mpc^3 (no 1/h^3).
   pknl_class(kovh) = cosmo.pk(kovh * h0, z) * h0^3
   pklin_class(kovh) = cosmo.pk_lin(kovh * h0, z) * h0^3
   # Spline interpolate in log(k)
   lnk = log(1e-4):0.05:log(100)
   pknl = Spline1D(lnk, pknl_class.(exp.(lnk)))
   pklin = Spline1D(lnk, pklin_class.(exp.(lnk)))
   # %% Compute Ωgrav from non-linear total matter P(k)
   Ωgpk[ired] = ograv_pk(x -> pknl(log(x)), z, Ωm)
   # %% Compute Ωgrav from linear total matter P(k)
   Ωglin[ired] = ograv_pk(x -> pklin(log(x)), z, Ωm)
end
Ωgrav_pklin = Spline1D(a, Ωglin)
Ωgrav_pknl = Spline1D(a, Ωgpk)

# %% Solve the Layzer-Irvine equation
tspan = (1 / (1 + zini), 1.0)
t0 = tspan[1]
## Non-linear P(k)
W(x) = 2 * Ωgrav_pknl(x)
dWda(x) = 2 * derivative(Ωgrav_pknl, x)
f(u, p, t) = -(2 * u + W(t)) / t - dWda(t)
u0 = KovW[1] * W(t0)
prob = ODEProblem(f, u0, tspan)
sol_pknl = solve(prob, Tsit5())
## Linear P(k): This result must agree with the analytical linear theory result
W(x) = 2 * Ωgrav_pklin(x)
dWda(x) = 2 * derivative(Ωgrav_pklin, x)
f(u, p, t) = -(2 * u + W(t)) / t - dWda(t)
u0 = KovW[1] * W(t0)
prob = ODEProblem(f, u0, tspan)
sol_pklin = solve(prob, Tsit5())

#%% Save table data to table2.csv
redshift = [0, 0.3, 0.5, 0.7, 1, 1.3, 1.5]
a = 1 ./ (1 .+ redshift)
t = Tables.table([redshift sol_pklin.(a) sol_pknl.(a)])
header = ["z", "Omega_K_pklin", "Omega_K_pknl"]
CSV.write("table2.csv", t, header = header)

# %% Clean CLASS (the equivalent of the struct_free() in the `main`
# of CLASS. This step is primordial when running in a loop over different
# cosmologies, as you will saturate your memory very fast if you ommit
# it.
cosmo.struct_cleanup()
# If you want to change completely the cosmology, you should also
# clean the arguments, otherwise, if you are simply running on a loop
# of different values for the same parameters, this step is not needed
cosmo.empty()


text:
# This file reproduces figures 2-3 of
# Identification of LTV Dynamical Models with Smooth or Discontinuous Time Evolution by means of Convex Optimization
# Fredrik Bagge Carlson, Anders Robertsson, Rolf Johansson
# https://arxiv.org/abs/1802.09794
# The script requires DynamicMovementPrimitives.jl to be installed

using OrdinaryDiffEq
using LTVModels

using ValueHistories, IterTools, MLDataUtils, OrdinaryDiffEq, Parameters, InteractNext

# Import robot simulator from DynamicMovementPrimitives
include(Pkg.dir("DynamicMovementPrimitives","src","two_link.jl"))
using TwoLink

@with_kw struct TwoLinkSys
    N  = 1000
    n  = 2
    ns = 2n
    h = 0.02
    σ0 = 0
    sind = 1:ns
    uind = ns+1:(ns+n)
    s1ind = (ns+n+1):(ns+n+ns)
end

# Problem parameters and trajectory
n,m,N,h,σ0 = 4,3,200,0.02,0.0
np         = n*(n+m)
sys        = TwoLinkSys(N=N, h=h, σ0 = σ0)#; warn("Added noise")
seed       = 1

ni         = N+2
traj       = connect_points([0 -0.5; 0.9 -0.1],ni)
q          = inverse_kin(traj[1],:up)
p          = forward_kin(q)
qd         = centraldiff(q')'/sys.h
qdd        = centraldiff(qd')'/sys.h
u0         = torque(q,qd,qdd) # Computed torque
u0       .+= 0.1randn(size(u0)) # Random noise on input

# Define callback to simulate the robot hitting a wall
const wall = -0.2
function affect!(i)
    q = i.u[1:2]
    qd = i.u[3:4]
    p1,p2 = forward_kin(q)[2]
    J = TwoLink.jacobian(q)
    v1,v2 = J*qd
    d = [0, wall-p2]
    δq = J\d
    i.u[3:4] = J\[v1, 0]
    i.u[1:2] .+= δq
end
condition(x,t,i) = forward_kin(x)[2][2] > wall
cb = DiscreteCallback(condition,affect!)

# Simulate robot using ODE-solver
t      = 0:h:N*h
q0     = [q[:,1]; qd[:,1]]
prob   = OrdinaryDiffEq.ODEProblem((x,p,t)->time_derivative(x, u0[:,floor(Int,t/h)+1]),q0,(t[[1,end]]...,))
sol    = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8, callback=cb)
x      = hcat(sol(t)...)
u      = u0[:,1:N]
u      = [u; ones(1,size(u,2))]
x,y    = x[:,1:N], x[:,2:N+1]

plot(forward_kin(x[1:2,:])', layout=2, subplot=1, title="End-effector positions")
plot!([1,N],[wall, wall], l=(:dash, :black), grid=false, subplot=1, lab="Wall")
plot!(u', subplot=2, title="Control torques")

# Normalize input data to mean 0 and unit variance
sx = std(x,dims=2)
su = std(u[1:2,:])
mx, mu = mean(x,dims=2), mean(u[1:2,:],dims=2)
xt = (x.-mx) ./ sx
ut = copy(u)
ut[1:2,:] = (ut[1:2,:].-mu) ./ su

# Fit model using ADMM
model = LTVModels.fit_statespace_admm(xt,ut,10^-0.5, extend=false, tol=1e-8, iters=20000, D=1, zeroinit=false)

signch = findfirst(u[2,:] .< 0) # Detect when control signal changes sign
contact = findfirst(forward_kin(x[1:2,:])[2,:] .> wall) # Detect first contact with wall

fig=plot(xt', label="States", show=false, layout=grid(3,2), subplot=1:4, size=(1800,1000), title="State trajectories", title=["\$q_1\$" "\$q_2\$" "\$\\dot{q}_1\$" "\$\\dot{q}_2\$"], c=:blue, linewidth=2)
vline!([signch].*ones(1,4), show=false, subplot=1:4, lab="Velocity sign change", l=(:black, :dash))
vline!([contact].*ones(1,4), show=false, subplot=1:4, lab="Stiff contact", l=(:black, :dash))
plot!(forward_kin(x[1:2,:])', subplot=5, title="End-effector positions", lab=["x" "y"], c=[:blue :red])
hline!([wall], l=(:dash, :black), grid=false, subplot=5, lab="Constraint")
vline!([contact], subplot=5, lab="Stiff contact", l=(:black, :dash))
plot!(ut[1:2,:]', subplot=6, title="Control torques", lab="", c=[:blue :red])
vline!([signch], subplot=6, lab="Velocity sign change", l=(:black, :dash))
plot!(LTVModels.simulate(model, xt[:,1], ut)', label="Simulation", show=false, subplot=1:4, l=(:red,), legend=[false false :left false])
gui()

# Validation data. Simulate system again with different input trajectory
uv   = torque(q,qd,qdd)
uv .+= 0.1randn(size(uv))
prob   = OrdinaryDiffEq.ODEProblem((x,p,t)->time_derivative(x, uv[:,floor(Int,t/h)+1]),q0,(t[[1,end]]...,))
sol   = solve(prob,Tsit5(),reltol = 1e-8,abstol = 1e-8, callback = cb)
xv    = hcat(sol(t)...)
uv    = uv[:,1:N]
uv    = [uv; ones(1,size(uv,2))] # Extend to estimate friction
xv,yv = xv[:,1:N], xv[:,2:N+1]

# Normalize data with same parameters as during training
xt2        = (xv.-mx) ./ sx
ut2        = copy(uv)
ut2[1:2,:] = (uv[1:2,:].-mu) ./ su

plot(xt2', label="States", show=false, layout=grid(2,2), subplot=1:4, size=(1800,1000), title="State trajectories", title=["\$q_1\$" "\$q_2\$" "\$\\dot{q}_1\$" "\$\\dot{q}_2\$"], c=:blue, linewidth=2, reuse=false, grid=false)
vline!([signch].*ones(1,4), show=false, subplot=1:4, lab="Original Velocity sign change", l=(:black, :dash))
vline!([contact].*ones(1,4), show=false, subplot=1:4, lab="Original Stiff contact", l=(:black, :dash))

plot!(LTVModels.simulate(model, xt2[:,1], ut2)', label="Simulation", show=false, subplot=1:4, l=(:red,), legend=[false false :left false])
gui()


text:
include("IOM/src/share/LogSystem.jl")
include("IOM/src/share/PolelikeCoordinate.jl")
include("IOM/src/models/EMOM/ENGINE_EMOM.jl")
include("IOM/src/driver/driver_working.jl")
include("IOM/src/share/CyclicData.jl")
include("DataLoader.jl")
include("func_reinitModel.jl")
include("func_data2SOM.jl")

using MPI
using CFTime, Dates
using ArgParse
using TOML
using DataStructures
using NCDatasets
using .PolelikeCoordinate
using .LogSystem
using .CyclicData
using .DataLoader

function parse_commandline()

    s = ArgParseSettings()

    @add_arg_table s begin

        "--config-file"
            help = "Configuration file."
            arg_type = String
            required = true

        "--hist-dir"
            help = "Annual forcing file. It should contain: TAUX, TAUY, SWFLX, NSWFLX, VSFLX"
            arg_type = String
            required = true

        "--year-rng"
            help = "The year range that the user wants to run. The begin date will be set to 12/31 of the year before the first year."
            nargs = 2
            required = true
            arg_type = Int64

        "--SOM"
            help = "If set then HMXL will be set to the deepest z_w_bot. Remember to double check Nz_bot file."
            arg_type = Bool
            default = false

    end

    return parse_args(s)
end

parsed = parse_commandline()

MPI.Init()
comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)
is_master = rank == 0

config = nothing
SOM_HMXL = nothing
if is_master

    year_rng = parsed["year-rng"]

    config = TOML.parsefile(parsed["config-file"])
    config["DRIVER"]["sync_thermo_state_before_stepping"] = true

    # We will strictly use a day
    Δt_float = 86400.0
    Δt = Dates.Second(Δt_float)
    read_restart = false

    cfgmc = config["MODEL_CORE"]
    cfgmm = config["MODEL_MISC"]

    cfgmc["Qflx"] = "off"
    cfgmc["weak_restoring"] = "off"
    cfgmc["transform_vector_field"] = false

    Nz = length(cfgmc["z_w"])-1

    first_run = true

    if parsed["SOM"]
        global SOM_HMXL = minimum(cfgmc["z_w"])
        cfgmc["Ks_V"] = 0.0

        println("`SOM` also set convective_adjustment = off and Ks_V = 0.0")
        println("HMXL will always be set to $(SOM_HMXL)")

    end

    dli = DataLoader.DataLoaderInfo(
        hist_dir = parsed["hist-dir"],
        casename = "paper2021_POP2_CTL",
        state_varnames = ["TEMP", "SALT", "HMXL",],
        forcing_varnames = ["SFWF", "TAUX", "TAUY", "SHF", "SHF_QSW",],
        year_rng = year_rng,
        layers = 1:Nz,
    )

    pred_steps = 5
    pred_cnt = 0

    t_start = DateTimeNoLeap(year_rng[1], 1, 1, 0, 0, 0) - pred_steps*Δt
    t_end   = DateTimeNoLeap(year_rng[2], 1, 1, 0, 0, 0)

end

global data = nothing
global QFLXT = nothing
global QFLXS = nothing
global Nz_bot = nothing

coupler_funcs = (

    master_before_model_init = function()

        return read_restart, t_start
    end,

    master_after_model_init! = function(OMMODULE, OMDATA)
            # setup forcing
        println("# Model beg time: $(string(t_start))")
        println("# Model end time: $(string(t_end))")

        if parsed["SOM"]
            println("SOM specific: Load Nz_bot for later usage.")
            global Nz_bot = OMDATA.mb.ev.topo.Nz_bot_sT
        end

    end,

    master_before_model_run! = function(OMMODULE, OMDATA)

        global first_run, data, SOM_HMXL

        if first_run

            first_run = false

            # only put the TEMP and SALT of
            # the day before start date into the model

            data_init = DataLoader.loadInitDataAndForcing(dli, OMDATA.clock.time)
            if parsed["SOM"]
                println("SOM would propagate SST and SSS to the whole column")
                data2SOM!(data_init["TEMP"], Nz_bot)
                data2SOM!(data_init["SALT"], Nz_bot)
            end

            reinitModel!(
                OMDATA,
                data_init;
                forcing = true,
                thermal = true,
                SOM_HMXL = SOM_HMXL,
            )

            global QFLXT = copy(data_init["TEMP"])
            global QFLXS = copy(QFLXT)
            QFLXT .= 0.0
            QFLXS .= 0.0

            global pred_reinit = false
            global pred_cnt = 0

        end

        write_restart = OMDATA.clock.time == t_end
        end_phase = OMDATA.clock.time > t_end

        if ! end_phase
            return_values = ( :RUN,  Δt, write_restart )
        else
            return_values = ( :END, 0.0, write_restart  )
        end

        return return_values
    end,

    master_after_model_run! = function(OMMODULE, OMDATA)

        global pred_reinit, pred_cnt, data

        fi = OMDATA.mb.fi

        data_next = DataLoader.loadInitDataAndForcing(dli, OMDATA.clock.time + Δt)
        if parsed["SOM"]
            println("SOM would propagate SST and SSS to the whole column")
            data2SOM!(data_next["TEMP"], Nz_bot)
            data2SOM!(data_next["SALT"], Nz_bot)
        end

        #_, m, d = DataLoader.getYMD(OMDATA.clock.time)

        pred_cnt += 1
        #if any(d .== compute_qflx_dates[m])
        if pred_cnt == pred_steps

            writeLog("# [Warning] Compare to OGCM, compute QFLX and reinitialize")

            # compute QFLX
            _Δt = Δt_float * pred_steps
#            @. fi.sv[:QFLXT] = (data_next["TEMP"] - fi.sv[:TEMP]) / _Δt
#            @. fi.sv[:QFLXS] = (data_next["SALT"] - fi.sv[:SALT]) / _Δt
            @. QFLXT = (data_next["TEMP"] - fi.sv[:TEMP]) / _Δt
            @. QFLXS = (data_next["SALT"] - fi.sv[:SALT]) / _Δt

            reinitModel!(OMDATA, data_next; forcing=true, thermal=true, SOM_HMXL=SOM_HMXL)

            pred_cnt = 0
        else
            reinitModel!(OMDATA, data_next; forcing=true, thermal=false, SOM_HMXL=SOM_HMXL)
            #fi._QFLXX_ .= 0.0
        end

        fi.sv[:QFLXT] .= QFLXT
        fi.sv[:QFLXS] .= QFLXS

    end,

    master_finalize! = function(OMMODULE, OMDATA)
        writeLog("[Coupler] Finalize")
    end,
)

runModel(
    ENGINE_EMOM,
    coupler_funcs,
    config,
)

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
module BlockSparseMatrices

using LinearAlgebra
using SparseArrays
# import Base: Order.Forward

# include Block type for indexing
import BlockArrays.Block
export Block

export SparseMatrixBSC, SparseMatrixCSC!, getCSCordering
# export blocksize, getCSCindex
# export nnzblocks, nblocks, nzblockrange

include("SparseMatrixBSC.jl")

end # module

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
isgpu(x) = false
ifgpufree(x) = nothing
function __init__()
    @require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
        gpu_or_cpu(x::CuArrays.CuArray) = CuArrays.CuArray
        gpu_or_cpu(x::Transpose{<:Any,<:CuArrays.CuArray}) = CuArrays.CuArray
        gpu_or_cpu(x::Adjoint{<:Any,<:CuArrays.CuArray}) = CuArrays.CuArray
        isgpu(::CuArrays.CuArray) = true
        isgpu(::Transpose{<:Any,<:CuArrays.CuArray}) = true
        isgpu(::Adjoint{<:Any,<:CuArrays.CuArray}) = true
        ifgpufree(x::CuArrays.CuArray) = CuArrays.unsafe_free!(x)
        ifgpufree(x::Transpose{<:Any,<:CuArrays.CuArray}) = CuArrays.unsafe_free!(x.parent)
        ifgpufree(x::Adjoint{<:Any,<:CuArrays.CuArray}) = CuArrays.unsafe_free!(x.parent)

        @require Tracker="9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" begin
            TrackedArray = Tracker.TrackedArray
            gpu_or_cpu(x::TrackedArray{<:Any,<:Any,<:CuArrays.CuArray}) = CuArrays.CuArray
            gpu_or_cpu(x::Adjoint{<:Any,TrackedArray{<:Any,<:Any,<:CuArrays.CuArray}}) = CuArrays.CuArray
            gpu_or_cpu(x::Transpose{<:Any,TrackedArray{<:Any,<:Any,<:CuArrays.CuArray}}) = CuArrays.CuArray
            isgpu(::TrackedArray{<:Any,<:Any,<:CuArrays.CuArray}) = true
            ifgpufree(x::TrackedArray{<:Any,<:Any,<:CuArrays.CuArray}) = CuArrays.unsafe_free!(x.data)
            isgpu(::Adjoint{<:Any,TrackedArray{<:Any,<:Any,<:CuArrays.CuArray}}) = true
            isgpu(::Transpose{<:Any,TrackedArray{<:Any,<:Any,<:CuArrays.CuArray}}) = true
            ifgpufree(x::Adjoint{<:Any,TrackedArray{<:Any,<:Any,<:CuArrays.CuArray}}) = CuArrays.unsafe_free!((x.data).parent)
            ifgpufree(x::Transpose{<:Any,TrackedArray{<:Any,<:Any,<:CuArrays.CuArray}}) = CuArrays.unsafe_free!((x.data).parent)
        end
    end

    @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" begin
        gpu_or_cpu(x::CUDA.CuArray) = CUDA.CuArray
        gpu_or_cpu(x::Transpose{<:Any,<:CUDA.CuArray}) = CUDA.CuArray
        gpu_or_cpu(x::Adjoint{<:Any,<:CUDA.CuArray}) = CUDA.CuArray
        isgpu(::CUDA.CuArray) = true
        isgpu(::Transpose{<:Any,<:CUDA.CuArray}) = true
        isgpu(::Adjoint{<:Any,<:CUDA.CuArray}) = true
        ifgpufree(x::CUDA.CuArray) = CUDA.unsafe_free!(x)
        ifgpufree(x::Transpose{<:Any,<:CUDA.CuArray}) = CUDA.unsafe_free!(x.parent)
        ifgpufree(x::Adjoint{<:Any,<:CUDA.CuArray}) = CUDA.unsafe_free!(x.parent)

        @require Tracker="9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" begin
            TrackedArray = Tracker.TrackedArray
            gpu_or_cpu(x::TrackedArray{<:Any,<:Any,<:CUDA.CuArray}) = CUDA.CuArray
            gpu_or_cpu(x::Adjoint{<:Any,TrackedArray{<:Any,<:Any,<:CUDA.CuArray}}) = CUDA.CuArray
            gpu_or_cpu(x::Transpose{<:Any,TrackedArray{<:Any,<:Any,<:CUDA.CuArray}}) = CUDA.CuArray
            isgpu(::Adjoint{<:Any,TrackedArray{<:Any,<:Any,<:CUDA.CuArray}}) = true
            isgpu(::TrackedArray{<:Any,<:Any,<:CUDA.CuArray}) = true
            isgpu(::Transpose{<:Any,TrackedArray{<:Any,<:Any,<:CUDA.CuArray}}) = true
            ifgpufree(x::TrackedArray{<:Any,<:Any,<:CUDA.CuArray}) = CUDA.unsafe_free!(x.data)
            ifgpufree(x::Adjoint{<:Any,TrackedArray{<:Any,<:Any,<:CUDA.CuArray}}) = CUDA.unsafe_free!((x.data).parent)
            ifgpufree(x::Transpose{<:Any,TrackedArray{<:Any,<:Any,<:CUDA.CuArray}}) = CUDA.unsafe_free!((x.data).parent)
        end
    end
end


text:
function Calctranpiration(VPD_z::T,Sₑ::T;
    E_cm_ref::T=1.812,s_VPD_z::T=3.121,s_Sₑ::T=18.342) where {T<:Float64}
    #Tor-Ngern et al. 2017 model for estimating canopy transpiration E_c (mm/d)
    #Default parameter values taken are form the original paper
    #for estimating E_c for Rosinedal Scots pine stand
    #---Input---
    #VPD_z day-length normalized vapor pressure deficit (VPD*Day length/24, Pa)
    #Sₑ effective saturation (-)
    E_cm = E_cm_ref*(1-exp(-s_VPD_z*VPD_z*10^-3))
    E_c = E_cm*(1-exp(-s_Sₑ*Sₑ))

    return E_c
end

function Est_E_C(daylight::T;
    env::CCPH.EnvironmentStruct=CCPH.EnvironmentStruct(),
    E_cm_ref::T=1.812,s_VPD_z::T=3.121,s_Sₑ::T=18.342) where {T<:Float64}
    #Estimating canopy transpiration E_C (mm/d) based on weather data,
    #using the Tor-Ngern et al. 2017 model for estimating canopy transpiration
    #---Input---
    # daylight lenght of day (s)
    VPD = env.VPD
    θₛ = env.θₛ
    VPD_z = VPD.*daylight/(24*3600)
    Sₑ = CCPH.CalcSₑ.(θₛ)
    E_C = Calctranpiration.(VPD_z,Sₑ;E_cm_ref=E_cm_ref,s_VPD_z=s_VPD_z,s_Sₑ=s_Sₑ) #mm/day

    return E_C
end

function Est_gₛ(LAI::T,daylight::T;cons::CCPH.Constants=CCPH.Constants(),
    env::CCPH.EnvironmentStruct=CCPH.EnvironmentStruct(),
    E_cm_ref::T=1.812,s_VPD_z::T=3.121,s_Sₑ::T=18.342) where {T<:Float64}
    #Estimating stomatal conductance gₛ (mol m⁻² leaf area s⁻¹) based on weather data,
    #using the Tor-Ngern et al. 2017 model for estimating canopy transpiration
    #---Input---
    # LAI leaf area index (m² leaf area m⁻² ground area)
    # daylight lenght of day (s)
    VPD = env.VPD
    θₛ = env.θₛ
    P = env.P
    VPD_z = VPD.*daylight/(24*3600)
    Sₑ = CCPH.CalcSₑ.(θₛ)
    E_c = Calctranpiration.(VPD_z,Sₑ;E_cm_ref=E_cm_ref,s_VPD_z=s_VPD_z,s_Sₑ=s_Sₑ) #mm/day
    E_c *= 10^-3*cons.ρ_H2O/cons.M_H2O/daylight #mol m⁻² s⁻¹
    gₛ = E_c*P/(VPD*cons.r*LAI)
    return gₛ
end

function Find_data_ind(i::Integer)
    if 1<=i<=21
        return 1
    elseif 22<=i<=44
        return 2
    elseif 45<=i<=64
        return 3
    else
        return 4
    end
end

mutable struct RoData{T<:Float64}
    H::Array{T,1}
    B::Array{T,1}
    GPP::Array{T,1}
    Wf::Array{T,1}
    N::Array{T,1}
    Hc::T
    gₛ::Array{T,1}
    E_C::Array{T,1}
end

function Create_RoData_C_F(GPP_data_F::Array{T,1},GPP_data_C::Array{T,1},
    weatherts_F::CCPH.WeatherTS,weatherts_C::CCPH.WeatherTS;treepar::TreePar=TreePar()) where {T<:Float64}

    H_data_F = [19.07,19.34,19.64,19.87]
    B_data_F = [0.0452,0.0467,0.0479,0.0487]
    Wf_data_F = [8.33,8.47,8.49,8.48]
    N_data_F = [850.0000,846.6667,846.6667,846.6667]/10000
    LAI_F = Wf_data_F.*N_data_F/treepar.LMA

    Hc_F = 10.89

    H_data_C = [20.86,21.01,21.19,21.36]
    B_data_C = [0.0348,0.0356,0.0362,0.0368]
    Wf_data_C = [5.025,5.110,5.141,5.165]
    N_data_C = [1010.0000,1010.0000,1006.6667,1006.6667]/10000
    Hc_C = 10.74
    LAI_C = Wf_data_C.*N_data_C/treepar.LMA

    gₛ_data_F = Float64[]
    gₛ_data_C = Float64[]
    E_C_data_F = Float64[]
    E_C_data_C = Float64[]

    for i = 1:length(weatherts_F.date)
        env_F = EnvironmentStruct(weatherts_F,i)
        env_C = EnvironmentStruct(weatherts_C,i)

        daylight_F = weatherts_F.daylight[i]/7
        daylight_C = weatherts_C.daylight[i]/7

        push!(gₛ_data_F,Est_gₛ(LAI_F[Find_data_ind(i)],daylight_F;env=env_F))
        push!(gₛ_data_C,Est_gₛ(LAI_C[Find_data_ind(i)],daylight_C;env=env_C))

        push!(E_C_data_F,Est_E_C(daylight_F;env=env_F))
        push!(E_C_data_C,Est_E_C(daylight_C;env=env_C))
    end

    data_F = RoData(H_data_F,B_data_F,GPP_data_F,Wf_data_F,N_data_F,Hc_F,gₛ_data_F,E_C_data_F)
    data_C = RoData(H_data_C,B_data_C,GPP_data_C,Wf_data_C,N_data_C,Hc_C,gₛ_data_C,E_C_data_C)

    return data_F,data_C
end

function calcR²(y_data::Array{Float64,1},f_model::Array{Float64,1})
    ymean = mean(y_data)
    SStot = sum((y_data.-ymean).^2)
    SSres = sum((y_data.-f_model).^2)
    return 1-SSres/SStot
end

function clac_GPP_R²_annual(GPP::Array{Float64,1},GPP_model::Array{Float64,1},dates::Array{Dates.DateTime,1})
    years = [Dates.year(dates[1])]
    growth_end = Int64[]
    for i = 2:length(dates)
        if Dates.year(dates[i]) != Dates.year(dates[i-1])
            push!(growth_end,i-1)
            push!(years,Dates.year(dates[i]))
        end
    end
    push!(growth_end,length(dates))

    @show growth_end
    @show years

    println("Year: ",years[1])
    println("R²: ",calcR²(GPP[1:growth_end[1]],GPP_model[1:growth_end[1]]))
    println("Cor: ",cor(GPP[1:growth_end[1]],GPP_model[1:growth_end[1]]))

    for i in 2:length(years)
        println("Year: ",years[i])
        println("R²: ",calcR²(GPP[growth_end[i-1]+1:growth_end[i]],GPP_model[growth_end[i-1]+1:growth_end[i]]))
        println("Cor: ",cor(GPP[growth_end[i-1]+1:growth_end[i]],GPP_model[growth_end[i-1]+1:growth_end[i]]))
    end
end

function Max_k_cost(model::CCPHStruct)
    ψ₅₀,b,g,ρ_H2O,θₛ = model.hydPar.ψ₅₀,model.hydPar.b,model.cons.g,model.cons.ρ_H2O,model.env.θₛ
    #Calculate tree height
    H = model.treesize.H
    #Caluclate soil water potential
    ψₛ = CCPH.θₛ2ψₛ(θₛ)
    #Soil water potential adjusted for gravitational pressure (MPa)
    ψₛ_g = ψₛ-H*ρ_H2O*g*10^-6

    k_cost_max = CCPH.Pfun(ψₛ_g,ψ₅₀,b)

    return k_cost_max
end

function OptimCCPH(growthlength::Float64,model::CCPH.CCPHStruct;K_cost_crit::Float64 = 0.12)

    gₛ_crit = CCPH.Calc_K_costᵢₙᵥ(K_cost_crit,model)

    gₛ_opt,Nₘ_f_opt = CCPH.CCPHTraitmodel(growthlength,model;gₛ_lim_hi=gₛ_crit)

    modeloutput = CCPH.CCPH_run(gₛ_opt,Nₘ_f_opt,growthlength,model)

    return modeloutput,gₛ_opt,Nₘ_f_opt
end
function OptimCCPH!(i::Integer,model::CCPH.CCPHStruct,
    weatherts::CCPH.WeatherTS,photo_kinetic::CCPH.PhotoKineticRates)
    growthlength,step_length = CCPH.Init_weather_par!(i,model,weatherts,photo_kinetic)
    modeloutput,gₛ_opt,Nₘ_f_opt = OptimCCPH(growthlength,model)

    return modeloutput,gₛ_opt,Nₘ_f_opt
end

function CalcE_C(gₛ::Float64,growthlength::Float64,model::CCPH.CCPHStruct)
    VPD = model.env.VPD
    P = model.env.P
    LAI = model.treesize.N*model.treesize.Wf/model.treepar.LMA
    daylight = growthlength/7

    E_c = model.cons.r*gₛ*VPD*LAI/P #mol m⁻² s⁻¹
    E_c *= model.cons.M_H2O*daylight*10^3/model.cons.ρ_H2O #mm day⁻¹

    return E_c
end

#Calc model GPP
function CalcGPP(ccphts::CCPHTS,weatherts::WeatherTS)
    GPP = ccphts.N.*ccphts.P.*weatherts.daylight./weatherts.tot_annual_daylight*1000/7 #g C day⁻¹ m⁻² ground area

    return GPP
end
function CalcGPP(step_length::Float64,model::CCPH.CCPHStruct,modeloutput::CCPH.CCPHOutput)
    GPP = model.treesize.N*modeloutput.P*step_length*1000/7 #g C day⁻¹ m⁻² ground area

    return GPP
end
function CalcGPP!(i::Integer,model::CCPH.CCPHStruct,
    weatherts::CCPH.WeatherTS,photo_kinetic::CCPH.PhotoKineticRates)
    growthlength,step_length = CCPH.Init_weather_par!(i,model,weatherts,photo_kinetic)
    modeloutput,gₛ_opt,Nₘ_f_opt = OptimCCPH(growthlength,model)
    GPP = CalcGPP(step_length,model,modeloutput)

    return GPP, gₛ_opt
end

function CalcModelOutput!(i::Integer,model::CCPH.CCPHStruct,
    weatherts::CCPH.WeatherTS,photo_kinetic::CCPH.PhotoKineticRates)
    growthlength,step_length = CCPH.Init_weather_par!(i,model,weatherts,photo_kinetic)

    modeloutput,gₛ_opt,Nₘ_f_opt = OptimCCPH(growthlength,model)
    GPP = CalcGPP(step_length,model,modeloutput)

    E_C = CalcE_C(gₛ_opt,growthlength*step_length,model)

    return GPP, E_C, modeloutput, gₛ_opt, Nₘ_f_opt
end

function clac_GPP_R²_annual(GPP::Array{Float64,1},GPP_model::Array{Float64,1},dates::Array{Dates.DateTime,1})
    years = [Dates.year(dates[1])]
    growth_end = Int64[]
    for i = 2:length(dates)
        if Dates.year(dates[i]) != Dates.year(dates[i-1])
            push!(growth_end,i-1)
            push!(years,Dates.year(dates[i]))
        end
    end
    push!(growth_end,length(dates))

    @show growth_end
    @show years

    println("Year: ",years[1])
    println("R²: ",calcR²(GPP[1:growth_end[1]],GPP_model[1:growth_end[1]]))
    println("Cor: ",cor(GPP[1:growth_end[1]],GPP_model[1:growth_end[1]]))

    for i in 2:length(years)
        println("Year: ",years[i])
        println("R²: ",calcR²(GPP[growth_end[i-1]+1:growth_end[i]],GPP_model[growth_end[i-1]+1:growth_end[i]]))
        println("Cor: ",cor(GPP[growth_end[i-1]+1:growth_end[i]],GPP_model[growth_end[i-1]+1:growth_end[i]]))
    end
end

function Est_C_assimilation(GPP_Data::T,LAI::T,N::T,weatherts::CCPH.WeatherTS,
    ind::Integer;
    treepar::CCPH.TreePar=CCPH.TreePar(),
    cons::CCPH.Constants=CCPH.Constants()) where {T<:Float64}
    #Estimate per leaf area carbon assimilation rate A (mol C m⁻² leaf area s⁻¹)
    #based on GPP data (g C m⁻² ground area day⁻¹)
    Xₜ = weatherts.acclimation_fac[ind]
    daylight = weatherts.daylight[ind]

    GPP = GPP_Data*10^-3/(daylight/7) #Kg C m⁻² ground area day⁻¹
    P = GPP/N #Kg C tree⁻¹ s⁻¹
    A = P*(N*treepar.k)/(1-exp(-treepar.k*LAI))/(Xₜ*cons.M_C) #mol C m⁻² leaf area s⁻¹

    return A
end