text:
# QP example
#
#    minimize 2 x^2 + y^2 + xy + x + y
#
#       s.t.  x, y >= 0
#             x + y = 1
#
#    solution: (0.25, 0.75), objv = 1.875
#

using CPLEX

env = CPLEX.Env()

model = CPLEX.Model(env, "qp_02")

CPLEX.add_vars!(model, [1., 1.], 0., Inf)

CPLEX.add_qpterms!(model, [1, 1, 2], [1, 2, 2], [2., 1., 1.])
CPLEX.add_constr!(model, [1., 1.], '=', 1.)

CPLEX.optimize!(model)

println("sol = $(CPLEX.get_solution(model))")
println("obj = $(CPLEX.get_objval(model))")


text:
using SparseArrays
using ModelingToolkit
import ModelingToolkit: bareiss!, find_pivot_col, bareiss_update!, swaprows!
import Base: swapcols!

function det_bareiss!(M)
    parity = 1
    _swaprows!(M, i, j) = (i != j && (parity = -parity); swaprows!(M, i, j))
    _swapcols!(M, i, j) = (i != j && (parity = -parity); swapcols!(M, i, j))
    # We only look at the last entry, so we don't care that the sub-diagonals are
    # garbage.
    zero!(M, i, j) = nothing
    rank = bareiss!(M, (_swapcols!, _swaprows!, bareiss_update!, zero!);
                    find_pivot = find_pivot_col)
    return parity * M[end, end]
end

@testset "bareiss tests" begin
# copy gives a dense matrix
@testset "bareiss tests: $T" for T in (copy, sparse)
    # matrix determinent pairs
    for (M, d) in ((BigInt[9 1 8 0; 0 0 8 7; 7 6 8 3; 2 9 7 7], -1),
                   (BigInt[1 big(2)^65+1; 3 4], 4 - 3 * (big(2)^65 + 1)))
        # test that the determinent was correctly computed
        @test det_bareiss!(T(M)) == d
    end
end end

ext: jl


text:
# # Creating a `System` representing the entire U.S.

# **Originally Contributed by**: Clayton Barrows

# ## Introduction

# This example demonstrates how to assemble a `System` representing the entire U.S. using
# [PowerSystems.jl](https://github.com/NREL-SIIP/powersystems.jl) and the data assembled by
# [Xu, et. al.](https://arxiv.org/abs/2002.06155). We'll use the same tabular data parsing
# capability [demonstrated on the RTS-GMLC dataset](../../notebook/2_PowerSystems_examples/parse_tabulardata.ipynb).

# ### Environment
# This notebook depends on the SIIPExamples.jl environment which is loaded by default

#nb using Pkg
#nb Pkg.status()

# ### Dependencies
using SIIPExamples
using PowerSystems
using TimeSeries
using Dates
using TimeZones
using DataFrames
using CSV

# ### Fetch Data
# PowerSystems.jl links to some test data that is suitable for this example.
# Let's download the test data
@info "downloading data..."
datadir = joinpath(dirname(dirname(pathof(SIIPExamples))), "US-System")
siip_data = joinpath(datadir, "SIIP")
if !isdir(datadir)
    mkdir(datadir)
    mkdir(siip_data)
    tempfilename =
        download("https://zenodo.org/record/3753177/files/USATestSystem.zip?download=1")
    SIIPExamples.unzip(SIIPExamples.os, tempfilename, datadir)
end

config_dir = joinpath(
    dirname(dirname(pathof(SIIPExamples))),
    "script",
    "2_PowerSystems_examples",
    "US_config",
)

# ### Data Formatting
# This is a big dataset. Typically one would only want to include one of the interconnects
# available. Lets use Texas to start. You can set `interconnect = nothing` if you want everything.
interconnect = "Texas"
timezone = FixedTimeZone("UTC-6")
initial_time = ZonedDateTime(DateTime("2016-01-01T00:00:00"), timezone)

# There are a few minor incompatibilities between the data and the supported tabular data
# format. We can resolve those here.
#
# First, PowerSystems.jl only supports parsing piecewise linear generator costs from tabular
# data. So, we can sample the quadratic polynomial cost curves and provide PWL points.
@info "formatting data ..."
!isnothing(interconnect) && @info "filtering data to include $interconnect ..."
gen = DataFrame(CSV.File(joinpath(datadir, "plant.csv")))
filter!(row -> row[:interconnect] == interconnect, gen)
gencost = DataFrame(CSV.File(joinpath(datadir, "gencost.csv")))
gen = innerjoin(gen, gencost, on = :plant_id, makeunique = true, validate = (false, false))

function make_pwl(gen::DataFrame, traunches = 2)
    output_pct_cols = ["output_point_" * string(i) for i in 0:traunches]
    hr_cols = ["heat_rate_incr_" * string(i) for i in 1:traunches]
    pushfirst!(hr_cols, "heat_rate_avg_0")
    pwl = DataFrame(repeat([Float64], 6), Symbol.(vcat(output_pct_cols, hr_cols)))
    for row in eachrow(gen)
        traunch_len = (1.0 - row.Pmin / row.Pmax) / traunches
        pct = [row.Pmin / row.Pmax + i * traunch_len for i in 0:traunches]
        #c(pct) = pct * row.Pmax * (row.GenIOB + row.GenIOC^2 + row.GenIOD^3)
        c(pct) = pct * row.Pmax * (row.c1 + row.c2^2) + row.c0 #this formats the "c" columns to hack the heat rate parser in PSY
        hr = [c(pct[1])]
        [push!(hr, c(pct[i + 1]) - hr[i]) for i in 1:traunches]
        push!(pwl, vcat(pct, hr))
    end
    return hcat(gen, pwl)
end

gen = make_pwl(gen);

gen[!, "fuel_price"] .= 1000.0;  #this formats the "c" columns to hack the heat rate parser in PSY

# There are some incomplete aspects of this dataset. Here, I've assigned some approximate
# minimum up/down times, and some minor adjustments to categories. There are better
# ways to do this, but this works for this script...
gen[:, :unit_type] .= "OT"
gen[:, :min_up_time] .= 0.0
gen[:, :min_down_time] .= 0.0
gen[:, :ramp_30] .= gen[:, :ramp_30] ./ 30.0 # we need ramp rates in MW/min
[
    gen[gen.type .== "wind", col] .= ["Wind", 0.0, 0.0][ix]
    for (ix, col) in enumerate([:unit_type, :min_up_time, :min_down_time])
]
[
    gen[gen.type .== "solar", col] .= ["PV", 0.0, 0.0][ix]
    for (ix, col) in enumerate([:unit_type, :min_up_time, :min_down_time])
]
[
    gen[gen.type .== "hydro", col] .= ["HY", 0.0, 0.0][ix]
    for (ix, col) in enumerate([:unit_type, :min_up_time, :min_down_time])
]
[
    gen[gen.type .== "ng", col] .= [4.5, 8][ix]
    for (ix, col) in enumerate([:min_up_time, :min_down_time])
]
[
    gen[gen.type .== "coal", col] .= [24, 48][ix]
    for (ix, col) in enumerate([:min_up_time, :min_down_time])
]
[
    gen[gen.type .== "nuclear", col] .= [72, 72][ix]
    for (ix, col) in enumerate([:min_up_time, :min_down_time])
]

# At the moment, PowerSimulations can't do unit commitment with generators that have Pmin = 0.0
idx_zero_pmin = [
    g.type in ["ng", "coal", "hydro", "nuclear"] && g.Pmin <= 0
    for g in eachrow(gen[:, [:type, :Pmin]])
]
gen[idx_zero_pmin, :Pmin] = gen[idx_zero_pmin, :Pmax] .* 0.05

gen[:, :name] = "gen" .* string.(gen.plant_id)
CSV.write(joinpath(siip_data, "gen.csv"), gen)

# Let's also merge the zone.csv with the bus.csv and identify bus types
bus = DataFrame(CSV.File(joinpath(datadir, "bus.csv")))
!isnothing(interconnect) && filter!(row -> row[:interconnect] == interconnect, bus)
zone = DataFrame(CSV.File(joinpath(datadir, "zone.csv")))
bus = leftjoin(bus, zone, on = :zone_id)
int2bustype(b) =
    replace(split(string(PowerSystems.BusTypes.BusType(b)), ".")[end], "]" => "")
bus.bustype = int2bustype.(bus.type)
bus.name = "bus" .* string.(bus.bus_id)
CSV.write(joinpath(siip_data, "bus.csv"), bus)

# We need branch names as strings
branch = DataFrame(CSV.File(joinpath(datadir, "branch.csv")))
branch = join(
    branch,
    DataFrames.rename!(bus[:, [:bus_id, :baseKV]], [:from_bus_id, :from_baseKV]),
    on = :from_bus_id,
    kind = :left,
)
branch = join(
    branch,
    DataFrames.rename!(bus[:, [:bus_id, :baseKV]], [:to_bus_id, :to_baseKV]),
    on = :to_bus_id,
    kind = :left,
)
!isnothing(interconnect) && filter!(row -> row[:interconnect] == interconnect, branch)
branch.name = "branch" .* string.(branch.branch_id)
branch.tr_ratio = branch.from_baseKV ./ branch.to_baseKV
CSV.write(joinpath(siip_data, "branch.csv"), branch)

# The PowerSystems parser expects the files to be named a certain way.
# And, we need a "control_mode" column in dc-line data
dcbranch = DataFrame(CSV.File(joinpath(datadir, "dcline.csv")))
!isnothing(interconnect) && filter!(row -> row[:from_bus_id] in bus.bus_id, dcbranch)
!isnothing(interconnect) && filter!(row -> row[:to_bus_id] in bus.bus_id, dcbranch)
dcbranch.name = "dcbranch" .* string.(dcbranch.dcline_id)
dcbranch[:, :control_mode] .= "Power"
CSV.write(joinpath(siip_data, "dc_branch.csv"), dcbranch)

# ### We need to create a reference for where to get timeseries data for each component.
timeseries = []
ts_csv = ["wind", "solar", "hydro", "demand"]
plant_ids = Symbol.(string.(gen.plant_id))
for f in ts_csv
    @info "formatting $f.csv ..."
    csvpath = joinpath(siip_data, f * ".csv")
    csv = DataFrame(CSV.File(joinpath(datadir, f * ".csv")))
    (category, name_prefix, label) =
        f == "demand" ? ("Area", "", "max_active_power") :
        ("Generator", "gen", "max_active_power")
    if !(:DateTime in names(csv))
        DataFrames.rename!(
            csv,
            (names(csv)[occursin.("UTC", String.(names(csv)))][1] => :DateTime),
        )
        #The timeseries data is in UTC, this converts it to a fixed UTC offset
        csv.DateTime =
            ZonedDateTime.(
                DateTime.(csv.DateTime, "yyyy-mm-dd HH:MM:SS"),
                timezone,
                from_utc = true,
            )
        delete!(csv, csv.DateTime .< initial_time)
        csv.DateTime = Dates.format.(csv.DateTime, "yyyy-mm-ddTHH:MM:SS")
    end
    device_names = f == "demand" ? unique(bus.zone_name) : gen.name
    for id in names(csv)
        colname = id
        if f == "demand"
            if Symbol(id) in Symbol.(zone.zone_id)
                colname = Symbol(zone[Symbol.(zone.zone_id) .== Symbol(id), :zone_name][1])
                DataFrames.rename!(csv, (id => colname))
            end
            sf = sum(bus[string.(bus.zone_id) .== id, :Pd])
        else
            if Symbol(id) in plant_ids
                colname = Symbol(gen[Symbol.(gen.plant_id) .== Symbol(id), :name][1])
                DataFrames.rename!(csv, (id => colname))
            end
            sf = maximum(csv[:, colname]) == 0.0 ? 1.0 : "Max"
        end
        if String(colname) in device_names
            push!(
                timeseries,
                Dict(
                    "simulation" => "DA",
                    "category" => category,
                    "module" => "InfrastructureSystems",
                    "type" => "SingleTimeSeries",
                    "component_name" => String(colname),
                    "name" => label,
                    "resolution" => 300,
                    "scaling_factor_multiplier" => "get_max_active_power",
                    "scaling_factor_multiplier_module" => "PowerSystems",
                    "normalization_factor" => sf,
                    "data_file" => csvpath,
                ),
            )
        end
    end
    CSV.write(csvpath, csv)
end

timeseries_pointers = joinpath(siip_data, "timeseries_pointers.json")
open(timeseries_pointers, "w") do io
    PowerSystems.InfrastructureSystems.JSON3.write(io, timeseries)
end

# ### The tabular data format relies on a folder containing `*.csv` files and `.yaml` files
# describing the column names of each file in PowerSystems terms, and the PowerSystems
# data type that should be created for each generator type. The respective "us_decriptors.yaml"
# and "US_generator_mapping.yaml" files have already been tailored to this dataset.
@info "parsing csv files..."
rawsys = PowerSystems.PowerSystemTableData(
    siip_data,
    100.0,
    joinpath(config_dir, "us_descriptors.yaml"),
    generator_mapping_file = joinpath(config_dir, "US_generator_mapping.yaml"),
)

# ### Create a `System`
# Next, we'll create a `System` from the `rawsys` data. Since a `System` is predicated on a
# time series resolution and the `rawsys` data includes both 5-minute and 1-hour resolution
# time series, we also need to specify which time series we want to include in the `System`.
# The `time_series_resolution` kwarg filters to only include time series with a matching resolution.

@info "creating System"
sys = System(rawsys; config_path = joinpath(config_dir, "us_system_validation.json"));
sys

# This all took reasonably long, so we can save our `System` using the serialization
# capability included with PowerSystems.jl:
to_json(sys, joinpath(siip_data, "sys.json"), force = true)


text:
using Documenter
using Peridynamics

DocMeta.setdocmeta!(Peridynamics.Discretization, :DocTestSetup, :(using Peridynamics.Discretization); recursive=true)

makedocs(
    sitename = "Peridynamics",
    format = Documenter.HTML(),
    modules = [Peridynamics]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/johntfoster/Peridynamics.jl.git"
)


text:
using Test, JSXGraph, JSExpr
const J = JSXGraph

include("test_utils.jl")

include("utils.jl")
include("jsfun.jl")
include("board.jl")

include("objects/controllers.jl")
include("objects/curves.jl")

# --- integration tests

include("examples/curves.jl")
include("examples/point.jl")

# include("objects/controllers.jl")


text:
function updateP!(P::Array{Any,1}, A::Array{Any,1}, DataInvMat::Array{Any,1}, DataMat::Array{Any,1}, τ::Float64)
# Update P by Eq. (10)

    for (i, Aᵢ) in enumerate(A)
        @inbounds tempMat::Matrix{Float64} = Aᵢ * DataMat[i]'
        scale!(tempMat, τ)
        @inbounds P[i] = tempMat * DataInvMat[i]
    end
end


text:
using NBody

function main()
  nb = NBody.read_nbody_json()
  NBody.set_ids(nb)
  NBody.write_snapshot(nb)
end

main()


text:
"""
$(TYPEDEF)

Statistics from the differential equation solver about the solution process.

## Fields

- nf: Number of function evaluations. If the differential equation is a split function,
  such as a `SplitFunction` for implicit-explicit (IMEX) integration, then `nf` is the
  number of function evaluations for the first function (the implicit function)
- nf2: If the differential equation is a split function, such as a `SplitFunction`
  for implicit-explicit (IMEX) integration, then `nf2` is the number of function
  evaluations for the second function, i.e. the function treated explicitly. Otherwise
  it is zero.
- nw: The number of W=I-gamma*J (or W=I/gamma-J) matrices constructed during the solving
  process.
- nsolve: The number of linear solves `W\b` required for the integration.
- njacs: Number of Jacobians calculated during the integration.
- nnonliniter: Total number of iterations for the nonlinear solvers.
- nnonlinconvfail: Number of nonlinear solver convergence failures.
- ncondition: Number of calls to the condition function for callbacks.
- naccept: Number of accepted steps.
- nreject: Number of rejected steps.
- maxeig: Maximum eigenvalue over the solution. This is only computed if the
  method is an auto-switching algorithm.
"""
mutable struct DEStats
  nf::Int
  nf2::Int
  nw::Int
  nsolve::Int
  njacs::Int
  nnonliniter::Int
  nnonlinconvfail::Int
  ncondition::Int
  naccept::Int
  nreject::Int
  maxeig::Float64
end

DEStats(x::Int = -1) = DEStats(x, x, x, x, x, x, x, x, x, x, 0.0)

function Base.show(io::IO, s::DEStats)
  println(io, summary(s))
  @printf io "%-50s %-d\n" "Number of function 1 evaluations:" s.nf
  @printf io "%-50s %-d\n" "Number of function 2 evaluations:" s.nf2
  @printf io "%-50s %-d\n" "Number of W matrix evaluations:" s.nw
  @printf io "%-50s %-d\n" "Number of linear solves:" s.nsolve
  @printf io "%-50s %-d\n" "Number of Jacobians created:" s.njacs
  @printf io "%-50s %-d\n" "Number of nonlinear solver iterations:" s.nnonliniter
  @printf io "%-50s %-d\n" "Number of nonlinear solver convergence failures:" s.nnonlinconvfail
  @printf io "%-50s %-d\n" "Number of rootfind condition calls:" s.ncondition
  @printf io "%-50s %-d\n" "Number of accepted steps:"  s.naccept
  @printf io "%-50s %-d" "Number of rejected steps:" s.nreject
  iszero(s.maxeig) || @printf io "\n%-50s %-d" "Maximum eigenvalue recorded:" s.maxeig
end


text:
module TestServices

import SimpleTasks.Services.Cache
import SimpleTasks.Services.FileSystemCache

export TEST_BASE_DIRECTORY, make_valid_file_system_cache
const TEST_BASE_DIRECTORY = "/var/tmp"
function make_valid_file_system_cache()
    return FileSystemCache.FileSystemCacheService(TEST_BASE_DIRECTORY)
end

end # module TestServices

max_forks_repo_name: seung-lab/SimpleTasks


text:
import PyPlot
import PlotlyBase, Plotly
using Dates
using CSV
using XLSX
using DataFrames
using Printf
using Gumbo
using HTTP

# Number of column containing country name
c_country=2

# Number of start column of time series
c_timeseries_start=5

#dsname=""
dsname=""

# Source file for data (for download)

# New version
jhu_infected="https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv"
jhu_deaths="https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv"
jhu_recovered="https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_recovered_global.csv"
rki_nowcast="https://www.rki.de/DE/Content/InfAZ/N/Neuartiges_Coronavirus/Projekte_RKI/Nowcasting_Zahlen.xlsx?__blob=publicationFile"

# Data URL (for data source in plot)
dataurl="https://github.com/CSSEGISandData/COVID-19/"

# Download and read data for confirmed infection cases
# This creates a dataframe
function load_jhu_infected(;download=true)
    if download
        Base.download(jhu_infected,"jhu_infected.csv")
    end
    CSV.read("jhu_infected.csv",DataFrame)
end

function load_jhu_recovered(;download=true)
    if download
        Base.download(jhu_recovered,"jhu_recovered.csv")
    end
    CSV.read("jhu_recovered.csv",DataFrame)
end
function load_jhu_deaths(;download=true)
    if download
        Base.download(jhu_deaths,"jhu_deaths.csv")
    end
    CSV.read("jhu_deaths.csv",DataFrame)
end

############################################################################
# web scraping for rki data

function scrape_rki_from_wikipedia()
    function find_table(element,text)
	if isa(element,HTMLText)
     	    return nothing
        end
	if length(element.children)>0 && element[1] isa Gumbo.HTMLElement{:caption} && occursin(text,String(element[1][1].text))
	    return element
	else
	    for i=1:length(element.children)
		result=find_table(element.children[i],text)
                if result!=nothing
		    return result
		end
	    end
	end
    end
    header(t)=push!([t[2][1][i][1][1].text for i=3:18],"DE")
    nrows(t)=length(t[2].children)-2

    function row(t,i)
	tr=t[2][i+1]
	start=2
	if haskey(tr[1].attributes,"rowspan")
	    start=3
	end
	[ isa(tr[i][1],HTMLText) ? replace(tr[i][1].text, r"\n|\t|\." => s"") : "" for i=start:length(tr.children)-1]
    end

    num(s)= try parse(Int64,s) catch e 0 end

    function dataframe(t,date_start)
        df=DataFrame(Datum=[date_start+Day(i) for i=1:nrows(t)])
        for icol=1:length(header(t))
	    df[!,Symbol(header(t)[icol])]=[num(row(t,irow)[icol]) for irow=1:nrows(t)]
        end
        df
    end

    rawpage=HTTP.get("https://de.wikipedia.org/wiki/COVID-19-Pandemie_in_Deutschland/Statistik");
    parsed_page=parsehtml(String(rawpage.body));
    t=find_table(parsed_page.root,"Infektionsfälle")

    df=dataframe(find_table(parsed_page.root,"Infektionsfälle"),Date(2020,2,23))
    println(df[!,:Datum][end])
    CSV.write("rki.csv",df)
    df=dataframe(find_table(parsed_page.root,"Todesfälle"),Date(2020,3,7))
    println(df[!,:Datum][end])
    CSV.write("rki_dead.csv",df)
end

############################################################
# Download and read data for confirmed infection cases
# This creates a dataframe
function load_rki()
    csv=CSV.read("rki.csv",DataFrame)
end

function load_nowcast(;download=false)
    if download
        Base.download(rki_nowcast,"Nowcasting_Zahlen.xlsx")
    end
    xf=XLSX.readxlsx("Nowcasting_Zahlen.xlsx")
    sheet=xf["Nowcast_R"]
    nowcast=[]
    i=1
    for row in XLSX.eachrow(sheet)
        if i>1
            push!(nowcast,row[2])
        end
        i=i+1
    end
    nowcast
end

function alldata_world(;download=false)
    (dead=load_jhu_deaths(download=download),
     infected=load_jhu_infected(download=download),
     recovered=load_jhu_recovered(download=download),
     #     https://data.worldbank.org/indicator/SP.POP.TOTL
     popdata=CSV.read("population.csv",DataFrame))
end

function alldata_blaender(;download=false)
    (
        dead=CSV.read("rki_dead.csv",DataFrame),
        infected=CSV.read("rki.csv",DataFrame),
        recovered=nothing,
        #     https://www.statistik-bw.de/VGRdL/tbls/tab.jsp?rev=RV2014&tbl=tab20&lang=de-DE
        popdata=CSV.read("einwohnerzahl-bundeslaender.csv",DataFrame)
    )
end

# Create a new dataframe wich contains only data for given countries
# For US and China, data are given state/province/county-wise
# For former colonial powers with overseas territories etc. there are multiple rows as well
function select_countries_rows(df_jhu,countries)
    countries_mask=BitArray(undef,length(countries),size(df_jhu,1))
    for i=1:length(countries)
        @show countries[i]
        countries_mask[i,:].=df_jhu[:,c_country].==countries[i]
    end
    mask=vec(reduce(|,countries_mask,dims=1))
    df_jhu[mask,:]
end

# Create timeseries for list of countries
function create_countries_timeseries(df_jhu,countries)
    crows=select_countries_rows(df_jhu,countries)
    data=sum(Matrix(crows[:,c_timeseries_start:end]),dims=1)'
end

# Calculate growth factor from growth rate (in %)
growth_factor(growth_rate)=growth_rate/100.0+1

# Calculate doubling time from growth factor
doubling_time(gfactor)= log(2.0)/log(gfactor)

growth_factor_from_doubling_time(dtime)=exp(log(2.0)/dtime)

# Calculate growth rate (in %) from growth factor
growth_rate(gfactor)=(gfactor-1)*100

# Plot data for country
function plotcountry(df,
                     countries,  # Array of countries
                     kind;       # Kind of plot
                     label="", # label
                     lw=2,       # plot line width
                     lt="-",     # plot line type
                     averaging_period=15,   # averaging period: averaging_period days
                     Nstart=100  # starting data for shifting curves
                     )
    rki=false
    if size(df,2) <20
        rki=true
    end

    # If label name is not given, take the first name from the array
    if label===""
        label=countries[1]
    end

    # Add 1 to data for allowing to  logarithm or division
    if kind=="log"
        logdiv_regularization=1.0e-10
    elseif  kind=="growthrate"
        logdiv_regularization=1.0
    elseif  kind=="doublingtime"
        logdiv_regularization=1.0
    else
        logdiv_regularization=0.0
    end

    # Create shifted time series

    # Shift timeseries to the left cutting of the days until Nstart infections occur

    if rki
        basedata=Array{Float64}(df[Symbol(label)])
    else
        basedata=create_countries_timeseries(df,countries).+logdiv_regularization
    end
    shift=1
    while basedata[shift]<Nstart
        shift=shift+1
    end
    data=basedata[shift:end]
    days=collect(0:length(data)-1)
    # print for control purposes (to verify  with the numbers given in the map app)
    println("$(label), $(maximum(data))")

    # Perform plots
    if kind=="abs"
        plot(days,data,label="$(label)",lt,linewidth=lw,markersize=6)
    end

    if kind=="log"
        semilogy(days,data,label="$(label)",lt,linewidth=lw,markersize=6)
    end

    if kind=="growthrate" || kind=="doublingtime"
        # Calculate daily growth factors
        gfactors=basedata[2:end]./basedata[1:end-1]
        if averaging_period > length(gfactors)-1
            println("skip: $(label)")
            return
        end

        # Calculate the average over averaging_period
        averaged_gfactors=ones(length(gfactors)-averaging_period)
        j=1
        for i=averaging_period+1:length(gfactors)
            for d=1:averaging_period
                averaged_gfactors[j]*=gfactors[i-d+1]
            end
            averaged_gfactors[j]=averaged_gfactors[j]^(1.0/averaging_period)
            j=j+1
        end

        # Calculate growth rates
        grates=growth_rate.(averaged_gfactors)
        dtimes=doubling_time.(averaged_gfactors)

        # Adjust starting day due to change of reporting on US data
        day0=1
        if rki
            day0=33
        end

        if kind=="growthrate"
            plot(day0:day0+length(grates)-1,grates[1:end],label="$(label)",lt,linewidth=lw,markersize=6)
        end
        if kind=="doublingtime"
            plot(day0:day0+length(grates)-1,dtimes[1:end],label="$(label)",lt,linewidth=lw,markersize=6)
        end
    end
end

# Plot data for all countries
function plotcountries(df_jhu,df_rki,
                       kind;       # Kind of plot
                       averaging_period=15,   # averaging period: averaging_period days
                       Nstart=500  # starting data for shifting curves
                       )
    # List of countries belonging to Europa
    # (omitting those with less than 100 cases)
    Europe=[
        "Austria",
        "Sweden",
        "Bulgaria",
        "Croatia",
        "Cyprus",
        "Czechia",
        "Denmark",
        "Estonia",
        "Finland",
        "France",
        "Germany",
        "Greece",
        "Hungary",
        "Ireland",
        "Italy",
        "Latvia",
        "Lithuania",
        "Luxembourg",
        "Malta",
        "Netherlands",
        "Poland",
        "Portugal",
        "Romania",
        "Turkey",
        "Slovenia",
        "Spain",
        "Sweden",
        "Norway",
        "Switzerland",
        "United Kingdom",
        "Serbia"
    ]

    plotcountry(df_jhu,["Italy"],kind,Nstart=Nstart, lt="o-", averaging_period=averaging_period)
    plotcountry(df_jhu,["France"],kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_jhu,["Spain"],kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_jhu,["Iran"],kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_jhu,["Korea, South"],kind,Nstart=Nstart,lt="o-", averaging_period=averaging_period)
    plotcountry(df_jhu,["China"],kind,Nstart=Nstart,lt="o-", averaging_period=averaging_period)
    plotcountry(df_jhu,["Switzerland"],kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_jhu,["Netherlands"],kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_jhu,["Austria"],kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_jhu,["Sweden"],kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_jhu,["Turkey"],kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_jhu,["Canada"],lt="y-",kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_jhu,["Russia"],lt="mo-",kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_jhu,["Brazil"],lt="g-",kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_jhu,Europe,label="Europe",kind,lt="b-o",Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_jhu,["Germany"],lw=3,lt="r-o",kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_jhu,["US"],kind,lt="k-",Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_jhu,["United Kingdom"],kind,lt="k-o", Nstart=Nstart, averaging_period=averaging_period)
end

# Plot data for all countries
function plotbundeslaender(df_jhu,df_rki,
                           kind;       # Kind of plot
                           averaging_period=15,   # averaging period: averaging_period days
                           Nstart=100  # starting data for shifting curves
                           )

    plotcountry(df_jhu,["Germany"], label="JHU",lw=3,lt="r-",kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_rki,["BW"],lt="y-o",kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_rki,["BY"],lt="k-o",kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_rki,["BB"],kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_rki,["HB"],kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_rki,["HH"],kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_rki,["HE"],kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_rki,["MV"],kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_rki,["NI"],kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_rki,["NW"],lt="b-o",kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_rki,["RP"],kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_rki,["SL"],kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_rki,["SN"],kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_rki,["ST"],kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_rki,["SH"],lt="y-",kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_rki,["TH"],lt="k-",kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_rki,["BE"],lw=3,lt="g-o",kind,Nstart=Nstart, averaging_period=averaging_period)
    plotcountry(df_rki,["Gesamt"],lw=3,lt="r-o",kind,Nstart=Nstart, averaging_period=averaging_period)
end

# Create the plots
function create_world_old(;averaging_periods=[7,15],Nstart=500,dtime=true)

    df_jhu=load_jhu_infected()
    df_rki=load_rki()
    N0=10000
    trailer="\nData source: JHU $(Dates.today())\nData processing: https://github.com/j-fu/coronaplot  License: CC-BY 2.0"

    # Plot absolute values (linear Y scale) to show exponential behavior
    fig = PyPlot.figure(1)
    fig = PyPlot.gcf()
    fig.set_size_inches(10,5)
    clf()
    title("Corona Virus Development in countries with more than $(N0) infections$(trailer)")
    plotcountries(df_jhu,df_rki,"abs",Nstart=Nstart)
#    PyPlot.ylim(1,350_000)
    PyPlot.ylim(1,200_000)
    PyPlot.xlim(0,60)
    PyPlot.grid()
    PyPlot.xlabel("Days since occurence of at least $(Nstart) infections")
    PyPlot.ylabel("Number of infections (linear scale)")
    PyPlot.legend(loc="upper left")
    PyPlot.savefig("../docs/infected-exp.png")
    PyPlot.savefig("../infected-exp.png")

    # Log plot
    fig = PyPlot.figure(2)
    fig = PyPlot.gcf()
    fig.set_size_inches(10,5)
    clf()
    title("Corona Virus Development in countries with more than  $(N0) infections$(trailer)")
    plotcountries(df_jhu,df_rki,"log",Nstart=Nstart)
    PyPlot.xlim(0,60)
    PyPlot.grid()
    PyPlot.xlabel("Days since occurence of at least $(Nstart) infections")
    PyPlot.ylabel("Number of infections (logarithmic scale)")
    PyPlot.legend(loc="upper left")
    PyPlot.savefig("../docs/infected.png")
    PyPlot.savefig("../infected.png")

    ifig=3
    for averaging_period in averaging_periods
        # Plot evolution of growth rate average
        fig = PyPlot.figure(ifig)
        fig = PyPlot.gcf()
        fig.set_size_inches(10,5)
        clf()
        if dtime
            title("$(averaging_period) day average of doubling time in countries with > $(N0) infections$(trailer)")
            plotcountries(df_jhu,df_rki,"doublingtime",averaging_period=averaging_period)
        else
            title("$(averaging_period) day average of daily growth rate of COVID-19 infections in countries with > $(N0) infections$(trailer)")
            plotcountries(df_jhu,df_rki,"growthrate",averaging_period=averaging_period)
        end
        PyPlot.xlim(10,90)
        if dtime
            PyPlot.ylim(0,80)
        else
            PyPlot.ylim(0,100)
        end

        PyPlot.grid()
        month="February"
        day=averaging_period-10
        if day<=0
            month="January"
            day=averaging_period+22
            PyPlot.xlim(16,95)
        end
        PyPlot.xlabel("Days since $(month) $(day), 2020")

        dtime_label="Doubling time/days"
        grate_label="$(averaging_period) day average of daily growth/%"
        if dtime
            PyPlot.ylabel(dtime_label)
        else
            PyPlot.ylabel(grate_label)
        end
        PyPlot.legend(loc="upper left")

        # Add second y axis
        # see  https://stackoverflow.com/a/10517481/8922290
        if dtime
            ax1 = PyPlot.gca()
            ax2 = ax1.twinx()
            ax2.set_ylim(ax1.get_ylim())
            dtimes=collect(0:5:60)
            ax2.set_yticks(dtimes)
            grates=growth_rate.(growth_factor_from_doubling_time.(dtimes))
            ax2.set_yticklabels([ @sprintf("%.2f",grates[i]) for i=1:length(dtimes)])
            ax2.set_ylabel(grate_label)
        else
            ax1 = PyPlot.gca()
            ax2 = ax1.twinx()
            ax2.set_ylim(ax1.get_ylim())
            growth_rates= collect(0:5:100)
            ax2.set_yticks(growth_rates)
            dtimes=doubling_time.(growth_factor.(growth_rates))
            ax2.set_yticklabels([ @sprintf("%.2f",dtimes[i]) for i=1:length(dtimes)])
            ax2.set_ylabel(dtime_label)
        end

        if averaging_period==15
            PyPlot.savefig("../docs/infected-growthrate.png")
            PyPlot.savefig("../infected-growthrate.png")
        end
        if averaging_period==7
            PyPlot.savefig("../docs/infected-growthrate-weeklyavg.png")
            PyPlot.savefig("../infected-growthrate-weeklyavg.png")
        end
        ifig=ifig+1
    end

end

# Create the plots
function create_blaender_old(;averaging_periods=[7,15],Nstart=100,dtime=true)

    df_jhu=load_jhu_infected()
    df_rki=load_rki()

    trailer="\nDatenquelle: RKI+ JHU $(Dates.today())\nDatenverarbeitung: https://github.com/j-fu/coronaplot  Lizenz: CC-BY 2.0"

    # Plot absolute values (linear Y scale) to show exponential behavior
    fig = PyPlot.figure(1)
    fig = PyPlot.gcf()
    fig.set_size_inches(10,5)
    clf()
    title("Entwicklung der COVID19-Infektionszahlen in Deutschland$(trailer)")
    plotbundeslaender(df_jhu,df_rki,"abs",Nstart=Nstart)
#    PyPlot.ylim(1,350_000)
    PyPlot.ylim(1,40_000)
    PyPlot.xlim(0,55)
    PyPlot.grid()
    PyPlot.xlabel("Tage seit dem Auftreten von $(Nstart) Infektionen")
    PyPlot.ylabel("Anzahl der Infektionen (lineare Skala)")
    PyPlot.legend(loc="upper left")
    PyPlot.savefig("../docs/de-infected-exp.png")

    # Log plot
    fig = PyPlot.figure(2)
    fig = PyPlot.gcf()
    fig.set_size_inches(10,5)
    clf()
    title("Entwicklung der COVID19-Infektionszahlen in Deutschland$(trailer)")
    plotbundeslaender(df_jhu,df_rki,"log",Nstart=Nstart)
    PyPlot.xlim(0,55)
    PyPlot.grid()
    PyPlot.xlabel("Tage seit dem Auftreten von $(Nstart) Infektionen")
    PyPlot.ylabel("Anzahl der Infektionen (logarithmische Skala)")
    PyPlot.legend(loc="lower right")
    PyPlot.savefig("../docs/de-infected.png")

    ifig=3
    for averaging_period in averaging_periods
        # Plot evolution of growth rate average
        fig = PyPlot.figure(ifig)
        fig = PyPlot.gcf()
        fig.set_size_inches(10,5)
        clf()
        if dtime
            title("$(averaging_period)-Tage-Mittel der Vertopplungszeiten$(trailer)")
            plotbundeslaender(df_jhu,df_rki,"doublingtime",averaging_period=averaging_period)
        else
            title("$(averaging_period)-Tage-Mittel der täglichen Wachstumsraten der COVID19-Infektionen in Deutschland$(trailer)")
            plotbundeslaender(df_jhu,df_rki,"growthrate",averaging_period=averaging_period)
        end
        if dtime
            PyPlot.ylim(0,100)
        else
            PyPlot.ylim(0,60)
        end
        PyPlot.xlim(50,90)
        PyPlot.grid()
        month="Februar"
        day=averaging_period-10
        if day<=0
            month="Januar"
            day=averaging_period+22
            PyPlot.xlim(50,95)
        end
        PyPlot.xlabel("Tage seit $(day). $(month)  2020")
        dtime_label="Vedopplungszeinten/d"
        grate_label="$(averaging_period)-tägiges Mittel der täglichen Wachstumsraten/%"

        if dtime
            PyPlot.ylabel(dtime_label)
        else
            PyPlot.ylabel(grate_label)
        end

        PyPlot.legend(loc="upper left")

        # Add second y axis with doubling time
        # see  https://stackoverflow.com/a/10517481/8922290
        if dtime
            ax1 = PyPlot.gca()
            ax2 = ax1.twinx()
            ax2.set_ylim(ax1.get_ylim())
            dtimes=collect(0:5:60)
            ax2.set_yticks(dtimes)
            grates=growth_rate.(growth_factor_from_doubling_time.(dtimes))
            ax2.set_yticklabels([ @sprintf("%.2f",grates[i]) for i=1:length(dtimes)])
            ax2.set_ylabel(grate_label)
        else
            ax1 = PyPlot.gca()
            ax2 = ax1.twinx()
            ax2.set_ylim(ax1.get_ylim())
            growth_rates= collect(0:2.5:60)
            ax2.set_yticks(growth_rates)
            dtimes=doubling_time.(growth_factor.(growth_rates))
            ax2.set_yticklabels([ @sprintf("%.2f",dtimes[i]) for i=1:length(dtimes)])
            ax2.set_ylabel("dtime_label")
        end
        if averaging_period==15
            PyPlot.savefig("../docs/de-infected-growthrate.png")
        end
        if averaging_period==7
            PyPlot.savefig("../docs/de-infected-growthrate-weeklyavg.png")
        end
        ifig=ifig+1
    end

end

function create_plots_old()
    create_world_old()
    create_blaender_old()
end
# publish on github
function publish(;msg="data update")
    run(`git commit -a -m $(msg)`)
    run(`git push`)
end

#####################################################################
#####################################################################
#####################################################################
#####################################################################
#####################################################################
# New plots for active cases and r0

"""
Moving average of timeseries over window
"""
function mvavg(ts,window)
    tsavg=[]
    for i=1:length(ts)-window
        push!(tsavg,sum(ts[i:i+window])/window)
    end
    tsavg
end

"""
Estimate of reproduction rate from time series.
"""
function r0(ts,window)
    r0=[]
    for i=window:length(ts)-window
        # Newly infected during the current period
        new=ts[i+1:i+window]
        # Newly infected during previous period
        old=ts[i-window+1:i]
#        @show length(old), length(new)
        push!(r0,sum(new)/sum(old))
    end
    r0
end

"""
Adjust length of times series to len, padding with values at the
beginning
"""
function leftpad(ts,len;pad=0.0)
    ts_padded=[]
    for i=1:len-length(ts)
        push!(ts_padded,pad)
    end
    for i=1:length(ts)
        push!(ts_padded,ts[i])
    end
    ts_padded
end

"""
Calculate results for given country
"""
function country_results(data, country;
                         world=true,   # do we have country or bundesland ?
                         avg_window=7, # Window for moving average of time series of infected people
                         active_period=15, # Period during which we assume an infection is active
                         infection_period=5, # Time it takes for an infected person to infect the next (RKI uses 4)
                         population_base=100_000.0
                         )

    # Access data frames via tuple
    (dead,infected,recovered,popdata)=data

    # Filter population data for country
    population=filter((row)-> row[1]==country, popdata)[1,2]

    popfac=Float64(population_base)/population

    # Create time series from data
    if world
        ts_dead=create_countries_timeseries(dead,[country]).*popfac
        ts_infected=create_countries_timeseries(infected,[country]).*popfac
        ts_recovered=create_countries_timeseries(recovered,[country]).*popfac
        mvavg_dead=mvavg(ts_dead,avg_window)
    else
        ts_dead=Array{Float64}(dead[!,Symbol(country)]).*popfac
        ts_infected=Array{Float64}(infected[!,Symbol(country)]).*popfac
        mvavg_dead=mvavg(ts_dead,avg_window)
    end

    mvavg_infected=mvavg(ts_infected,avg_window)

    if world
        mvavg_recovered=mvavg(ts_recovered,avg_window)
        mvavg_active=mvavg_infected-mvavg_recovered-mvavg_dead
    else
        mvavg_recovered=nothing
        mvavg_active=nothing
    end

    # Estimate number of active cases. We just assume that reported
    # cases are active during the active period. If this number is chosen
    # large enough, we seem to get some crude, but conservative estimate.
    #
    # Just subtract the number from `active_period` ago from th ecurrent one
    # We base this on the moving average calculated before
    est_active=leftpad(mvavg_infected[active_period+1:end]-mvavg_infected[1:end-active_period],length(mvavg_infected))

    # Calculate time series of new cases from moving average of infected people
    ts_new=ts_infected[2:end]-ts_infected[1:end-1]
    weekly_new=[]
    for i=8:length(ts_new)
        push!(weekly_new,sum(ts_new[i-6:i]))
    end
    # This is based on the description given in
    # https://www.heise.de/newsticker/meldung/Corona-Pandemie-Die-Mathematik-hinter-den-Reproduktionszahlen-R-4712676.html
    # describing the RKI method for estimating R0.
    #
    # There are two main differences here: we possibly assume a longer infection period,
    # and we base the results on the moving average of the timeseries of infected instead of the "nowcast"
    ts_new=mvavg_infected[2:end]-mvavg_infected[1:end-1]
    est_r0=leftpad(r0(ts_new,infection_period),length(mvavg_infected))

    # Create a time series of dates to get the x axis right
    if world
        date_start=Date(2020,1,22)+Day(avg_window)
    else
        date_start=Date(2020,2,23)+Day(avg_window)
    end
    dates=[date_start+Day(i) for i=1:length(mvavg_infected)]

    (mvavg_dead=mvavg_dead,
     mvavg_infected=mvavg_infected,
     mvavg_recovered=mvavg_recovered,
     mvavg_active=mvavg_active,
     est_active=est_active,
     est_r0=est_r0,
     population=population,
     mvavg_new=leftpad(weekly_new,length(mvavg_infected)),
     dates=dates)
end

function nowcast_results(nowcast;
                         active_period=15,   # Period during which we assume an infection is active
                         infection_period=5, # Time it takes for an infected person to infect the next (RKI uses 4)
                         population_base=100_000)
    pop=82927922.0
    popfac=population_base/pop
    newly=nowcast*popfac
    date_start=Date(2020,3,2)
    dates=[date_start+Day(i) for i=1:length(nowcast)]
    active=[]
    for i=active_period+1:length(nowcast)
        push!(active,sum(newly[i-active_period:i]))
    end
    active=leftpad(active,length(nowcast))
    est_r0=leftpad(r0(newly,infection_period),length(nowcast))
    est_new=[]
    for i=8:length(nowcast)
        push!(est_new,sum(newly[i-7:i]))
    end
    est_new=leftpad(est_new,length(nowcast))

    (dates=dates,
     est_new=est_new,
     est_active=active,
     est_r0=est_r0)
end

function plot_active_r0(country;download=false,world=true,avg_window=7)

    if world
        data=alldata_world(download=download)
    else
        data=alldata_blaender()
    end

    population_base=100_000
    results=country_results(data,country,
                            world=world,
                            avg_window=avg_window, # Window for moving average of time series of infected people
                            active_period=15, # Period during which we assume an infection is active
                            infection_period=5, # Time it takes for an infected person to infect the next (RKI uses 4)
                            population_base=population_base)

    fig = PyPlot.figure(1)
    fig = PyPlot.gcf()
    fig.set_size_inches(10,5)
    clf()
    title(@sprintf("Country: %s, population: %.1f million",country,results.population/1.0e6))

    d0=1
    PyPlot.plot_date(results.dates[d0:end],results.mvavg_infected[d0:end],label="infected","b-")
    PyPlot.plot_date(results.dates[d0:end],results.mvavg_new[d0:end],label="new","y-")
    PyPlot.plot_date(results.dates[d0:end],results.est_active[d0:end],label="active (est)","rx")
    if world
        PyPlot.plot_date(results.dates[d0:end],results.mvavg_active[d0:end],label="active (jhu)","r-")
        PyPlot.plot_date(results.dates[d0:end],results.mvavg_recovered[d0:end],label="recovered","g-")
        PyPlot.plot_date(results.dates[d0:end],results.mvavg_dead[d0:end],label="dead","k-")
    end
    PyPlot.xlabel("Date")
    PyPlot.ylabel("Cases per $(population_base) inhabitants")
    PyPlot.legend(loc="upper left")
    PyPlot.grid()

    ax1 = PyPlot.gca()
    ax2 = ax1.twinx()
    ax2.set_ylim(0,4)
    ax2.plot_date(results.dates[d0:end],results.est_r0[d0:end],label="R0","m--")
    ax2.plot_date(results.dates[d0:end],[1.0 for i=d0:length(results.dates)],label="R0=1","m-",linewidth=2)
    ax2.set_ylabel("R0")
    ax2.legend(loc="upper right")

    PyPlot.show()
end

function countrylist()
   [
    ["Italy", "o-"],
    ["France", "o-"],
    ["Spain", "o-"],
    ["Iran", "-"],
    ["Korea, South","o-"],
    ["China","o-"],
    ["Israel", "-"],
    ["India", "-"],
    ["Czechia", "-"],
    ["Sweden", "o-"],
    ["Turkey", "-"],
    ["Canada","y-"],
    ["Russia","mo-"],
    ["Brazil","g-"],
    ["Germany","r-o"],
    ["US","k-o"],
    ["United Kingdom","k-"],
    ["Poland", "-"],
    ["Austria", "-"],
   ]
end

function blaenderlist()
[
["BW","y-o" ],
["BY","k-o" ],
["BB","-"   ],
["HB","-"   ],
["HH","-"   ],
["HE","-"   ],
["MV","-"   ],
["NI","-"   ],
["NW","b-o" ],
["RP","-"   ],
["SL","-"   ],
["SN","-"   ],
["ST","-"   ],
["SH","y-"  ],
["TH","k-"  ],
["DE","r-o" ],
["BE","g-o" ],
]
end

function plot_active_r0(;download=false, world=true, infection_period=5,avg_window=7)
    if world
        data=alldata_world(download=download)
        countries=countrylist()
        trailer="Data source: JHU $(Dates.today())\nData processing: https://github.com/j-fu/coronaplot  License: CC-BY 2.0"
        prefix="world"
    else
        if download
            scrape_rki_from_wikipedia()
        end
        data=alldata_blaender()
        countries=blaenderlist()
        trailer="Datenquelle: RKI via de.wikipedia.org/wiki/COVID-19-Pandemie_in_Deutschland/Statistik $(Dates.today())\nDatenverarbeitung: https://github.com/j-fu/coronaplot  Lizenz: CC-BY 2.0"
        prefix="de"
        nowcast=load_nowcast(download=false)
    end
    fig = PyPlot.figure(1)
    fig = PyPlot.gcf()
    PyPlot.clf()
    fig.set_size_inches(10,5)
    if world
        d0=40
    else
        d0=10
    end
    population_base=100_000

    if world
        PyPlot.title("Estimated number of infectious persons in selected countries\n$(trailer)")
    else
        PyPlot.title("Schätzung der Zahl der infektiösen Personen\n$(trailer)")
    end

    for country in countries
        results=country_results(data,country[1],
                                world=world,
                                avg_window=avg_window, # Window for moving average of time series of infected people
                                active_period=15, # Period during which we assume an infection is active
                                infection_period=infection_period, # Time it takes for an infected person to infect the next (RKI uses 4)
                                population_base=population_base)
        PyPlot.plot_date(results.dates[d0:end],results.est_active[d0:end],label=country[1],country[2])
    end
    if !world
        results=nowcast_results(nowcast,active_period=15,infection_period=infection_period,population_base=population_base)
        PyPlot.plot_date(results.dates[d0:end],results.est_active[d0:end],label="nowcast","r*-",linewidth=2)
    end

    if world
        PyPlot.xlabel("Date")
        PyPlot.ylabel("Cases per $(population_base) inhabitants")
    else
        PyPlot.xlabel("Datum")
        PyPlot.ylabel("Fälle pro $(population_base) Einwohner")
    end
    PyPlot.legend(loc="upper left")
    PyPlot.grid()
    PyPlot.show()
    PyPlot.savefig("../docs/$(prefix)-active.png")

    ###########################################################################################################################
    fig = PyPlot.figure(2)
    fig = PyPlot.gcf()
    PyPlot.clf()
    if world
        PyPlot.title("Estimated reproduction number of SARS CoV-2 infections in selected countries\n$(trailer)")
    else
        PyPlot.title("Schätzung der Reproduktionszahl für die SARS CoV-2 Pandemie in den Bundesländern\n$(trailer)")
    end
    fig.set_size_inches(10,5)
    PyPlot.ylim(0,2)
    for country in countries
        @show country[1]
        results=country_results(data,country[1],world=world,infection_period=infection_period,avg_window=avg_window)
        PyPlot.plot_date(results.dates[d0:end],results.est_r0[d0:end],label=country[1],country[2])
        if country[1]=="Italy" || country[1]=="BY"
            PyPlot.plot_date(results.dates[d0:end],[1.0 for i=d0:length(results.dates)],"k--",linewidth=2)
        end
    end
    if !world
        results=nowcast_results(nowcast,active_period=15,infection_period=infection_period,population_base=population_base)
        PyPlot.plot_date(results.dates[d0:end],results.est_r0[d0:end],label="nowcast","r*-")
    end
    if world
        PyPlot.xlabel("Date")
        PyPlot.ylabel("Reproduction number")
    else
        PyPlot.xlabel("Datum")
        PyPlot.ylabel("Reproduktionszahl")
    end
    PyPlot.legend(loc="upper left")
    PyPlot.grid()
    PyPlot.show()
    PyPlot.savefig("../docs/$(prefix)-repro.png")

    ###########################################################################################################################
    fig = PyPlot.figure(3)
    fig = PyPlot.gcf()
    PyPlot.clf()
    if world
        PyPlot.title("Number of newly observed SARS-CoV2 infections in the last 7 days\n$(trailer)")
    else
        PyPlot.title("Anzahl der Neuinfektionen in den zurückliegenden 7 Tagen\n$(trailer)")
    end
    fig.set_size_inches(10,5)
    for country in countries
        results=country_results(data,country[1],world=world,infection_period=infection_period,avg_window=avg_window)
        PyPlot.plot_date(results.dates[d0:end],results.mvavg_new[d0:end],label=country[1],country[2])
    end
    if !world
        results=nowcast_results(nowcast,active_period=15,infection_period=infection_period,population_base=population_base)
        PyPlot.plot_date(results.dates[d0:end],results.est_new[d0:end],label="nowcast","r*-")
    end
    if world
        PyPlot.xlabel("Date")
        PyPlot.ylabel("Cases per $(population_base) inhabitants")
    else
        PyPlot.xlabel("Datum")
        PyPlot.ylabel("Fälle pro $(population_base) Einwohner")
    end
    PyPlot.legend(loc="upper left")
    PyPlot.grid()
    PyPlot.show()
    PyPlot.savefig("../docs/$(prefix)-new.png")

    ###########################################################################################################################
    fig = PyPlot.figure(4)
    fig = PyPlot.gcf()
    PyPlot.clf()
    if world
        PyPlot.title("Number of COVID-19 deaths\n$(trailer)")
    else
        PyPlot.title("Anzahl der COVID-19-Toten\n$(trailer)")
    end
    fig.set_size_inches(10,5)
    for country in countries
        results=country_results(data,country[1],world=world,infection_period=infection_period,avg_window=avg_window)
        ndead=length(results.mvavg_dead[d0:end])
        PyPlot.plot_date(results.dates[d0:d0+ndead-1],results.mvavg_dead[d0:end],label=country[1],country[2])
    end
    if world
        PyPlot.xlabel("Date")
        PyPlot.ylabel("Cases per $(population_base) inhabitants")
    else
        PyPlot.xlabel("Datum")
        PyPlot.ylabel("Fälle pro $(population_base) Einwohner")
    end
    PyPlot.legend(loc="upper left")
    PyPlot.grid()
    PyPlot.show()
    PyPlot.savefig("../docs/$(prefix)-dead.png")
end

function create_plots()
    plot_active_r0(download=true, world=true )
    plot_active_r0(download=false, world=false )
end

# publish on github
function publish(;msg="data update")
    run(`git commit -a -m $(msg)`)
    run(`git push`)
end

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

function countrylist2()
   [
       ["Austria","g-"],
       ["Brazil","g-"],
       ["Canada","y-"],
       ["China","o-"],
       ["Czechia", "-"],
       ["France", "o-"],
       ["Germany","r-o"],
       ["India", "-"],
       ["Italy", "o-"],
       ["Iran", "-"],
       ["Israel", "-"],
       ["Korea, South","o-"],
       ["Poland","-"],
       ["Russia","mo-"],
       ["Slovakia", "-"],
       ["Spain", "o-"],
       ["Sweden", "o-"],
       ["Turkey", "-"],
       ["United Kingdom","k-"],
       ["US","k-o"],
   ]
end

function blaenderlist2()
    [
        ["BE","g-o" ],
        ["BW","y-o" ],
        ["BY","k-o" ],
        ["BB","-"   ],
        ["HB","-"   ],
        ["HH","-"   ],
        ["HE","-"   ],
        ["MV","-"   ],
        ["NI","-"   ],
        ["NW","b-o" ],
        ["RP","-"   ],
        ["SL","-"   ],
        ["SN","-"   ],
        ["ST","-"   ],
        ["SH","y-"  ],
        ["TH","k-"  ],
        ["DE","r-o" ],
    ]
end

function plot_plotly(;download=false, world=true, infection_period=5,avg_window=7)
    if world
        data=alldata_world(download=download)
        countries=countrylist2()
        trailer="<br>Data source: JHU $(Dates.today()) <br> Data processing: https://github.com/j-fu/coronaplot  License: CC-BY 2.0"
        prefix="world"
    else
        if download
            scrape_rki_from_wikipedia()
        end
        data=alldata_blaender()
        countries=blaenderlist2()
        trailer="<br>Datenquelle: RKI via de.wikipedia.org/wiki/COVID-19-Pandemie_in_Deutschland/Statistik $(Dates.today())<br> Datenverarbeitung: https://github.com/j-fu/coronaplot  Lizenz: CC-BY 2.0"
        prefix="de"
        nowcast=load_nowcast(download=false)
    end
    d0=10
    if world
        d0=40
    end

    population_base=100_000

    ###########################################################################################################################
    title="Anzahl der Neuinfektionen in den zurückliegenden 7 Tagen\n$(trailer)"
    xlabel="Datum"
    ylabel="Fälle pro $(population_base) Einwohner"
    if world
        title="Number of newly observed SARS-CoV2 infections in the last 7 days\n$(trailer)"
        xlabel="Date"
        ylabel="Cases per $(population_base) inhabitants"
    end
    traces=Plotly.GenericTrace{Dict{Symbol,Any}}[]
    for country in countries
        results=country_results(data,country[1],world=world,infection_period=infection_period,avg_window=avg_window)
        trace=Plotly.scatter(x=results.dates[d0:end],y=results.mvavg_new[d0:end], mode="lines",name=country[1])
        push!(traces,trace)
    end
    if !world
        results=nowcast_results(nowcast,active_period=15,infection_period=infection_period,population_base=population_base)
        push!(traces,Plotly.scatter(x=results.dates[d0:end],y=results.est_new[d0:end],mode="markers",name="nowcast"))
    end

    layout=Plotly.Layout(;title=title, yaxis_title=ylabel)
    p=Plotly.plot(traces,layout)
    Plotly.savehtml(p,"../docs/$(prefix)-new.html",:remote)

    ################################################################################
    title="Schätzung der Zahl der infektiösen Personen\n$(trailer)"
    ylabel="Fälle pro $(population_base) Einwohner"
    if world
        title="Estimated number of infectious persons in selected countries\n$(trailer)"
        ylabel="Cases per $(population_base) inhabitants"
    end

    traces=Plotly.GenericTrace{Dict{Symbol,Any}}[]
    for country in countries
        results=country_results(data,country[1],
                                world=world,
                                avg_window=avg_window, # Window for moving average of time series of infected people
                                active_period=15, # Period during which we assume an infection is active
                                infection_period=infection_period, # Time it takes for an infected person to infect the next (RKI uses 4)
                                population_base=population_base)
        trace=Plotly.scatter(x=results.dates[d0:end],y=results.est_active[d0:end], mode="lines",name=country[1])
        push!(traces,trace)
    end
    if !world
        results=nowcast_results(nowcast,active_period=15,infection_period=infection_period,population_base=population_base)
        push!(traces,Plotly.scatter(x=results.dates[d0:end],y=results.est_new[d0:end],mode="markers",name="nowcast"))
    end

    layout=Plotly.Layout(;title=title, yaxis_title=ylabel)
    p=Plotly.plot(traces,layout)
    Plotly.savehtml(p,"../docs/$(prefix)-active.html",:remote)

    ###########################################################################################################################
    title="Schätzung der Reproduktionszahl für die SARS CoV-2 Pandemie in den Bundesländern\n$(trailer)"
    ylabel="Reproduktionszahl"
    if world
        title="Estimated reproduction number of SARS CoV-2 infections in selected countries\n$(trailer)"
        ylabel="Reproduction number"
    end
    traces=Plotly.GenericTrace{Dict{Symbol,Any}}[]
    for country in countries
        results=country_results(data,country[1],world=world,infection_period=infection_period,avg_window=avg_window)
        trace=Plotly.scatter(x=results.dates[d0:end],y=results.est_r0[d0:end], mode="lines",name=country[1])
        push!(traces,trace)
    end
    if !world
        results=nowcast_results(nowcast,active_period=15,infection_period=infection_period,population_base=population_base)
        push!(traces,Plotly.scatter(x=results.dates[d0:end],y=results.est_r0[d0:end],mode="markers",name="nowcast"))
    end

    layout=Plotly.Layout(;title=title, yaxis_title=ylabel, yaxis_range=[0,3])
    p=Plotly.plot(traces,layout)
    Plotly.savehtml(p,"../docs/$(prefix)-repro.html",:remote)

    ######################################################################
    title="Anzahl der COVID-19-Toten\n$(trailer)"
    ylabel="Fälle pro $(population_base) Einwohner"

    if world
        title="Number of COVID-19 deaths\n$(trailer)"
        ylabel="Cases per $(population_base) inhabitants"
    end
    traces=Plotly.GenericTrace{Dict{Symbol,Any}}[]
    for country in countries
        results=country_results(data,country[1],world=world,infection_period=infection_period,avg_window=avg_window)
        ndead=length(results.mvavg_dead[d0:end])
        trace=Plotly.scatter(x=results.dates[d0:d0+ndead-1],y=results.mvavg_dead[d0:end], mode="lines",name=country[1])
        push!(traces,trace)
    end
    layout=Plotly.Layout(;title=title, yaxis_title=ylabel)
    p=Plotly.plot(traces,layout)
    Plotly.savehtml(p,"../docs/$(prefix)-dead.html",:remote)
end

function create_plots()
    plot_active_r0(download=true, world=true )
    plot_active_r0(download=true, world=false )
    plot_plotly(download=false, world=true )
    plot_plotly(download=false, world=false )
end

# publish on github
function publish(;msg="data update")
    run(`git commit -a -m $(msg)`)
    run(`git push`)
end

max_stars_repo_name: j-fu/coronaplot

max_issues_repo_name: j-fu/coronaplot

max_forks_repo_name: j-fu/coronaplot


text:
#!/usr/bin/env julia

#
#   Weather update client
#   Connects SUB socket to tcp://localhost:5556
#   Collects weather updates and finds avg temp in zipcode
#

using ZMQ

context = Context()
socket = Socket(context, SUB)

println("Collecting updates from weather server...")
ZMQ.connect(socket, "tcp://localhost:5556")

# Subscribe to zipcode, default is NYC, 10001
zip_filter = length(ARGS) > 0 ? int(ARGS[1]) : 10001

ZMQ.set_subscribe(socket, string(zip_filter))

# Process 5 updates
update_nbr = 5

total_temp = 0
for update in [1:update_nbr]
    message = unsafe_string(ZMQ.recv(socket))
    zipcode, temperature, relhumidity = split(message)
    total_temp += parse(temperature)
end

avg_temp = total_temp / update_nbr

println("Average temperature for zipcode $zip_filter was $(avg_temp)F")

# Making a clean exit.
ZMQ.close(socket)
ZMQ.close(context)

max_stars_repo_name: JailbreakFox/LightWeightRepository

max_issues_repo_name: JailbreakFox/LightWeightRepository

max_forks_repo_name: JailbreakFox/LightWeightRepository


text:
# create the mutable struct and test it
@testset "Test mutable struct" begin
    p = Pet("asdf", 0)
    @test p.hunger == 0
    go_for_a_walk!(p)
    @test p.hunger == 1
    p.hunger = -1
    go_for_a_walk!(p)
    @test p.hunger == 0
end


text:
#This file contains a list of the currently available mean functions

import Base.show

abstract Mean

function show(io::IO, m::Mean, depth::Int = 0)
    pad = repeat(" ", 2*depth)
    print(io, "$(pad)Type: $(typeof(m)), Params: ")
    show(io, get_params(m))
    print(io, "\n")
end

# Calculates the stack [dm / dθᵢ] of mean matrix gradients
function grad_stack(m::Mean, X::Matrix{Float64})
    n = num_params(m)
    d, nobsv = size(X)
    mat = Array(Float64, nobsv, n)
    for i in 1:nobsv
        @inbounds mat[i,:] = grad_mean(m, X[:,i])
    end
    return mat
end

include("mConst.jl")         # Constant mean function, which also contains the zero mean function
include("mLin.jl")           # Linear mean function
include("mPoly.jl")          # Polynomial mean function
include("sum_mean.jl")       # Sum mean functions
include("prod_mean.jl")      # Product of mean functions


text:
##################################################
# Visualising TR Bilevel and data-collecting iteration tools
##################################################

module BilevelVisualise

using Printf
using FileIO
using Setfield
using ColorTypes: Gray
using ImageContrastAdjustment
import GR

using AlgTools.Util
using AlgTools.StructTools
using AlgTools.LinkedLists
using AlgTools.Comms

using VariationalImaging.GradientOps

##############
# Our exports
##############

export BilevelLogEntry,
       bg_bilevel_visualise,
       bilevel_visualise,
       clip,
       #grayimg,
       secs_ns,
       iterate_bilevel_visualise,
       initialise_bilevel_visualisation,
       finalise_bilevel_visualisation

##################
# Data structures
##################

struct BilevelLogEntry <: IterableStruct
    iter :: Int
    time :: Float64
    function_value :: Float64
    gradient_value :: Float64
    radius_value :: Float64
    stopping_criteria :: Float64
end

struct BilevelState
    vis :: Union{Channel,Bool,Nothing}
    visproc :: Union{Nothing,Task}
    start_time :: Union{Real,Nothing}
    wasted_time :: Real
    log :: LinkedList{BilevelLogEntry}
end

##################
# Helper routines
##################

@inline function secs_ns()
    return convert(Float64, time_ns())*1e-9
end

clip = x -> min(max(x, 0.0), 1.0)
grayimg = im -> Gray.(clip.(im))

################
# Visualisation
################

function process_bilevel_channel(fn, rc)
    while true
        d=take!(rc)
        # Take only the latest image to visualise
        while isready(rc)
            d=take!(rc)
        end
        # We're done if we were fed nothing
        if isnothing(d)
            break
        end
        try
            fn(d)
        catch ex
            error("Exception in process_channel handler. Terminating.\n")
            rethrow(ex)
        end
    end
end

function bg_bilevel_visualise(rc)
    process_bilevel_channel(do_visualise, rc)
end

# function do_visualise(imgs)
#     plt = im -> plot(grayimg(im), showaxis=false, grid=false, aspect_ratio=:equal, margin=2mm)
#     display(plot([plt(imgs[i]) for i =1:length(imgs)]..., reuse=true, margin=0mm))
# end

grayGR = x -> begin
    y = round(UInt32, 0xff*clip(x))
    return 0x010101*y + 0xff000000
end

function fill_viewport(vp, c)
    GR.savestate()
    GR.selntran(0)
    GR.setscale(0)
    GR.setfillintstyle(GR.INTSTYLE_SOLID)
    GR.setfillcolorind(c)
    GR.fillrect(vp...)
    GR.selntran(1)
    GR.restorestate()
end

function do_visualise(imgs; refresh=true, fullscreen=false)
    n = length(imgs)

    # Get device dimensions in metres and pixels
    scrw, scrh, pw, ph = GR.inqdspsize()
    imgaspect = n > 0 ? float(size(imgs[1], 1))/float(size(imgs[1], 2)) : 1
    # Scaling to maximum size window
    sc=0.7
    # Set up window and transformations
    GR.clearws()
    GR.setscale(0);
    GR.selntran(1)
    # - First OS window size
    if fullscreen
        w, h = float(scrw), float(scrh)
    elseif scrw/n>scrh
        w, h = float(sc*scrh*n), float(sc*scrh*imgaspect)
    else
        w, h = float(sc*scrw), float(sc*scrw/n*imgaspect)
    end
    GR.setwsviewport(0, w, 0, h)
    # NDC to device
    if w>h
        canvas=[0, 1, 0, h/w]
    else
        canvas=[0, w/h, 0, 1]
    end
    GR.setwswindow(canvas...)
    fill_viewport(canvas, 1)
    # World coordinates to NDC
    if imgaspect/n<h/w
        y0 = (canvas[3]+canvas[4])/2
        ww = (canvas[2]-canvas[1])/2
        y1 = y0-ww*(imgaspect/n)
        y2 = y0+ww*(imgaspect/n)
        GR.setviewport(canvas[1], canvas[2], y1, y2)
    else
        x0 = (canvas[1]+canvas[2])/2
        hh = (canvas[4]-canvas[3])/2
        x1 = x0-hh*(n/imgaspect)
        x2 = x0+hh*(n/imgaspect)
        GR.setviewport(x1, x2, canvas[3], canvas[4])
    end
    GR.setwindow(0, n, 0, 1)
    # Clear background
    # Plot images
    for i=1:n
        im = imgs[i]'
        sz = size(im)
        GR.drawimage(i-1, i, 0, 1, sz[1], sz[2], grayGR.(im))
    end
    if refresh
        GR.updatews()
    end

end

function bilevel_visualise(channel_or_toggle, data)
    if isa(channel_or_toggle, Channel)
        put_onlylatest!(channel_or_toggle, data)
    elseif isa(channel_or_toggle, Bool) && channel_or_toggle
        do_visualise(data)
    end
end

######################################################
# Iterator that does visualisation and log collection
######################################################

function iterate_bilevel_visualise(st :: BilevelState,
                           step :: Function,
                           params :: NamedTuple) where DisplacementT
    try
        stop_flag = false
        for iter=1:params.maxiter
            st = step() do calc_objective
                if isnothing(st.start_time)
                    # The Julia precompiler is a miserable joke, apparently not crossing module
                    # boundaries, so only start timing after the first iteration.
                    st = @set st.start_time=secs_ns()
                end

                verb = params.verbose_iter!=0 && mod(iter, params.verbose_iter) == 0

                if verb || iter ≤ 20 || (iter ≤ 200 && mod(iter, 10) == 0)
                    verb_start = secs_ns()
                    tm = verb_start - st.start_time - st.wasted_time
                    par, x, value, g, Δ, stopping_criteria = calc_objective()

                    entry = BilevelLogEntry(iter, tm, value, g, Δ, stopping_criteria)

                    # (**) Collect a singly-linked list of log to avoid array resizing
                    # while iterating
                    st = @set st.log=LinkedListEntry(entry, st.log)

                    if verb
                        @printf("%d/%d x=%e, f=%.3e, g=%.4e, Δ=%.3e, stop=%.3e\n", iter, params.maxiter, norm₂(par), value, g, Δ, stopping_criteria)
                        if par isa Union{Real,AbstractVector}
                            bilevel_visualise(st.vis, (x,))
                        elseif par isa AbstractArray{Float64,3}
                            pOp = PatchOp(par,x)
                            par_ = Gray.(pOp(par))
                            for i = 1:3
                                if abs(maximum(par_[:,:,i])-minimum(par)) < sqrt(eps())
                                    par_[:,:,i] = (par_[:,:,i] .- minimum(par_[:,:,i])) ./ (maximum(par_[:,:,i]))
                                else
                                    par_[:,:,i] = (par_[:,:,i] .- minimum(par_[:,:,i])) ./ (maximum(par_[:,:,i])-minimum(par_[:,:,i]))
                                end
                            end
                            bilevel_visualise(st.vis, (x,par_[:,:,1],par_[:,:,2],par_[:,:,3],))
                        else
                            pOp = PatchOp(par,x)
                            par_ = Gray.(pOp(par))
                            if abs(maximum(par_)-minimum(par)) < sqrt(eps())
                                par_ = (par_ .- minimum(par_)) ./ (maximum(par_))
                            else
                                par_ = (par_ .- minimum(par_)) ./ (maximum(par_)-minimum(par_))
                            end
                            bilevel_visualise(st.vis, (x,par_,))
                        end
                    end

                    if params.save_iterations
                        fn = t -> "$(params.save_prefix)_$(t)_iter$(iter).png"
                        save(File(format"PNG", fn("reco")), grayimg(x))
                    end

                    st = @set st.wasted_time += (secs_ns() - verb_start)

                    # stop criteria
                    if Δ < params.tol
                        stop_flag = true
                    end
                end

                return st
            end
            if stop_flag
                break
            end
        end
    catch ex
        if isa(ex, InterruptException)
            # If SIGINT is received (user pressed ^C), terminate computations,
            # returning current status. Effectively, we do not call `step()` again,
            # ending the iterations, but letting the algorithm finish up.
            # Assuming (**) above occurs atomically, `st.log` should be valid, but
            # any results returned by the algorithm itself may be partial, as for
            # reasons of efficiency we do *not* store results of an iteration until
            # the next iteration is finished.
            printstyled("\rUser interrupt—finishing up.\n", bold=true, color=202)
        else
            throw(ex)
        end
    end

    return st
end

####################
# Launcher routines
####################

function initialise_bilevel_visualisation(visualise; iterator=iterate_bilevel_visualise)
    # Create visualisation
    if visualise
        rc = Channel(1)
        visproc = Threads.@spawn bg_bilevel_visualise(rc)
        bind(rc, visproc)
        vis = rc
    else
        vis = false
        visproc = nothing
    end

    st = BilevelState(vis, visproc, nothing, 0.0, nothing)
    iterate = curry(iterate_bilevel_visualise, st)

    return st, iterate
end

function finalise_bilevel_visualisation(st)
    if isa(st.vis, Channel)
        # Tell subprocess to finish, and wait
        put!(st.vis, nothing)
        close(st.vis)
        wait(st.visproc)
    end
end

end # Module

max_stars_repo_name: dvillacis/BPLDenoising

max_issues_repo_name: dvillacis/BPLDenoising

max_forks_repo_name: dvillacis/BPLDenoising


text:
module MittagLeffler

export mittleff
export mittlefferr
export mittleffderiv

include("mittag_leffler.jl")

end # module


text:
@testset "Plots" begin
    pomdp = DoubleCovidPOMDP()
    sol = CovidPOMDPs.ProportionalControlSolver()
    planner = solve(sol, pomdp)

    b0 = initialstate(pomdp, Distributions.Uniform(1,10_000))
    s0 = rand(b0)
    a = rand(actions(pomdp))

    h = simulate(
        pomdp,
        b0,
        planner;
        T = 50,
        upd = BootstrapFilter,
        n_p = 1_000,
        progress = false
    )

    @test plot(h) isa Figure;
    @test plot(h; kind=:line) isa Figure;
    @test_throws DomainError plot(h; kind=:wrong);

    @test CovidPOMDPs.plot_inf_belief(h) isa Figure;
    @test CovidPOMDPs.plot_inf_belief(h; particle_samples=100) isa Figure;
end


text:
## 调包
using CSV
using Interpolations
using PyCall
using PyPlot
using DataFrames

## 调 python 包
@pyimport numpy as np
pygui(true)

## 原始数据
"""
原始数据输入
(塔顶产品组成, 塔釜组成) 的质量分数 (ω_d, ω_w)
"""

# (塔顶产品组成, 塔釜组成) 的质量分数
ω_d = 93.5107e-2
ω_w = 12.4097e-2

# 质量分数 -> 摩尔分数 换算函数
f(ω_i) = (ω_i / 46.07) / (ω_i / 46.07 + (1 - ω_i) / 18.02)

# (塔顶产品组成, 塔釜组成) 摩尔分数
x_d = f(ω_d)
x_w = f(ω_w)

## 对角线
x = y = np.linspace(0, 1, 1000)

## 平衡线
# EtOH-H2O 相平衡数据
EtOH = CSV.read("EtOH-H2O.csv", DataFrame)
# EtOH-H2O 相平衡组成插值
interp_y = Interpolations.LinearInterpolation(EtOH.x, EtOH.y)
## 相平衡线
yy = interp_y(x)

## 绘图
# 对角线
PyPlot.plot(x, y, "-")
# 相平衡线
PyPlot.plot(x, yy, "-")
# x_w 垂直线
PyPlot.plot([x_w, x_w], [0, x_w], "b-")
# x_d 垂直线
PyPlot.plot([x_d, x_d], [0, x_d], "b-")

# 绘图设置
PyPlot.xlabel("\$ x \$")
PyPlot.ylabel("\$ y \$")
PyPlot.xlim([0, 1])
PyPlot.ylim([0, 1])

## 塔板数计算
# 迭代点
point_1 = [x_d, x_d]
point_2 = zeros(2)
# 塔板数
num = 0
# 迭代记数
n = 1
# 迭代计算 & 绘图
println("开始迭代")
while true
    global point_1, point_2, num, n
    if point_1[1] < x_w
        break
    end
    println("迭代第 $(n) 次")
    n += 1
    for x_i in np.linspace(0, 1, 10000000)
        if abs(interp_y(x_i) - point_1[2]) < 0.00001 * point_1[2]
            pointTemp = [x_i, point_1[2]]
            point_2 = [x_i, x_i]
            num += 1
            PyPlot.plot([point_1[1], pointTemp[1]], [point_1[2], pointTemp[2]], "b-")
            PyPlot.plot([point_2[1], pointTemp[1]], [point_2[2], pointTemp[2]], "b-")
            # 迭代后新点
            point_1 = point_2
            point_2 = zeros(2)
            println(point_1[1])
            break
        end
    end
end

# 理论塔板数结果
println("理论塔板数 = $(num)")

max_stars_repo_path: Experiment 07/total_reflux.jl

max_issues_repo_path: Experiment 07/total_reflux.jl

max_forks_repo_path: Experiment 07/total_reflux.jl


text:
module BioinformaticsBISC195

export normalizeDNA,
        composition,
        gc_content,
        complement,
        reverse_complement,
        parse_fasta,
        unique_kmers,
        kmerset_distance

# # uncomment the following line if you intend to use BioSequences types
# using BioSequences

"""
    normalizeDNA(::AbstractString)

Ensures that a sequence only contains valid bases
(or `'N'` for unknown bases).
Returns a String.
"""
function normalizeDNA(seq)
    seq = uppercase(string(seq))
    for base in seq
        # note: `N` indicates an unknown base
        occursin(base, "AGCTNYRWMKSHVDB") || error("invalid base $base")
    end
    return seq # change to `return LongDNASeq(seq)` if you want to try to use BioSequences types
end

# Your code here.
# Don't forget to export your functions!

function composition(seq)
    seq = normalizeDNA(seq) # make uppercase string, check invalid bases

    base_dict = Dict('A'=>  0, 'C' => 0, 'G' => 0, 'T'=> 0, 'N' => 0) #create dictionary w/ keys A,C,T,G,N
    for base in seq
         if base == 'A'
            base_dict['A'] = base_dict['A'] +1
         elseif base == 'C'
            base_dict['C'] = base_dict['C'] +1
         elseif base == 'G'
            base_dict['G'] = base_dict['G']+1
        elseif base == 'T'
            base_dict['T'] = base_dict['T']+1
         elseif base == 'N'
            base_dict['N'] = base_dict['N']+1
         end
        end
    return base_dict

end

function gc_content(seq)
    base_dict=composition(seq)
    c=base_dict['C']
    g=base_dict['G']

    total=base_dict['N'] + base_dict['G'] + base_dict['T'] + base_dict['C'] + base_dict['A']   #same as total=content[1]+content[2]+content[3]+content[4]
    #c in second position of the tuple
    #g in third position
    #total is sum of all positions
    return (c+g)/total
end

 function complement(seq::String)
    seq= normalizeDNA(seq)
    complements= Dict('A' => 'T', 'T' => 'A', 'G' => 'C', 'C' => 'G', 'N' => 'N')
    my_complement=[]
    for base in seq
        if occursin(base, "AGTCN")
        push!(my_complement, complements[base])
        end
    end

return join(my_complement)
end

function reverse_complement(seq)
        sequence = reverse(uppercase(seq))
        arr=[]
        for i in sequence
        push!(arr, string(i))  #pushes each base into the array as a string (b/c complement only takes strings)
        end

        comp = map(complement, arr)
        join(comp) #turns array back into a string
end

function parse_fasta(path)

    header_vector = []
    combined_seq_vector = []
    currentvector = [] #for sequences
    for line in eachline(path)
        if startswith(line, ">")
            header = lstrip(line, ['>'])
            push!(header_vector, header) #pushes the header string into the vector
            joined_string = join(currentvector)
            push!(combined_seq_vector, joined_string)
            currentvector = []
        else
            line = strip(line)
            normalizeDNA(line) #only allow valid DNA sequences including N to enter seq vector
            push!(currentvector, line)

        end
    end
    joined_string = join(currentvector)
    push!(combined_seq_vector, joined_string)#need to deal with last sequence
    return header_vector, combined_seq_vector[2:end] #dont return the empty string from the first header
end

end # module BioinformaticsBISC195

function count_kmers(sequence, k) #k is an integer
    1 <= k <= length(sequence) || error("k must be a positive integer less than the length of the sequence")
    kmers = Dict() # initialize dictionary
    stopindex = lastindex(sequence)-(k-1)
    for i in 1:stopindex
        kmer = sequence[i:i+(k-1)] # Change to index the sequence from i to i+k-1
        kmer = normalizeDNA(kmer)
       if haskey(kmers, kmer) # if this kmer is a key the dictionary
       kmers[kmer] = kmers[kmer] + 1 #   add 1 to the value referenced by that kmer
       else # otherwise
       kmers[kmer] = 1 #  make a new entry in the dictionary with a value of 1
    end
end
return kmers #return all unique kmers of length k
end

function unique_kmers(sequence, k) #makes a vector of all of the kmers of length k
    1 <= k <= length(sequence) || error("k must be a positive integer less than the length of the sequence")
    kmers = [] # initialize vector
    stopindex = lastindex(sequence)-(k-1)
    for i in 1:stopindex
        push!(kmers, sequence[i:i+(k-1)]) # Change to index the sequence from i to i+k-1
        kmers[i] = normalizeDNA(kmers[i])
    end
return kmers #return all unique kmers of length k
end

function kmerset_distance(kmerset1, kmerset2) #takes two Sets of kmers (Set is a type, which is an unordered collection of unique elements) and calulates the distance between them
   distance = 1- (length(intersect(kmerset1, kmerset2)) / length(union(kmerset2, kmerset1))) #a distance of 0 is identical and a distance of 1 is no similarity
   return distance
end


text:
## --- General conversions

    # Europium anomalies
    @test eustar(6.5433, 5.9037) ≈ 2.0329978601003864
    @test eustar(34.7773, 6.5433, 5.9037, 0.8904) ≈ 2.0825737578695205

    # Iron oxide conversions
    @test feoconversion(3.5, NaN, NaN, NaN) == 3.5
    @test feoconversion(3.5, NaN, 7.5, NaN) == 7.5
    @test feoconversion(3.5, NaN, 7.5, 10) == 7.5
    @test feoconversion(3.5, 4.4, NaN, NaN) ≈ 7.45916511675966
    @test feoconversion(NaN, 4.4, NaN, NaN) ≈ 3.9591651167596607

    # Other oxide conversion
    D = elementify(["Fe" "Mg" "Ca" "P"; 10000 10000 10000 10000; 10000 10000 10000 10000], importas=:Dict)
    D = oxideconversion(D)
    @test all(D["FeOT"] .≈ (molarmass["Fe"]+molarmass["O"])/molarmass["Fe"])
    @test all(D["MgO"] .≈ (molarmass["Mg"]+molarmass["O"])/molarmass["Mg"])
    @test all(D["CaO"] .≈ (molarmass["Ca"]+molarmass["O"])/molarmass["Ca"])
    @test all(D["P2O5"] .≈ (molarmass["P"]+2.5*molarmass["O"])/molarmass["P"])

    # Weathering indices
    @test CIA(14.8577, 4.5611, 3.29641, 2.3992) ≈ 47.66582778067264
    @test WIP(3.2964, 4.5611, 2.3992, 5.9121) ≈ 78.40320264846837

## -- Perplex name abbreviations

    abbreviations = ("ak", "alm", "and", "andr", "chum", "cz", "crd", "ep", "fa", "fctd", "fcrd", "fep", "fosm", "fst", "fo", "geh", "gr", "hcrd", "tpz", "ky", "larn", "law", "merw", "mctd", "mst", "mnctd", "mncrd", "mnst", "mont", "osm1", "osm2", "phA", "pump", "py", "rnk", "sill", "spss", "sph", "spu", "teph", "ty", "vsv", "zrc", "zo", "acm", "cats", "di", "en", "fs", "hed", "jd", "mgts", "pswo", "pxmn", "rhod", "wo", "anth", "cumm", "fanth", "fgl", "ftr", "ged", "gl", "grun", "parg", "rieb", "tr", "ts", "deer", "fcar", "fspr", "mcar", "spr4", "spr7", "ann", "cel", "east", "fcel", "ma", "mnbi", "mu", "naph", "pa", "phl", "afchl", "ames", "clin", "daph", "fsud", "mnchl", "sud", "atg", "chr", "fta", "kao", "pre", "prl", "ta", "tats", "ab", "anl", "an", "coe", "crst", "heu", "abh", "kals", "lmt", "lc", "me", "mic", "ne", "q", "san", "stlb", "stv", "trd", "wrk", "bdy", "cor", "geik", "hem", "herc", "ilm", "oilm", "lime", "mft", "mt", "mang", "bunsn", "per", "pnt", "ru", "sp", "usp", "br", "dsp", "gth", "ank", "arag", "cc", "dol", "mag", "rhc", "sid", "diam", "gph", "iron", "Ni", "CO2", "CO", "H2", "CH4", "O2", "H2O", "abL", "anL", "diL", "enL", "faL", "fliq", "foL", "h2oL", "hliq", "kspL", "mliq", "qL", "silL", "H+", "Cl-", "OH-", "Na+", "K+", "Ca++", "Mg++", "Fe++", "Al+++", "CO3", "AlOH3", "AlOH4-", "KOH", "HCL", "KCL", "NaCl", "CaCl2", "CaCl+", "MgCl2", "MgCl", "FeCl2", "aqSi", "Augite(G)", "Cpx(JH)", "Cpx(l)", "Cpx(h)", "Cpx(stx)", "Cpx(stx7)", "Omph(HP)", "Cpx(HP)", "Cpx(m)", "Cpx(stx8)", "Omph(GHP)", "cAmph(G)", "Cumm", "Gl", "Tr", "GlTrTsPg", "Amph(DHP)", "Amph(DPW)", "Ca-Amph(D)", "Na-Amph(D)", "Act(M)", "GlTrTsMr", "cAmph(DP)", "melt(G)", "melt(W)", "melt(HP)", "melt(HGP)", "pMELTS(G)", "mMELTS(G)", "LIQ(NK)", "LIQ(EF)", "Chl(W)", "Chl(HP)", "Chl(LWV)", "O(JH)", "O(SG)", "O(HP)", "O(HPK)", "O(stx)", "O(stx7)", "Ol(m)", "O(stx8)", "Sp(JH)", "GaHcSp", "Sp(JR)", "Sp(GS)", "Sp(HP)", "Sp(stx)", "CrSp", "Sp(stx7)", "Sp(WPC)", "Sp(stx8)", "Pl(JH)", "Pl(h)", "Pl(stx8)", "Kf", "San", "San(TH)", "Grt(JH)", "Gt(W)", "CrGt", "Gt(MPF)", "Gt(B)", "Gt(GCT)", "Gt(HP)", "Gt(EWHP)", "Gt(WPH)", "Gt(stx)", "Gt(stx8)", "Gt(WPPH)", "ZrGt(KP)", "Maj", "Opx(JH)", "Opx(W)", "Opx(HP)", "CrOpx(HP)", "Opx(stx)", "Opx(stx8)", "Mica(W)", "Pheng(HP)", "MaPa", "Mica(CF)", "Mica(CHA1)", "Mica(CHA)", "Mica+(CHA)", "Mica(M)", "Mica(SGH)", "Ctd(W)", "Ctd(HP)", "Ctd(SGH)", "St(W)", "St(HP)", "Bi(W)", "Bio(TCC)", "Bio(WPH)", "Bio(HP)", "Crd(W)", "hCrd", "Sa(WP)", "Sapp(HP)", "Sapp(KWP)", "Sapp(TP)", "Osm(HP)", "F", "F(salt)", "COH-Fluid", "Aq_solven0", "WADDAH", "T", "Scap", "Carp", "Carp(M)", "Carp(SGH)", "Sud(Livi)", "Sud", "Sud(M)", "Anth", "o-Amph", "oAmph(DP)", "feldspar", "feldspar_B", "Pl(I1,HP)", "Fsp(C1)", "Do(HP)", "M(HP)", "Do(AE)", "Cc(AE)", "oCcM(HP)", "Carb(M)", "oCcM(EF)", "dis(EF)", "IlHm(A)", "IlGkPy", "Ilm(WPH)", "Ilm(WPH0)", "Neph(FB)", "Chum", "Atg(PN)", "B", "Pu(M)", "Stlp(M)", "Wus",)
    common_names = ("akermanite", "almandine", "andalusite", "andradite", "clinohumite", "clinozoisite", "cordierite", "epidote", "fayalite", "Fe-chloritoid", "Fe-cordierite", "Fe-epidote", "Fe-osumilite", "Fe-staurolite", "forsterite", "gehlenite", "grossular", "hydrous cordierite", "hydroxy-topaz", "kyanite", "larnite", "lawsonite", "merwinite", "Mg-chloritoid", "Mg-staurolite", "Mn-chloritoid", "Mn-cordierite", "Mn-staurolite", "monticellite", "osumilite(1)", "osumilite(2)", "phase A", "pumpellyite", "pyrope", "rankinite", "sillimanite", "spessartine", "sphene", "spurrite", "tephroite", "tilleyite", "vesuvianite", "zircon", "zoisite", "acmite", "Ca-tschermakite", "diopside", "enstatite", "ferrosilite", "hedenbergite", "jadeite", "Mg-tschermakite", "pseudowollastonite", "pyroxmangite", "rhodonite", "wollastonite", "anthophyllite", "cummingtonite", "Fe-anthophyllite", "Fe-glaucophane", "ferroactinolite", "gedrite", "glaucophane", "grunerite", "pargasite", "riebeckite", "tremolite", "tschermakite", "deerite", "Fe-carpholite", "Fe-sapphirine(793)", "Mg-carpholite", "sapphirine(442)", "sapphirine(793)", "annite", "celadonite", "eastonite", "Fe-celadonite", "margarite", "Mn-biotite", "muscovite", "Na-phlogopite", "paragonite", "phlogopite", "Al-free chlorite", "amesite", "clinochlore", "daphnite", "Fe-sudoite", "Mn-chlorite", "sudoite", "antigorite", "chrysotile", "Fe-talc", "kaolinite", "prehnite", "pyrophyllite", "talc", "tschermak-talc", "albite", "analcite", "anorthite", "coesite", "cristobalite", "heulandite", "highalbite", "kalsilite", "laumontite", "leucite", "meionite", "microcline", "nepheline", "quartz", "sanidine", "stilbite", "stishovite", "tridymite", "wairakite", "baddeleyite", "corundum", "geikielite", "hematite", "hercynite", "ilmenite", "ilmenite(ordered)", "lime", "magnesioferrite", "magnetite", "manganosite", "nickel oxide", "periclase", "pyrophanite", "rutile", "spinel", "ulvospinel", "brucite", "diaspore", "goethite", "ankerite", "aragonite", "calcite", "dolomite", "magnesite", "rhodochrosite", "siderite", "diamond", "graphite", "iron", "nickel", "carbon dioxide", "carbon monoxide", "hydrogen", "methane", "oxygen", "water fluid", "albite liquid", "anorthite liquid", "diopside liquid", "enstatite liquid", "fayalite liquid", "Fe-liquid (in KFMASH)", "forsterite liquid", "H2O liquid", "H2O liquid (in KFMASH)", "K-feldspar liquid", "Mg liquid (in KFMASH)", "Silica liquid", "Sillimanite liquid", "H+(aq)", "Cl(aq)", "OH(aq)", "Na+(aq)", "K+(aq)", "Ca2+(aq)", "Mg2+(aq)", "Fe2+(aq)", "Al3+(aq)", "CO3--(aq)", "Al(OH)3(aq)", "Al(OH)4----(aq)", "KOH(aq)", "HCl(aq)", "KCl(aq)", "NaCl(aq)", "CaCl(aq)", "CaCl+(aq)", "MgCl2(aq)", "MgCl+(aq)", "FeCl(aq)", "Aqueous silica", "clinopyroxene", "clinopyroxene", "clinopyroxene", "clinopyroxene", "clinopyroxene", "clinopyroxene", "clinopyroxene", "clinopyroxene", "clinopyroxene", "clinopyroxene", "clinopyroxene", "clinoamphibole", "clinoamphibole", "clinoamphibole", "clinoamphibole", "clinoamphibole", "clinoamphibole", "clinoamphibole", "clinoamphibole", "clinoamphibole", "clinoamphibole", "clinoamphibole", "clinoamphibole", "melt", "melt", "melt", "melt", "melt", "melt", "melt", "melt", "chlorite", "chlorite", "chlorite", "olivine", "olivine", "olivine", "olivine", "olivine", "olivine", "olivine", "olivine", "spinel", "spinel", "spinel", "spinel", "spinel", "spinel", "spinel", "spinel", "spinel", "spinel", "plagioclase", "plagioclase", "plagioclase", "k-feldspar", "k-feldspar", "k-feldspar", "garnet", "garnet", "garnet", "garnet", "garnet", "garnet", "garnet", "garnet", "garnet", "garnet", "garnet", "garnet", "garnet", "garnet", "orthopyroxene", "orthopyroxene", "orthopyroxene", "orthopyroxene", "orthopyroxene", "orthopyroxene", "white mica", "white mica", "white mica", "white mica", "white mica", "white mica", "white mica", "white mica", "white mica", "chloritoid", "chloritoid", "chloritoid", "staurolite", "staurolite", "biotite", "biotite", "biotite", "biotite", "cordierite", "cordierite", "sapphirine", "sapphirine", "sapphirine", "sapphirine", "osumilite", "fluid", "fluid", "fluid", "fluid", "fluid", "talc", "scapolite", "carpholite", "carpholite", "carpholite", "sudoite", "sudoite", "sudoite", "orthoamphibole", "orthoamphibole", "orthoamphibole", "ternary feldspar", "ternary feldspar", "ternary feldspar", "ternary feldspar", "calcite", "calcite", "calcite", "calcite", "calcite", "calcite", "calcite", "calcite", "ilmenite", "ilmenite", "ilmenite", "ilmenite", "nepheline", "clinohumite", "serpentine", "brucite", "pumpellyite", "stilpnomelane", "wüstite",)

    @test perplex_common_name.(abbreviations) == common_names

    abbreviations = ("ak", "alm", "and", "andr", "chum", "cz", "crd", "ep", "fa", "fctd", "fcrd", "fep", "fosm", "fst", "fo", "geh", "gr", "hcrd", "tpz", "ky", "larn", "law", "merw", "mctd", "mst", "mnctd", "mncrd", "mnst", "mont", "osm1", "osm2", "phA", "pump", "py", "rnk", "sill", "spss", "sph", "spu", "teph", "ty", "vsv", "zrc", "zo", "acm", "cats", "di", "en", "fs", "hed", "jd", "mgts", "pswo", "pxmn", "rhod", "wo", "anth", "cumm", "fanth", "fgl", "ftr", "ged", "gl", "grun", "parg", "rieb", "tr", "ts", "deer", "fcar", "fspr", "mcar", "spr4", "spr7", "ann", "cel", "east", "fcel", "ma", "mnbi", "mu", "naph", "pa", "phl", "afchl", "ames", "clin", "daph", "fsud", "mnchl", "sud", "atg", "chr", "fta", "kao", "pre", "prl", "ta", "tats", "ab", "anl", "an", "coe", "crst", "heu", "abh", "kals", "lmt", "lc", "me", "mic", "ne", "q", "san", "stlb", "stv", "trd", "wrk", "bdy", "cor", "geik", "hem", "herc", "ilm","oilm","lime", "mft", "mt", "mang", "bunsn", "per", "pnt", "ru", "sp", "usp", "br", "dsp", "gth", "ank", "arag", "cc", "dol", "mag", "rhc", "sid", "diam", "gph", "iron", "Ni", "CO2", "CO", "H2", "CH4", "O2", "H2O", "abL", "anL", "diL", "enL", "faL", "fliq", "foL", "h2oL", "hliq", "kspL", "mliq", "qL", "silL", "H+", "Cl-", "OH-", "Na+", "K+", "Ca++", "Mg++", "Fe++", "Al+++", "CO3", "AlOH3", "AlOH4-", "KOH", "HCL", "KCL", "NaCl", "CaCl2", "CaCl+", "MgCl2", "MgCl", "FeCl2", "aqSi",)
    full_names = ("akermanite", "almandine", "andalusite", "andradite", "clinohumite", "clinozoisite", "cordierite", "epidote(ordered)", "fayalite", "Fe-chloritoid", "Fe-cordierite", "Fe-epidote", "Fe-osumilite", "Fe-staurolite", "forsterite", "gehlenite", "grossular", "hydrous cordierite", "hydroxy-topaz", "kyanite", "larnite-bredigite", "lawsonite", "merwinite", "Mg-chloritoid", "Mg-staurolite", "Mn-chloritoid", "Mn-cordierite", "Mn-staurolite", "monticellite", "osumilite(1)", "osumilite(2)", "phase A", "pumpellyite", "pyrope", "rankinite", "sillimanite", "spessartine", "sphene", "spurrite", "tephroite", "tilleyite", "vesuvianite", "zircon", "zoisite", "acmite", "Ca-tschermaks pyroxene", "Diopside", "enstatite", "ferrosilite", "hedenbergite", "jadeite", "mg-tschermak", "pseudowollastonite", "pyroxmangite", "rhodonite", "wollastonite", "anthophyllite", "cummingtonite", "Fe-anthophyllite", "Fe-glaucophane", "ferroactinolite", "gedrite(Na-free)", "glaucophane", "grunerite", "pargasite", "riebeckite", "tremolite", "tschermakite", "deerite", "fe-carpholite", "fe-sapphirine(793)", "mg-carpholite", "sapphirine(442)", "sapphirine(793)", "annite", "celadonite", "eastonite", "Fe-celadonite", "margarite", "Mn-biotite", "muscovite", "Na-phlogopite", "paragonite", "phlogopite", "Al-free chlorite", "amesite(14Ang)", "clinochlore(ordered)", "daphnite", "Fe-sudoite", "Mn-chlorite", "Sudoite", "antigorite", "chrysotile", "Fe-talc", "Kaolinite", "prehnite", "pyrophyllite", "talc", "tschermak-talc", "albite", "analcite", "anorthite", "coesite", "cristobalite", "heulandite", "highalbite", "kalsilite", "laumontite", "leucite", "meionite", "microcline", "nepheline", "quartz", "sanidine", "stilbite", "stishovite", "tridymite", "wairakite", "baddeleyite", "corundum", "geikielite", "hematite", "hercynite", "ilmenite", "ilmenite(ordered)","lime", "magnesioferrite", "magnetite", "manganosite", "nickel oxide", "periclase", "pyrophanite", "rutile", "spinel", "ulvospinel", "brucite", "diaspore", "goethite", "ankerite", "aragonite", "calcite", "dolomite", "magnesite", "rhodochrosite", "siderite", "diamond", "graphite", "iron", "nickel", "carbon dioxide", "carbon monoxide", "hydrogen", "methane", "oxygen", "water fluid", "albite liquid", "anorthite liquid", "diopside liquid", "enstatite liquid", "fayalite liquid", "Fe-liquid (in KFMASH)", "Forsterite liquid", "H2O liquid", "H2O liquid (in KFMASH)", "K-feldspar liquid", "Mg liquid (in KFMASH)", "Silica liquid", "Sillimanite liquid", "H+(aq)", "Cl(aq)", "OH(aq)", "Na+(aq)", "K+(aq)", "Ca2+(aq)", "Mg2+(aq)", "Fe2+(aq)", "Al3+(aq)", "CO3--(aq)", "Al(OH)3(aq)", "Al(OH)4----(aq)", "KOH(aq)", "HCl(aq)", "KCl(aq)", "NaCl(aq)", "CaCl(aq)", "CaCl+(aq)", "MgCl2(aq)", "MgCl+(aq)", "FeCl(aq)", "Aqueous silica",)

    @test perplex_expand_name.(abbreviations) == full_names
    @test perplex_abbreviate_name.(full_names) == abbreviations
    @test perplex_phase_is_solid.(("melt(HGP)", "q", "diL", "andr", "T(K)")) == (false, true, false, true, false)

    @test findall(germ_perplex_name_matches.(germ_kd["minerals"], germ_kd["minerals"])) == [3, 12, 18]

## --- Saturation models

    #          SiO2,  TiO2,  Al2O3,  FeOT,   MnO,   MgO,   CaO,   Na2O,  K2O,  P2O5
    majors = [58.509, 1.022, 14.858, 4.371, 0.141, 4.561, 5.912, 3.296, 2.399, 0.279]
    @test tzirc(majors..., 100) ≈ 602.8489762809595
    @test tzircZr(majors..., 800) ≈ 832.9689080567883
    @test all(tzircM((repeat([m],2) for m in majors)...,) .≈ 2.335918319204001)

    @test StatGeochem.Ayers_tsphene(majors...) ≈ 637.139776663209
    @test StatGeochem.Ayers_tspheneTiO2(majors..., 800) ≈ 2.3545537746637324
    @test all(StatGeochem.Ayers_tspheneC.((repeat([m],2) for m in majors)...,) .≈ 2.4338232746497326)

    #                SiO2,   TiO2,  Al2O3,  FeOT,   MgO,  CaO,   Na2O,  K2O,   Li2O, H2O
    montel_elems = [58.509, 1.022, 14.858, 4.371, 4.561, 5.912, 3.296, 2.399, 0.01, 4.0]
    @test StatGeochem.Montel_tmonaziteREE(montel_elems..., 750.0) ≈ 12.03834338792398
    @test StatGeochem.Montel_tmonazite(montel_elems..., 100,100,100,0,0,0) ≈ 630.4499586271999

    @test StatGeochem.Rusiecka_tmonaziteREE(200, 750) ≈ 0.27430570654873154
    @test StatGeochem.Rusiecka_txenotimeY(200, 750) ≈ 41.9312030248943

    @test StatGeochem.Harrison_tapatiteP2O5(58.509, 14.858, 5.912, 3.296, 2.399, 750.) ≈ 0.10142278764336987
    @test StatGeochem.Harrison_tapatiteP(58.509, 14.858, 5.912, 3.296, 2.399, 750.) ≈ 442.6362451135793
    @test StatGeochem.Harrison_tapatiteP2O5(58.509, 750.) ≈ 0.10142278764336987
    @test StatGeochem.Harrison_tapatite(58.509, 0.1) ≈ 748.6127179814277

    #          SiO2,  TiO2,  Al2O3,  FeOT,   MgO,   CaO,   Na2O,  K2O,  P2O5
    majors = [58.509, 1.022, 14.858, 4.371, 4.561, 5.912, 3.296, 2.399, 0.279]
    @test StatGeochem.Tollari_tapatite(majors...) ≈ 521.0594433599132
    @test StatGeochem.Tollari_tapatiteP2O5(58.509,5.912,750.) ≈ 0.5011681927262436

## --- End of File