text:
import TimeSeries: TimeArray

@testset "Time resolution" begin
    twomins = TimeArray([DateTime(today()) + Dates.Minute(i * 2) for i in 1:5], ones(5))
    oneday = TimeArray([DateTime(today()) + Dates.Day(i) for i in 1:5], ones(5))
    onesec = TimeArray([DateTime(today()) + Dates.Second(i) for i in 1:5], ones(5))
    onehour = TimeArray([DateTime(today()) + Dates.Hour(i) for i in 1:5], ones(5))

    @test PowerSystems.getresolution(twomins) == Dates.Minute(2)
    @test PowerSystems.getresolution(oneday) == Dates.Day(1)
    @test PowerSystems.getresolution(onesec) == Dates.Second(1)
    @test PowerSystems.getresolution(onehour) == Dates.Hour(1)
end

@testset "Angle limits" begin
    nodes5 = [
        Bus(1, "nodeA", PowerSystems.PV::BusType, 0, 1.0, (min = 0.9, max=1.05), 230),
        Bus(2, "nodeB", PowerSystems.PQ::BusType, 0, 1.0, (min = 0.9, max=1.05), 230),
        Bus(3, "nodeC", PowerSystems.PV::BusType, 0, 1.0, (min = 0.9, max=1.05), 230),
        Bus(4, "nodeD", PowerSystems.REF::BusType, 0, 1.0, (min = 0.9, max=1.05), 230),
        Bus(5, "nodeE", PowerSystems.PV::BusType, 0, 1.0, (min = 0.9, max=1.05), 230),
    ]

    branches_test = [
        Line("1", true, 0.0, 0.0, Arc(from=nodes5[1], to=nodes5[2]), 0.00281, 0.0281,
             (from=0.00356, to=0.00356), 400.0, (min=-360.0, max=360.0)),
        Line("2", true, 0.0, 0.0, Arc(from=nodes5[1], to=nodes5[4]), 0.00304, 0.0304,
             (from=0.00329, to=0.00329), 3960.0, (min=-360.0, max=75.0)),
        Line("3", true, 0.0, 0.0, Arc(from=nodes5[1], to=nodes5[5]), 0.00064, 0.0064,
             (from=0.01563, to=0.01563), 18812.0, (min=-75.0, max=360.0)),
        Line("4", true, 0.0, 0.0, Arc(from=nodes5[2], to=nodes5[3]), 0.00108, 0.0108,
             (from=0.00926, to=0.00926), 11148.0, (min=0.0, max=0.0)),
        Line("5", true, 0.0, 0.0, Arc(from=nodes5[3], to=nodes5[4]), 0.00297, 0.0297,
             (from=0.00337, to=0.00337), 4053.0, (min=-1.2, max=60.0)),
        Line("6", true, 0.0, 0.0, Arc(from=nodes5[4], to=nodes5[5]), 0.00297, 0.0297,
             (from=0.00337, to=00.00337), 240.0, (min=-1.17, max=1.17))
    ]

    foreach(x -> PowerSystems.check_angle_limits!(x), branches_test)

    @test branches_test[1].anglelimits == (min=-pi/2, max=pi/2)
    @test branches_test[2].anglelimits == (min=-pi/2, max=75.0 * (π / 180))
    @test branches_test[3].anglelimits == (min=-75.0 * (π / 180), max=pi/2)
    @test branches_test[4].anglelimits == (min=-pi/2, max=pi/2)
    @test branches_test[5].anglelimits == (min=-1.2, max=60.0 * (π / 180))
    @test branches_test[6].anglelimits == (min=-1.17, max=1.17)

    bad_angle_limits = Line("1", true, 0.0, 0.0, Arc(from=nodes5[1], to=nodes5[2]), 0.00281, 0.0281,
             (from=0.00356, to=0.00356), 400.0, (min=360.0, max=-360.0))

    @test_throws(PowerSystems.DataFormatError,
                 PowerSystems.check_angle_limits!(bad_angle_limits))
end


text:
# taken from https://github.com/mapbox/wgs84/blob/master/index.js

const RADIUS = 6378137
const FLATTENING_DENOM = 298.257223563
const FLATTENING = 1.0 / FLATTENING_DENOM
const POLAR_RADIUS = RADIUS * (1.0 - FLATTENING)

ext: jl


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	write_emissions(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)

Function for reporting time-dependent CO$_2$ emissions by zone.

"""
function write_emissions(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
	dfGen = inputs["dfGen"]
	G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)
	T = inputs["T"]     # Number of time steps (hours)
	Z = inputs["Z"]     # Number of zones
	L = inputs["L"]     # Number of transmission lines
	W = inputs["REP_PERIOD"]     # Number of subperiods
    SEG = inputs["SEG"] # Number of load curtailment segments

	if (setup["WriteShadowPrices"]==1 || setup["UCommit"]==0 || (setup["UCommit"]==2 && (setup["Reserves"]==0 || (setup["Reserves"]>0 && inputs["pDynamic_Contingency"]==0)))) # fully linear model
		# CO2 emissions by zone

		if setup["CO2Cap"]>=1
			# Dual variable of CO2 constraint = shadow price of CO2
			tempCO2Price = zeros(Z,inputs["NCO2Cap"])
			if has_duals(EP) == 1
				for cap in 1:inputs["NCO2Cap"]
					for z in findall(x->x==1, inputs["dfCO2CapZones"][:,cap])
						tempCO2Price[z,cap] = dual.(EP[:cCO2Emissions_systemwide])[cap]
						# when scaled, The objective function is in unit of Million US$/kton, thus k$/ton, to get $/ton, multiply 1000
						if setup["ParameterScale"] ==1
							tempCO2Price[z,cap] = tempCO2Price[z,cap]* ModelScalingFactor
						end
					end
				end
			end
			dfEmissions = hcat(DataFrame(Zone = 1:Z), DataFrame(tempCO2Price, :auto), DataFrame(AnnualSum = Array{Union{Missing,Float64}}(undef, Z)))
			auxNew_Names=[Symbol("Zone"); [Symbol("CO2_Price_$cap") for cap in 1:inputs["NCO2Cap"]]; Symbol("AnnualSum")]
			rename!(dfEmissions,auxNew_Names)
		else
			dfEmissions = DataFrame(Zone = 1:Z, AnnualSum = Array{Union{Missing,Float32}}(undef, Z))
		end

		for i in 1:Z
			if setup["ParameterScale"]==1
				dfEmissions[!,:AnnualSum][i] = sum(inputs["omega"].*value.(EP[:eEmissionsByZone])[i,:])*ModelScalingFactor
			else
				dfEmissions[!,:AnnualSum][i] = sum(inputs["omega"].*value.(EP[:eEmissionsByZone])[i,:])/ModelScalingFactor
			end
		end

		if setup["ParameterScale"]==1
			dfEmissions = hcat(dfEmissions, DataFrame(value.(EP[:eEmissionsByZone])*ModelScalingFactor, :auto))
		else
			dfEmissions = hcat(dfEmissions, DataFrame(value.(EP[:eEmissionsByZone])/ModelScalingFactor, :auto))
		end

		if setup["CO2Cap"]>=1
			auxNew_Names=[Symbol("Zone");[Symbol("CO2_Price_$cap") for cap in 1:inputs["NCO2Cap"]];Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
			rename!(dfEmissions,auxNew_Names)
			total = DataFrame(["Total" zeros(1,inputs["NCO2Cap"]) sum(dfEmissions[!,:AnnualSum]) fill(0.0, (1,T))], :auto)
			for t in 1:T
				total[:,t+inputs["NCO2Cap"]+2] .= sum(dfEmissions[:,Symbol("t$t")][1:Z])
			end
			rename!(total,auxNew_Names)
			dfEmissions = vcat(dfEmissions, total)
		else
			auxNew_Names=[Symbol("Zone"); Symbol("AnnualSum"); [Symbol("t$t") for t in 1:T]]
			rename!(dfEmissions,auxNew_Names)
			total = DataFrame(["Total" sum(dfEmissions[!,:AnnualSum]) fill(0.0, (1,T))], :auto)
			for t in 1:T
				total[:,t+2] .= sum(dfEmissions[:,Symbol("t$t")][1:Z])
			end
			rename!(total,auxNew_Names)
			dfEmissions = vcat(dfEmissions, total)
		end

## Aaron - Combined elseif setup["Dual_MIP"]==1 block with the first block since they were identical. Why do we have this third case? What is different about it?
	else
		# CO2 emissions by zone
		dfEmissions = hcat(DataFrame(Zone = 1:Z), DataFrame(AnnualSum = Array{Union{Missing,Float64}}(undef, Z)))
		for i in 1:Z
			if setup["ParameterScale"]==1
				dfEmissions[!,:AnnualSum][i] = sum(inputs["omega"].*value.(EP[:eEmissionsByZone])[i,:]) *ModelScalingFactor
			else
				dfEmissions[!,:AnnualSum][i] = sum(inputs["omega"].*value.(EP[:eEmissionsByZone])[i,:])/ModelScalingFactor
			end
		end
		if setup["ParameterScale"]==1
			dfEmissions = hcat(dfEmissions, DataFrame(value.(EP[:eEmissionsByZone])*ModelScalingFactor, :auto))
		else
			dfEmissions = hcat(dfEmissions, DataFrame(value.(EP[:eEmissionsByZone])/ModelScalingFactor, :auto))
		end
		auxNew_Names=[Symbol("Zone");Symbol("AnnualSum");[Symbol("t$t") for t in 1:T]]
		rename!(dfEmissions,auxNew_Names)
		total = DataFrame(["Total" sum(dfEmissions[!,:AnnualSum]) fill(0.0, (1,T))], :auto)
		for t in 1:T
			total[:,t+2] .= sum(dfEmissions[:,Symbol("t$t")][1:Z])
		end
		rename!(total,auxNew_Names)
		dfEmissions = vcat(dfEmissions, total)
	end
	CSV.write(joinpath(path, "emissions.csv"), dftranspose(dfEmissions, false), writeheader=false)
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
include("../test_fem_sa.jl")

nparts = (2,2)
prun(test_fem_sa,sequential,nparts)

nparts = 4
prun(test_fem_sa,sequential,nparts)


text:
module Execs

using ..PFB: targets

using Printf
using Statistics

using CSV
using DataFrames

# Load runs from execs.csv
function load(path)
    df = DataFrame(CSV.File(path))
    disallowmissing!(df)
    return df
end

# Agggregate execs accross runs with mean
aggexecs(df) = combine(groupby(df, [:subject, :fuzzer]),
                       :execs => mean => :execs_mean,
                       :execs => std => :execs_std)

function multiple(dir::String)
    options, choices = targets(dir, false)
    length(choices) > 0 || return

    df = DataFrame()
    for choice in choices
        target = options[choice]
        append!(df, load(target))
    end

    return df
end

function mktable(dir::String)
    df = multiple(dir)
    df == nothing && return
    agg = aggexecs(df)
    table = DataFrame()
    for gd in groupby(agg, :subject)
        row = Dict{String, Any}("subject" => gd.subject[1])
        mx = -Inf
        mx_f = nothing
        for fuzzer in unique(gd.fuzzer)
            x = @view gd[gd.fuzzer .== fuzzer, :]
            m = x[1, :execs_mean]
            if m >= mx
                mx = m
                mx_f = fuzzer
            end
        end
        for fuzzer in unique(gd.fuzzer)
            x = @view gd[gd.fuzzer .== fuzzer, :]
            m = x[1, :execs_mean]
            s = x[1, :execs_std]
            if fuzzer == mx_f
                row[fuzzer] = @sprintf("\$\\mathbf{%.1f} \\pm %.1f\$", m, s)
            else
                row[fuzzer] = @sprintf("\$%.1f \\pm %.1f\$", m, s)
            end
        end
        if !("aflpp" in keys(row))
            row["aflpp"] = "-"
        end
        push!(table, row, cols=:union)
    end
    select!(table, "subject" => "target", "aflnet", "aflnet-no-state", "aflnwe", "aflpp",
           "nyx", "nyx-balanced", "nyx-aggressive")
    sort!(table, :target)
    for row in eachrow(table)
        println(join(Tuple(row), " & "), " \\\\")
    end
    return table
end

end # module

max_stars_repo_name: cyruscyliu/nyx-net-profuzzbench

max_issues_repo_name: cyruscyliu/nyx-net-profuzzbench

max_forks_repo_name: cyruscyliu/nyx-net-profuzzbench


text:
using Images, ImageSegmentation, Statistics, ProgressMeter, Printf, FileIO
using CoordinateTransformations, Rotations, OffsetArrays, Optim, FFTW
using REALMSupport
kernel = Kernel.gaussian((5, 5, 5))
FFTW.set_num_threads(16)

function beadstest(img, filename)
    x_width = []
    y_width = []
    z_width = []
    p_axis = []
    #peak = []

    (x, y, z) = size(img)
    filtered = zeros(x, y, z)

    M = [
        -0.140292 -0.0912973 1.13447
        -0.0273906 -0.940779 -0.0213417
        -1.00816 0.0545783 0.161258
    ]
    v = [0.0, 0.0, 0.0]
    rot = AffineMap(M, v)
    #tilt = LinearMap(RotY(-0.19))
    img_r = warp(img, rot)
    #img_r = warp(img,tilt)
    filtered = warp(filtered, rot)

    #ro = imfilter(img_r,kernel)

    replace!(img_r, NaN => 0)

    beads = findlocalmaxima(img_r)
    select = img_r[beads] .> mean(img_r) + 2 * std(img_r)
    @show(size(beads[select]), mean(img_r) + 1 * std(img_r))

    @showprogress @sprintf("Filtering of Record %s...", filename) for coord in beads[select]
        checkbounds(Bool, img_r, coord[1] - 10, coord[2] - 30, coord[3] - 10) || continue
        checkbounds(Bool, img_r, coord[1] + 10, coord[2] + 30, coord[3] + 10) || continue
        img_window =
            img_r[coord[1]-10:coord[1]+10, coord[2]-30:coord[2]+30, coord[3]-10:coord[3]+10]
        img_x = OffsetArray(vec(mean(mean(img_window, dims = 3), dims = 2)), -10:10)
        x_res = gauss_line_fit(img_x; g_abstol = 1e-14)

        img_y = OffsetArray(vec(mean(mean(img_window, dims = 3), dims = 1)), -30:30)
        y_res = gauss_line_fit(img_y; g_abstol = 1e-14)

        img_z = OffsetArray(vec(mean(mean(img_window, dims = 2), dims = 1)), -10:10)
        z_res = gauss_line_fit(img_z; g_abstol = 1e-14)

        x_params = Optim.minimizer(x_res)
        y_params = Optim.minimizer(y_res)
        z_params = Optim.minimizer(z_res)

        a = x_params[2]
        b = y_params[2]
        c = z_params[2]
        if a > 50 || b > 50 || c > 50 || a < 0 || b < 0 || c < 0
            continue
        end
        #if (x_params[1]+y_params[1]+z_params[1])/3 < mean(img_r) + 3*std(img_r)
        #	continue
        #end
        push!(x_width, a)
        push!(y_width, b)
        push!(z_width, c)
        push!(p_axis, coord[2])
        #push!(peak,(x_params[1]+y_params[1])/2)
        @show(x_params[2], y_params[2], z_params[2], coord)

        x_ = Int(floor(x_params[2]))
        y_ = Int(floor(y_params[2]))
        z_ = Int(floor(z_params[2]))

        checkbounds(Bool, img_r, coord[1] - x_, coord[2] - y_, coord[3] - z_) || continue
        checkbounds(Bool, img_r, coord[1] + x_, coord[2] + y_, coord[3] + z_) || continue

        filtered[
            coord[1]-x_:coord[1]+x_,
            coord[2]-y_:coord[2]+y_,
            coord[3]-z_:coord[3]+z_,
        ] = img_r[coord[1]-x_:coord[1]+x_, coord[2]-y_:coord[2]+y_, coord[3]-z_:coord[3]+z_]

    end
    replace!(filtered, NaN => 0)
    img3 = Gray.(convert(Array{N0f16}, OffsetArrays.no_offset_view(filtered)))
    img4 = Gray.(convert.(Normed{UInt16,16}, img3))

    img_save(img4, "/home/jchang/image/result/", @sprintf("%s-fi.tif", filename))
    @show(mean(x_width), mean(y_width), mean(z_width))
    return x_width, y_width, z_width, p_axis
end

function supportfunc(x_width, y_width, p_axis, img)

    boxrange = 50
    # data convert for plotting
    x_ = convert(Array{Float32,1}, x_width)
    y_ = convert(Array{Float32,1}, y_width)
    y_box = [Vector{Float64}() for _ = 1:boxrange]
    p_ = [Vector{Float64}() for _ = 1:maximum(p_axis)]

    # group x data
    for i = 1:size(x_, 1)
        push!(p_[p_axis[i]], x_[i])
        push!(y_box[Int(ceil(x_[i]))], y_[i])
    end

    xy_ = zeros(boxrange, 1)
    for i = 1:boxrange
        if size(y_box, 1) != 0
            xy_[i] = mean(y_box[i])
        end
    end

    scene_xy = Scene()
    #GLMakie.scatter!(scene_xy,x_,y_)
    scene_xy = Plots.plot(x_, y_, seriestype = :scatter)
    scene_box = Plots.plot(1:boxrange, xy_)

    xp_ = zeros(maximum(p_axis), 1)
    for i = 1:maximum(p_axis)
        if size(p_[i], 1) != 0
            xp_[i] = Int(round(mean(p_[i])))
        end
    end

    scene_xp = Plots.plot(
        1:maximum(p_axis),
        vec(xp_),
        axis = (limits = (1, maximum(p_axis), 1, 20),),
    )
    return scene_xy, scene_box, scene_xp
end


text:
########################################################
# Bilevel Parameter Learning via Nonsmooth Trust Region (Vectorized Mode)
########################################################

using SparseArrays

using AlgTools
using AlgTools.Util
import AlgTools.Iterate

using VariationalImaging.GradientOps
using VariationalImaging.Util

using LinearAlgebra
using LinearOperators,Krylov

export bilevel_learn

#############
# Data types
#############

ImageSize = Tuple{Integer,Integer}
Image = Array{Float64,2}
Primal = Image
Dual = Array{Float64,3}
Parameter = Union{Real,AbstractArray}
Dataset = Tuple{Array{Float64,3},Array{Float64,3}}

#########################
# Iterate initialisation
#########################

function init_rest(x::Real,learning_function::Function,Δ,ds)
    x̄ = copy(x)
    u,fx,gx = learning_function(x,ds,Δ)
    ū = copy(u)
    fx̄ = copy(fx)
    gx̄ = copy(gx)
    B = 0.1
    return x, x̄, u, ū, fx, gx, fx̄, gx̄, Δ, B
end

function init_rest(x::Union{AbstractArray{T,3},AbstractArray{T,2},AbstractArray{T,1}},learning_function::Function,Δ,ds) where T
    x̄ = copy(x)
    u,fx,gx = learning_function(x,ds,Δ)
    ū = copy(u)
    fx̄ = copy(fx)
    gx̄ = copy(gx)
    B = LBFGSOperator(length(x[:]))
    return x, x̄, u, ū, fx, gx, fx̄, gx̄, Δ, B
end

############
# Auxiliary Functions
############

# Dogleg step calculation using l_\infty ball
function dogleg_box(x::Real,gx,B,Δ)
    lb,ub = get_bounds(x,Δ)
    #@info lb,ub
    pn = B\gx
    if in_bounds(lb,Δ,pn)
        return pn
    end
    p = -(norm(gx)^2/(gx'*(B*gx)))*gx # Cauchy Step
    if in_bounds(lb,Δ,p) == false
        #@info "Scaled"
        t = step_to_bound(p/norm(p),lb,Δ)
        return (p/norm(p))*t
    end
    #@info "Dogleg"
    t = step_to_bound(pn-p,lb,Δ)
    return p + t * (pn-p)
end

function dogbox(x::Real,gx,B,Δ)
    lb,ub = get_bounds(x,Δ)
    #@info lb,ub
    pn = B\gx
    if in_bounds(lb,Δ,pn)
        return pn
    end
    p = -(norm(gx)^2/(gx'*(B*gx)))*gx # Cauchy Step
    if in_bounds(lb,Δ,p) == false
        @info "Dogbox"
        t = step_to_bound(p,lb,Δ)
        psc = p*t
        t2 = step_to_bound(psc-pn,lb,Δ)
        return psc + t2 * (psc-pn)
    end
    #@info "Dogleg"
    t = step_to_bound(pn-p,lb,Δ)
    return p + t * (pn-p)
end

# Dogleg step calculation using l_\infty ball
function dogleg_box(x::Union{AbstractArray{T,3},AbstractArray{T,2},AbstractArray{T,1}},gx,B,Δ) where T
    lb,ub = get_bounds(x,Δ)
    pn = newton_step(B,gx)
    if in_bounds(lb,Δ,pn)
        return pn
    end
    p = cauchy_step(B,gx)
    if in_bounds(lb,Δ,p) == false
        @info "Scaled"
        t = step_to_bound(p/norm₂(p),lb,Δ)
        return (p/norm₂(p)) .*t
    end
    @info "Dogleg"
    t = step_to_bound(pn-p,lb,Δ)
    return p + t .* (pn-p)
end

function dogbox(x::Union{AbstractArray{T,3},AbstractArray{T,2},AbstractArray{T,1}},gx,B,Δ) where T
    lb,ub = get_bounds(x,Δ)
    pn = newton_step(B,gx)
    if in_bounds(lb,Δ,pn)
        return pn
    end
    p = cauchy_step(B,gx)
    if in_bounds(lb,Δ,p) == false
        @info "Dogbox"
        t = step_to_bound(p,lb,Δ)
        psc = t .* p
        t2 = step_to_bound(pn-psc,lb,Δ)
        return psc + t2 .* (pn-psc)
    end
    @info "Dogleg"
    t = step_to_bound(pn-p,lb,Δ)
    return p + t .* (pn-p)
end

function newton_step(B::LinearOperators.LBFGSOperator,gx::AbstractArray)
    pn, ks = Krylov.cg_lanczos(B,-gx[:])
    if ks.solved == false
        @error ks
    end
    return reshape(pn,size(gx)...)
end

function cauchy_step(B::LinearOperators.LBFGSOperator,gx::AbstractArray)
    p = -(norm(gx[:])^2/(gx[:]'*(B*gx[:])))*gx[:]
    return reshape(p,size(gx)...)
end

# Distance to l_\infty bound
function step_to_bound(p,lb,ub)
    dist = max.(lb ./ p, ub ./ p)
    return dist
end

# Test if vector is within the l_\infty ball
function in_bounds(lb,ub,x)
    return all(x .>= lb) && all(x .<= ub)
end

# Get bounds for the l_\infty ball intersected with positive quadrant
function get_bounds(x,Δ)
    lb = max.(-Δ,eps() .-x)
    ub = Δ * ones(size(x))
    return lb,ub
end

function pred(B::LinearOperators.LBFGSOperator,p,gx)
    return -p[:]'*gx[:] -0.5 * p[:]'*(B*p[:])
end

function pred(B::Real,p,gx)
    return -p*gx - 0.5*p*B*p
end

function updateBFGS!(B::LinearOperators.LBFGSOperator,y,s)
    if y[:]'*(B*y[:]) > 0
        push!(B,y[:],s[:])
    end
    return B
end

function updateBFGS!(B::Real,y,s)
    if y'*(B*y) > 0
        B = B + (y*y)/(y*s) - (B*s*s*B)/(s*B*s)
    end
    return B
end

############
# Algorithm
############

function bilevel_learn(ds :: Dataset,
    learning_function::Function;
    xinit :: Parameter,
    iterate = Iterate.simple_iterate,
    params::NamedTuple)

    ################################
    # Extract and set up parameters
    ################################

    η₁, η₂ = params.η₁, params.η₂
    β₁, β₂ =  params.β₁, params.β₂
    Δ₀ = params.Δ₀
    parsize = size(xinit)

    ######################
    # Initialise iterates
    ######################
    x, x̄, u, ū, fx, gx, fx̄, gx̄, Δ, B = init_rest(xinit,learning_function,Δ₀,ds)
    residual = x-x̄

    ####################
    # Run the algorithm
    ####################

    v = iterate(params) do verbose :: Function

        #println("x=$(norm(x,1)/length(x)), Δ=$Δ, gx=$gx")

        p = dogleg_box(x,gx,B,Δ) # solve tr subproblem
        #p = dogbox(x,gx,B,Δ)

        x̄ = x + p  # test new point
        #println(minimum(x̄))

        ū,fx̄,gx̄ = learning_function(x̄,ds,Δ)

        predf = pred(B,p,gx)
        ρ = (fx-fx̄)/predf # ared/pred
        if predf == 0
            @error "Problems with step calculated"
        end
        #println("$fx, $fx̄, $(fx-fx̄), $predf, $p")
        #println("ρ=$ρ, ared = $(fx-fx̄), pred = $(predf)")

        updateBFGS!(B,gx̄-gx,p)

        if ρ < η₁               # radius update
            Δ = β₁*Δ
        elseif ρ > η₂
            if norm₂(p) > 0.8*Δ
                Δ = β₂*Δ
            end
        end

        if predf < 0
            Δ = β₁*Δ
        end

        if ρ > 0
            residual = x-x̄
            x = x̄
            u = ū
            fx = fx̄
            gx = gx̄
        end

        ################################
        # Give function value if needed
        ################################
        v = verbose() do
            x, u[:,:,1], fx, norm₂(gx), Δ, norm₂(residual) # just show the first image on the dataset
            #fx,u[:,:,1]
        end

        v
    end

    return x, u, v
end

max_stars_repo_name: dvillacis/BPLDenoising

max_issues_repo_name: dvillacis/BPLDenoising

max_forks_repo_name: dvillacis/BPLDenoising


text:
using GeoIP
using Documenter

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

makedocs(;
    modules=[GeoIP],
    authors = "Andrey Oskin, Seth Bromberger, contributors: https://github.com/JuliaWeb/GeoIP.jl/graphs/contributors",
    repo="https://github.com/JuliaWeb/GeoIP.jl/blob/{commit}{path}#{line}",
    sitename="GeoIP.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaWeb.github.io/GeoIP.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaWeb/GeoIP.jl",
)


text:
module MocosSimLauncher

using ArgParse
using Base.Threads
using DataFrames
using FileIO
using JLD2
using JSON
using ProgressMeter
using Random

import MocosSim
import MocosSim: ContactKind, NoContact, contactkind, time

const OptTimePoint = Union{Missing, MocosSim.TimePoint}
optreal2float32(optreal::Union{Missing,T} where T<:Real) = ismissing(optreal) ? NaN32 : Float32(optreal)

include("cmd_parsing.jl")
include("callback.jl")
include("load_params.jl")
include("outputs.jl")

export launch

function launch()
  @info "Stated" nthreads()
  if nthreads() == 1
    @warn "using single thread, set more threads by passing --threads agrument to julia or setting JULIA_NUM_THREADS environment variable"
  end

  cmd_args = parse_commandline()
  @info "Parsed args" cmd_args
  json = JSON.parsefile(cmd_args["JSON"])

  max_num_infected = json["stop_simulation_threshold"] |> Int
  num_trajectories = json["num_trajectories"] |> Int
  num_initial_infected = json["initial_conditions"]["cardinalities"]["infectious"] |> Int
  params_seed = get(json, "params_seed", 0)

  @info "loading population and setting up parameters" params_seed
  rng = MersenneTwister(params_seed)
  params = read_params(json, rng)
  num_individuals =  MocosSim.numindividuals(params)

  states = [MocosSim.SimState(num_individuals) for _ in 1:nthreads()]
  callbacks = [DetectionCallback(num_individuals, max_num_infected) for _ in 1:nthreads()]

  outputs = make_outputs(cmd_args, num_trajectories)
  foreach(o->beforetrajectories(o, params), outputs)

  @info "starting simulation" num_trajectories
  writelock = ReentrantLock()
  progress = ProgressMeter.Progress(num_trajectories)
  GC.gc()
  @threads for trajectory_id in 1:num_trajectories
    state = states[threadid()]
    MocosSim.reset!(state, trajectory_id)
    MocosSim.initialfeed!(state, num_initial_infected)

    callback = callbacks[threadid()]
    reset!(callback)
    try
      MocosSim.simulate!(state, params, callback)
      foreach(o->pushtrajectory!(o, trajectory_id, writelock, state, params, callback), outputs)
    catch err
      println(stderr, "Failed on thread ", threadid(), " iteration ", trajectory_id, " failed: ", err)
      foreach(x -> println(stderr, x), stacktrace(catch_backtrace()))
    end

    ProgressMeter.next!(progress) # is thread-safe
  end
  foreach(o->aftertrajectories(o, params), outputs)
end

end

max_stars_repo_name: Novirmor/MocosSimLauncher

max_issues_repo_name: Novirmor/MocosSimLauncher

max_forks_repo_name: Novirmor/MocosSimLauncher


text:
# Copyright 2018 IBM Corp.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

""" Command-line interface for raw and semantic flow graphs.
"""
module CLI
export main, invoke, parse

using ArgParse
import DefaultApplication
using Requires
import JSON
import Serd

using Catlab.WiringDiagrams, Catlab.Graphics
import Catlab.Graphics: Graphviz
using ..RawFlowGraphs
using ..Ontology, ..SemanticEnrichment
using ..Serialization

# CLI arguments
###############

const settings = ArgParseSettings()

@add_arg_table settings begin
  "record"
    help = "record code as raw flow graph"
    action = :command
  "enrich"
    help = "convert raw flow graph to semantic flow graph"
    action = :command
  "visualize"
    help = "visualize flow graph"
    action = :command
  "ontology"
    help = "export ontology"
    action = :command
end

@add_arg_table settings["record"] begin
  "path"
    help = "input Python or R file, or directory"
    required = true
  "-o", "--out"
    help = "output raw flow graph or directory"
  "-f", "--from"
    help = "input format (one of: \"python\", \"r\") (default: from file name)"
    default = nothing
  "-t", "--to"
    help = "output format (one of: \"graphml\", \"json\")"
    default = "graphml"
  "--graph-outputs"
    help = "whether and how to retain outputs of raw flow graph (Python only)"
    default = "none"
    metavar = "GRAPH-OUT"
  "--annotations"
    help = "load annotations from JSON file (default: load from remote database)"
    metavar = "FILE"
end

@add_arg_table settings["enrich"] begin
  "path"
    help = "input raw flow graph or directory"
    required = true
  "-o", "--out"
    help = "output semantic flow graph or directory"
  "-f", "--from"
    help = "input format (one of: \"graphml\", \"json\") (default: from file name)"
    default = nothing
  "-t", "--to"
    help = "output format (one of: \"graphml\", \"json\")"
    default = "graphml"
  "--ontology"
    help = "load ontology from JSON file (default: load from remote database)"
    metavar = "FILE"
end

@add_arg_table settings["visualize"] begin
  "path"
    help = "input raw or semantic flow graph, or directory"
    required = true
  "-o", "--out"
    help = "output file or directory"
  "-f", "--from"
    help = "input format (one of: \"raw-graphml\", \"raw-json\", \"semantic-graphml\", \"semantic-json\") (default: from file name)"
    default = nothing
  "-t", "--to"
    help = "Graphviz output format (default: Graphviz input only)"
  "--open"
    help = "open output using OS default application"
    action = :store_true
end

@add_arg_table settings["ontology"] begin
  "json"
    help = "export ontology as JSON"
    action = :command
  "rdf"
    help = "export ontology as RDF/OWL"
    action = :command
end

@add_arg_table settings["ontology"]["json"] begin
  "-o", "--out"
    help = "output file (default: stdout)"
  "--indent"
    help = "number of spaces to indent (default: compact output)"
    arg_type = Int
    default = nothing
  "--no-concepts"
    help = "exclude concepts from export"
    dest_name = "concepts"
    action = :store_false
  "--no-annotations"
    help = "exclude annotations from export"
    dest_name = "annotations"
    action = :store_false
end

@add_arg_table settings["ontology"]["rdf"] begin
  "-o", "--out"
    help = "output file (default: stdout)"
  "-t", "--to"
    help = "output format (one of: \"turtle\", \"ntriples\", \"nquads\", \"trig\")"
    default = "turtle"
  "--no-concepts"
    help = "omit concepts"
    dest_name = "concepts"
    action = :store_false
  "--no-annotations"
    help = "omit annotations"
    dest_name = "annotations"
    action = :store_false
  "--no-schema"
    help = "omit preamble defining OWL schema for concepts and annotations"
    dest_name = "schema"
    action = :store_false
  "--no-provenance"
    help = "omit interoperability with PROV Ontology (PROV-O)"
    dest_name = "provenance"
    action = :store_false
  "--no-wiring-diagrams"
    help = "omit wiring diagrams in concepts and annotations"
    dest_name = "wiring"
    action = :store_false
end

# Record
########

function record(args::Dict)
  paths = parse_io_args(args["path"], args["out"],
    format = args["from"],
    formats = Dict(
      ".py" => "python",
      ".r" => "r",
      ".R" => "r",
    ),
    out_ext = format -> ".raw." * args["to"]
  )
  for (inpath, lang, outpath) in paths
    if args["to"] == "graphml"
      record_file(abspath(inpath), abspath(outpath), args, Val(Symbol(lang)))
    else
      # Both the Python and R packages emit GraphML.
      diagram = mktemp() do temp_path, io
        record_file(abspath(inpath), temp_path, args, Val(Symbol(lang)))
        read_graph_file(temp_path, format="graphml")
      end
      write_graph_file(diagram, outpath, format=args["to"])
    end
  end
end

function record_file(inpath::String, outpath::String, args::Dict,
                     ::Val{lang}) where lang
  if lang == :python
    error("PyCall.jl has not been imported")
  elseif lang == :r
    error("RCall.jl has not been imported")
  else
    error("Unsupported language: $lang")
  end
end

# Enrich
########

function enrich(args::Dict)
  paths = parse_io_args(args["path"], args["out"],
    format = args["from"],
    formats = Dict(
      ".raw.graphml" => "graphml",
      ".raw.json" => "json",
    ),
    out_ext = format -> ".semantic." * args["to"],
  )

  # Load ontology from local or remote source.
  db = OntologyDB()
  if isnothing(args["ontology"])
    # Only load concepts from remote. Annotations will be loaded on-the-fly.
    load_concepts(db)
  else
    load_ontology_file(db, args["ontology"])
  end

  # Run semantic enrichment.
  for (inpath, format, outpath) in paths
    raw = read_graph_file(inpath, kind="raw", format=format)
    raw = rem_literals!(raw)
    semantic = to_semantic_graph(db, raw)
    write_graph_file(semantic, outpath, format=args["to"])
  end
end

# Visualize
###########

function visualize(args::Dict)
  format = isnothing(args["from"]) ? nothing : Tuple(split(args["from"],"-",1))
  out_format = isnothing(args["to"]) ? "dot" : args["to"]
  paths = parse_io_args(args["path"], args["out"],
    format = format,
    formats = Dict(
      ".raw.graphml" => ("raw", "graphml"),
      ".raw.json" => ("raw", "json"),
      ".semantic.graphml" => ("semantic", "graphml"),
      ".semantic.json" => ("semantic", "json"),
    ),
    out_ext = format -> "." * first(format) * "." * out_format,
  )
  for (inpath, (kind, format), outpath) in paths
    # Read flow graph and convert to Graphviz AST.
    diagram = read_graph_file(inpath, kind=kind, format=format)
    graphviz = if kind == "raw"
      raw_graph_to_graphviz(diagram)
    else
      semantic_graph_to_graphviz(diagram)
    end

    # Pretty-print Graphviz AST to output file.
    if isnothing(args["to"])
      # Default: no output format, yield Graphviz dot input.
      open(outpath, "w") do f
        Graphviz.pprint(f, graphviz)
      end
    else
      # Run Graphviz with given output format.
      open(`dot -T$out_format -o $outpath`, "w", stdout) do f
        Graphviz.pprint(f, graphviz)
      end
    end

    if args["open"]
      DefaultApplication.open(outpath)
    end
  end
end

function raw_graph_to_graphviz(diagram::WiringDiagram)
  to_graphviz(rem_unused_ports(diagram);
    graph_name = "raw_flow_graph",
    labels = true,
    label_attr = :xlabel,
    graph_attrs = Graphviz.Attributes(
      :fontname => "Courier",
    ),
    node_attrs = Graphviz.Attributes(
      :fontname => "Courier",
    ),
    edge_attrs = Graphviz.Attributes(
      :fontname => "Courier",
      :arrowsize => "0.5",
    )
  )
end

function semantic_graph_to_graphviz(diagram::WiringDiagram)
  to_graphviz(diagram;
    graph_name = "semantic_flow_graph",
    labels = true,
    label_attr = :xlabel,
    graph_attrs=Graphviz.Attributes(
      :fontname => "Helvetica",
    ),
    node_attrs=Graphviz.Attributes(
      :fontname => "Helvetica",
    ),
    edge_attrs=Graphviz.Attributes(
      :fontname => "Helvetica",
      :arrowsize => "0.5",
    ),
    cell_attrs=Graphviz.Attributes(
      :style => "rounded",
    )
  )
end

# Ontology
##########

function ontology_as_json(args::Dict)
  docs = AbstractDict[]
  db = OntologyDB()
  if args["concepts"]
    append!(docs, OntologyDBs.api_get(db, "/concepts"))
  end
  if args["annotations"]
    append!(docs, OntologyDBs.api_get(db, "/annotations"))
  end
  if args["out"] != nothing
    open(args["out"], "w") do out
      JSON.print(out, docs, args["indent"])
    end
  else
    JSON.print(stdout, docs, args["indent"])
  end
end

function ontology_as_rdf(args::Dict)
  # Load ontology data from remote database.
  db = OntologyDB()
  if args["concepts"]
    load_concepts(db)
  end
  if args["annotations"]
    load_annotations(db)
  end

  # Load ontology schemas from filesystem.
  stmts = Serd.RDF.Statement[]
  if args["schema"]
    append!(stmts, [
      read_ontology_rdf_schema("list.ttl");
      args["concepts"] ? read_ontology_rdf_schema("concept.ttl") : [];
      args["annotations"] ? read_ontology_rdf_schema("annotation.ttl") : [];
      args["wiring"] ? read_ontology_rdf_schema("wiring.ttl") : [];
    ])
  end

  # Convert to RDF.
  prefix = Serd.RDF.Prefix("dso", "https://www.datascienceontology.org/ns/dso/")
  append!(stmts, ontology_to_rdf(db, prefix,
    include_provenance=args["provenance"],
    include_wiring_diagrams=args["wiring"]))

  # Serialize RDF to file or stdout.
  syntax = args["to"]
  if args["out"] != nothing
    open(args["out"], "w") do out
      Serd.write_rdf(out, stmts, syntax=syntax)
    end
  else
    Serd.write_rdf(stdout, stmts, syntax=syntax)
  end
end

function read_ontology_rdf_schema(name::String)
  Serd.read_rdf_file(joinpath(ontology_rdf_schema_dir, name))
end

const ontology_rdf_schema_dir = joinpath(@__DIR__, "ontology", "rdf", "schema")

# CLI main
##########

function main(args)
  invoke(parse(args)...)
end

function parse(args)
  cmds = String[]
  parsed_args = parse_args(args, settings)
  while haskey(parsed_args, "%COMMAND%")
    cmd = parsed_args["%COMMAND%"]
    parsed_args = parsed_args[cmd]
    push!(cmds, cmd)
  end
  return (cmds, parsed_args)
end

function invoke(cmds, cmd_args)
  try
    cmd_fun = command_table
    for cmd in cmds; cmd_fun = cmd_fun[cmd] end
    cmd_fun(cmd_args)
  catch err
    # Handle further "parsing" errors ala ArgParse.jl.
    isa(err, ArgParseError) || rethrow()
    settings.exc_handler(settings, err)
  end
end

const command_table = Dict(
  "record" => record,
  "enrich" => enrich,
  "visualize" => visualize,
  "ontology" => Dict(
    "json" => ontology_as_json,
    "rdf" => ontology_as_rdf,
  ),
)

# CLI extras
############

function __init__()
  @require PyCall="438e738f-606a-5dbb-bf0a-cddfbfd45ab0" include("extras/CLI-Python.jl")
  @require RCall="6f49c342-dc21-5d91-9882-a32aef131414" include("extras/CLI-R.jl")
end

# CLI utilities
###############

""" Map CLI input/output arguments to pairs of input/output files.
"""
function parse_io_args(input::String, output::Union{String,Nothing};
                       format::Any=nothing, formats::AbstractDict=Dict(),
                       out_ext::Function=format->"")
  function get_format_and_output(input, output=nothing)
    if isnothing(output) || isnothing(format)
      name, ext = match_ext(input, keys(formats))
      inferred_format = formats[ext]
      inferred_output = name * out_ext(inferred_format)
    end
    (isnothing(format) ? inferred_format : format,
     isnothing(output) ? inferred_output : output)
  end

  if isdir(input)
    inexts = collect(keys(formats))
    names = filter(name -> any(endswith.(name, inexts)), readdir(input))
    outdir = if isnothing(output); input
      elseif isdir(output); output
      else; throw(ArgParseError(
        "Output must be directory when input is directory"))
      end
    map(names) do name
      format, output = get_format_and_output(name)
      (joinpath(input, name), format, joinpath(outdir, output))
    end
  elseif isfile(input)
    format, output = get_format_and_output(input, output)
    [ (input, format, output) ]
  else
    throw(ArgParseError("Input must be file or directory"))
  end
end

function match_ext(name::String, exts)
  # Don't use splitext because we allow extensions with multiple dots.
  for ext in exts
    if endswith(name, ext)
      return (name[1:end-length(ext)], ext)
    end
  end
  throw(ArgParseError("Cannot match extension in filename: \"$name\""))
end

function read_graph_file(filename::String;
    kind::Union{String,Nothing}=nothing, format::String="graphml")
  reader = get(graph_readers, (kind, format)) do
    error("Invalid graph kind \"$kind\" or format \"$format\"")
  end
  reader(filename)
end

function write_graph_file(diagram::WiringDiagram, filename::String;
    format::String="graphml")
  writer = get(graph_writers, format) do
    error("Invalid graph format \"$format\"")
  end
  writer(diagram, filename)
end

const graph_readers = Dict(
  (nothing, "graphml") => x -> read_graphml(Dict, Dict, Dict, x),
  (nothing, "json") => x -> read_json_graph(Dict, Dict, Dict, x),
  ("raw", "graphml") => read_raw_graphml,
  ("raw", "json") => read_raw_graph_json,
  ("semantic", "graphml") => read_semantic_graphml,
  ("semantic", "json") => read_semantic_graph_json,
)

const graph_writers = Dict(
  "graphml" => write_graphml,
  "json" => (x, filename) -> write_json_graph(x, filename, indent=2),
)

end

max_stars_repo_name: cclauss/semanticflowgraph

max_issues_repo_name: cclauss/semanticflowgraph

max_forks_repo_name: cclauss/semanticflowgraph


text:
# BioSequence
# ===========
#
# A general purpose biological sequence representation.
#
# It is cheap to create a subsequence from a sequence because sequences can
# share underlying data: creating a subsequence is just allocating a new
# sequence object defining a part of the original sequence without copy.
# Destructive operations create a copy of the underlying data if and only if it
# is shared between (sub)sequences. This is often called as copy-on-write
# strategy in computer science and should be transparent to the user.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

# About Internals
# ---------------
#
# The `data` field of a `BioSequence{A}` object contains binary representation
# of a biological character sequence. Each character is encoded with an encoder
# corresponding to the alphabet `A` and compactly packed into `data`. To extract
# a character from a sequence, you should decode this binary sequence with a
# decoder that is a pair of the encoder. The length of encoded binary bits is
# fixed, and hence a character at arbitrary position can be extracted in a
# constant time. To know the exact location of a character at a position, you
# can use the `bitsof(seq, i)` function, which returns a pair of element's index
# containing binary bits and bits' offset. As a whole, character extraction
# `seq[i]` can be written as:
#
#     j = bitindex(seq, i)
#     decode(A, (seq.data[index(j)] >> offset(j)) & mask(A))
#
#  index :        index(j) - 1       index(j)       index(j) + 1
#   data : ....|xxxxx...........|xxXxxxxxxxxxxxxx|............xxxx|....
# offset :                          |<-offset(j)-|
#  width :      |<---- 64 ---->| |<---- 64 ---->| |<---- 64 ---->|
#
#  * '.' : unused (4 bits/char)
#  * 'x' : used
#  * 'X' : used and pointed by index `i`

"""
Biological sequence data structure indexed by an alphabet type `A`.
"""
mutable struct BioSequence{A<:Alphabet} <: Sequence
    data::Vector{UInt64}  # encoded character sequence data
    part::UnitRange{Int}  # interval within `data` defining the (sub)sequence
    shared::Bool          # true if and only if `data` is shared between sequences

    function BioSequence{A}(data::Vector{UInt64},
                            part::UnitRange{Int},
                            shared::Bool) where A
        return new(data, part, shared)
    end
end

const DNASequence       = BioSequence{DNAAlphabet{4}}
const RNASequence       = BioSequence{RNAAlphabet{4}}
const AminoAcidSequence = BioSequence{AminoAcidAlphabet}
const CharSequence      = BioSequence{CharAlphabet}

"Gets the alphabet encoding of a given BioSequence."
BioSymbols.alphabet(::Type{BioSequence{A}}) where {A} = alphabet(A)

Base.length(seq::BioSequence) = length(seq.part)
Base.eltype(::Type{BioSequence{A}}) where {A} = eltype(A)

function seq_data_len(::Type{A}, len::Integer) where {A}
    return cld(len, div(64, bitsof(A)))
end

# Replace a BioSequence's data with a copy, copying only what's needed.
# The user should never need to call this, as it has no outward effect on the
# sequence.
function orphan!(seq::BioSequence{A}, size::Integer=length(seq), force::Bool=false) where {A}
    if !seq.shared && !force
        return seq
    end

    j, r = bitindex(seq, 1)
    data = Vector{UInt64}(undef, seq_data_len(A, size))

    if !isempty(seq) && !isempty(data)
        x = seq.data[j] >> r
        m = index(bitindex(seq, lastindex(seq))) - j + 1
        l = min(lastindex(data), m)
        @inbounds @simd for i in 1:l-1
            y = seq.data[j + i]
            data[i] = x | y << (64 - r)
            x = y >> r
        end
        if m <= l
            data[l] = x
        else
            y = seq.data[j + l]
            data[l] = x | y << (64 - r)
        end
    end

    seq.data = data
    seq.part = 1:length(seq)
    seq.shared = false
    return seq
end

function enc64(::BioSequence{A}, x) where {A}
    return UInt64(encode(A, convert(eltype(A), x)))
end

# Summaries
# ---------

Base.summary(seq::BioSequence{<:DNAAlphabet}) = string(length(seq), "nt ", "DNA Sequence")
Base.summary(seq::BioSequence{<:RNAAlphabet}) = string(length(seq), "nt ", "RNA Sequence")
Base.summary(seq::AminoAcidSequence) = string(length(seq), "aa ", "Amino Acid Sequence")
Base.summary(seq::CharSequence) = string(length(seq), "char ", "Char Sequence")

include("indexing.jl")
include("constructors.jl")
include("iteration.jl")
include("copying.jl")
include("conversion.jl")
include("stringliterals.jl")
include("transformations.jl")
include("operators.jl")


text:
module NQCDynamics

using Reexport: @reexport

@reexport using NQCBase
@reexport using NQCModels

include("RingPolymers/RingPolymers.jl")
@reexport using .RingPolymers
@reexport using RingPolymerArrays

include("Calculators/Calculators.jl")

include("simulations.jl")
export Simulation,
       RingPolymerSimulation,
       natoms,
       masses

@reexport using NQCDistributions

include("DynamicsUtils/DynamicsUtils.jl")
@reexport using .DynamicsUtils: get_positions, get_velocities
export DynamicsUtils

include("Estimators.jl")
export Estimators

include("DynamicsOutputs.jl")

include("TimeCorrelationFunctions.jl")
export TimeCorrelationFunctions

include("DynamicsMethods/DynamicsMethods.jl")
@reexport using .DynamicsMethods

include("InitialConditions/InitialConditions.jl")
export InitialConditions

include("Ensembles/Ensembles.jl")
export Ensembles
@reexport using .Ensembles

end # module


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

@userplot VarPlane

@recipe function f(vp::VarPlane;
                   theta=0, phi=90,
                   ptol=1e-6, dtol=0.5,
                   nangs=50, nlags=20, maxlag=nothing,
                   showrange=true, varmodel=GaussianVariogram)
  # sanity checks
  @assert 0 ≤ theta ≤ 360 "theta must lie in [0,360]"
  @assert -90 ≤ phi ≤ 90 "phi must lie in [-90,90]"

  # get inputs
  spatialdata = vp.args[1]
  var₁ = vp.args[2]
  var₂ = length(vp.args) == 3 ? vp.args[3] : var₁
  N = ndims(spatialdata)

  # basis for variogram plane
  if N == 2
    planes = [spatialdata]
    u, v = (1.,0.), (0.,1.)
  elseif N == 3
    θ = deg2rad(theta); ϕ = deg2rad(phi)
    n = (cos(ϕ)cos(θ), cos(ϕ)sin(θ), sin(ϕ))

    planes = partition(spatialdata, PlanePartitioner(n, tol=ptol))
    u, v = plane_basis(n)
  else
    @error "variogram plane only supported in 2D or 3D"
  end

  # polar coordinates and variogram values
  rs = Vector{Float64}(undef, nlags)
  θs = range(0, stop=π, length=nangs)
  zs = Matrix{Float64}(undef, nlags, nangs)

  # ranges for each angle
  ls = Vector{Float64}(undef, nangs)

  # loop over half of the plane
  for (j, θ) in Iterators.enumerate(θs)
    dir = ntuple(i -> cos(θ)*u[i] + sin(θ)*v[i], N)

    dpart = DirectionPartitioner(dir, tol=dtol)

    # compute directional variogram across planes
    plane, _ = iterate(planes)
    γ = EmpiricalVariogram(partition(plane, dpart), var₁, var₂;
                           nlags=nlags, maxlag=maxlag)
    for plane in Iterators.drop(planes, 1)
      γplane = EmpiricalVariogram(partition(plane, dpart), var₁, var₂;
                                  nlags=nlags, maxlag=maxlag)
      merge!(γ, γplane)
    end

    rs, ys, _ = values(γ)

    # clean NaN values (i.e. empty bins)
    isnan(ys[1]) && (ys[1] = 0)
    for i in 2:length(ys)
      isnan(ys[i]) && (ys[i] = ys[i-1])
    end

    # save variogram for this angle
    zs[:,j] = ys

    if showrange
      γtheo = fit(varmodel, γ)
      ls[j] = range(γtheo)
    end
  end

  # exploit symmetry
  θs = range(0, stop=2π, length=2nangs)
  zs = hcat(zs, zs)

  # show variables in plot title
  ptitle = var₁ == var₂ ? "$var₁" : "$var₁-$var₂"

  projection --> :polar
  ylim --> (0, Inf)
  title --> ptitle

  # plot variogram plane
  @series begin
    seriestype --> :heatmap
    color --> :curl

    θs, rs, zs
  end

  # plot ranges
  if showrange
    @series begin
      seriestype --> :path
      linecolor --> :black
      label --> "range"

      θs, ls
    end
  end
end

function plane_basis(normal::NTuple{3,T}) where {T<:Real}
  # normalize input
  n = normal ./ sqrt(sum(normal[i]^2 for i in 1:3))

  # find last non-zero component
  idx = -1
  for (i, c) in enumerate(reverse(n))
    if c != 0
      idx = length(n) - i + 1
      break
    end
  end

  @assert idx > 0 "invalid normal vector"

  # first basis vector (perturb and subtract projection)
  u = ntuple(i -> i == idx%3 + 1 ? n[i] + one(T) : n[i], 3)
  l = sum(u[i]*n[i] for i in 1:3)
  u = ntuple(i -> u[i] - l*n[i], 3)

  # second basis vector (cross product)
  nx, ny, nz = n
  ux, uy, uz = u
  v = (ny*uz - nz*uz, nz*ux - nx*uz, nx*uy - ny*ux)

  # normalize output
  u = u ./ sqrt(sum(u[i]^2 for i in 1:3))
  v = v ./ sqrt(sum(v[i]^2 for i in 1:3))

  u, v
end


text:
 # script to plot model predictions
#
# The script can be called in multiple ways:
#
# (i) discover the model parameters
#
# plotmodel.jl [model]
#
# (ii) plot model fits for previously fitted model
#
# plotmodel.jl [model] [objective]
#
# (iii) plot model fits for given parameter values
#
# plotmodel.jl [model] [p1] [p2] ...
#
# The plots are written to figs/[model]_[objective]_[categ/ident].pdf

using Printf, Random, CSV
import Cairo, Fontconfig

include("common.jl")
include("simmodel.jl")
include("behaviorstats.jl")
include("models.jl")

const nrats = 1
const t = Theme(minor_label_font_size=6pt, major_label_font_size=8pt)
const fitpath = "nrats1_0.1tndvar_fits"
const datapath = "../fitdata/$(fitpath)"
const figpath = "../fitdata/$(fitpath)"

# returns model for given name
function createmodel(modelname)
    modeldict = Dict("opt" => OptimalModel,
                     "optlapse" => OptimalLapseModel,
                     "optcollapse" => OptimalColLapseModel,
                     "optbiaslapse" => OptimalBiasLapseModel,
                     "optcolbiaslapse" => OptimalColBiasLapseModel,
                     "adflapse" => ADFLapseModel,
                     "adfcollapse" => ADFColLapseModel,
                     "adfbiaslapse" => ADFBiasLapseModel,
                     "adfcolbiaslapse" => ADFColBiasLapseModel,
                     "randomcolbiaslapse" => RandomColBiasLapseModel,
                     "randomcollapse" => RandomColLapseModel,
                     "gammalapse" => GammaLapseModel,
                     "gammacollapse" => GammaColLapseModel,
                     "gammabiaslapse" => GammaBiasLapseModel,
                     "gammacolbiaslapse" => GammaColBiasLapseModel,
                     "deltalapse" => DeltaLapseModel,
                     "deltacollapse" => DeltaColLapseModel,
                     "deltabiaslapse" => DeltaBiasLapseModel,
                     "deltacolbiaslapse" => DeltaColBiasLapseModel,
                     "deltaexpcolbiaslapse" => DeltaExpColBiasLapseModel,
                     "lcalapse" => LCALapseModel,
                     "lcacollapse" => LCAColLapseModel)
    if modelname in keys(modeldict)
        return modeldict[modelname]()
    else
        error("Uknown model identifier $modelname")
    end
end

# plots given dataset
function plotdataset(dataset::Symbol, m_behav, stimx, magctrst,
                     stimlabel, ctrstlabel, ttitle,
                     plotfilename)
    d = load_data(dataset)
    # compute stats
    d_bystim, d_bymag = getpsychchronstats(d, nrats)
    d_seq, d_seqfit, d_biasid, d_biasmag = getseqstats(d, nrats)
    m_bystim, m_bymag = getpsychchronstats(m_behav)
    m_seq, m_seqfit, m_biasid, m_biasmag = getseqstats(m_behav)

    # generate plots
    pchronid, ppsychid = plotchronpsychidm(d_bystim, m_bystim, stimx, stimlabel, ttitle, t)
    pchronmag, ppsychmag = plotchronpsychmagm(d_bymag, m_bymag, magctrst, ctrstlabel, ttitle, t)
    pcondpsych = plotcondpsychm(d_seq, m_seq, m_seqfit, stimx, "prev. $stimlabel", ttitle, t)
    pbiasid = plotbiasidm(d_biasid, m_biasid, stimx, "prev. $stimlabel", ttitle, t)
    pbiasmag = plotbiasmagm(d_biasmag, m_biasmag, magctrst, "prev. $ctrstlabel", ttitle, t)
    # plot all of those onto single "sheet"
    p = vstack(
        hstack(vstack(pchronid, ppsychid), vstack(pchronmag, ppsychmag)),
        vstack(
            hstack(pcondpsych[1], pcondpsych[2]),
            hstack(pcondpsych[3], pcondpsych[4])
        ),
        hstack(pbiasid, pbiasmag)
    )
    # output to file
    draw(PDF(plotfilename, 8.5inch, 11inch), p)

end

# process arguments
if isempty(ARGS)
    error("Plotting script needs to be called with at least on argument")
end
const modelname = lowercase(ARGS[1])
const m = createmodel(modelname)
const ϕns = ϕnames(m)
const ϕn = length(ϕns)
if length(ARGS) == 1
    println("Model $modelname has $ϕn parameters:")
    println(foldl((a, b) -> a * ", " * b, [string(i) for i in ϕns]))
    quit()

elseif length(ARGS) == 2
    basefilename = "$(modelname)_$(ARGS[2])"
    fitfilename = "$(datapath)/fit_$(basefilename).csv"
    if !isfile(fitfilename)
        error("Fitting data $fitfilename not found")
    end
    println("Loading fitting data from $fitfilename")
    d = CSV.read(fitfilename)
    (maxllh, maxi) = findmax(d[:llh])
    if maxi != size(d,1) && maxllh > d[end,end]
        warn("Maximum llh $maxllh found in line $maxi instead of last line")
    end
    ϕ = Vector(d[maxi,1:end-1])
else
    basefilename = modelname
    if length(ARGS) != ϕn+1
        error("Expected $(ϕn+1) parameters for model $modelname")
    end
    # process and print arguments
    ϕ = [parse(Float64, s) for s in ARGS[2:end]]
end
println("Plotting model $modelname with parameters")
println(foldl((a, b) -> a * ", " * b,
        [@sprintf("%s=%.2f", ϕns[i], ϕ[i]) for i in 1:ϕn]))

# simulate model
Random.seed!(0)
const evidence, ids = catevidence(ident_cs, categ_cs, 10, 10^5);
const m_ident, m_categ = simmodel(m, ϕ, evidence, ids)

# plot both conditions
const categ_x = [0, 0.2, 0.32, 0.44, 0.56, 0.68, 0.8, 1]
const categ_ctrst = [0, 0, 0, 0, 0.12, 0.36, 0.6, 1.0]
plotdataset(:categ, m_categ, categ_x, categ_ctrst,
            "odor A (frac)", "concentration contrast", "categorization",
            "$(figpath)/$(basefilename)_categ.pdf")
println("Plot written to $(figpath)/$(basefilename)_categ.pdf")

const ident_x = [-4, -3, -2, -1, 1, 2, 3, 4]
const ident_ctrst = [-5.0, -5, -5, -5, -4, -3, -2, -1]
plotdataset(:ident, m_ident, ident_x, ident_ctrst,
            "concentration (10^x)", "concentration contrast", "identification",
            "$(figpath)/$(basefilename)_ident.pdf")
println("Plot written to $(figpath)/$(basefilename)_ident.pdf")


text:
using Media
import REPL

using Logging: with_logger
using .Progress: JunoProgressLogger

ends_with_semicolon(x) = REPL.ends_with_semicolon(split(x,'\n',keepempty = false)[end])

withpath(f, path) =
  CodeTools.withpath(f, path == nothing || isuntitled(path) ? nothing : path)

const eval_channel_in = Channel(1)
const eval_channel_out = Channel(1)
const eval_backend_task = Ref{Any}(nothing)
const is_backend_working = Ref{Bool}(false)

function is_evaling()
  return is_backend_working[] || inREPL[]
end

function run_with_backend(f, args...)
  put!(eval_channel_in, (f, args))
  r = take!(eval_channel_out)
  if r isa EvalError || r isa ErrorException
    render′(Editor(), r)
  else
    r
  end
end

function start_eval_backend()
  global eval_backend_task[] = @async begin
    Base.sigatomic_begin()
    while true
      try
        f, args = take!(eval_channel_in)
        Base.sigatomic_end()
        is_backend_working[] = true
        res = @errs Base.invokelatest(f, args...)
        is_backend_working[] = false
        Base.sigatomic_begin()
        put!(eval_channel_out, res)
      catch err
        put!(eval_channel_out, err)
      finally
        is_backend_working[] = false
      end
    end
    Base.sigatomic_end()
  end
end

handle("evalshow") do data
  @destruct [
    text,
    line,
    path,
    mod
  ] = data

  run_with_backend(evalshow, text, line, path, mod)
  nothing
end

function evalshow(text, line, path, mod)
  fixjunodisplays()
  @dynamic let Media.input = Editor()
    mod = getmodule(mod)

    result = hideprompt() do
      with_logger(JunoProgressLogger()) do
        withpath(path) do
          res = @errs include_string(mod, text, path, line)

          Base.invokelatest() do
            if res isa EvalError
              Base.showerror(IOContext(stderr, :limit => true), res)
            elseif res !== nothing && !ends_with_semicolon(text)
              display(res)
            end
          end

          res
        end
      end
    end

    Base.invokelatest() do
      display = Media.getdisplay(typeof(result), Media.pool(Editor()), default = Editor())
      !isa(result, EvalError) && ends_with_semicolon(text) && (result = nothing)
      if display ≠ Editor() && result !== nothing
        @ierrs render(display, result)
      end
    end
  end
end

handle("eval") do data
  @destruct [
    text,
    line,
    path,
    mod,
    errorinrepl = :errorInRepl || false
  ] = data
  run_with_backend(eval, text, line, path, mod, errorinrepl)
end

function eval(text, line, path, mod, errorinrepl = false)
  fixjunodisplays()
  @dynamic let Media.input = Editor()
    mod = getmodule(mod)

    result = hideprompt() do
      with_logger(JunoProgressLogger()) do
        withpath(path) do
          res = @errs include_string(mod, text, path, line)
          if errorinrepl && res isa EvalError
            Base.invokelatest() do
              try
                Base.showerror(IOContext(stderr, :limit => true), res)
              catch err
                show(stderr, err)
              end
            end
          end
          return res
        end
      end
    end

    Base.invokelatest() do
      !isa(result, EvalError) && ends_with_semicolon(text) && (result = nothing)
      @ierrs displayandrender(result)
    end
  end
end

# dummy handler for Revise compat
handle("evalrepl") do data
  @warn("""
    Juno's evalrepl handler is deprecated.
    """, _id="evalrepl", maxlog=1)
end

handle("evalall") do data
  @destruct [
    code,
    mod = :module || nothing,
    path || "untitled"
  ] = data
  run_with_backend(evalall, code, mod, path)
  nothing
end

function evalall(code, mod = nothing, path = "untitled")
  fixjunodisplays()
  @dynamic let Media.input = Editor()
    mod = if mod !== nothing
       getmodule(mod)
    elseif isabspath(path)
      getmodule(CodeTools.filemodule(path))
    else
      Main
    end

    hideprompt() do
      with_logger(JunoProgressLogger()) do
        withpath(path) do
          result = nothing
          try
            result = include_string(mod, code, path)
          catch err
            ee = EvalError(err, stacktrace(catch_backtrace()))

            # show error in REPL:
            Base.invokelatest() do
              Base.showerror(IOContext(stderr, :limit => true), ee)
              # show notification (if enabled in Atom):
              @msg error(d(:msg => "Error evaluating $(basename(path))",
                           :detail => string(ee),
                           :dismissable => true))
            end
          end

          Base.invokelatest() do
            @ierrs displayandrender(result)
          end
        end
      end
    end
  end
end

handle("docs") do data
  @destruct [
    word,
    mod || "Main"
  ] = data
  docs(word, mod)
end

function docs(word, mod = "Main")
  docstring = @errs getdocs(mod, word)
  docstring isa EvalError && return Dict(:error => true)

  mtable = try getmethods(mod, word)
    catch e
      []
    end

  Dict(
    :error    => false,
    :type     => :dom,
    :tag      => :div,
    :contents =>  map(x -> render(Inline(), x), [docstring; mtable])
  )
end


text:
#! /usr/local/bin

# comments are allowed before frontmatter

# empty lines are allowed, too

#---
#title: Custom Title
#id: custom_id
#---

# This is the content


text:
#See Philippe, *Simulation of right and left truncated gamma
# distributions by mixtures* and Dagpunar *Sampling of variates from
# a truncated gamma distribution*
# for pdf, cdf and quantile only right truncation is supported, left bound is 0.0
# for rand both left and right truncation is supported

TruncatedGamma(α::Float64, θ::Float64, lower::Float64, upper::Float64) =
    Distributions.Truncated(Distributions.Gamma(α, θ), lower, upper)
TruncatedGamma(α::Float64, θ::Float64, lower::Real, upper::Real) =
    TruncatedGamma(Float64(α), Float64(θ), lower, Float64(upper))

# ### statistics

function Distributions.minimum(
    d::Distributions.Truncated{Distributions.Gamma{T},Distributions.Continuous},
) where {T<:Real}
    min(d.lower, Distributions.minimum(d.untruncated))
end
function Distributions.maximum(
    d::Distributions.Truncated{Distributions.Gamma{T},Distributions.Continuous},
) where {T<:Real}
    max(d.lower, Distributions.maximum(d.untruncated))
end

# density of the right truncated gamma distribution
function Distributions.pdf(
    d::Distributions.Truncated{Distributions.Gamma{T},Distributions.Continuous},
    x::Real,
) where {T<:Real}
    if minimum(d.untruncated) == d.lower && maximum(d.untruncated) == d.upper
        Distributions.pdf(Distributions.Gamma(d.untruncated.α, d.untruncated.θ), x)
    elseif minimum(d.untruncated) == d.lower
        # right truncated
        return Distributions.pdf(Distributions.Gamma(d.untruncated.α, d.untruncated.θ), x) /
               Distributions.cdf(
            Distributions.Gamma(d.untruncated.α, d.untruncated.θ),
            d.upper,
        )
    else
        # left truncated
        println("left truncation still not supported.")
        return nothing
    end
end

function Distributions.cdf(
    d::Distributions.Truncated{Distributions.Gamma{T},Distributions.Continuous},
    x::Real,
) where {T<:Real}
    if minimum(d.untruncated) == d.lower && maximum(d.untruncated) == d.upper
        Distributions.cdf(Distributions.Gamma(d.untruncated.α, d.untruncated.θ), x)
    elseif minimum(d.untruncated) == d.lower
        # right truncated
        if (x > d.upper)
            x = d.upper
        end
        return Distributions.cdf(d.untruncated, x) /
               Distributions.cdf(d.untruncated, d.upper)
    else
        # left truncated
        println("left truncation still not supported.")
        return nothing
    end
end

function Distributions.quantile(
    d::Distributions.Truncated{Distributions.Gamma{T},Distributions.Continuous},
    x::Real,
) where {T<:Real}
    if minimum(d.untruncated) == d.lower && maximum(d.untruncated) == d.upper
        Distributions.quantile(Distributions.Gamma(d.untruncated.α, d.untruncated.θ), x)
    elseif minimum(d.untruncated) == d.lower
        # right truncated
        if (d.untruncated.α == 0.0)
            return 0.0
        end

        val = x * Distributions.cdf(d.untruncated, d.upper)
        return Distributions.quantile(d.untruncated, val)
    else
        # left truncated
        println("left truncation still not supported.")
        return nothing
    end
end

#Sampling

# optimal number of components for p = 0.95 fixed
function ncomp_optimal(θ::Float64)
    q = 1.644853626951
    ans = 0.25 * (q * sqrt(q * q + 4.0 * θ)^2)
    return trunc(Int64, ans)::Int64
end

#random number generation from the gamma with right truncation point t = 1,
# i.e. TG^-(α,θ,0.0,1.0)
function rng_right_tgamma_t1(
    rng::Distributions.AbstractRNG,
    d::Distributions.Truncated{Distributions.Gamma{T},Distributions.Continuous},
) where {T<:Real}
    n = ncomp_optimal(d.untruncated.θ)
    wl = zeros(Float64, n + 2)
    wlc = similar(wl)

    wl[1] = 1.0
    wlc[1] = 1.0
    for i = 2:(n+1)
        wl[i] = wl[i-1] * d.untruncated.θ / (d.untruncated.α + i)
        wlc[i] = wlc[i-1] + wl[i]
    end
    wlc = wlc ./ wlc[end-1]
    y = 1.0
    yy = 1.0
    for i = 2:(n+1)
        yy *= d.untruncated.θ / i
        y += yy
    end
    cont = true
    z = 0.0
    while (cont)
        u = Distributions.rand(rng)
        j = 1
        while (u > wlc[j])
            j += 1
            x = Distributions.rand(rng, Distributions.Beta(d.untruncated.α, j))
            u = Distributions.rand(rng)
            z = 1.0
            zz = 1.0
            for i = 2:(n+1)
                zz *= (1 - x) * d.untruncated.θ / i
                z += zz
            end
            z = _exp_c(-d.untruncated.θ * x) * y / z
            if (u <= z)
                cont = false
            end
        end
    end
    return z::Float64
end

function ltgamma(
    rng::Distributions.AbstractRNG,
    d::Distributions.Truncated{Distributions.Gamma{T},Distributions.Continuous},
) where {T<:Real}
    trunc = d.lower
    shape = d.untruncated.α
    b = 1 / d.untruncated.θ * trunc

    if (shape == 1)
        return Random.randexp(rng) / (1 / d.untruncated.θ) + trunc
    end

    d1 = b - shape
    d3 = shape - 1
    c0 = 0.5 * (d1 + sqrt(d1 * d1 + 4 * b)) / b

    x = 0.0
    accept = false
    while (!accept)
        x = b + Random.randexp(rng) / c0
        u = Distributions.rand(rng)
        l_rho = d3 * _log_c(x) - x * (1 - c0)
        l_M = d3 * _log_c(d3 / (1 - c0)) - d3
        accept = _log_c(u) <= (l_rho - l_M)
    end

    return trunc * (x / b)
end

# random number generation from the gamma with right truncation point#
# i.e. TG^-(α,θ,0.0,upper)
function Distributions.rand(
    rng::Distributions.AbstractRNG,
    d::Distributions.Truncated{Distributions.Gamma{T},Distributions.Continuous},
) where {T<:Real}
    if minimum(d.untruncated) == d.lower && maximum(d.untruncated) == d.upper
        Distributions.rand(rng, Distributions.Gamma(d.untruncated.α, d.untruncated.θ))
    elseif minimum(d.untruncated) == d.lower
        # right truncated
        d1 = TruncatedGamma(d.untruncate.α, d.untruncated.θ * d.upper, 0.0, d.upper)
        return rng_right_tgamma_t1(rng, d1) * d.upper
    else
        # left truncated
        return ltgamma(rng, d)
    end
end


text:
#Unscented Kalman Inversion: specific structures and function definitions

"""
    Unscented{FT<:AbstractFloat, IT<:Int} <: Process

An unscented Kalman Inversion process.
"""
mutable struct Unscented{FT <: AbstractFloat, IT <: Int} <: Process
    "a vector of arrays of size `N_parameters` containing the mean of the parameters (in each `uki` iteration a new array of mean is added)"
    u_mean::Vector{Vector{FT}}
    "a vector of arrays of size (`N_parameters x N_parameters`) containing the covariance of the parameters (in each `uki` iteration a new array of `cov` is added)"
    uu_cov::Vector{Matrix{FT}}
    "a vector of arrays of size `N_y` containing the predicted observation (in each `uki` iteration a new array of predicted observation is added)"
    obs_pred::Vector{Vector{FT}}
    "weights in UKI"
    c_weights::Vector{FT}
    mean_weights::Vector{FT}
    cov_weights::Vector{FT}
    "covariance of the artificial evolution error"
    Σ_ω::Matrix{FT}
    "covariance of the artificial observation error"
    Σ_ν_scale::FT
    "regularization parameter"
    α_reg::FT
    "regularization vector"
    r::Vector{FT}
    "update frequency"
    update_freq::IT
    "current iteration number"
    iter::IT
end

# outer constructors
function EnsembleKalmanProcess(
    obs_mean::Vector{FT},
    obs_noise_cov::Matrix{FT},
    process::Unscented{FT, IT};
    Δt = FT(1),
) where {FT <: AbstractFloat, IT <: Int}

    #initial parameters stored as columns
    init_params = [DataContainer(update_ensemble_prediction!(process), data_are_columns = true)]
    # Number of parameter
    N_u = length(process.u_mean[1])
    # ensemble size
    N_ens = 2N_u + 1 #stored with data as columns
    #store for model evaluations
    g = []
    # error store
    err = FT[]
    # timestep store
    Δt = Array([Δt])

    EnsembleKalmanProcess{FT, IT, Unscented}(init_params, obs_mean, obs_noise_cov, N_ens, g, err, Δt, process)
end

"""
    Unscented(
        u0_mean::Vector{FT},
        uu0_cov::Matrix{FT},
        α_reg::FT = 1.0;
        update_freq::IT = 1,
        modified_unscented_transform::Bool = true,
        prior_mean::Union{Vector{FT}, Nothing} = nothing,
    ) where {FT <: AbstractFloat, IT <: Int}

Construct an Unscented Inversion EnsembleKalmanProcess.

Arguments
=========

  - `u0_mean`: Mean at initialization.
  - `uu0_cov`: Covariance at initialization.
  - `α_reg`: Hyperparameter controlling regularization toward the prior mean (0 < `α_reg` ≤ 1),
  default should be 1, without regulariazion.
  - `update_freq`: Set to 0 when the inverse problem is not identifiable,
  namely the inverse problem has multiple solutions, the covariance matrix
  will represent only the sensitivity of the parameters, instead of
  posterior covariance information; set to 1 (or anything > 0) when
  the inverse problem is identifiable, and the covariance matrix will
  converge to a good approximation of the posterior covariance with an
  uninformative prior.
  - `modified_unscented_transform`: Modification of the UKI quadrature given
    in Huang et al (2021).
  - `prior_mean`: Prior mean used for regularization.
"""
function Unscented(
    u0_mean::Vector{FT},
    uu0_cov::Matrix{FT};
    α_reg::FT = 1.0,
    update_freq::IT = 1,
    modified_unscented_transform::Bool = true,
    prior_mean::Union{Vector{FT}, Nothing} = nothing,
) where {FT <: AbstractFloat, IT <: Int}

    N_u = size(u0_mean, 1)
    # ensemble size
    N_ens = 2 * N_u + 1

    c_weights = zeros(FT, N_u)
    mean_weights = zeros(FT, N_ens)
    cov_weights = zeros(FT, N_ens)

    # set parameters λ, α
    α = min(sqrt(4 / N_u), 1.0)
    λ = α^2 * N_u - N_u

    c_weights[1:N_u] .= sqrt(N_u + λ)
    mean_weights[1] = λ / (N_u + λ)
    mean_weights[2:N_ens] .= 1 / (2 * (N_u + λ))
    cov_weights[1] = λ / (N_u + λ) + 1 - α^2 + 2.0
    cov_weights[2:N_ens] .= 1 / (2 * (N_u + λ))

    if modified_unscented_transform
        mean_weights[1] = 1.0
        mean_weights[2:N_ens] .= 0.0
    end

    u_mean = Vector{FT}[]  # array of Vector{FT}'s
    push!(u_mean, u0_mean) # insert parameters at end of array (in this case just 1st entry)
    uu_cov = Matrix{FT}[]  # array of Matrix{FT}'s
    push!(uu_cov, uu0_cov) # insert parameters at end of array (in this case just 1st entry)

    obs_pred = Vector{FT}[]  # array of Vector{FT}'s

    Σ_ω = (2 - α_reg^2) * uu0_cov
    Σ_ν_scale = 2.0

    r = isnothing(prior_mean) ? u0_mean : prior_mean
    iter = 0

    Unscented(
        u_mean,
        uu_cov,
        obs_pred,
        c_weights,
        mean_weights,
        cov_weights,
        Σ_ω,
        Σ_ν_scale,
        α_reg,
        r,
        update_freq,
        iter,
    )
end

"""
    construct_sigma_ensemble(
        process::Unscented,
        x_mean::Array{FT},
        x_cov::Matrix{FT},
    ) where {FT <: AbstractFloat, IT <: Int}

Construct the sigma ensemble based on the mean `x_mean` and covariance `x_cov`.
"""
function construct_sigma_ensemble(
    process::Unscented,
    x_mean::Array{FT},
    x_cov::Matrix{FT},
) where {FT <: AbstractFloat, IT <: Int}

    N_x = size(x_mean, 1)

    c_weights = process.c_weights

    chol_xx_cov = cholesky(Hermitian(x_cov)).L

    x = zeros(FT, N_x, 2 * N_x + 1)
    x[:, 1] = x_mean
    for i in 1:N_x
        x[:, i + 1] = x_mean + c_weights[i] * chol_xx_cov[:, i]
        x[:, i + 1 + N_x] = x_mean - c_weights[i] * chol_xx_cov[:, i]
    end

    return x
end

"""
construct_mean `x_mean` from ensemble `x`.
"""
function construct_mean(
    uki::EnsembleKalmanProcess{FT, IT, Unscented},
    x::Matrix{FT},
) where {FT <: AbstractFloat, IT <: Int}
    N_x, N_ens = size(x)

    @assert(uki.N_ens == N_ens)

    x_mean = zeros(FT, N_x)

    mean_weights = uki.process.mean_weights

    for i in 1:N_ens
        x_mean += mean_weights[i] * x[:, i]
    end

    return x_mean
end

"""
construct_cov `xx_cov` from ensemble `x` and mean `x_mean`.
"""
function construct_cov(
    uki::EnsembleKalmanProcess{FT, IT, Unscented},
    x::Matrix{FT},
    x_mean::Array{FT},
) where {FT <: AbstractFloat, IT <: Int}
    N_ens, N_x = uki.N_ens, size(x_mean, 1)

    cov_weights = uki.process.cov_weights

    xx_cov = zeros(FT, N_x, N_x)

    for i in 1:N_ens
        xx_cov .+= cov_weights[i] * (x[:, i] - x_mean) * (x[:, i] - x_mean)'
    end

    return xx_cov
end

"""
construct_cov `xy_cov` from ensemble x and mean `x_mean`, ensemble `obs_mean` and mean `y_mean`.
"""
function construct_cov(
    uki::EnsembleKalmanProcess{FT, IT, Unscented},
    x::Matrix{FT},
    x_mean::Array{FT},
    obs_mean::Matrix{FT},
    y_mean::Array{FT},
) where {FT <: AbstractFloat, IT <: Int, P <: Process}
    N_ens, N_x, N_y = uki.N_ens, size(x_mean, 1), size(y_mean, 1)

    cov_weights = uki.process.cov_weights

    xy_cov = zeros(FT, N_x, N_y)

    for i in 1:N_ens
        xy_cov .+= cov_weights[i] * (x[:, i] - x_mean) * (obs_mean[:, i] - y_mean)'
    end

    return xy_cov
end

"""
uki prediction step : generate sigma points
"""
function update_ensemble_prediction!(process::Unscented) where {FT <: AbstractFloat, IT <: Int}

    process.iter += 1
    # update evolution covariance matrix
    if process.update_freq > 0 && process.iter % process.update_freq == 0
        process.Σ_ω = (2 - process.α_reg^2) * process.uu_cov[end]
    end

    u_mean = process.u_mean[end]
    uu_cov = process.uu_cov[end]

    α_reg = process.α_reg
    r = process.r
    Σ_ω = process.Σ_ω

    N_u = length(process.u_mean[1])
    ############# Prediction step:

    u_p_mean = α_reg * u_mean + (1 - α_reg) * r
    uu_p_cov = α_reg^2 * uu_cov + Σ_ω

    ############ Generate sigma points
    u_p = construct_sigma_ensemble(process, u_p_mean, uu_p_cov)
    return u_p
end

"""
uki analysis step
g is the predicted observations  `Ny x N_ens` matrix
"""
function update_ensemble_analysis!(
    uki::EnsembleKalmanProcess{FT, IT, Unscented},
    u_p::Matrix{FT},
    g::Matrix{FT},
) where {FT <: AbstractFloat, IT <: Int}

    obs_mean = uki.obs_mean
    Σ_ν = uki.process.Σ_ν_scale * uki.obs_noise_cov

    ############# Prediction step:

    u_p_mean = construct_mean(uki, u_p)
    uu_p_cov = construct_cov(uki, u_p, u_p_mean)

    ###########  Analysis step

    g_mean = construct_mean(uki, g)
    gg_cov = construct_cov(uki, g, g_mean) + Σ_ν
    ug_cov = construct_cov(uki, u_p, u_p_mean, g, g_mean)

    tmp = ug_cov / gg_cov

    u_mean = u_p_mean + tmp * (obs_mean - g_mean)
    uu_cov = uu_p_cov - tmp * ug_cov'

    ########### Save results
    push!(uki.process.obs_pred, g_mean) # N_ens x N_data
    push!(uki.process.u_mean, u_mean) # N_ens x N_params
    push!(uki.process.uu_cov, uu_cov) # N_ens x N_data

    push!(uki.g, DataContainer(g, data_are_columns = true))

    compute_error!(uki)
end

function update_ensemble!(
    uki::EnsembleKalmanProcess{FT, IT, Unscented},
    g_in::Matrix{FT},
) where {FT <: AbstractFloat, IT <: Int}
    #catch works when g_in non-square
    if !(size(g_in)[2] == uki.N_ens)
        throw(DimensionMismatch("ensemble size in EnsembleKalmanProcess and g_in do not match, try transposing or check ensemble size"))
    end

    u_p_old = get_u_final(uki)

    #perform analysis on the model runs
    update_ensemble_analysis!(uki, u_p_old, g_in)
    #perform new prediction output to model parameters u_p
    u_p = update_ensemble_prediction!(uki.process)

    push!(uki.u, DataContainer(u_p, data_are_columns = true))

    return u_p
end

function get_u_mean_final(uki::EnsembleKalmanProcess{FT, IT, Unscented}) where {FT <: AbstractFloat, IT <: Int}
    return uki.process.u_mean[end]
end

function compute_error!(uki::EnsembleKalmanProcess{FT, IT, Unscented}) where {FT <: AbstractFloat, IT <: Int}
    mean_g = uki.process.obs_pred[end]
    diff = uki.obs_mean - mean_g
    X = uki.obs_noise_cov \ diff # diff: column vector
    newerr = dot(diff, X)
    push!(uki.err, newerr)
end

function get_error(uki::EnsembleKalmanProcess{FT, IT, Unscented}) where {FT <: AbstractFloat, IT <: Int}
    return uki.err
end

function Gaussian_2d(
    u_mean::Vector{FT},
    uu_cov::Matrix{FT},
    Nx::IT,
    Ny::IT;
    xx = nothing,
    yy = nothing,
) where {FT <: AbstractFloat, IT <: Int}
    # 2d Gaussian plot
    u_range = [min(5 * sqrt(uu_cov[1, 1]), 5); min(5 * sqrt(uu_cov[2, 2]), 5)]

    if xx === nothing
        xx = Array(LinRange(u_mean[1] - u_range[1], u_mean[1] + u_range[1], Nx))
    end
    if yy == nothing
        yy = Array(LinRange(u_mean[2] - u_range[2], u_mean[2] + u_range[2], Ny))
    end
    X, Y = repeat(xx, 1, Ny), Array(repeat(yy, 1, Nx)')
    Z = zeros(FT, Nx, Ny)

    det_uu_cov = det(uu_cov)

    for ix in 1:Nx
        for iy in 1:Ny
            Δxy = [xx[ix] - u_mean[1]; yy[iy] - u_mean[2]]
            Z[ix, iy] = exp(-0.5 * (Δxy' / uu_cov * Δxy)) / (2 * pi * sqrt(det_uu_cov))
        end
    end

    return xx, yy, Z
end


text:
#=
delete_node:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-31
=#

module NodeDelete

    # node class
    mutable struct Node
        data
        p_next
    end

    # head node
    head = nothing

    # method to print contents of linked list
    function print_list()
        tmp = head
        while tmp != nothing
            println(tmp.data)
            tmp = tmp.p_next
        end
    end

    function delete_node(data)
        tmp = head

        if tmp.data == data
            return head.p_next
        end

        while tmp.p_next != nothing
            if tmp.p_next.data == data
                tmp.p_next = tmp.p_next.p_next
                return head
            end
            tmp = tmp.p_next
        end

        return head
    end

    function main()
        # create simple linked list
        # next pointer of each node is initialized as Null
        global head = Node(1, nothing)
        second = Node(2, nothing)
        third = Node(3, nothing)
        fourth = Node(4, nothing)
        fifth = Node(5, nothing)

        head.p_next = second # link first node with second
        second.p_next = third # link second node with third
        third.p_next = fourth # link third node with fourth
        fourth.p_next = fifth # link forth node with fifth

        println("Original Linked List")
        print_list()
        println("")
        println("Delete 3 from Linked List")
        head = delete_node(3)
        print_list()
    end

end

if abspath(PROGRAM_FILE) == @__FILE__
    using .NodeDelete
    NodeDelete.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
"""
A web client for the Hydrologic Reference Stations of the Australian Bureau of Meteorology
at http://www.bom.gov.au/water/hrs.
"""
module HydroRefStations

import Base.show
using Dates: Date
import HTTP
using CSV: File
using DataStructures: OrderedDict
using DataFrames: DataFrame, first, rename, rename!, vcat, sort!, select!

export HRS, get_data, close!

const HRS_URL = "http://www.bom.gov.au/water/hrs/"
const SITES_URL = HRS_URL * "content/hrs_station_details.csv"
const DATA_URL = HRS_URL * "content/data/"

const HEADER_DELIM = "#"

const URL_SUFFIXES = OrderedDict(
    "daily data" => "daily_ts",
    "daily flow duration curve" => "daily_fdc_linear",
    "event frequency analysis" => "daily_event_histogram_frequency",
    "event volume analysis" => "daily_event_histogram_volume",

    "monthly data" => "",
    "january data" => "monthly_total_01",
    "february data" => "monthly_total_02",
    "march data" => "monthly_total_03",
    "april data" => "monthly_total_04",
    "may data" => "monthly_total_05",
    "june data" => "monthly_total_06",
    "july data" => "monthly_total_07",
    "august data" => "monthly_total_08",
    "september data" => "monthly_total_09",
    "october data" => "monthly_total_10",
    "november data" => "monthly_total_11",
    "december data" => "monthly_total_12",
    "january anomaly" => "monthly_anomaly_01",
    "february anomaly" => "monthly_anomaly_02",
    "march anomaly" => "monthly_anomaly_03",
    "april anomaly" => "monthly_anomaly_04",
    "may anomaly" => "monthly_anomaly_05",
    "june anomaly" => "monthly_anomaly_06",
    "july anomaly" => "monthly_anomaly_07",
    "august anomaly" => "monthly_anomaly_08",
    "september anomaly" => "monthly_anomaly_09",
    "october anomaly" => "monthly_anomaly_10",
    "november anomaly" => "monthly_anomaly_11",
    "december anomaly" => "monthly_anomaly_12",
    "monthly boxplot" => "monthly_boxplot",

    "seasonal data" => "",
    "summer data" => "seasonal_total_Summer",
    "autumn data" => "seasonal_total_Autumn",
    "winter data" => "seasonal_total_Winter",
    "spring data" => "seasonal_total_Spring",
    "summer anomaly" => "seasonal_anomaly_Summer",
    "autumn anomaly" => "seasonal_anomaly_Autumn",
    "winter anomaly" => "seasonal_anomaly_Winter",
    "spring anomaly" => "seasonal_anomaly_Spring",
    "seasonal boxplot" => "seasonal_boxplot",

    "annual data" => "annual_total",
    "cease to flow" => "annual_total_cease_to_flow",
    "annual anomaly" => "annual_anomaly",
    "3 year moving average" => "annual_anomaly_3MA",
    "5 year moving average" => "annual_anomaly_5MA",
    "11 year moving average" => "annual_anomaly_11MA"
)

const DATA_TYPES = Dict(
    "day" => [
        "daily data",
        "daily flow duration curve",
        "event frequency analysis",
        "event volume analysis"],
    "month"=> [
        "monthly data",
        "january data",
        "february data",
        "march data",
        "april data",
        "may data",
        "june data",
        "july data",
        "august data",
        "september data",
        "october data",
        "november data",
        "december data",
        "january anomaly",
        "february anomaly",
        "march anomaly",
        "april anomaly",
        "may anomaly",
        "june anomaly",
        "july anomaly",
        "august anomaly",
        "september anomaly",
        "october anomaly",
        "november anomaly",
        "december anomaly",
        "monthly boxplot"],
    "season" => [
        "seasonal data",
        "summer data",
        "autumn data",
        "winter data",
        "spring data",
        "summer anomaly",
        "autumn anomaly",
        "winter anomaly",
        "spring anomaly",
        "seasonal boxplot"],
    "year" => [
        "annual data",
        "cease to flow",
        "annual anomaly",
        "3 year moving average",
        "5 year moving average",
        "11 year moving average"])

const COMPOSITE_DATA_ATTRS = Dict(
    "seasonal data" => Dict(
        "data types" => ["autumn data", "winter data", "spring data", "summer data"],
        "months" => 3:3:12,
        "column name" =>  "Seasonal streamflow (ML/season)"),
    "monthly data" => Dict(
        "data types" => [
            "january data", "february data", "march data", "april data", "may data",
            "june data", "july data", "august data", "september data", "october data",
            "november data", "december data"],
        "months" => 1:12,
        "column name" =>  "Monthly streamflow (ML/season)")
)

"""
    hrs = HRS()

Open the connection and download the information of HRS service sites e.g. AWRC ID, name and location.

# Fields
* `sites`: Site information table
* `header`: Header of site information table
* `data_types`: Types of data available at the web site

# Examples
```julia
julia> hrs = HRS();
julia> hrs.sites
467×8 DataFrame. Omitted printing of 5 columns
│ Row │ AWRC Station Number │ Station Name                                 │ Latitude │
│     │ String              │ String                                       │ Float64  │
├─────┼─────────────────────┼──────────────────────────────────────────────┼──────────┤
│ 1   │ 410713              │ Paddy's River at Riverlea                    │ -35.3843 │
│ 2   │ 410730              │ Cotter River at Gingera                      │ -35.5917 │
...

julia> hrs.header
6-element Array{String,1}:
 "Australian Bureau of Meteorology"
 "Hydrologic Reference Stations"
...

julia> hrs.data_types
Dict{String,Array{String,1}} with 4 entries:
  "day"    => ["daily data", "daily flow duration curve", "event frequency analysis", "event vo…
  "month"  => ["monthly data", "january data", "february data", "march data", "april data", "ma…
  "year"   => ["annual data", "cease to flow", "annual anomaly", "3 year moving average", "5 ye…
  "season" => ["seasonal data", "summer data", "autumn data", "winter data", "spring data", "su…

julia> hrs.data_types["year"]
6-element Array{String,1}:
 "annual data"
 "cease to flow"
 "annual anomaly"
 "3 year moving average"
 "5 year moving average"
 "11 year moving average"

```
"""
struct HRS
    sites::DataFrame
    header::Array{String,1}
    data_types::Dict{String,Array{String,1}}

    function HRS()
        sites, header = get_sites()

        new(sites, header, DATA_TYPES)
    end
end

function show(io::IO, hrs::HRS)
    isempty(hrs.sites) && return

    for line in hrs.header[1:3]
        println(io, line)
    end
    println(io)
    show(io, first(hrs.sites, 6))
    println(io)
    println(io, "...")
end

"""
    close!(hrs::HRS)

Close the connection and reset the HRS object.
"""
function close!(hrs::HRS)
    select!(hrs.sites, Int[])
    empty!(hrs.header)
    empty!(hrs.data_types)
    return
end

"""
    data, header = get_data(hrs::HRS,
        awrc_id::AbstractString, data_type::AbstractString)

Return the data of a site.

# Arguments
* `hrs`: `HRS` object
* `awrc_id`: AWRC ID of the site. The ID can found in `hrs.sites`
* `data_type`: Type of the data.  Data type strings can be found in `hrs.data_types`

# Examples
```julia
julia> hrs = HRS();
julia> data, header = get_data(hrs, "410730", "annual data");
julia> data
55×2 DataFrame
│ Row │ Water Year (March to February) │ Annual streamflow (GL/water year) │
│     │ Int64                          │ Float64                           │
├─────┼────────────────────────────────┼───────────────────────────────────┤
│ 1   │ 1964                           │ 80.3924                           │
│ 2   │ 1965                           │ 19.7936                           │
...
```
"""
function get_data(hrs::HRS,
                  awrc_id::AbstractString,
                  data_type::AbstractString)::Tuple{DataFrame,Array{String,1}}

    if findfirst(isequal(awrc_id), hrs.sites[!,"AWRC Station Number"]) === nothing
        throw(ArgumentError("cannot find a site with AWRC ID, $(awrc_id)."))
    end

    data_type = lowercase(data_type)

    url_suffix = nothing
    try
        url_suffix = URL_SUFFIXES[data_type]
    catch e
        if isa(e, KeyError)
            throw(ArgumentError("unsupported data type, $(data_type)."))
        else
            rethrow()
        end
    end

    if length(url_suffix) > 0
        return get_raw_data(awrc_id, data_type)
    else
        return get_composite_data(awrc_id, data_type)
    end
end

"""
    sites, header = get_sites()

Download site information e.g. AWRC ID and descriptions.
"""
function get_sites()::Tuple{DataFrame,Array{String,1}}
    r = HTTP.get(SITES_URL)

    body_buf = IOBuffer(String(r.body))
    header = extract_header!(body_buf, HEADER_DELIM)
    new_header = prune_header(header, HEADER_DELIM)

    body_buf = seek(body_buf, 0)
    sites = DataFrame(File(body_buf, comment=HEADER_DELIM))

    return sites, new_header
end

"""
    data, header = get_raw_data(
        awrc_id::AbstractString, data_type::AbstractString)

Download multiple datasets, combine them to a time series.
"""
function get_raw_data(awrc_id::AbstractString,
                      data_type::AbstractString)::Tuple{DataFrame,Array{String,1}}

    url = get_url(awrc_id, data_type)
    r = HTTP.get(url)

    body_buf = IOBuffer(String(r.body))
    header = extract_header!(body_buf, HEADER_DELIM)
    new_header = prune_header(header, HEADER_DELIM)
    body_buf = seek(body_buf, 0)
    data = DataFrame(File(body_buf, comment=HEADER_DELIM))
    return data, new_header
end

"""
    data, header = get_composite_data(awrc_id::AbstractString, data_type::AbstractString)

Download multiple datasets and combine them to a time series.
"""
function get_composite_data(awrc_id::AbstractString,
                            data_type::AbstractString)::Tuple{DataFrame,Array{String,1}}
    attrs = COMPOSITE_DATA_ATTRS[data_type]
    data_list = []
    header = [""]
    for i in 1:length(attrs["data types"])
        dtype = attrs["data types"][i]
        data, header = get_raw_data(awrc_id, dtype)
        rename!(data, [:Year, :Q])
        data[!,:Date] = [Date(year, attrs["months"][i], 1) for year in data.Year]
        push!(data_list, data)
    end

    composite_data = vcat(data_list...)
    sort!(composite_data, :Date)
    composite_data = rename(composite_data[:,[:Date, :Q]],
                            ["Start Date", attrs["column name"]])

    return composite_data, header
end

"""
    url = get_url(awrc_id::AbstractString, data_type::AbstractString)

Return the URL to download data of a site.
"""
function get_url(awrc_id::AbstractString, data_type::AbstractString)::AbstractString
    url_suffix = URL_SUFFIXES[data_type]
    url = "$(DATA_URL)/$(awrc_id)/$(awrc_id)_$(url_suffix).csv"
    return url
end

"""
    header = extract_header!(body_buf::Base.GenericIOBuffer{Array{UInt8,1}}, delim::AbstractString)

Return the header document of the data file. Note that it moves the position of body_buf.
"""
function extract_header!(body_buf::Base.GenericIOBuffer{Array{UInt8,1}},
                         delim::AbstractString)::Array{String,1}
    header = String[]
    for line in eachline(body_buf)
        startswith(line, delim) ? push!(header, line) : break
    end
    return header
end

"""
    new_header = prune_header(header::AbstractString, delim::AbstractString)

Prune the extracted header document. It drops blank lines and double quotations from the raw document.
"""
function prune_header(header::Array{String,1}, delim::AbstractString)::Array{String,1}
    new_header = String[]
    for line in header
        #TODO: Give more information about the format error.
        startswith(line, delim) || throw(IOError("A header line does not start with $(delim)."))

        the_line = rstrip(line[2:end])
        isempty(the_line) && continue
        istart = 1
        i = findfirst('\"', the_line)
        (i != nothing) && (istart = i+1)
        iend = length(the_line)
        (the_line[end] == '\"') && (iend -= 1)

        push!(new_header, the_line[istart:iend])
    end
    return new_header
end

end


text:
## Vorticity-streamfunction formulation
#===========================================================================#

# Somewhat unnecessary, as the point is transformed to coordinates with the first point of the panel as the origin
function transform(x, n, x1, x2)
    s1, s2 = x - x1, x - x2

    # println("ses: $((s1, s2))")

    r1, r2 = √(s1^2 + n^2), √(s2^2 + n^2)
    β1, β2 = atan(s1, n), atan(s2, n)
    # r1, r2, β1, β2 = ifelse(r1 == 0 || r2 == 0, (1., 1., 0., 0.), (r1, r2, β1, β2))

    s1, s2, r1, r2, β1, β2
end

function vortex_stream_plus(str, x, n, x1, x2)
    s1, s2, r1, r2, β1, β2 = transform(x, n, x1, x2)
    str / 4π * (s1 * log(r1) - s2 * log(r2) + s2 - s1 + n * (β1 - β2))
end

function vortex_stream_minus(str, x, n, x1, x2)
    s1, s2, r1, r2, β1, β2 = transform(x, n, x1, x2)
    ψ_plus = vortex_stream_plus(1., x, n, x1, x2)

    str / (4π * (s1 - s2)) * ((s1 + s2) * ψ_plus + r2^2 * log(r2) - r1^2 * log(r1) + 0.5 * (s1^2 - s2^2))
end

function source_stream_trailing_edge(str, x, n, x1, x2)
    s1, s2, r1, r2, β1, β2 = transform(x, n, x1, x2)
    θ1, θ2 = atan(n, s1), atan(n, s2)

    str / 4π * ((s2 * θ2 - s1 * θ1) + n * log(r1 / r2))
end

# two_piece_source(λ1, λ2, r1, r2, r) = (r - r1) / (r2 - r1) * λ1 + (r2 - r) / (r2 - r1) * λ2

# function source_stream_plus(str, x, n, x_p, x_q)
#   s_p, s_q, r_p, r_q, β_p, β_q = transform(x, n, x_p, x_q)
#   str / 4π * (x_p * β_p - x_q * β_q + n * log(r_p / r_q))
# end

# function source_stream_minus(str, x, n, x_p, x_q)
#   s_p, s_q, r_p, r_q, β_p, β_q = transform(x, n, x_p, x_q)
#   λ_plus = source_stream_1(1., x, n, x_p, x_q)

#   str / (4π * (x_p - x_q)) * ((x_p + x_q) * λ_plus + r_p^2 * β_p - r_q^2 * β_q + n * (x_q - x_p))
# end


text:
function linear_compartment_model(graph, sinks)
    """
    Input:
        - graph - graph of the model network represented via adjacency lists
        - sinks - the indices of nodes having a sink
    Output: the corresponding ODE object where each parameter a_ij is replaced
            with a_ij + b_ij * x_0, where x_0 is a constant input encoded as a constant output
    """
    n = length(graph)
    x_vars_names = ["x$i" for i in 0:n]
    edges_vars_names = Array{String, 1}()
    for i in 1:n
        for j in graph[i]
            push!(edges_vars_names, "a_$(i)_$(j)")
            push!(edges_vars_names, "b_$(i)_$(j)")
        end
    end
    for s in sinks
        push!(edges_vars_names, "a_$(s)_0")
        push!(edges_vars_names, "b_$(s)_0")
    end
    R, vars = PolynomialRing(QQ, vcat(x_vars_names, edges_vars_names))
    x_vars = vars[2:(n + 1)]
    x0 = vars[1]
    equations = Dict{fmpq_mpoly, Union{fmpq_mpoly, Generic.Frac{fmpq_mpoly}}}(x => R(0) for x in x_vars)
    equations[x0] = R(0)
    for i in 1:n
        for j in graph[i]
            rate = str_to_var("a_$(i)_$(j)", R) + str_to_var("b_$(i)_$(j)", R) * x0

            if i != j
                equations[x_vars[j]] += x_vars[i] * rate
                equations[x_vars[i]] -= x_vars[i] * rate
            else
                equations[x_vars[i]] -= x_vars[i] * rate
            end
        end
        if i in sinks
            rate = str_to_var("a_$(i)_0", R) + str_to_var("b_$(i)_0", R) * x0
            equations[x_vars[i]] += -x_vars[i] * rate
        end
    end
    return ODE{fmpq_mpoly}(equations)
end

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

function bicycle(n)
    """
    Generates a bidirected cycle of length n
    """
    graph = []
    for i in 1:n
        prev = (i == 1) ? n : (i - 1)
        next = (i == n) ? 1 : i + 1
        push!(graph, [prev, next])
    end
    return graph
end

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

function cycle(n)
    """
    Single directed cycle
    """
    graph = [[(i == n) ? 1 : (i + 1)] for i in 1:n]
    return graph
end

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

function catenary(n)
    """
    Bidirected chain from 1 to n
    """
    graph = [[] for i in 1:n]
    for i in 1:n
        if i != 1
            push!(graph[i], i - 1)
        end
        if i != n
            push!(graph[i], i + 1)
        end
    end
    return graph
end

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

function mammilary(n)
    """
    Bidirected 'star' with center at 1 and rays to 2, ..., n
    """
    graph = []
    push!(graph, [i for i in 2:n])
    for i in 2:n
        push!(graph, [1])
    end
    return graph
end

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

max_stars_repo_name: pogudingleb/ExperimentsBound

max_issues_repo_name: pogudingleb/ExperimentsBound

max_forks_repo_name: pogudingleb/ExperimentsBound