text:
struct JSString
    s::String
end

JSString(s::JSString) = s
jsstring(xs::JSString...) = JSString(string([x.s for x in xs]...))
jsstring(xs...) = jsstring(JSString.(xs)...)

# Required to allow JSStrings to interpolate into normal strings correctly.
Base.print(io::IO, x::JSString) = print(io, x.s)

Base.string(s::JSString) = s.s
Base.:(==)(x::JSString, y::JSString) = x.s==y.s

JSON.lower(x::JSString) = JSON.lower(x.s)

"""
    _Interpolator(s)

An iterator that yields _interpolation chunks_ (_i.e.,_ a string if the chunk
was a string literal and an `Expr` if the chunk was an interpolation).
"""
struct _Interpolator{S <: AbstractString}
    s::S
end

function Base.iterate(interp::_Interpolator, args...)
    return _iterate_interpolations(interp.s, args...)
end
Base.IteratorSize(::_Interpolator) = Base.SizeUnknown()

"""
    _iterate_interpolations(s[, i])

Advance the iterator to obtain the next string or expression to be interpolated. If no parts
remain (where a part is either a string literal or an interpolation), return nothing; otherwise,
return a tuple with either a string or an `Expr` and the new iteration state.

# Examples
See `test/syntax.jl`.
"""
function _iterate_interpolations(s, i=firstindex(s))
    # Notation:
    # i0 - start index
    # l - last index (length)
    # c - the most recently iterated character
    # i - current index (while iterating character-by-character)
    #     note: this is actually the index of the character **after** c
    # j - the end of the current "chunk" of string literal (when set)
    # escape_next - true if we should escape the next character
    # escape_current - true if c was preceded by an escape

    i0 = i
    l = lastindex(s)

    # Return nothing if we're past the end.
    if i0 > l
        return nothing
    end

    # Get the initial character (we have to handle the first character specially).
    c, i = iterate(s, i0)

    # Special case: when the first character is a $, it's an interpolation, so we return the
    # `Expr` of the interpolated _thing_.
    # For example, if `s` is `$foo</script>` (and `i` is 1), we should return `:foo`.
    if c == '$'
        # Note: at this point, `i` is after the `$`, so we parse `foo</script>` and return `foo`
        return Meta.parse(s, i, greedy=false, raise=false)
    end

    escape_next = (c == '\\')
    escape_current = false

    # Iterate over characters, stopping when we hit a `$` or run iterate through end of the string.
    while (c != '$') && i <= l
        c, i = iterate(s, i)
        if escape_next
            escape_current = true
            break
        elseif c == '\\'
            escape_next = true
        end
    end

    # If we have `\$`, we **don't** want to interpolate.
    if c == '$' && escape_current
        # Suppose the string is `this.\$refs = $myvar`. We don't want to interpolate a variable
        # named refs, but we do want to interpolate `myvar`. We set up indices j and k as follows.
        # t h i s . \ $ r e f s ␣ = ␣ \ $ m y v a r
        # i0      j     i           k             l
        # We form the string as s[i0:j] * "$" * s[i:k].
        # We get k by continuing to iterate, starting at position i (which is past the $ and thus
        # treated like a literal string).
        j = prevind(s, i, 3)
        # If we're at the end of the string, `iterate_interpolations` returns `nothing`.
        next_s, k = i < l ? _iterate_interpolations(s, i) : ("", i)
        return (s[i0:j] * "\$" * next_s), k
    end

    if escape_current
        # Julia's string literal escaping logic gets wonky when quotes are involved.
        # For example, `console.log("\\")` gets passed to us as `console.log("\")`
        # whereas `foo = '\\\\'` gets passed to us exactly as written. The issue seems to be that
        # whenever quotes are involved, it eagerly consumes the backslashes in a weird way; but the
        # upshot of this is that if a user ever writes `\"`, it gets passed to as as just `"`, so
        # if we ever see `\"`, we know that the user actually wrote `\\"`.
        if c == '"'
            # j is the index of the quote (_i.e._ s[j] = '"')
            j = prevind(s, i)
            # We get the next chunk of string...
            next_s, k = i < l ? _iterate_interpolations(s, i) : ("", i)
            # And then concatenate it.
            return s[i0:j] * next_s, k
        end
        j = prevind(s, i, 3)
        return s[i0:j] * "$c", i
    end

    # We don't have a preceeding backslash, so we **do** want to interpolate the **next** piece.
    # We do this by returning the current chunk of string literal and setting the iteration pointer
    # to the `$` so that the next call to `iterate_interpolations` knows to return an `Expr`.
    if c == '$'
        # Suppose s is `alert($thing)`.
        # Indices:           j i
        # We get j by shifting back twice (once so that s[i] == '$', and again so that s[i] == '(').
        # We also shift the position of the iterator back one so that the next call will have the
        # string beginning with the $ and will thus start interpolating.
        j = prevind(s, i, 2)
        return s[i0:j], prevind(s, i)
    end

    # If we get here, we're at the end of the string and we didn't find any interpolations or
    # escapes; so just return the rest of the string.
    return s[i0:prevind(s, i)], i
end

"""
    tojs(x)

Returns a JSString object that constructs the same object as `x`
"""
tojs(x) = x

"""
    showjs(io, x)

Print Javascript code to `io` that constructs the equivalent of `x`.
"""
showjs(io, x::Any) = JSON.show_json(io, JSEvalSerialization(), x)
showjs(io, x::AbstractString) = write(io, JSON.json(x))

"""
    @js_str(s)

Create a `JSString` using a string literal and perform interpolations from Julia.

# Examples
```julia-repl
julia> mystr = "foo"; mydict = Dict("foo" => "bar", "spam" => "eggs");
julia> println(js"const myStr = \$mystr; const myObject = \$mydict;")
const myStr = "foo"; const myObject = {"spam":"eggs","foo":"bar"};
```
"""
macro js_str(s)
    writes = map(_Interpolator(s)) do x
        if isa(x, AbstractString)
            # If x is a string, it was specified in the js"..." literal so let it
            # through as-is.
            :(write(io, $(x)))
        else
            # Otherwise, it's some kind of interpolation so we need to generate a
            # JavaScript representation of whatever it is/whatever it evaluates to.
            :(showjs(io, tojs($(esc(x)))))
       end
   end

   :(JSString(sprint(io->(begin; $(writes...) end))))
end

const JSONContext = JSON.Writer.StructuralContext
const JSONSerialization = JSON.Serializations.CommonSerialization

struct JSEvalSerialization <: JSONSerialization end

# adapted (very slightly) from JSON.jl test/serializer.jl
function JSON.show_json(io::JSONContext, ::JSEvalSerialization, x::JSString)
    Base.print(io, x.s)
end


text:
# --------------------------------------------------------------------
# ProcessObservations: all functions to read and post-process
# observations from altimetry and tide gauges
#   Run a virtual-station average for US tide-gauge observations
#   Read altimetry observations
#   Return annual-mean reconstructions
# These routines are called from LocalProjections and
# RegionalProjections
# --------------------------------------------------------------------

module ProcessObservations
using NCDatasets
using NetCDF
using Dates
using XLSX
using DelimitedFiles
using Statistics
using MAT
using Interpolations
using CSV
dir_code = homedir()*"/Projects/2021_NOAA_TR/Code/"
include(dir_code*"Masks.jl")
include(dir_code*"Hector.jl")

function ReadLocalObs(settings)
    # -------------------------------------------------------------------------------------------
    # Read tide-gauge data and return list of local annual-mean tide-gauge data,
    # the station names and coordinates
    #
    # To compute annual-mean sea level from monthly-mean sea level:
    # - First remove mean seasonal cycle, which would give trouble for years with missing months
    # - For years with at least 10 months with data, compute annual mean. Otherwise, NaN32
    # -------------------------------------------------------------------------------------------
    println("  Reading tide-gauge data...")
    # Read data file
    tg_data_raw = XLSX.readxlsx(settings["fn_tg_data"])["monthly MSL"][:]

    # Station name, lon,lat coords. Longitude in [0 360] interval
    station_names = Vector{String}(tg_data_raw[1,2:122])
    station_coords = reverse!(Array{Float32}(tg_data_raw[2:3,2:122]'),dims=2)
    @. station_coords[:,1] = mod(station_coords[:,1],360)

    # Time steps
    tsteps_raw = tg_data_raw[4:end,1]
    tvec = zeros(Int32,length(tsteps_raw),2)
    for t in 1:length(tsteps_raw)
        tvec[t,1] = Year(tsteps_raw[t]).value
        tvec[t,2] = Month(tsteps_raw[t]).value
    end
    tval = @. tvec[:,1] + (tvec[:,2]/12 - 1/24)

    # Monthly-mean RSL
    rsl_raw = @views tg_data_raw[4:end,2:122]
     @. rsl_raw[rsl_raw=="NaN"] = NaN
    rsl_monthly = Array{Float32}(rsl_raw')
    rsl_monthly .*= 1000.0f0 # From meters to millimeters

    # Remove station Seward because of large vertical datum discontinuity
    stat_acc = ones(Bool,size(rsl_monthly,1))
    stat_acc[107] = false
    station_names = station_names[stat_acc]
    station_coords = station_coords[stat_acc,:]
    rsl_monthly = rsl_monthly[stat_acc,:]

    # From monthly to annual-mean sea level
    local_obs = Dict()
    local_obs["name"]  = station_names
    local_obs["coords"] = station_coords
    local_obs["rsl"]            = zeros(Float32,size(rsl_monthly,1),length(settings["years_tg"])) .* NaN32
    for stat in 1:length(station_names)
        # 1. Remove mean seasonal cycle before computing annual means
        acc_mnth = @. isfinite(rsl_monthly[stat,:])
        amat = ones(Float32,sum(acc_mnth),6)
        amat[:,2] = tval[acc_mnth] .- mean(tval[acc_mnth])
        amat[:,3] = @. sin(2*π*tval[acc_mnth])
        amat[:,4] = @. cos(2*π*tval[acc_mnth])
        amat[:,5] = @. sin(4*π*tval[acc_mnth])
        amat[:,6] = @. cos(4*π*tval[acc_mnth])
        sol = amat\rsl_monthly[stat,acc_mnth]
        sol[2] = 0
        rsl_monthly[stat,acc_mnth]  -= amat*sol
        rsl_rfmt = reshape(rsl_monthly[stat,:],(12,length(settings["years_tg"])))
        for yr in 1:length(settings["years_tg"])
            sum(isnan.(rsl_rfmt[:,yr])) < 3 ? local_obs["rsl"][stat,yr] = mean(filter(isfinite,rsl_rfmt[:,yr])) : nothing
        end
    end
    local_obs["region"] = MapToRegions(local_obs,settings)

    # Find nearest PSMSL ID (For Nga)
    psmsl_data = readdlm(settings["dir_project"] * "Data/filelist_psmsl.txt" ,';')
    lon_psmsl = mod.(convert.(Float32,psmsl_data[:,3]),360)
    lat_psmsl = convert.(Float32,psmsl_data[:,2])
    # name_psmsl = psmsl_data[:,4]
    id_psmsl = convert.(Int32,psmsl_data[:,1])
    local_obs["psmsl_id"] = zeros(Int,length(local_obs["name"]))
    for tg in 1:length(local_obs["name"])
        dst_arr = (local_obs["coords"][tg,1].-lon_psmsl).^2 + (local_obs["coords"][tg,2].-lat_psmsl).^2
        if minimum(dst_arr) .< 0.001
            local_obs["psmsl_id"][tg] = id_psmsl[argmin(dst_arr)]
            # println(rpad(local_obs["name"][tg],40," ")*" "*name_psmsl[argmin(dst_arr)])
        else
            local_obs["psmsl_id"][tg] = -1
            # println("NOMATCH: "*rpad(local_obs["name"][tg],40," ")*" "*name_psmsl[argmin(dst_arr)])
        end
    end

    # Flag Chesapeake bridge and La Jolla because of possible data problems
    local_obs["qc_flag"] = zeros(Bool,length(local_obs["name"]))
    local_obs["qc_flag"][findfirst(==(256),local_obs["psmsl_id"])] = true
    local_obs["qc_flag"][findfirst(==(1635),local_obs["psmsl_id"])] = true

    return local_obs
end

function MapToRegions(local_obs, settings)
    # --------------------------------------------------------------
    # Map the individual tide-gauge records to one of the 11 regions
    # Return a list with region numbers
    # --------------------------------------------------------------
    println("  Mapping tide-gauge stations onto regions...")
    mask = Masks.ReadMask(settings)
    region_num = zeros(Bool,length(settings["regions"]),size(local_obs["coords"],1))

    # Find Alaska stations
    in_alaska = zeros(Bool,size(local_obs["coords"],1))
    for stat ∈ 1:size(local_obs["coords"],1)
        strip(local_obs["name"][stat][end-2:end]) == "AK" ? in_alaska[stat] = true : nothing
    end
    NAL_list = ["Unalaska, AK","Prudhoe Bay, AK","Nome, AK","Adak Island, AK ","Unalaska, AK"]
    # Determine region
    for (region_idx,region) ∈ enumerate(settings["regions"])
        if (region!="NAL") & (region!="SAL")
            for stat ∈ 1:size(local_obs["coords"],1)
                ϕ_idx = argmin(@. abs(local_obs["coords"][stat,1] - mask["ϕ"]))
                θ_idx = argmin(@. abs(local_obs["coords"][stat,2] - mask["θ"]))
                if mask[region*"_unc"][ϕ_idx,θ_idx]
                    region_num[region_idx,stat] = true
                end
            end
        elseif region=="NAL"
            for stat ∈ 1:size(local_obs["coords"],1)
                local_obs["name"][stat] in NAL_list ? region_num[region_idx,stat] = true : nothing
            end
        elseif region=="SAL"
            for stat ∈ 1:size(local_obs["coords"],1)
                in_alaska[stat] & !region_num[region_idx-1,stat] ? region_num[region_idx,stat] = true : nothing
            end
        end
    end
    return region_num
end

function ComputeRegionObs(local_obs,settings)
    # -------------------------------------------------------------
    # Read tide-gauge observations and merge using virtual stations
    # Read altimetry observations using region mask
    # Return both curves per region
    # -------------------------------------------------------------
    println("  Computing regional observations...")
    region_obs = Dict()
    [region_obs[region] = Dict() for region in settings["regions"]]

    ComputeVirstat!(local_obs,region_obs,settings)
    ReadAltimetry!(region_obs,settings) # Altimetry

    # Contiguous USA curve is average of the six regions that circumference the Contiguous US
    region_obs["USA"]["rsl"] = (region_obs["NE"]["rsl"] + region_obs["SE"]["rsl"] + region_obs["WGOM"]["rsl"] + region_obs["EGOM"]["rsl"] + region_obs["SW"]["rsl"] + region_obs["NW"]["rsl"]) ./6

    # Remove natural variability
    clim_indices = ReadClimIndices(settings)
    for region in settings["regions"]
        region_obs[region]["rsl"] = RemoveVariability(settings["years_tg"],region_obs[region]["rsl"],clim_indices)
        region_obs[region]["alt"] = RemoveVariability(settings["years_tg"],region_obs[region]["alt"],clim_indices)
    end
    return region_obs
end

function ComputeVirstat!(local_obs,region_obs,settings)
    # ------------------------------------------------------------------------------------------------
    # Merge the tide gauges in each region to an average curve using the
    # Dangendorf et al. (2017) virtual station method:
    # Iteratively find the two nearest stations, subtract local mean sea level
    # over the years where both stations have data from both stations. Then take
    # the average RSL of both stations. Remove these two stations from the list,
    # but add a new one with the averaged sea-level curve halfway both stations.
    # Repeat this procedure until only one virtual station is left. The resulting
    # curve is the estimate of regional sea level.
    #
    # For more details, see:
    #  Dangendorf, S., Marcos, M., Wöppelmann, G., Conrad, C. P., Frederikse, T., & Riva, R. (2017).
    #  Reassessment of 20th century global mean sea level rise. Proceedings of the National Academy
    #  of Sciences, 114(23), 5946–5951. https://doi.org/10.1073/pnas.1616007114
    # ------------------------------------------------------------------------------------------------
    for (region_idx, region) in enumerate(settings["regions"])
        # List with locations and rsl data
        rsl_reg   = local_obs["rsl"][local_obs["region"][region_idx,:],:]
        coord_reg = local_obs["coords"][local_obs["region"][region_idx,:],:]
        while size(coord_reg,1) > 1
            # Find two nearest stations
            merge_idx, ϕ_mid, θ_mid = FindStatsToMerge(rsl_reg,coord_reg)
            merge_data = rsl_reg[[merge_idx[1],merge_idx[2]],:]
            # Subtract the mean over the overlap period
            ovl = @. isfinite(merge_data[1,:]) & isfinite(merge_data[2,:])
            merge_data[1,:] .-= mean(merge_data[1,ovl])
            merge_data[2,:] .-= mean(merge_data[2,ovl])
            merged_data = MergeTwoStats(merge_data)
            # Replace stations by virtaul station
            keep_idx = ones(Bool,size(rsl_reg,1))
            keep_idx[[merge_idx[1],merge_idx[2]]] .= false
            rsl_reg = rsl_reg[keep_idx,:]
            coord_reg = coord_reg[keep_idx,:]
            rsl_reg = vcat(rsl_reg, merged_data')
            coord_reg = vcat(coord_reg, [ϕ_mid, θ_mid]')
        end
        region_obs[region]["rsl"] = rsl_reg[:]
    end
    return nothing
end

function ReadAltimetry!(region_obs,settings)
    # -----------------------------------------------------------------------------------
    # Read satellite altimetry and add GIA and GRD corrections to obtain
    # relative sea level, which is comparable to tide-gauge observations
    #
    # Altimetry observations ahve been obtained from Copernicus Climate Data Store
    # https://cds.climate.copernicus.eu/cdsapp#!/dataset/satellite-sea-level-global
    # This dataset contains modified Copernicus Climate Change Service information [2020]
    # -----------------------------------------------------------------------------------
    println("  Reading altimetry...")
    mask = Masks.ReadMask(settings)
    area = ComputeGridArea(mask["ϕ"],mask["θ"])

    # Read altimetry data
    ssh_grid = convert.(Float32,ncread(settings["fn_altimetry"],"ssh"))
    @. ssh_grid[ssh_grid > 31000] = 0
    @. ssh_grid .*= 0.05
    alt_time = convert.(Float32,ncread(settings["fn_altimetry"],"time"))
    alt_year = [1993:2020...]
    alt_acc  = findall(in(alt_year),settings["years_tg"])

    # From high frquency data to annual-means
    ssh_grid_year = zeros(Float32,size(ssh_grid,1),size(ssh_grid,2),length(alt_year));
    for (idx,yr) ∈ enumerate(alt_year)
        acc_idx = floor.(alt_time) .== yr
        ssh_grid_year[:,:,idx] = mean(ssh_grid[:,:,acc_idx],dims=3)
    end

    # Take regional average
    for region ∈ settings["regions"]
        region_obs[region]["alt"] = zeros(Float32,length(settings["years_tg"])) .* NaN32
        region_obs[region]["alt"][alt_acc] = sum( @.(ssh_grid_year * mask[region] * area),dims=(1,2))[:] ./ sum(mask[region] .* area)
    end

    # 3 Add GIA and contemporary GRD
    GIA = ComputeGIA(mask,area,settings) # GIA from Caron et al. (2019)
    GRD = ComputeGRD(mask,area,alt_year,settings) # Contemporary GRD effects from Frederikse e al. (2020)
    for (region_idx,region) ∈ enumerate(settings["regions"])
        region_obs[region]["alt"][alt_acc] = region_obs[region]["alt"][alt_acc] - (GIA[region_idx] .* (alt_year .- mean(alt_year))) - GRD[:,region_idx]
    end
    return nothing
end

function ComputeGIA(mask,area,settings)
    # --------------------------------------------------------------------------------------
    # Read GIA solid-Earth deformation from:
    # Caron, L., Ivins, E. R., Larour, E., Adhikari, S., Nilsson, J., & Blewitt, G. (2018).
    # GIA Model Statistics for GRACE Hydrology, Cryosphere, and Ocean Science.
    # Geophysical Research Letters, 45(5), 2203–2212. https://doi.org/10.1002/2017GL076644
    # --------------------------------------------------------------------------------------
    println("  Reading GIA...")
    GIA = zeros(Float32,length(settings["regions"]))
    GIA_rad = ncread(settings["fn_GIA"],"rad_mean")
    for (region_idx,region) ∈ enumerate(settings["regions"])
        GIA[region_idx] = sum(@.(GIA_rad * mask[region] * area)) ./ sum(mask[region] .* area)
    end
    return GIA
end

function ComputeGRD(mask,area,years,settings)
    # --------------------------------------------------------------------------------------
    # Read GRD solid-Earth deformation from:
    # Frederikse, T., Landerer, F., Caron, L., Adhikari, S., Parkes, D., Humphrey, V. W.,
    # Dangendorf, S., Hogarth, P., Zanna, L., Cheng, L., & Wu, Y.-H. (2020).
    # The causes of sea-level rise since 1900.
    # Nature, 584(7821), 393–397. https://doi.org/10.1038/s41586-020-2591-3
    # --------------------------------------------------------------------------------------
    println("  Reading contemporary GRD...")
    GRD = zeros(Float32,length(years),length(settings["regions"]))
    GRD_rad = ncread(settings["fn_GRD"],"rad") .* 0.05
    GRD_years = ncread(settings["fn_GRD"],"time")
    for (region_idx,region) ∈ enumerate(settings["regions"])
        GRD_lcl = sum( @.(GRD_rad * mask[region] * area),dims=(1,2))[:] ./ sum(mask[region] .* area)
        GRD[:,region_idx] = LinearInterpolation(GRD_years, GRD_lcl,extrapolation_bc=Line())(years)
    end
    return GRD
end

function RemoveVariability(years,tseries,clim_indices)
    # ---------------------------------------------------
    # Remove the effects of climate indices (NAO/PDO/MEI)
    # from the sea level curves using OLS.
    # ---------------------------------------------------
    acc_idx = @. isfinite(tseries)
    amat = ones(Float32,sum(acc_idx),5);
    amat[:,2] = years[acc_idx].- mean(years[acc_idx]);
    amat[:,3] = LinearInterpolation(clim_indices[1]["years"], clim_indices[1]["index"],extrapolation_bc=Line())(years[acc_idx])
    amat[:,4] = LinearInterpolation(clim_indices[2]["years"], clim_indices[2]["index"],extrapolation_bc=Line())(years[acc_idx])
    amat[:,5] = LinearInterpolation(clim_indices[3]["years"], clim_indices[3]["index"],extrapolation_bc=Line())(years[acc_idx])
    sol = amat\tseries[acc_idx]
    sol[2] = 0
    tseries_var = zeros(Float32,length(tseries)) .* NaN32
    tseries_var[acc_idx] = amat*sol
    tseries_novar = tseries - tseries_var
    return tseries_novar
end

function ReadClimIndices(settings)
    # ---------------------------------------------------------------
    # Read climate indices used for removing natural variability
    # All indices come from NOAA Physical Sciences Laboratory (PSL)
    # https://psl.noaa.gov/data/climateindices/
    # and NOAA Climate Prediction Centre (CPC)
    # https://www.cpc.ncep.noaa.gov/data/teledoc/telecontents.shtml
    # ---------------------------------------------------------------
    clim_indices = Array{Dict}(undef,3)
    [clim_indices[i] = Dict() for i ∈ 1:3]

    # NAO
    clim_indices[1]["name"] = ["NAO"]
    clim_indices[1]["years"] = [1899:2020...]
    clim_indices[1]["index"] =  mean(reshape(readdlm(settings["fn_NAO"],skipstart=1)[1:122,2:end]'[:],(12,:)),dims=1)[1,:]

    # PDO
    clim_indices[2]["name"] = ["PDO"]
    clim_indices[2]["years"] = [1854:2020...]
    clim_indices[2]["index"] =  mean(reshape(readdlm(settings["fn_PDO"],skipstart=2)[1:167,2:end]'[:],(12,:)),dims=1)[1,:]

    # MEI
    MEI_1 = readdlm(settings["fn_MEI_1"])[1:end,2:13]'[:]
    MEI_1_time = [1871+1/24:1/12:2005+23/24...]
    MEI_2 = readdlm(settings["fn_MEI_2"],skipstart=1)[1:42,2:13]'[:]
    MEI_2_time = [1979+1/24:1/12:2020+23/24...]
    MEI_time = union(MEI_1_time,MEI_2_time)
    MEI = zeros(length(MEI_time))
    MEI[end-503:end] = MEI_2
    idx_1 = findall(in(MEI_2_time),MEI_1_time)
    idx_2 = findall(in(MEI_1_time),MEI_2_time)
    idx_3 = findall(!in(MEI_2_time),MEI_1_time)
    MEI_1 = MEI_1 .- mean(MEI_1[idx_1]) .+ mean(MEI_2[idx_2])
    MEI[1:end-504] = MEI_1[idx_3]
    clim_indices[3]["name"] = ["MEI"]
    clim_indices[3]["years"] = [1871:2020...]
    clim_indices[3]["index"] =  mean(reshape(MEI,(12,:)),dims=1)[1,:]
    return clim_indices
end

### General functions
function MergeTwoStats(merge_data)
    merged_data = zeros(Float32,size(merge_data,2))
    for t ∈ 1:length(merged_data)
        acc = isfinite.(merge_data[:,t])
        merged_data[t] = ifelse(sum(acc)>0,mean(merge_data[acc,t]),NaN)
    end
    return merged_data
end

function FindStatsToMerge(rsl_reg,loc_reg)
    dst_arr = ones(Float32,size(loc_reg,1),size(loc_reg,1)) * 100000
    for i ∈ 1:size(loc_reg,1), j ∈ i+1:size(loc_reg,1)
        if sum(@. isfinite(rsl_reg[i,:] * rsl_reg[j,:])) > 10
            dst_arr[i,j] = ComputeDistance(loc_reg[i,1],loc_reg[i,2],loc_reg[j,1],loc_reg[j,2])
        end
    end
    merge_idx = argmin(dst_arr)
    ϕ_mid, θ_mid = ComputeMidpoint(loc_reg[merge_idx[1],1],loc_reg[merge_idx[1],2],loc_reg[merge_idx[2],1],loc_reg[merge_idx[2],2])
    return merge_idx, ϕ_mid, θ_mid
end

function ComputeDistance(ϕ1,θ1,ϕ2,θ2)
    dist = 2*6371*asin(sqrt(sind(0.5*(θ2-θ1))^2+cosd(θ1)*cosd(θ2)*sind(0.5*(ϕ1-ϕ2))^2))
    return dist
end

function ComputeMidpoint(ϕ1,θ1,ϕ2,θ2)
    Bx = cos(deg2rad(θ2)) * cos(deg2rad(ϕ2 - ϕ1))
    By = cos(deg2rad(θ2)) * sin(deg2rad(ϕ2 - ϕ1))
    ϕ_mid = rad2deg(deg2rad(ϕ1)+atan(By, cos(deg2rad(θ1)) + Bx))
    θ_mid = rad2deg(atan(sin(deg2rad(θ1))+sin(deg2rad(θ2)), sqrt((cos(deg2rad(θ1))+Bx)^2 + By^2)))
    return ϕ_mid, θ_mid
end

function ComputeGridArea(lon,lat)
    gridsize = abs(lat[2]-lat[1])
    area = @. deg2rad(gridsize) * (sind(lat+gridsize/2)-sind(lat-gridsize/2)) * 6371000^2
    return repeat(area',size(lon,1))
end
end

ext: jl


text:
function tag_package(api::GitHub.GitHubAPI, rname, ver::VersionNumber, mcs, auth; tag_name = "v$ver")
    tagger = Dict("name" => CONFIG["github"]["user"],
                  "email" => CONFIG["github"]["email"],
                  "date" => Dates.format(now(), dateformat"YYYY-mm-ddTHH:MM:SSZ"))
    create_tag(api, rname; auth=auth,
               params=Dict("tag" => tag_name,
                           "message" => "Release: $tag_name",
                           "object" => mcs,
                           "type" => "commit",
                           "tagger" => tagger))
end

function get_metadata_from_pr_body(api::GitHub.GitHubAPI, rp::RequestParams, auth)
    reg_name = rp.reponame
    reg_prid = rp.trigger_src.prid

    pr = pull_request(api, reg_name, reg_prid; auth=auth)

    mstart = match(r"<!--", pr.body)
    mend = match(r"-->", pr.body)

    key = CONFIG["enc_key"]
    try
        enc_meta = strip(pr.body[mstart.offset+4:mend.offset-1])
        meta = String(decrypt(MbedTLS.CIPHER_AES_128_CBC, key, hex2bytes(enc_meta), key))
        return JSON.parse(meta)
    catch ex
        @debug("Exception occured while parsing PR body", get_backtrace(ex))
    end

    nothing
end

function handle_approval(api::GitHub.GitHubAPI, rp::RequestParams{ApprovalTrigger})
    auth = get_access_token(api, rp.evt)
    d = get_metadata_from_pr_body(api, rp, auth)

    if d === nothing
        return "Unable to get registration metdata for this PR"
    end

    reg_name = rp.reponame
    reg_prid = rp.trigger_src.prid
    reponame = d["pkg_repo_name"]
    ver = VersionNumber(d["version"])
    tree_sha = d["tree_sha"]
    trigger_id = d["trigger_id"]
    request_type = d["request_type"]
    subdir = d["subdir"]

    if request_type == "pull_request"
        pr = pull_request(api, reponame, trigger_id; auth=auth)
        tree_sha = pr.merge_commit_sha
        if pr.state == "open"
            @debug("Merging pull request on package repo", reponame, trigger_id)
            merge_pull_request(api, reponame, trigger_id; auth=auth,
                               params=Dict("merge_method" => "squash"))
        else
            @debug("Pull request already merged", reponame, trigger_id)
        end
    end

    tag_exists = false
    tag = tag_name(ver, subdir)
    # Get tags in a try-catch block as GitHub.jl error if no tag exists
    try
        ts = tags(reponame; auth=auth, page_limit=1, params=Dict("per_page" => 15))[1]
        for t in ts
            if split(t.url.path, "/")[end] == tag
                if t.object["sha"] != tree_sha
                    return "Tag with name `$tag` already exists and points to a different commit"
                end
                tag_exists = true
                @debug("Tag already exists", reponame, ver, tree_sha)
                break
            end
        end
    catch e
        if occursin("Status Code: 404", e.msg)
            @debug("No tag exists", reponame)
        else
            rethrow(e)
        end
    end

    if !tag_exists
        @debug("Creating new tag", reponame, ver, tree_sha)
        tag_package(api, reponame, ver, tree_sha, auth; tag_name = tag)
    end

    release_exists = false
    if tag_exists
        # Look for release in last 15 releases
        rs = releases(reponame; auth=auth, page_limit=1, params=Dict("per_page"=>15))[1]
        for r in rs
            if r.name == tag
                release_exists = true
                @debug("Release already exists", r.name)
                break
            end
        end
    end

    if !release_exists
        @debug("Creating new release", ver)
        create_release(reponame; auth=auth,
                       params=Dict("tag_name" => tag, "name" => tag))
    end

    if request_type == "issue"
        iss = issue(reponame, Issue(trigger_id); auth=auth)
        if iss.state == "open"
            @debug("Closing issue", reponame, trigger_id)
            edit_issue(reponame, trigger_id; auth=auth, params=Dict("state"=>"closed"))
        else
            @debug("Issue already closed", reponame, trigger_id)
        end
    end

    reg_pr = pull_request(api, reg_name, reg_prid; auth=auth)
    if reg_pr.state == "open"
        @debug("Merging pull request on registry", reg_name, reg_prid)
        merge_pull_request(api, reg_name, reg_prid; auth=auth)
    else
        @debug("Pull request on registry already merged", reg_name, reg_prid)
    end
    nothing
end

function print_entry_log(rp::RequestParams{ApprovalTrigger})
    @info "Approving Pull request" reponame=rp.reponame prid=rp.trigger_src.prid
end

function action(api::GitHub.GitHubAPI, rp::RequestParams{ApprovalTrigger}, zsock)
    @info("Processing approval event", reponame=rp.reponame, rp.trigger_src.prid)
    try
        err = handle_approval(api, rp)
        if err !== nothing
            @debug(err)
            make_comment(api, rp.evt, "Error in approval process: $err")
        end
    catch ex
        bt = get_backtrace(ex)
        @info("Unexpected error: $bt")
        raise_issue(api, rp.evt, rp.phrase, bt)
    end
    @info("Done processing approval event", reponame=rp.reponame, rp.trigger_src.prid)
end

string(::RequestParams{ApprovalTrigger}) = "approval"


text:
# Testing:
#
# - type hierarchy

using DiffModels
using Test

@test ConstDrift <: AbstractDrift
@test VarDrift <: AbstractDrift

@test ConstBound <: AbstractBound
@test VarBound <: AbstractBound
@test LinearBound <: AbstractBound

@test SymBounds <: AbstractBounds
@test AsymBounds <: AbstractBounds
@test ConstSymBounds <: AbstractBounds
@test ConstSymBounds <: ConstBounds
@test VarSymBounds <: AbstractBounds
@test ConstAsymBounds <: AbstractBounds
@test ConstAsymBounds <: ConstBounds
@test VarAsymBounds <: AbstractBounds


text:
module DbIdNameMapping

function allGeneReferenceProduct(dbidfile)
    genes = []
    header = true
    for line in readlines(dbidfile)
        if header == false
            lineparts = split(chomp(line), "\t")
            node = String(lineparts[2])
            if contains(lineparts[6], "Reference") == true
                push!(genes, node)
            end
        end
        header = false
    end

    return genes
end

function geneToNodes(dbidfile)
    genetonodes = Dict()
    header = true
    for line in readlines(dbidfile)
        if header == false
            lineparts = split(chomp(line), "\t")
            node = String(lineparts[2])
            gene = String(lineparts[3])
            if haskey(genetonodes, gene)
                nodes = genetonodes[gene]
                push!(nodes, node)
            else
                genetonodes[gene] = [node]
            end
        end
        header = false
    end

    return genetonodes
end

function nodeToGene(dbidfile)
   genetonodes = geneToNodes(dbidfile)

   nodetogene = Dict()
   for (gene, nodes) in genetonodes
       for node in nodes
           nodetogene[node] = gene
       end
   end

   return nodetogene
end

function geneToRootNodes(pinodes, dbidfile)
    genetonodes = Dict()
    header = true
    for line in readlines(dbidfile)
        if header == false
            lineparts = split(chomp(line), "\t")
            node = String(lineparts[2])
            gene = String(lineparts[3])
            if haskey(pinodes, node) && pinodes[node].relation == "ROOT"
    #        if match(r"Reference", lineparts[6]) != nothing && haskey(pinodes, node) && pinodes[node].relation == "ROOT"
                if haskey(genetonodes, gene)
                    nodes = genetonodes[gene]
                    push!(nodes, node)
                else
                    genetonodes[gene] = [node]
                end
            end
        end
        header = false
    end

    return genetonodes
end

end

max_stars_repo_name: OICR/mp-biopath

max_issues_repo_name: OICR/mp-biopath

max_forks_repo_name: OICR/mp-biopath


text:
const Hdwf = Cint

@enum(ENUMFILTER,
  ENUMFILTERALL=0,
  ENUMFILTEREEXPLORER=1,
  ENUMFILTERDISCOVERY=2)

@enum(DEVID,
  DEVIDEEXPLORER=0,
  DEVIDDISCOVERY=1,
  DEVIDDISCOVERY2=2)

@enum(DEVVER,
  DEVVEREEXPLORERCDISCOVERYB=2,
  DEVVEREEXPLORERE=4,
  DEVVEREEXPLORERF=5,
  DEVVERDISCOVERYA=1,
  DEVVERDISCOVERYC=3)

@enum(DWFENUMCONFIGINFO,
  DECIANALOGINCHANNELCOUNT=1,
  DECIANALOGOUTCHANNELCOUNT=2,
  DECIANALOGIOCHANNELCOUNT=3,
  DECIDIGITALINCHANNELCOUNT=4,
  DECIDIGITALOUTCHANNELCOUNT=5,
  DECIDIGITALIOCHANNELCOUNT=6,
  DECIANALOGINBUFFERSIZE=7,
  DECIANALOGOUTBUFFERSIZE=8,
  DECIDIGITALINBUFFERSIZE=9,
  DECIDIGITALOUTBUFFERSIZE=0)

@enum(DWFSTATE,
  DWFSTATEREADY=0,
  DWFSTATECONFIG=4,
  DWFSTATEPREFILL=5,
  DWFSTATEARMED=1,
  DWFSTATEWAIT=7,
  DWFSTATETRIGGEREDRUNNING=3,
  DWFSTATEDONE=2)

@enum(ACQMODE,
ACQMODESINGLE=0,
ACQMODESCANSHIFT=1,
ACQMODESCANSCREEN=2,
ACQMODERECORD=3)

@enum(FILTER,
FILTERDECIMATE=0,
FILTERAVERAGE=1,
FILTERMINMAX=2)

@enum(TRIGSRC,
TRIGSRCNONE=0,
TRIGSRCPC=1,
TRIGSRCDETECTORANALOGIN=2,
TRIGSRCDETECTORDIGITALIN=3,
TRIGSRCANALOGIN=4,
TRIGSRCDIGITALIN=5,
TRIGSRCDIGITALOUT=6,
TRIGSRCANALOGOUT1=7,
TRIGSRCANALOGOUT2=8,
TRIGSRCANALOGOUT3=9,
TRIGSRCANALOGOUT4=10,
TRIGSRCEXTERNAL1=11,
TRIGSRCEXTERNAL2=12,
TRIGSRCEXTERNAL3=13,
TRIGSRCEXTERNAL4=14)

@enum(DWFANALOGOUTMODE,
DWFANALOGOUTMODEVOLTAGE=0,
DWFANALOGOUTMODECURRENT=1)

@enum(DWFANALOGOUTIDLE,
DWFANALOGOUTIDLEDISABLE=0,
DWFANALOGOUTIDLEOFFSET=1,
DWFANALOGOUTIDLEINITIAL=2)

@enum(ANALOGOUTNODE,
ANALOGOUTNODECARRIER=0,
ANALOGOUTNODEFM=1,
ANALOGOUTNODEAM=2)

@enum(FUNC,
FUNCDC=0,
FUNCSINE=1,
FUNCSQUARE=2,
FUNCTRIANGLE=3,
FUNCRAMPUP=4,
FUNCRAMPDOWN=5,
FUNCNOISE=6,
FUNCCUSTOM=30,
FUNCPLAY=31)


text:
using Documenter, RecipesBase, Plots

makedocs(
    sitename = "RecipesBase.jl",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    pages = [
        "Home" => "index.md",
        "Recipe Syntax" => "syntax.md",
        "Recipe Types" => "types.md",
        "Internals" => "internals.md",
        "Library" => "api.md"
    ],
    modules = [RecipesBase]
)

deploydocs(
    repo = "github.com/JuliaPlots/RecipesBase.jl.git",
    push_preview = true,
)


text:
"""
    function getS(m::SBML.Model; zeros=spzeros)::Tuple{Vector{String},Vector{String},AbstractMatrix{Float64}}

Extract the vector of species (aka metabolite) identifiers, vector of reaction
identifiers, and the (dense) stoichiometry matrix from an existing `SBML.Model`.
Returns a tuple with these values.

The matrix is sparse by default (initially constructed by
`SparseArrays.spzeros`). You can fill in a custom empty matrix constructed to
argument `zeros`; e.g. running with `zeros=zeros` will produce a dense matrix.
"""
function getS(
    m::SBML.Model;
    zeros = spzeros,
)::Tuple{Vector{String},Vector{String},AbstractMatrix{Float64}}
    rows = [k for k in keys(m.species)]
    cols = [k for k in keys(m.reactions)]
    rowsd = Dict(k => i for (i, k) in enumerate(rows))
    S = zeros(Float64, length(rows), length(cols))
    for col = 1:length(cols)
        stoi = m.reactions[cols[col]].stoichiometry
        S[getindex.(Ref(rowsd), keys(stoi)), col] .= values(stoi)
    end
    return rows, cols, S
end

"""
    getLBs(m::SBML.Model)::Vector{Tuple{Float64,String}}

Extract a vector of lower bounds of reaction rates from the model. All bounds
are accompanied with the unit of the corresponding value (this behavior is
based on SBML specification).
"""
getLBs(m::SBML.Model)::Vector{Tuple{Float64,String}} =
    broadcast(x -> x.lb, values(m.reactions))

"""
    getUBs(m::SBML.Model)::Vector{Tuple{Float64,String}}

Likewise to [`getLBs`](@ref), extract the upper bounds.
"""
getUBs(m::SBML.Model)::Vector{Tuple{Float64,String}} =
    broadcast(x -> x.ub, values(m.reactions))

"""
    getOCs(m::SBML.Model)::Vector{Float64}

Extract the vector of objective coefficients of each reaction.
"""
getOCs(m::SBML.Model)::Vector{Float64} = broadcast(x -> x.oc, values(m.reactions))

"""
    initial_amounts(m::SBML.Model; convert_concentrations = false)

Return initial amounts for each species as a generator of pairs
`species_name => initial_amount`; the amount is set to `nothing` if not
available. If `convert_concentrations` is true and there is information about
initial concentration available together with compartment size, the result is
computed from the species' initial concentration.

In the current version, units of the measurements are completely ignored.

# Example
```
# get the initial amounts as dictionary
Dict(initial_amounts(model, convert_concentrations = true))

# remove the empty entries
Dict(k => v for (k,v) in initial_amounts(model) if !isnothing(v))
```
"""
initial_amounts(m::SBML.Model; convert_concentrations = false) = (
    k => if !isnothing(s.initial_amount)
        s.initial_amount[1]
    elseif convert_concentrations &&
           !isnothing(s.initial_concentration) &&
           haskey(m.compartments, s.compartment) &&
           !isnothing(m.compartments[s.compartment].size)
        s.initial_concentration[1] * m.compartments[s.compartment].size
    else
        nothing
    end for (k, s) in m.species
)

"""
    initial_concentrations(m::SBML.Model; convert_amounts = false)

Return initial concentrations of the species in the model. Refer to work-alike
[`initial_amounts`](@ref) for details.
"""
initial_concentrations(m::SBML.Model; convert_amounts = false) = (
    k => if !isnothing(s.initial_concentration)
        s.initial_concentration[1]
    elseif convert_amounts &&
           !isnothing(s.initial_amount) &&
           haskey(m.compartments, s.compartment) &&
           !isnothing(m.compartments[s.compartment].size)
        s.initial_amount[1] / m.compartments[s.compartment].size
    else
        nothing
    end for (k, s) in m.species
)

"""
    function extensive_kinetic_math(
        m::SBML.Model,
        formula::SBML.Math;
        handle_empty_compartment_size = (id::String) -> throw(
            DomainError(
                "Non-substance-only-unit reference to species `\$id' in an unsized compartment `\$(m.species[id].compartment)'",
            ),
        ),
    )

Convert a SBML math `formula` to "extensive" kinetic laws, where the references
to species that are marked as not having only substance units are converted
from amounts to concentrations.

If the data is missing, you can supply a function that adds them. A common way
to handle errors is to assume that unsized compartments have volume 1.0 (of
whatever units), you can specify that behavior by supplying
`handle_empty_compartment_size = _ -> 1.0`.

Handling of units in the conversion process is ignored in this version.
"""
function extensive_kinetic_math(
    m::SBML.Model,
    formula::SBML.Math;
    handle_empty_compartment_size = (id::String) -> throw(
        DomainError(
            "Non-substance-only-unit reference to species `$id' in an unsized compartment `$(m.species[id].compartment)'",
        ),
    ),
)
    conv(x::SBML.MathIdent) = begin
        haskey(m.species, x.id) || return x
        sp = m.species[x.id]
        sp.only_substance_units && return x
        sz = m.compartments[sp.compartment].size
        isnothing(sz) && (sz = handle_empty_compartment_size(x.id))
        SBML.MathApply("/", [x, SBML.MathVal(sz)])
    end
    conv(x::SBML.MathApply) = SBML.MathApply(x.fn, conv.(x.args))
    conv(x::SBML.Math) = x

    conv(formula)
end

"""
    get_error_messages(doc::VPtr, error::Exception, report_severities)

Show the error messages reported by SBML in the `doc` document and throw the
`error` if they are more than 1.

`report_severities` switches the reporting of certain error types defined by
libsbml; you can choose from `["Fatal", "Error", "Warning", "Informational"]`.
"""
function get_error_messages(doc::VPtr, error::Exception, report_severities)
    n_errs = ccall(sbml(:SBMLDocument_getNumErrors), Cuint, (VPtr,), doc)
    do_throw = false
    for i = 1:n_errs
        err = ccall(sbml(:SBMLDocument_getError), VPtr, (VPtr, Cuint), doc, i - 1)
        msg = string(strip(get_string(err, :XMLError_getMessage)))
        sev = string(strip(get_string(err, :XMLError_getSeverityAsString)))
        # keywords from `libsbml/src/sbml/xml/XMLError.cpp` xmlSeverityStringTable:
        if sev == "Fatal"
            sev in report_severities && @error "SBML reported fatal error: $(msg)"
            do_throw = true
        elseif sev == "Error"
            sev in report_severities && @error "SBML reported error: $(msg)"
            do_throw = true
        elseif sev == "Warning"
            sev in report_severities && @warn "SBML reported warning: $(msg)"
        elseif sev == "Informational"
            sev in report_severities && @info "SBML reported: $(msg)"
        end
    end
    do_throw && throw(error)
    nothing
end

"""
    check_errors(
        success::Integer,
        doc::VPtr,
        error::Exception,
        report_severities = ["Fatal", "Error"],
    )

If success is a 0-valued `Integer` (a logical `false`), then call
[`get_error_messages`](@ref) to show the error messages reported by SBML in the
`doc` document and throw the `error` if they are more than 1.  `success` is
typically the value returned by an SBML C function operating on `doc` which
returns a boolean flag to signal a successful operation.
"""
check_errors(
    success::Integer,
    doc::VPtr,
    error::Exception,
    report_severities = ["Fatal", "Error"],
) = Bool(success) || get_error_messages(doc, error, report_severities)

# NOTE: this mapping is valid for Level 3/Version 2, it *may* not be valid for
# other versions.  See
# https://github.com/sbmlteam/libsbml/blob/d4bc12abc4e72e451a0a0f2be4b0b6101ac94160/src/sbml/UnitKind.c#L46-L85
const UNIT_KIND_STRINGS_TO_UNIT = Dict(
    "ampere" => 1.0 * u"A", # UNIT_KIND_AMPERE
    "avogadro" => ustrip(u"mol^-1", Unitful.Na), # UNIT_KIND_AVOGADRO
    "becquerel" => 1.0 * u"Bq", # UNIT_KIND_BECQUEREL
    "candela" => 1.0 * u"cd", # UNIT_KIND_CANDELA
    "Celsius" => 1.0 * u"°C", # UNIT_KIND_CELSIUS
    "coulomb" => 1.0 * u"C", # UNIT_KIND_COULOMB
    "dimensionless" => 1, # UNIT_KIND_DIMENSIONLESS
    "farad" => 1.0 * u"F", # UNIT_KIND_FARAD
    "gram" => 1.0 * u"g", # UNIT_KIND_GRAM
    "gray" => 1.0 * u"Gy", # UNIT_KIND_GRAY
    "henry" => 1.0 * u"H", # UNIT_KIND_HENRY
    "hertz" => 1.0 * u"Hz", # UNIT_KIND_HERTZ
    "item" => 1, # UNIT_KIND_ITEM
    "joule" => 1.0 * u"J", # UNIT_KIND_JOULE
    "katal" => 1.0 * u"kat", # UNIT_KIND_KATAL
    "kelvin" => 1.0 * u"K", # UNIT_KIND_KELVIN
    "kilogram" => 1.0 * u"kg", # UNIT_KIND_KILOGRAM
    "liter" => 1.0 * u"L", # UNIT_KIND_LITER
    "litre" => 1.0 * u"L", # UNIT_KIND_LITRE
    "lumen" => 1.0 * u"lm", # UNIT_KIND_LUMEN
    "lux" => 1.0 * u"lx", # UNIT_KIND_LUX
    "meter" => 1.0 * u"m", # UNIT_KIND_METER
    "metre" => 1.0 * u"m", # UNIT_KIND_METRE
    "mole" => 1.0 * u"mol", # UNIT_KIND_MOLE
    "newton" => 1.0 * u"N", # UNIT_KIND_NEWTON
    "ohm" => 1.0 * u"Ω", # UNIT_KIND_OHM
    "pascal" => 1.0 * u"Pa", # UNIT_KIND_PASCAL
    "radian" => 1.0 * u"rad", # UNIT_KIND_RADIAN
    "second" => 1.0 * u"s", # UNIT_KIND_SECOND
    "siemens" => 1.0 * u"S", # UNIT_KIND_SIEMENS
    "sievert" => 1.0 * u"Sv", # UNIT_KIND_SIEVERT
    "steradian" => 1.0 * u"sr", # UNIT_KIND_STERADIAN
    "tesla" => 1.0 * u"T", # UNIT_KIND_TESLA
    "volt" => 1.0 * u"V", # UNIT_KIND_VOLT
    "watt" => 1.0 * u"W", # UNIT_KIND_WATT
    "weber" => 1.0 * u"W", # UNIT_KIND_WEBER
    "(Invalid UnitKind)" => 1, # UNIT_KIND_INVALID (let's treat is as a dimensionless quantity)
)

# Get a `Unitful` quantity out of a `Unit_t`.
get_unit(u::VPtr) =
    UNIT_KIND_STRINGS_TO_UNIT[unsafe_string(
        ccall(
            sbml(:UnitKind_toString),
            Cstring,
            (Cint,),
            ccall(sbml(:Unit_getKind), Cint, (VPtr,), u)
        ))] ^
    ccall(sbml(:Unit_getExponent), Cint, (VPtr,), u) *
    exp10(ccall(sbml(:Unit_getScale), Cint, (VPtr,), u)) *
    ccall(sbml(:Unit_getMultiplier), Cdouble, (VPtr,), u)

# Get `Unitful` quantity out of a `UnitDefinition_t`.
get_units(ud::VPtr) =
    prod(get_unit(ccall(sbml(:UnitDefinition_getUnit), VPtr, (VPtr, Cuint), ud, j - 1))
         for j = 1:ccall(sbml(:UnitDefinition_getNumUnits), Cuint, (VPtr,), ud))


text:
module MatDeforLinearElasticModule

__precompile__(true)

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtoolsDeforLinear.DeforModelRedModule: AbstractDeforModelRed
import FinEtoolsDeforLinear.MatDeforModule: AbstractMatDefor

"""
    AbstractMatDeforLinearElastic <: AbstractMatDefor

Abstract Linear Elasticity  material.

"""
abstract type AbstractMatDeforLinearElastic <: AbstractMatDefor; end

"""
    tangentmoduli!(self::AbstractMatDeforLinearElastic,  D::FFltMat,  t::FFlt, dt::FFlt, loc::FFltMat, label::FInt)

Calculate the material stiffness matrix.

- `D` = matrix of tangent moduli, supplied as a buffer and overwritten. Returned
as output.
"""
function tangentmoduli!(self::AbstractMatDeforLinearElastic,  D::FFltMat,  t::FFlt, dt::FFlt, loc::FFltMat, label::FInt)
    return self.tangentmoduli!(self, D, t, dt, loc, label)
end

"""
    update!(self::AbstractMatDeforLinearElastic,  stress::FFltVec, output::FFltVec,  strain::FFltVec, thstrain::FFltVec=zeros(6), t::FFlt= 0.0, dt::FFlt= 0.0,  loc::FFltMat=zeros(3,1), label::FInt=0, quantity=:nothing)

Update material state.

- `strain` = strain vector,
- `thstrain` = thermal strain vector,
- `t` = current time,
- `dt` = current time step,
- `loc` = location of the quadrature point in global Cartesian coordinates,
- `label` = label of the finite element in which the quadrature point is found.

# Output
- `stress` = stress vector, allocated by the caller with a size of the number of stress and
strain components, `nstressstrain`. The components of the stress vector are
calculated and stored in the `stress` vector.
- `output` =  array which is (if necessary) allocated  in an appropriate size, filled
  with the output quantity, and returned.
"""
function update!(self::AbstractMatDeforLinearElastic,  stress::FFltVec, output::FFltVec,  strain::FFltVec, thstrain::FFltVec=zeros(6), t::FFlt= 0.0, dt::FFlt= 0.0,  loc::FFltMat=zeros(3,1), label::FInt=0, quantity=:nothing)
    return self.update!(self, stress, output, strain, thstrain, t, dt, loc, label, quantity)
end

"""
    thermalstrain!(self::AbstractMatDeforLinearElastic, thstrain::FFltVec, dT= 0.0)

Compute thermal strain from the supplied temperature increment.

- `thstrain` = thermal strain vector, supplied as buffer, returned as output.
"""
function thermalstrain!(self::AbstractMatDeforLinearElastic, thstrain::FFltVec, dT= 0.0)
    return self.thermalstrain!(self, thstrain, dT)
end

end


text:
function operatorspawn(command,
    logicset::LogicalCombo;
    returnlogical=false,
    prefix=">>> ",
    verbose=true,
    printall=false)

    logicsetcopy = deepcopy(logicset)

    tempcommand = command
    m = eachmatch(r"(\{\{.*?\}\})", tempcommand)
    matches = [replace(x[1], r"\{|\}"=>"") for x in collect(m)] |> unique

    # Check if a {{j}}. structure or {{J}}. structure exists
    m_dot = eachmatch(r"(\{\{[jJ]\}\})(\.[a-zA-Z0-9_.])", tempcommand)
    matches_dot = string.([x.captures[2] for x in collect(m_dot)]) |> unique

    if occursin(r"^[0-9]+,[0-9]+$", matches[end])
        countrange = (x -> x[1]:x[2])(integer.(split(matches[end], ",")))
        tempcommand = replace(tempcommand, "{{$(matches[end])}}"=>"") |> strip
        constraintprint = "{{" * matches[end]  * "}}"
        matches = matches[1:(end-1)]
    elseif occursin(r"^[0-9]+,$", matches[end])
        countrange = integer(matches[end][1:(end-1)]):size(logicset,2)^2
        tempcommand = replace(tempcommand, "{{$(matches[end])}}"=>"") |> strip
        constraintprint = "{{" * replace(matches[end], r",$"=>",∞") * "}}"
        matches = matches[1:(end-1)]
    elseif occursin(r"^,[0-9]+$", matches[end])
        countrange = 0:integer(matches[end][2:end])
        tempcommand = replace(tempcommand, "{{$(matches[end])}}"=>"") |> strip
        constraintprint = "{{" *  replace(matches[end], r"^,"=>"0,")  * "}}"
        matches = matches[1:(end-1)]
    elseif occursin(r"^[0-9]+$", matches[end])
        countrange = (x -> x[1]:x[1])(integer(matches[end]))
        tempcommand = replace(tempcommand, "{{$(matches[end])}}"=>"") |> strip
        constraintprint = "{{" *  matches[end]  * "}}"
        matches = matches[1:(end-1)]
    else
        countrange = missing
        constraintprint = ""
    end

    ijJ = map(x -> x[1], eachmatch(r"(?:\b)(j|J|i)(?:\b)", join(matches, " "))) |>
              unique .|> string

    mykeys   = keys(logicset)
    # Reduce the keyset if `j` or `J` are used
    # Remove variables that have attributes
    (length(ijJ)>0) && (ijJ[1] == "J") && (mykeys = [v for v in mykeys if !occursin(".", string(v))])
    # Remove attributes but keep variable names
    (length(ijJ)>0) && (ijJ[1] == "j") && (mykeys = unique([Symbol(replace(v, r"\..*$"=>"")) for v in string.(mykeys)]))

    rangematches = matches[occursin.(r">|<|!.+$" , matches)]
    imatches = matches[.!occursin.(r">|<|!.+$" , matches)]
    imatches = [replace(v , "!"=>"") for v in imatches]

   function firstmatch(x::String, y::Array{String})
     for i in 1:length(y)
         (x == y[i]) && return i
     end
   end

   rangematches0 = copy(rangematches)

    for v in 1:length(rangematches)
      m = eachmatch(r"(:[a-zA-Z][a-zA-Z0-9._]*)", rangematches[v])
      for vv in collect(m)
        rangematches[v] =
        replace(rangematches[v], string(vv[1]) => firstmatch(string(vv[1])[2:end], string.(mykeys)))
      end
    end

    (length(rangematches)>3) && throw("Only three !i, >i, >=i, <=i, <i wildcard allowed with i (or j/J)")
    (length(ijJ)>1)     && throw("Only one type i or j or J wildcard allowed")
    (length(matches)>6) && throw("No more than 6 wildcard values allowed")

    mydomain = 1:length(mykeys)

    collection = []

    txtoutarray = String[]

    iset = Symbol[]

    for i in mydomain, j in mydomain, k in mydomain, l in mydomain
      (length(imatches) == 0) && (i>1) && continue
      (length(rangematches) < 1 ) && (j>1) && continue
      (length(rangematches) < 2 ) && (k>1) && continue
      (length(rangematches) < 3 ) && (l>1) && continue

      txtcmd = tempcommand |> x -> replace(x, "!}}"=>"}}")

        for v in imatches
          txtcmd = subout(txtcmd, i, v, mykeys)
        end

        if length(rangematches)  >= 1
            splitthenskipcheck(rangematches[1], i, j) && continue
            txtcmd = replace(txtcmd, "{{" * rangematches0[1] * "}}" => String(mykeys[j]))
        end
        if length(rangematches)  >= 2
            splitthenskipcheck(rangematches[2], i, k) && continue
            txtcmd = replace(txtcmd, "{{" * rangematches0[2] * "}}" => String(mykeys[k]))
        end
        if length(rangematches)  >= 3
            splitthenskipcheck(rangematches[3], i, l) && continue
            txtcmd = replace(txtcmd, "{{" * rangematches0[3] * "}}" => String(mykeys[l]))
        end

        occursin("~~OUTOFBOUNDS~~", txtcmd) & !occursin("!}}", tempcommand) && continue

        if occursin("~~OUTOFBOUNDS~~", txtcmd) && occursin("!}}", tempcommand)
           txtcmd = replace(txtcmd, "~~OUTOFBOUNDS~~"=> "999" )
           # push!(collection, fill(false, length(logicset[:])))
        end
        # verbose &&  println(prefix * "$txtcmd")

        ℧∇ = logicset[:]

        try
           verbose && push!(txtoutarray, "`$prefix` $txtcmd")
           ℧∇ = superoperator(txtcmd, logicset, verbose=verbose)[:]
           # verbose && (txtoutarray[end] *= " ✔")
        catch
           verbose && (txtoutarray[end] *= " X")
        end

       push!(collection, ℧∇)
       (length(matches_dot)==0) && push!(iset, mykeys[i])
       (length(matches_dot)>0)  && push!(iset, Symbol(string(mykeys[i]) * matches_dot[1]))
    end

    !printall && (length(txtoutarray)>4) &&
      (txtoutarray = [txtoutarray[1:2]..., " ... ", txtoutarray[end .- (1:-1:0)]...])

    verbose && push!(txtoutarray, constraintprint)

    verbose && (printmarkdown(join(txtoutarray, " ")); println())

    collector = hcat(collection...)

    returnlogical && return [collector, iset]

    if (countrange === missing)
      ℧Δ = logicset[:] .& [all(collector[i,:]) for i in 1:size(collector)[1]]
    else
      ℧Δ = logicset[:] .& [sum(collector[i,:]) ∈ countrange for i in 1:size(collector)[1]]
    end

    logicsetcopy[:] = ℧Δ
    logicsetcopy
end


text:
function createSampleLTISystem(; uses::Bool = false, usek::Bool = false, usee::Bool = false, usef::Bool =false )
    A = [1.0 1.0;
         0.0 1.0]
    B = [0.0;
         1.0]
    C = [1.0 0.0;
         0.0 1.0]
    D = [0.0;
         0.0]
    Ts = 0.1
    sys = StateSpace( A, B, C, D, Ts )

    Q = [1.0 1.0;
         1.0 1.0]
    R = [1.0]

    E = Array{Float64}(undef, 0, 2)
    F = Array{Float64}(undef, 0, 1)
    g = Array{Float64}(undef, 0)
    K = Array{Float64}(undef, 0, 0)

    if( uses )
      S = [1.0;
           0.0]
    else
      S = [ 0.0;
            0.0]
    end

    if( usek )
        # Define a simple controller that has both poles at 0.5
        K = [ 0.25 1.0 ]
    else
        K = [ 0.0 0.0 ]
    end

    if( usee )
        # Define a set of bound constraints between -1 <= u <= 1
        E = [ 0.0 0.0;
              0.0 0.0 ]
        F = [ 1.0;
             -1.0 ]
        g = [ g;
              1.0;
              1.0]
    end

    if( usef )
        # Define a set of bounds and an affine constraint
        E = [ E;
              1.0  0.0;
             -1.0  0.0;
              0.0  1.0;
              0.0 -1.0;
              1.0  1.0 ]
        F = [ F;
              0.0;
              0.0;
              0.0;
              0.0;
              0.0 ]
        g = [ g;
              1.0;
              1.0;
              1.0;
              1.0;
              4.0 ]
    end

    if( usef && usee )
        # Create a coupling constraint between the states and inputs
        F = [ F;
              1.0 ]
        E = [ E;
              1.0 0.0 ]
        g = [ g;
              1.0 ]
    end

    return ConstrainedTimeInvariantLQR( sys, 10, Q, R, :Q, S=S, K=K, E=E, F=F, g=g )

end


text:
################################################################################
# Squencial processing of multiple counts matrices
################################################################################

"""
Multiple Counts Distributed Multinomial Regression by running independent poisson gamma lasso regression to each column of counts,
picks the minimum AICc segement of each path, and returns a coefficient matrix (coefs) representing point estimates
for the entire multinomial (includes the intercept if one was included).
Unlike dmr(), mcdmr() takes a vector of counts matrices and sequencially
  1. regresses each counts matrix on the covars genrating an srproj Z (includes total counts m)
  2. accumulated srproj z into Z
  3. uses [covars Z] for subsequent regressions
Returns the accumulated Z matrix and a vector of coefficient matrices
"""
function mcdmr(covars::AbstractMatrix{T},multicounts::Vector,projdir::Int;
    verbose=true, kwargs...) where {T<:AbstractFloat}
  L = length(multicounts)
  n,p = size(covars)
  verbose && @info("fitting mcdmr to $L counts matrices")
  multicoefs = Vector{DMRCoefs}(undef,L)
  Z = Array{T}(undef,n,0)
	for l=1:L
		verbose && @info("fitting dmr to counts matrix #$l on $p covars + $(size(Z,2)) previous SR projections ...")
		# get l'th counts matrix
		counts = multicounts[l]
		# fit dmr
		multicoefs[l] = dmr([covars Z],counts; verbose=verbose, kwargs...)
		# collapse counts into low dimensional SR projection
		z = srproj(multicoefs[l],counts,projdir)
		# append z to Z
	  Z = [Z z]
	end
  Z, multicoefs
end


text:
using Waterfall

# include("")
include(joinpath(WATERFALL_DIR,"src","includes.jl"))
include(joinpath(WATERFALL_DIR,"bin","io.jl"))

# List of permutations.
perms = vcat([[
    ii,
    swapat!(copy(ii), 1, 2),
] for ii in [collect(1:4)]]...)

# !!!! Add support for defining a legend to violin plots.
for T in [Horizontal, Vertical, Parallel]
    for nsample in [1,10,50]
        for colorcycle in [true,false]
            for perm in perms

                (length(perm)>length(COLORCYCLE) && colorcycle) && continue

                # Define keyword arguments specific to this iteration.
                # By default, permute and correlate kwargs = true.
                global locals = (
                    # What modifications should be made to the input data?
                    nsample = nsample,
                    permutation = perm,
                    # How should the plot be formatted?
                    colorcycle = colorcycle,
                    ylabel = "Efficiency (%)",
                )

                global cascade = define_from(Cascade{Data}, df; locals..., kwargs...)
                global plot = define_from(Plot{T}, copy(cascade); locals..., kwargs...)

                Luxor.@png begin
                    Luxor.fontface("Gill Sans")
                    Luxor.fontsize(FONTSIZE)
                    Luxor.setline(1.0)
                    Luxor.setmatrix([1 0 0 1 LEFT_BORDER TOP_BORDER])

                    draw(plot)

                end WIDTH+LEFT_BORDER+RIGHT_BORDER HEIGHT+TOP_BORDER+BOTTOM_BORDER+padding(plot.axes[1]) plot.path
            end
        end
    end
end


text:
function save_csv(file_path::String, group_path, data; kwargs...)
    if group_path != ""
        @warn "Group path ignored when saving to .csv file."
    end
    path, file = splitdir(file_path)
    mkpath(path)
    CSV.write(file_path, data; kwargs...)
end

function load_csv(file_path::String, group_path; kwargs...)
    if group_path != ""
        @warn "Group path ignored when loading .csv file."
    end
    CSV.File(file_path; kwargs...) |> DataFrame
end


text:
# Julia wrapper for header: /Volumes/BOOTCAMP/programs/compa_libs/thelib/dll/laszip_dll.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

using Libdl

@static Sys.iswindows() ? (Sys.WORD_SIZE == 64 ? (const suffixed_name = "laszip_w64") : (const suffixed_name = "laszip_w32")) : (const suffixed_name = "laszip")
const thelib = Libdl.find_library(["laszip", suffixed_name])
#const thelib = "V:/laszip/build/bin/laszip_w64"

if isempty(thelib)
    error("could not find thelib library")
end

function laszip_get_version(version_major::Ptr{UInt8},version_minor::Ptr{UInt8},version_revision::Ptr{UInt16},version_build::Ptr{UInt32})
    ccall((:laszip_get_version,thelib),Cint,(Ptr{UInt8},Ptr{UInt8},Ptr{UInt16},Ptr{UInt32}),version_major,version_minor,version_revision,version_build)
end

function laszip_create(pointer::Ptr{Ptr{Cvoid}})
    ccall((:laszip_create,thelib),Cint,(Ptr{Ptr{Cvoid}},),pointer)
end

function laszip_get_error(pointer::Ptr{Cvoid}, error::Ptr{Ptr{UInt8}})
    ccall((:laszip_get_error,thelib),Cint,(Ptr{Cvoid},Ptr{Ptr{UInt8}}),pointer,error)
end
function laszip_get_error(pointer::Ptr{Cvoid}, error::Ptr{String})
    ccall((:laszip_get_error,thelib),Cint,(Ptr{Cvoid},Ptr{Ptr{UInt8}}),pointer,error)
end

function laszip_get_warning(pointer::Ptr{Cvoid},warning::Ptr{Ptr{UInt8}})
    ccall((:laszip_get_warning,thelib),Cint,(Ptr{Cvoid},Ptr{Ptr{UInt8}}),pointer,warning)
end

function laszip_clean(pointer::Ptr{Cvoid})
    ccall((:laszip_clean,thelib),Cint,(Ptr{Cvoid},),pointer)
end

function laszip_destroy(pointer::Ptr{Cvoid})
    ccall((:laszip_destroy,thelib),Cint,(Ptr{Cvoid},),pointer)
end

function laszip_get_header_pointer(pointer::Ptr{Cvoid},header_pointer::Ptr{Ptr{laszip_header}})
    ccall((:laszip_get_header_pointer,thelib),Cint,(Ptr{Cvoid},Ptr{Ptr{laszip_header}}),pointer,header_pointer)
end

function laszip_get_point_pointer(pointer::Ptr{Cvoid},point_pointer::Ptr{Ptr{laszip_point}})
    ccall((:laszip_get_point_pointer,thelib),Cint,(Ptr{Cvoid},Ptr{Ptr{laszip_point}}),pointer,point_pointer)
end

function laszip_get_point_count(pointer::Ptr{Cvoid},count::Ptr{Clonglong})
    ccall((:laszip_get_point_count,thelib),Cint,(Ptr{Cvoid},Ptr{Clonglong}),pointer,count)
end

function laszip_set_header(pointer::Ptr{Cvoid},header::Ptr{laszip_header})
    ccall((:laszip_set_header,thelib),Cint,(Ptr{Cvoid},Ptr{laszip_header}),pointer,header)
end

function laszip_check_for_integer_overflow(pointer::Ptr{Cvoid})
    ccall((:laszip_check_for_integer_overflow,thelib),Cint,(Ptr{Cvoid},),pointer)
end

function laszip_auto_offset(pointer::Ptr{Cvoid})
    ccall((:laszip_auto_offset,thelib),Cint,(Ptr{Cvoid},),pointer)
end

function laszip_set_point(pointer::Ptr{Cvoid},point::Ptr{laszip_point})
    ccall((:laszip_set_point,thelib),Cint,(Ptr{Cvoid},Ptr{laszip_point}),pointer,point)
end

function laszip_set_coordinates(pointer::Ptr{Cvoid},coordinates::Ptr{Cdouble})
    ccall((:laszip_set_coordinates,thelib),Cint,(Ptr{Cvoid},Ptr{Cdouble}),pointer,coordinates)
end

function laszip_get_coordinates(pointer::Ptr{Cvoid},coordinates::Ptr{Cdouble})
    ccall((:laszip_get_coordinates,thelib),Cint,(Ptr{Cvoid},Ptr{Cdouble}),pointer,coordinates)
end

function laszip_set_geokeys(pointer::Ptr{Cvoid},number::UInt32,key_entries::Ptr{laszip_geokey})
    ccall((:laszip_set_geokeys,thelib),Cint,(Ptr{Cvoid},UInt32,Ptr{laszip_geokey}),pointer,number,key_entries)
end

function laszip_set_geodouble_params(pointer::Ptr{Cvoid},number::UInt32,geodouble_params::Ptr{Cdouble})
    ccall((:laszip_set_geodouble_params,thelib),Cint,(Ptr{Cvoid},UInt32,Ptr{Cdouble}),pointer,number,geodouble_params)
end

function laszip_set_geoascii_params(pointer::Ptr{Cvoid},number::UInt32,geoascii_params::Ptr{UInt8})
    ccall((:laszip_set_geoascii_params,thelib),Cint,(Ptr{Cvoid},UInt32,Ptr{UInt8}),pointer,number,geoascii_params)
end

function laszip_add_vlr(pointer::Ptr{Cvoid},vlr::Ptr{laszip_vlr})
    ccall((:laszip_add_vlr,thelib),Cint,(Ptr{Cvoid},Ptr{laszip_vlr}),pointer,vlr)
end

function laszip_open_writer(pointer::Ptr{Cvoid},file_name::Ptr{UInt8},compress::Int)
    ccall((:laszip_open_writer,thelib),Cint,(Ptr{Cvoid},Ptr{UInt8},Cint),pointer,file_name,compress)
end
function laszip_open_writer(pointer::Ptr{Cvoid}, file_name::AbstractString, compress::Int)
    ccall((:laszip_open_writer,thelib),Cint,(Ptr{Cvoid},Ptr{UInt8},Cint),pointer,file_name,compress)
end

function laszip_write_point(pointer::Ptr{Cvoid})
    ccall((:laszip_write_point,thelib),Cint,(Ptr{Cvoid},),pointer)
end

function laszip_close_writer(pointer::Ptr{Cvoid})
    ccall((:laszip_close_writer,thelib),Cint,(Ptr{Cvoid},),pointer)
end

function laszip_open_reader(pointer::Ptr{Cvoid},file_name::Ptr{UInt8},is_compressed::Ptr{Cint})
    ccall((:laszip_open_reader,thelib),Cint,(Ptr{Cvoid},Ptr{UInt8},Ptr{Cint}),pointer,file_name,is_compressed)
end
function laszip_open_reader(pointer::Ptr{Cvoid},file_name::AbstractString, is_compressed::Ptr{Cint})
    ccall((:laszip_open_reader,thelib),Cint,(Ptr{Cvoid},Ptr{UInt8},Ptr{Cint}),pointer,file_name,is_compressed)
end

function laszip_seek_point(pointer::Ptr{Cvoid},index::Clonglong)
    ccall((:laszip_seek_point,thelib),Cint,(Ptr{Cvoid},Clonglong),pointer,index)
end

function laszip_read_point(pointer::Ptr{Cvoid})
    ccall((:laszip_read_point,thelib),Cint,(Ptr{Cvoid},),pointer)
end

function laszip_close_reader(pointer::Ptr{Cvoid})
    ccall((:laszip_close_reader,thelib),Cint,(Ptr{Cvoid},),pointer)
end

function create_empty_header()
    ccall((:create_empty_header,thelib),laszip_header,())
end
function create_empty_point()
    ccall((:create_empty_point,thelib),laszip_point,())
end
function create_empty_vlr()
    ccall((:create_empty_vlr,thelib),laszip_vlr,())
end

max_stars_repo_name: joa-quim/Laszip

max_issues_repo_name: joa-quim/Laszip

max_forks_repo_name: joa-quim/Laszip


text:
using Test, SpecialFunctions
import TruncatedNormal as TruncNorm

@test TruncNorm.tnmean(-Inf, Inf) == 0
@test TruncNorm.tnmean(0, +Inf) ≈ +√(2/π)
@test TruncNorm.tnmean(-Inf, 0) ≈ -√(2/π)
@test TruncNorm.tnmean(-Inf, -Inf) == -Inf
@test TruncNorm.tnmean(+Inf, +Inf) == +Inf

for x = -10:10
    @test TruncNorm.tnmean(x, +Inf) ≈ +√(2/π) / erfcx(+x/√2)
    @test TruncNorm.tnmean(-Inf, x) ≈ -√(2/π) / erfcx(-x/√2)
end

for x = exp.(-100:100)
    @test TruncNorm.tnmean(-x, x) ≈ 0
    @test TruncNorm.tnmean(0, x) ≈ -√(2/π) * expm1(-x^2 / 2) / erf(x / √2)
end

@test TruncNorm.tnmean(1e-44, 1e-43) ≈ 5.4999999999999999999999999999999999999999e-44

@test TruncNorm.tnmean(-Inf, Inf, 0, 1) == TruncNorm.tnmean(-Inf, Inf)
for x = -100:100
    @test TruncNorm.tnmean(x, +Inf, 0, 1) ≈ TruncNorm.tnmean(x, +Inf)
    @test TruncNorm.tnmean(-Inf, x, 0, 1) ≈ TruncNorm.tnmean(-Inf, x)
end

@test TruncNorm.tnmean(-1e4, 1e4, 0, 1) == TruncNorm.tnmean(-1e4, 1e4)

for x = -100:100, y = x + 1 : 100
    @test TruncNorm.tnmean(x, y, 0, 1) == TruncNorm.tnmean(x, y)
end

@test TruncNorm.tnmean(100, 115) ≈ 100.00999800099926070518490239457545847490332879043
@test TruncNorm.tnmean(-1e6, -999000) ≈ -999000.00000100100100099899498898098
@test TruncNorm.tnmean(+1e6, +Inf) ≈ +1.00000000000099999999999800000e6
@test TruncNorm.tnmean(-Inf, -1e6) ≈ -1.00000000000099999999999800000e6

@test TruncNorm.tnmean(-1e200, 1e200) ≈ 0
@test TruncNorm.tnmean(0, +1e200) ≈ +0.797884560802865355879892119869
@test TruncNorm.tnmean(-1e200, 0) ≈ -0.797884560802865355879892119869

@test TruncNorm.tnmean(50, 70, -2, 3) ≈ 50.171943499898757645751683644632860837133138152489
@test TruncNorm.tnmean(-100.0, 0.0, 0.0, 2.0986317998643735) ≈ -1.6744659119217125058885983754999713622460154892645
@test TruncNorm.tnmean(0.0, 0.9, 0.0, 0.07132755843183151) ≈ 0.056911157632522598806524588414964004271754161737065
@test TruncNorm.tnmean(-100.0, 100.0, 0.0, 17.185261847875548) ≈ 0
@test TruncNorm.tnmean(-100.0, 0.5, 0.0, 0.47383322897860064) ≈ -0.1267981330521791493635176736743283314399
@test TruncNorm.tnmean(-100.0, 100.0, 0.0, 17.185261847875548) ≈ 0

for x = 1e1 .^ (-400:100:400)
    @test TruncNorm.tnmean(-x, x) == TruncNorm.tnmean(-x, x) || isnan(TruncNorm.tnmean(-x, x)) && isnan(TruncNorm.tnmean(-x, x))
end

for a = exp.(-10:10), b = exp.(-10:10)
    a ≤ b || continue
    @test a ≤ TruncNorm.tnmean(a, b) ≤ b
end

# https://github.com/JuliaStats/Distributions.jl/issues/827
@test TruncNorm.tnmean(0, 1000, 1000000, 1) ≈ 999.99999899899899900100501101901899090472046236710608108591983


text:
struct FilterDef
    name ::Symbol
    id ::UInt8
    description ::String
end

struct Filter
    name ::Symbol
    id ::UInt8
    meta ::UInt8
end

max_filters_count() = Lib.BLOSC2_MAX_FILTERS

struct FilterPipeline
    filters ::NTuple{max_filters_count(), Filter}
    function FilterPipeline(filters::Vector)
        full = fill(make_filter(:nofilter), max_filters_count())
        copyto!(full, filters)
        return new((full...,))
    end
end

"""
    filter_pipeline(args...)

Make filter pipeline. Each argument can be filter name (`Symbol`) or pair name=>meta ('Pair{Symbol, Integer}').
Each argument represents a filter in a chain of filters that are sequentially applied during compression
Maximum number of filters is `max_filters_count()`

# Examples

```
filter_pipeline(:trunc_prec=>23, :shuffle)
```
"""
filter_pipeline(args...) = FilterPipeline([filter_element.(args)...,])

filter_element(name::Symbol) = make_filter(name)
filter_element(name_meta::Pair{Symbol, <:Integer}) = make_filter(name_meta[1], name_meta[2])

add_available_filter(name::Symbol, id::Integer, description::AbstractString = "") =
         _available_filters[name] = FilterDef(name, id, description)

function fill_buildin_filters()
    add_available_filter(:nofilter, Lib.BLOSC_NOFILTER, "No filter.")
    add_available_filter(:shuffle, Lib.BLOSC_SHUFFLE, "Byte-wise shuffle.")
    add_available_filter(:bitshuffle, Lib.BLOSC_BITSHUFFLE, "Bit-wise shuffle.")
    add_available_filter(:delta, Lib.BLOSC_DELTA, "Delta filter.")
    add_available_filter(:trunc_prec, Lib.BLOSC_TRUNC_PREC, "Truncate precision filter.")
end

"""
    has_filter(name::Symbol)

Check if filter with name `name` is exists
"""
has_filter(name::Symbol) = haskey(_available_filters, name)

"""
    make_filter(name::Symbol, [meta::UInt8 = 0])

Make filter
"""
function make_filter(name::Symbol, meta::Integer = UInt8(0))
    def = _available_filters[name]
    return Filter(def.name, def.id, meta)
end

"""
    available_filter_names()

Get vector with names of available filters
"""
available_filter_names() = collect(keys(_available_filters))

"""
    filter_description(name::Symbol)

Get desctription of filter `name`
"""
filter_description(name::Symbol) = _available_filters[name].description


text:
struct GitTreeRoot
    path::String
    read::Bool
    write::Bool
end

function GitTreeRoot(path::AbstractString; write=false, read=true)
    path = abspath(path)
    if !isdir(path)
        throw(ArgumentError("$(repr(path)) must be a directory"))
    end
    GitTreeRoot(path, read, write)
end

sys_abspath(root::GitTreeRoot) = root.path
sys_abspath(ap::AbsPath{GitTreeRoot}) = joinpath(sys_abspath(ap.root), sys_joinpath(ap.path))

function Base.open(f::Function, root::GitTreeRoot)
    git(subcmd) = setenv(`git $subcmd`, dir=root.path)
    s = read(git(`status --porcelain`), String)
    isempty(s) || error("Git working copy is dirty")
    result = f(BlobTree(root))
    # FIXME: From the point of view of this code, it seems unnatural to attach
    # `write` to GitTreeRoot.
    if root.write
        run(pipeline(git(`add $(root.path)`), stdout=devnull))
        run(pipeline(git(`commit --allow-empty $(root.path) -m "New version"`), stdout=devnull))
    end
    result
end

#-------------------------------------------------------------------------------
# FIXME: Factor together with BlobTreeRoot

function Base.haskey(tree::BlobTree{GitTreeRoot}, name::AbstractString)
    ispath(sys_abspath(joinpath(tree,name)))
end

function Base.open(func::Function, f::Blob{GitTreeRoot}; write=false, read=!write)
    if !f.root.write && write
        error("Error writing file at read-only path $f")
    end
    open(func, sys_abspath(f); read=read, write=write)
end

function Base.open(func::Function, p::AbsPath{GitTreeRoot}; write=false, read=!write)
    if !p.root.write && write
        error("Error writing file at read-only path $p")
    end
    open(func, sys_abspath(p); read=read, write=write)
end

function Base.mkdir(p::AbsPath{GitTreeRoot}, args...)
    if !p.root.write
        error("Cannot make directory in read-only tree root at $(sys_abspath(p.root))")
    end
    mkdir(sys_abspath(p), args...)
    return BlobTree(p.root, p.path)
end


text:
using UnicodePlots
using Base.Test

# write your own tests here
include("tst_plots.jl")


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Documenter
import TOML

makedocs(
    modules = [TOML],
    sitename = "TOML",
    checkdocs = :strict,
    doctest = true,
    pages = Any[
        "TOML" => "index.md"
        ]
    )

deploydocs(repo = "github.com/KristofferC/TOML.jl.git")

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
@model SDTmodel(data) = begin
    d ~ Normal(0, 1/sqrt(.5))
    c ~ Normal(0, 1/sqrt(2))
    data ~ SDT(d, c)
end

AHMC_SDT(hits, fas, Nd) = SDTmodel((hits, fas, Nd))

AHMCconfig = Turing.NUTS(1000, .85)

CmdStan_SDT = "
data {
  int<lower=1> Nd;
  int<lower=0> hits;
  int<lower=0> fas;
}
parameters {
  real d;
  real c;
}

model {
  real thetah;
  real thetaf;

  d ~ normal(0, inv_sqrt(.5));
  c ~ normal(0, inv_sqrt(2));

  thetah = Phi(d/2-c);
  thetaf = Phi(-d/2-c);

  // Observed counts
  hits ~ binomial(Nd, thetah);
  fas ~ binomial(Nd, thetaf);
}
"

CmdStanConfig = Stanmodel(name="CmdStan_SDT", model=CmdStan_SDT, nchains=1, output_format=:mcmcchains,
   Sample(num_samples=1000, num_warmup=1000, adapt=CmdStan.Adapt(delta=0.8),
   save_warmup=true))

  struct SDTProblem
      hits::Int64
      fas::Int64
      Nd::Int64
  end

  function (problem::SDTProblem)(θ)
      @unpack hits,fas,Nd=problem   # extract the data
      @unpack d,c=θ
      logpdf(SDT(d, c),[hits, fas, Nd]) + logpdf(Normal(0, 1/sqrt(.5)), d) +
      logpdf(Normal(0, 1/sqrt(2)), c)
  end

  # Define problem with data and inits.
  function sampleDHMC(hits, fas, Nd, nsamples, autodiff)
    p = SDTProblem(hits, fas, Nd)
    p((d=2.0, c=.0))
    # Write a function to return properly dimensioned transformation.
    trans = as((d=asℝ, c=asℝ))
    P = TransformedLogDensity(trans, p)
    ∇P = ADgradient(autodiff, P)
    # Sample from the posterior.
    results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, nsamples; reporter = NoProgressReport())
    # Undo the transformation to obtain the posterior from the chain.
    posterior = TransformVariables.transform.(trans, results.chain)
    chns = nptochain(results, posterior)
    return chns
  end

  function simulateSDT(;d=2., c=0., Nd, kwargs...)
    θhit  = cdf(Normal(0, 1), d/2-c)
    θfa = cdf(Normal(0, 1), -d/2-c)
    hits = rand(Binomial(Nd, θhit))
    fas = rand(Binomial(Nd, θfa))
    return (hits=hits, fas=fas, Nd=Nd)
  end


text:
struct LazyCoulomb{T,B<:AbstractQuasiMatrix,I<:Integer} <: AbstractQuasiMatrix{T}
    R::B
    ℓ::I
end

LazyCoulomb(R::B, ℓ::I) where {T,B<:AbstractQuasiMatrix{T}, I<:Integer} =
    LazyCoulomb{T,B,I}(R,ℓ)

function Base.axes(LC::LazyCoulomb)
    a = axes(LC.R,1)
    a,a
end
Base.axes(LC::LazyCoulomb,i) = axes(LC.R,1)

function Base.size(LC::LazyCoulomb)
    m = size(LC.R,1)
    m,m
end
Base.size(LC::LazyCoulomb, i) = size(LC.R,1)

function Base.getindex(LC::LazyCoulomb{T,B,I},i::Integer,j::Integer) where {T,B,I}
    i ≠ j && return zero(T)
    R.x[i]
end

Base.replace_in_print_matrix(::LazyCoulomb, i::Integer, j::Integer, s::AbstractString) =
    i == j ? s : Base.replace_with_centered_mark(s)

const CoulombIntegral{T} = Applied{<:Any, typeof(*),<:Tuple{
    Applied{<:Any, typeof(*), <:Tuple{
        <:Adjoint{T,<:AbstractVector},
        <:AbstractQuasiMatrix}},
    <:LazyCoulomb{T,<:AbstractQuasiMatrix},
    Applied{<:Any, typeof(*), <:Tuple{
        <:AbstractQuasiMatrix,
        <:AbstractVector}}}}

K(ℓ::I,r::T,r̃::T) where {I<:Integer,T} = min(r,r̃)^ℓ/max(r,r̃)^(ℓ+1)

function Base.copyto!(dest::AbstractVector{T}, M::CoulombIntegral{T}) where T
    u,LC,v = M.args
    R = LC.R
    r = locs(R)
    k = Diagonal(similar(r))

    Rv = R'v
    tmp = similar(v)

    tmp.args[2] .= 0

    kv = applied(*, k, Rv)

    for i in eachindex(r)
        k.diag .= K.(LC.ℓ,r,r[i])
        tmp.args[2] .= kv
        # u is already conjugated, hence we don't use the
        # transposition operator to calculate the dot product.
        dest[i:i] .= applied(*, u, tmp)
    end

    dest
end

function Base.similar(M::CoulombIntegral, ::Type{T}) where T
    LC = M.args[2]
    R = LC.R
    Vector{T}(undef, size(R,2))
end
LazyArrays.materialize(M::CoulombIntegral) = copyto!(similar(M, eltype(M)), M)

export LazyCoulomb


text:
using Ferrite
using BlockArrays, SparseArrays, LinearAlgebra

function create_cook_grid(nx, ny)
    corners = [Vec{2}((0.0,   0.0)),
               Vec{2}((48.0, 44.0)),
               Vec{2}((48.0, 60.0)),
               Vec{2}((0.0,  44.0))]
    grid = generate_grid(Triangle, (nx, ny), corners);
    # facesets for boundary conditions
    addfaceset!(grid, "clamped", x -> norm(x[1]) ≈ 0.0);
    addfaceset!(grid, "traction", x -> norm(x[1]) ≈ 48.0);
    return grid
end;

function create_values(interpolation_u, interpolation_p)
    # quadrature rules
    qr      = QuadratureRule{2,RefTetrahedron}(3)
    face_qr = QuadratureRule{1,RefTetrahedron}(3)

    # geometric interpolation
    interpolation_geom = Lagrange{2,RefTetrahedron,1}()

    # cell and facevalues for u
    cellvalues_u = CellVectorValues(qr, interpolation_u, interpolation_geom)
    facevalues_u = FaceVectorValues(face_qr, interpolation_u, interpolation_geom)

    # cellvalues for p
    cellvalues_p = CellScalarValues(qr, interpolation_p, interpolation_geom)

    return cellvalues_u, cellvalues_p, facevalues_u
end;

function create_dofhandler(grid, ipu, ipp)
    dh = DofHandler(grid)
    push!(dh, :u, 2, ipu) # displacement
    push!(dh, :p, 1, ipp) # pressure
    close!(dh)
    return dh
end;

function create_bc(dh)
    dbc = ConstraintHandler(dh)
    add!(dbc, Dirichlet(:u, getfaceset(dh.grid, "clamped"), (x,t) -> zero(Vec{2}), [1,2]))
    close!(dbc)
    t = 0.0
    update!(dbc, t)
    return dbc
end;

struct LinearElasticity{T}
    G::T
    K::T
end

function doassemble(cellvalues_u::CellVectorValues{dim}, cellvalues_p::CellScalarValues{dim},
                    facevalues_u::FaceVectorValues{dim}, K::SparseMatrixCSC, grid::Grid,
                    dh::DofHandler, mp::LinearElasticity) where {dim}

    f = zeros(ndofs(dh))
    assembler = start_assemble(K, f)
    nu = getnbasefunctions(cellvalues_u)
    np = getnbasefunctions(cellvalues_p)

    fe = PseudoBlockArray(zeros(nu + np), [nu, np]) # local force vector
    ke = PseudoBlockArray(zeros(nu + np, nu + np), [nu, np], [nu, np]) # local stiffness matrix

    # traction vector
    t = Vec{2}((0.0, 1/16))
    # cache ɛdev outside the element routine to avoid some unnecessary allocations
    ɛdev = [zero(SymmetricTensor{2, dim}) for i in 1:getnbasefunctions(cellvalues_u)]

    for cell in CellIterator(dh)
        fill!(ke, 0)
        fill!(fe, 0)
        assemble_up!(ke, fe, cell, cellvalues_u, cellvalues_p, facevalues_u, grid, mp, ɛdev, t)
        assemble!(assembler, celldofs(cell), fe, ke)
    end

    return K, f
end;

function assemble_up!(Ke, fe, cell, cellvalues_u, cellvalues_p, facevalues_u, grid, mp, ɛdev, t)

    n_basefuncs_u = getnbasefunctions(cellvalues_u)
    n_basefuncs_p = getnbasefunctions(cellvalues_p)
    u▄, p▄ = 1, 2
    reinit!(cellvalues_u, cell)
    reinit!(cellvalues_p, cell)

    # We only assemble lower half triangle of the stiffness matrix and then symmetrize it.
    @inbounds for q_point in 1:getnquadpoints(cellvalues_u)
        for i in 1:n_basefuncs_u
            ɛdev[i] = dev(symmetric(shape_gradient(cellvalues_u, q_point, i)))
        end
        dΩ = getdetJdV(cellvalues_u, q_point)
        for i in 1:n_basefuncs_u
            divδu = shape_divergence(cellvalues_u, q_point, i)
            δu = shape_value(cellvalues_u, q_point, i)
            for j in 1:i
                Ke[BlockIndex((u▄, u▄), (i, j))] += 2 * mp.G * ɛdev[i] ⊡ ɛdev[j] * dΩ
            end
        end

        for i in 1:n_basefuncs_p
            δp = shape_value(cellvalues_p, q_point, i)
            for j in 1:n_basefuncs_u
                divδu = shape_divergence(cellvalues_u, q_point, j)
                Ke[BlockIndex((p▄, u▄), (i, j))] += -δp * divδu * dΩ
            end
            for j in 1:i
                p = shape_value(cellvalues_p, q_point, j)
                Ke[BlockIndex((p▄, p▄), (i, j))] += - 1/mp.K * δp * p * dΩ
            end

        end
    end

    symmetrize_lower!(Ke)

    # We integrate the Neumann boundary using the facevalues.
    # We loop over all the faces in the cell, then check if the face
    # is in our `"traction"` faceset.
    @inbounds for face in 1:nfaces(cell)
        if onboundary(cell, face) && (cellid(cell), face) ∈ getfaceset(grid, "traction")
            reinit!(facevalues_u, cell, face)
            for q_point in 1:getnquadpoints(facevalues_u)
                dΓ = getdetJdV(facevalues_u, q_point)
                for i in 1:n_basefuncs_u
                    δu = shape_value(facevalues_u, q_point, i)
                    fe[i] += (δu ⋅ t) * dΓ
                end
            end
        end
    end
end

function symmetrize_lower!(K)
    for i in 1:size(K,1)
        for j in i+1:size(K,1)
            K[i,j] = K[j,i]
        end
    end
end;

function solve(ν, interpolation_u, interpolation_p)
    # material
    Emod = 1.
    Gmod = Emod / 2(1 + ν)
    Kmod = Emod * ν / ((1+ν) * (1-2ν))
    mp = LinearElasticity(Gmod, Kmod)

    # grid, dofhandler, boundary condition
    n = 50
    grid = create_cook_grid(n, n)
    dh = create_dofhandler(grid, interpolation_u, interpolation_p)
    dbc = create_bc(dh)

    # cellvalues
    cellvalues_u, cellvalues_p, facevalues_u = create_values(interpolation_u, interpolation_p)

    # assembly and solve
    K = create_sparsity_pattern(dh);
    K, f = doassemble(cellvalues_u, cellvalues_p, facevalues_u, K, grid, dh, mp);
    apply!(K, f, dbc)
    u = Symmetric(K) \ f;

    # export
    filename = "cook_" * (isa(interpolation_u, Lagrange{2,RefTetrahedron,1}) ? "linear" : "quadratic") *
                         "_linear"
    vtk_grid(filename, dh) do vtkfile
        vtk_point_data(vtkfile, dh, u)
    end
    return u,dh
end

linear    = Lagrange{2,RefTetrahedron,1}()
quadratic = Lagrange{2,RefTetrahedron,2}()

u,dh = solve(0.4999999, quadratic, linear);


text:
using DelimitedFiles
using Distributions
using StatsPlots
using Plots.PlotMeasures

#function that reads in data and return log normal parameters + plot
function GetLogNormal(path, titleArg)
    tmp = readdlm(path)
    params = fit(LogNormal, tmp)
    mu = params.μ
    sigma = params.σ
    if titleArg == "Real"
        p1 = histogram(tmp, normalize = :pdf, xlims = (0, 4), ylims = (0,10),
        title = titleArg, xlabel = "Mean Isoforms per Gene per Cell",
        axis = (font(18)), titlefontsize=18, yticks = 0:5:15, left_margin=1.0mm,
        legend = false)
    elseif titleArg == "3 Isoforms" && path != "figures/figure5_data/Weibull_H1_24_scnorm_counts.txt_3_isoforms_array.txt"
        p1 = histogram(tmp, normalize = :pdf, xlims = (0, 4), ylims = (0,10),
        title = titleArg, ylabel = "Normalised Frequency", axis = (font(18)),
        titlefontsize=18, yticks = 0:5:15, left_margin=1.0mm, legend = false)
    else
        p1 = histogram(tmp, normalize = :pdf, xlims = (0, 4), ylims = (0,10),
        title = titleArg, axis = (font(18)), titlefontsize=18, yticks = 0:5:15,
        left_margin=1.0mm, legend = false)
    end
    plot!(LogNormal(mu, sigma), linewidth = 2)

    return [p1, mu, sigma, tmp]
end

################################################################################
# Implementation of expectation maximisation - we are only looking for k
function find_relationships(mu_arr, sigma_arr, k_arr, k, x)
    top = k_arr[k] * pdf(LogNormal(mu_arr[k], sigma_arr[k]), x)

    bottom = 0.0

    for i in 1:length(mu_arr)
        bottom += k_arr[i] * pdf(LogNormal(mu_arr[i], sigma_arr[i]), x)
    end

    #@show top/bottom
    return top/bottom
end

function expectation!(relationships, dim1, dim2, mu_arr, sigma_arr, k_arr, real_arr)
    for point in 1:dim1
        for k in 1:dim2
            #@show point
            #@show k
            relationships[point, k] = find_relationships(mu_arr, sigma_arr,
            k_arr, k, real_arr[point])
        end
    end
end

function maximisation_k!(relationships, k_arr)

    for i in 1:length(k_arr)
        k_arr[i] = sum(relationships[:, i]) / length(relationships[:,i])
    end
end

################################################################################
#These functions are not for actual use but are convenient for testing
function maximisation_mu!(relationships, mu_arr, x)

    for i in 1:length(mu_arr)
        mu_arr[i] = (relationships[:, i] .* x) ./ (relationships[:, i])
    end
end

function maximisation_sigma!(relationships, sigma_arr, mu_arr, x)

    for i in 1:length(sigma_arr)
        sigma_arr[i] = (relationships[:, i] .* (x .- mu_arr[i]).^2) ./
        relationships[:, i]
    end
end

max_stars_repo_name: jenni-westoby/Obstacles

max_issues_repo_name: jenni-westoby/Obstacles

max_forks_repo_name: jenni-westoby/Obstacles


text:
using Tree
using Base.Test

# write your own tests here
tree(".")


text:
using Test, BenchmarkTools
include("../src/dqmc_framework.jl") # to be replaced by using DQMC or similar

@testset "All Speed Tests" begin

    # Linear algebra
    include("speedtests_linalg.jl")

end

nothing

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
using Revise
using JLD2
# using ViewerGL
using LarSurf
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

# created by surface_extraction_parallel_ircad01.jl
# @load "liver01.jld2" V FV

# @JLD2.load "ircad_$(mask_label).jld2" V FV
# @JLD2.load "ircad_$(mask_label)_tri.jld2" V FVtri
mask_labels=["liver", "portalvein"]
# mask_labels=["liver"]
# pth = Pio3d.datasets_join_path("medical/orig/3Dircadb1.$data_id/MASKS_DICOM/liver")

for mask_label in mask_labels

	@JLD2.load "ircad_$(mask_label)_sm.jld2" Vs FVtri
	println("size FVtri $(size(FVtri))")

	objlines = LarSurf.Lar.lar2obj(Vs, FVtri, "ircad_$(mask_label)_taubin.obj")
	@info "length of obj $(size(objlines))"

	# save also non triangulated version
	@info "save data withouth smoothing"
	@JLD2.load "ircad_$(mask_label).jld2" V FV
	objlines = LarSurf.Lar.lar2obj(V, FVtri, "ircad_$(mask_label).obj")
end
# print(obj)


text:
function _wc_helper(x::AbstractArray{Float64,2})
    view(x,:,1) ./ sum(x, dims = 2)
end

#= dont think we need the full version
function _wc_full(data::PopData)
    loc_names = loci(data)
    n_loci = length(loc_names)
    n_pops = length(unique(data.meta.population))
    # get genotype counts
    per_locpop = groupby(data.loci, [:locus, :population])
    nonmissing_counts = DataFrames.combine(per_locpop, :genotype => nonmissing => :n)
    # reshape as a matrix of loci x pop (row x col)
    # collect?
    n_per_locpop = reshape(nonmissing_counts.n', (n_pops, n_loci))'
    n_total = reduce(+, eachcol(n_per_locpop))
    # screen for completely absent loci
    missing_loc = findall(iszero, n_total)
    present_loc = findall(!iszero, n_total)
    if !isempty(missing_loc)
        # index locus names by the missing indices
        tmp_data = exclude(data, locus = loc_names[missing_loc])
        loc_names = loci(tmp_data)
        n_loci = length(loc_names)
        per_locpop = groupby(tmp_data.loci, [:locus, :population])
        nonmissing_counts = DataFrames.combine(per_locpop, :genotype => nonmissing => :n)
        # reshape as a matrix of loci x pop (row x col)
        # collect?
        n_per_locpop = reshape(nonmissing_counts.n', (n_pops, length(present_loci)))'
        n_total = reduce(+, eachcol(n_per_locpop))
    else
        tmp_data = data
    end
    n_pop_per_loc = map(count_nonzeros, eachrow(n_per_locpop))
    per_loc = groupby(tmp_data.loci, :locus)
    allele_counts = DataFrames.combine(
        per_loc,
        :genotype => allele_count => :allele_count,  # allele counts
        :genotype => (i -> allele_freq(i)) => :freqs,          # global allele freqs
    )
    # allele freqs per locus per population
    alle_freqs = DataFrames.combine(per_locpop, :genotype => allele_freq => :freqs)
    # expand out the n matrix to be same dimensions as unique_alleles x pop
    n_expanded = reduce(hcat, repeat.(eachrow(n_per_locpop), 1, allele_counts.allele_count)) |> permutedims
    # expand n_total matrix to be same dimensions as unique_alleles x pop
    n_tot_expanded = reduce(vcat, repeat.(eachrow(n_total), allele_counts.allele_count))
    # calculate corrected n per locus
    corr_n_per_loc = (n_total .- (sum(n_per_locpop .^2, dims = 2) ./ n_total)) ./ (n_pop_per_loc .- 1)
    # expand corr_n matrix to be same dimensions as unique_alleles x pop
    corr_n_per_loc_exp = reduce(vcat, repeat.(eachrow(corr_n_per_loc), allele_counts.allele_count))
    # list of alleles in each locus
    _alleles_perloc = map(j -> collect(keys(j)), allele_counts.freqs)
    # extremely convoluted, creates a big matrix of allele freqs per locus per population
    # TODO are there too many reshapes going on?
    #TODO move into its own function? This seems like it could be a recurring thing
    afreq_tmp = permutedims(reshape(alle_freqs.freqs, n_pops, :))
    allele_freq_pop = reshape(
        reduce(vcat,
            map(zip(eachrow(afreq_tmp), _alleles_perloc)) do (_freqs_p, _alle)
                reduce(hcat,
                    map(_freqs_p) do _freqs
                        [get(_freqs, i, 0.) for i in _alle]     # query the dict of allele freqs
                    end
                )
            end
            ),
       :, n_pops  # reshape by these dimensions
    )

    #_alleles = map(keys, allele_counts.freqs) |> Base.Iterators.flatten |> collect
    _freqs = map(values, allele_counts.freqs) |> Base.Iterators.flatten |> collect

    # # heterozygotes per allele per locus per population

    # gets emptied from the popfirst! calls below(!)
    genos_vec = [i.genotype for i in per_locpop]

    # create matrix of heterozygote occurrences per allele per pop
    het_mtx = reduce(vcat,     # vcat will concatenate the returned matrices into a single matrix
        map(_alleles_perloc) do _alleles          # map across the vector of alleles for each locus
            reduce(hcat,
                # each element in x is locus × population, so use a comprehension to
                # do counthet() as many times as there are populations, popfirst!'ing
                # the first element of x until it's ultimately empty
                # then concatenate it into a matrix
                [counthet(popfirst!(genos_vec), _alleles) for pop_n in 1: n_pops]
            )
        end
    )
    μ_het =  (2 * n_expanded .* allele_freq_pop - het_mtx) / 2
    SSG = sum(n_expanded .* allele_freq_pop - μ_het, dims = 2)
    #SSG = map(sum, eachrow(n_expanded .* allele_freq_pop - μ_het))
    SSi = sum(n_expanded .* (allele_freq_pop - 2 * allele_freq_pop .^ 2) + μ_het, dims = 2)
    #SSi = map(sum, eachrow(n_expanded .* (allele_freq_pop - 2 * allele_freq_pop .^ 2) + μ_het))

    SSP = 2 .* sum(n_expanded .* reduce(hcat, map(i -> (i .- _freqs) .^ 2, eachcol(allele_freq_pop))), dims = 2)
    #SSP = 2 .* map(sum, eachrow(n_expanded .* reduce(hcat, map(i -> (i .- _freqs) .^ 2, eachcol(allele_freq_pop)))))
    n_correction = reduce(vcat, fill.(n_pop_per_loc, allele_counts.allele_count))
    #return n_tot_expanded
    MSG = SSG ./ n_tot_expanded
    MSP = SSP ./ (n_correction .- 1)
    MSI = SSi ./ (n_tot_expanded - n_correction)
    σ_w = MSG
    σ_b = 0.5 * (MSI - MSG)
    σ_a = 0.5 ./ corr_n_per_loc_exp .* (MSP - MSI)
    σ = hcat(σ_a, σ_b, σ_w)

    fst_allele = _wc_helper(σ)
    fis_allele = _wc_helper(@view σ[:,2:end])
    #return fst_allele, fis_allele
    σ_loci = mapreduce(i -> sum(i, dims = 1), vcat, partitionarray(σ, allele_counts.allele_count))
    σ_total = sum(σ, dims = 1)
    fst_total = @views σ_total[1] / sum(σ_total)
    fis_total = @views σ_total[2] / sum(σ_total[:,2:end])
    return fst_total, fis_total

end
=#

function weircockerham1984(data::AbstractDataFrame)
    loc_names = unique(data.locus)
    n_loci = length(loc_names)
    n_pops = length(unique(data.population))
    # get genotype counts
    per_locpop = groupby(data, [:locus, :population])
    nonmissing_counts = DataFrames.combine(per_locpop, :genotype => nonmissing => :n)
    # reshape as a matrix of loci x pop (row x col)
    # collect?
    n_per_locpop = reshape(nonmissing_counts.n', (n_pops, n_loci))'
    n_total = reduce(+, eachcol(n_per_locpop))
    # screen for completely absent loci
    missing_loc = findall(iszero, n_total)
    if !isempty(missing_loc)
        present_loc = n_loci - length(missing_loc)
        # index locus names by the missing indices
        tmp_data = filter(:locus => loc -> loc ∉ loc_names[missing_loc], data, view = true)
        #tmp_data = exclude(data, locus = loc_names[missing_loc])
        loc_names = unique(tmp_data.locus)
        n_loci = length(loc_names)
        per_locpop = groupby(tmp_data, [:locus, :population])
        nonmissing_counts = DataFrames.combine(per_locpop, :genotype => nonmissing => :n)
        # reshape as a matrix of loci x pop (row x col)
        # collect?
        n_per_locpop = reshape(nonmissing_counts.n', (n_pops, present_loci))'
        n_total = reduce(+, eachcol(n_per_locpop))
    else
        tmp_data = data
    end
    n_pop_per_loc = map(count_nonzeros, eachrow(n_per_locpop))
    per_loc = groupby(tmp_data, :locus)
    allele_counts = DataFrames.combine(
        per_loc,
        :genotype => allele_count => :allele_count,  # allele counts
        :genotype => allele_freq => :freqs,          # global allele freqs
    )
    # allele freqs per locus per population
    alle_freqs = DataFrames.combine(per_locpop, :genotype => allele_freq => :freqs)
    # expand out the n matrix to be same dimensions as unique_alleles x pop
    n_expanded = reduce(hcat, repeat.(eachrow(n_per_locpop), 1, allele_counts.allele_count)) |> permutedims
    # expand n_total matrix to be same dimensions as unique_alleles x pop
    n_tot_expanded = reduce(vcat, repeat.(eachrow(n_total), allele_counts.allele_count))
    # calculate corrected n per locus
    corr_n_per_loc = (n_total .- (sum(n_per_locpop .^2, dims = 2) ./ n_total)) ./ (n_pop_per_loc .- 1)
    # expand corr_n matrix to be same dimensions as unique_alleles x pop
    corr_n_per_loc_exp = reduce(vcat, repeat.(eachrow(corr_n_per_loc), allele_counts.allele_count))
    # list of alleles in each locus
    _alleles_perloc = map(j -> collect(keys(j)), allele_counts.freqs)
    # extremely convoluted, creates a big matrix of allele freqs per locus per population
    # TODO are there too many reshapes going on?
    #TODO move into its own function? This seems like it could be a recurring thing
    afreq_tmp = permutedims(reshape(alle_freqs.freqs, n_pops, :))
    allele_freq_pop = reshape(
        reduce(vcat,
            map(zip(eachrow(afreq_tmp), _alleles_perloc)) do (_freqs_p, _alle)
                reduce(hcat,
                    map(_freqs_p) do _freqs
                        [get(_freqs, i, 0.) for i in _alle]     # query the dict of allele freqs
                    end
                )
            end
            ),
       :, n_pops  # reshape by these dimensions
    )

    #_alleles = map(keys, allele_counts.freqs) |> Base.Iterators.flatten |> collect
    _freqs = map(values, allele_counts.freqs) |> Base.Iterators.flatten |> collect

    #heterozygotes per allele per locus per population
    # gets emptied from the popfirst! calls below(!)
    genos_vec = [i.genotype for i in per_locpop]

    # create matrix of heterozygote occurrences per allele per pop
    het_mtx = reduce(vcat,     # vcat will concatenate the returned matrices into a single matrix
        map(_alleles_perloc) do _alleles          # map across the vector of alleles for each locus
            reduce(hcat,
                # each element in x is locus × population, so use a comprehension to
                # do counthet() as many times as there are populations, popfirst!'ing
                # the first element of x until it's ultimately empty
                # then concatenate it into a matrix
                [counthet(popfirst!(genos_vec), _alleles) for pop_n in 1: n_pops]
            )
        end
    )
    μ_het =  (2 * n_expanded .* allele_freq_pop - het_mtx) / 2
    SSG = sum(n_expanded .* allele_freq_pop - μ_het, dims = 2)
    SSi = sum(n_expanded .* (allele_freq_pop - 2 * allele_freq_pop .^ 2) + μ_het, dims = 2)
    SSP = 2 .* sum(n_expanded .* reduce(hcat, map(i -> (i .- _freqs) .^ 2, eachcol(allele_freq_pop))), dims = 2)
    n_correction = reduce(vcat, fill.(n_pop_per_loc, allele_counts.allele_count))

    MSG = SSG ./ n_tot_expanded
    MSP = SSP ./ (n_correction .- 1)
    MSI = SSi ./ (n_tot_expanded - n_correction)
    σ_w = MSG
    σ_b = 0.5 * (MSI - MSG)
    σ_a = 0.5 ./ corr_n_per_loc_exp .* (MSP - MSI)
    σ = hcat(σ_a, σ_b, σ_w)
    σ_total = sum(σ, dims = 1)
    σ_total
    fst_total = σ_total[1] / sum(σ_total)
end


text:
include("compare_with_experimets.jl")
include("prepare_metnet.jl")


text:
struct sCKSAFTParam <: EoSParam
    segment::SingleParam{Float64}
    sigma::PairParam{Float64}
    epsilon::PairParam{Float64}
    epsilon_assoc::AssocParam{Float64}
    bondvol::AssocParam{Float64}
end

abstract type sCKSAFTModel <: CKSAFTModel end
@newmodel sCKSAFT sCKSAFTModel sCKSAFTParam
export sCKSAFT
function sCKSAFT(components;
    idealmodel=BasicIdeal,
    userlocations=String[],
    ideal_userlocations=String[],
    verbose=false,
    assoc_options = AssocOptions())

    params,sites = getparams(components, ["SAFT/sCKSAFT","properties/molarmass.csv"]; userlocations=userlocations, verbose=verbose)
    segment = params["m"]
    k = params["k"]
    sigma = params["vol"]
    sigma.values .*= 6*0.74048/N_A/1e6/π
    sigma.values .^= 1/3
    sigma = sigma_LorentzBerthelot(sigma)
    epsilon = epsilon_LorentzBerthelot(params["epsilon"], k)
    epsilon_assoc = params["epsilon_assoc"]
    bondvol = params["bondvol"]

    packagedparams = sCKSAFTParam(segment, sigma, epsilon, epsilon_assoc, bondvol)
    references = ["TODO sCKSAFT", "TODO sCKSAFT"]

    model = sCKSAFT(packagedparams, sites, idealmodel; ideal_userlocations, references, verbose, assoc_options)
    return model
end

function a_disp(model::sCKSAFTModel, V, T, z)
    ∑z = ∑(z)
    m = model.params.segment.values
    m̄ = dot(z, m)/∑z
    vs = V/(N_A*∑z*m̄)
    it = @comps
    v̄Ȳ = ∑(z[i]*z[j]*m[i]*m[j]*(@f(d,i,j)^3/√(2))*(exp(@f(u,i,j)/T/2)-1) for i ∈ it for j ∈ it)/∑(z[i]*z[j]*m[i]*m[j] for i ∈ it for j ∈ it)
    #res = 36*log(vs/(vs+v̄Ȳ))
    #this formulation is prone to horrible floating point issues.
    #writing this term in log1p form solves the problem.
    res = -36*log1p(v̄Ȳ/vs)
    return res
end

function d(model::sCKSAFTModel, V, T, z, i)
    ϵ = model.params.epsilon.diagvalues
    σ = model.params.sigma.diagvalues
    res = σ[i] * (1 - 0.333exp(-3ϵ[i]/T))
    return res
end

function u(model::sCKSAFTModel, V, T, z, i, j)
    ϵ0 = model.params.epsilon.values[i,j]
    res = ϵ0*(1-10/T)
    return res
end

function Δ(model::sCKSAFTModel, V, T, z, i, j, a, b)
    ϵ_assoc = model.params.epsilon_assoc.values[i,j][a,b]
    κ = model.params.bondvol.values[i,j][a,b]
    g = @f(g_hsij,i,j)
    res = g*@f(d,i,j)^3*(exp(ϵ_assoc/T)-1)*κ
    return res
end