text:
# Julia code
# Copyright (c) 2022 the SLICOTMath.jl developers
# Portions extracted from SLICOT-Reference distribution:
# Copyright (c) 2002-2020 NICONET e.V.

function run_mb03lf(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    NMAX = 50
    LDB = NMAX÷2
    LDFG = NMAX÷2
    LDQ = 2*NMAX
    LDU = NMAX
    LDZ  = NMAX
    LDWORK = 10*NMAX*NMAX + max( NMAX*NMAX + max( NMAX÷2 + 252, 432 ), max( 8*NMAX +  48, 171 ) )
    LIWORK = max( NMAX + 18, NMAX÷2 + 48, 5*NMAX÷2 + 1 )
    Z = Array{Float64,2}(undef, LDZ,NMAX)
    B = Array{Float64,2}(undef, LDB,NMAX÷2)
    FG = Array{Float64,2}(undef, LDFG,NMAX÷2+1)
    Q = Array{Float64,2}(undef, LDQ,2*NMAX)
    U = Array{Float64,2}(undef, LDU,2*NMAX)
    ALPHAR = Array{Float64,1}(undef, NMAX÷2)
    ALPHAI = Array{Float64,1}(undef, NMAX÷2)
    BETA = Array{Float64,1}(undef, NMAX÷2)
    BWORK = Array{BlasBool,1}(undef, NMAX÷2)
    IWORK = Array{BlasInt,1}(undef, LIWORK)
    DWORK = Array{Float64,1}(undef, LDWORK)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    COMPQ =  vs[1][1]
    COMPU =  vs[2][1]
    ORTH =  vs[3][1]
    N = parse(BlasInt, vs[4])
if ( N<0 || N>NMAX || mod( N, 2 )!=0 )
else
    M = N÷2

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       Z[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (M,M)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       B[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (M,M+1)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       FG[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
# interp call 1

    NEIG, INFO, IWARN = SLICOT.mb03lf!(COMPQ, COMPU, ORTH, N, Z, B, FG, Q, U, ALPHAR, ALPHAI, BETA, LIWORK)
    println(io, "NEIG = $NEIG")
    @test INFO == 0
    INFO == 0 || return
    println(io, "IWARN = $IWARN")

if ( INFO!=0 )
else
# interp output 1
    println(io,"Z:")
    _nc = N
    _nr = N
    show(io,"text/plain",Z[1:_nr,1:_nc])
    println(io,)

# interp output 2
    println(io,"ALPHAR:")
    _nr = M
    show(io,"text/plain",ALPHAR[1:_nr])
    println(io,)

# interp output 3
    println(io,"ALPHAI:")
    _nr = M
    show(io,"text/plain",ALPHAI[1:_nr])
    println(io,)

# interp output 4
    println(io,"BETA:")
    _nr = M
    show(io,"text/plain",BETA[1:_nr])
    println(io,)

if ( LSAME( COMPQ, 'C' ) && NEIG>0 )
# interp output 5
    println(io,"Q:")
    _nc = NEIG
    _nr = N
    show(io,"text/plain",Q[1:_nr,1:_nc])
    println(io,)

end # if
if ( LSAME( COMPU, 'C' ) && NEIG>0 )
# interp output 6
    println(io,"U:")
    _nc = NEIG
    _nr = N
    show(io,"text/plain",U[1:_nr,1:_nc])
    println(io,)

end # if
end # if
end # if
    close(f)
end # run_X()


text:
module test_julia_functions

using Test

struct G
    a
end

f(g::G, c=g.a) = c
@test f(G(2)) == 2

g(;a) = a+1
@test invoke(Core.kwfunc(g), Tuple{NamedTuple{(:a,), Tuple{Int}}, typeof(g)}, (a=2,), g) == 3

function execute(nts::Vector{<:NamedTuple})
    1
end

@test execute([(a=1, b=2)]) == 1
@test execute([(c=1, d=2)]) == 1

function execute(nts::Vector{<:NamedTuple{(:a,:b)}})
    2
end

@test execute([(a=1, b=2)]) == 2
@test execute([(c=1, d=2)]) == 1

str(a; b="$a") = b
@test str(42) == "42"

@test str(42; b=3) == 3

b = 3
@test str(42; b) == 3

h(; kwargs...) = kwargs
@test (; h(a=1, b="")...).a == 1

end # module test_julia_functions

module test_julia_functions_return

using Test

function f()
    a::Float32 = 1
    a
end

function g()
    a::Float32 = 0.0
    a = 1
    a
end

function h()
    a::Float32 = 1
end

@test f() isa Float32
@test g() isa Float32
@test h() isa Int

end # module test_julia_functions_return

module test_julia_functions_inner

using Test

function f(val=true)
    function g()
        val = false
    end
    g()
    val
end

@test f() == false

end # module test_julia_functions_inner

module test_julia_functions_anonymous

using Test

@test (() -> 1)()   == 1
@test ((x) -> x)(1) == 1

end # module test_julia_functions_anonymous

module test_julia_functions_propertynames

using Test

function builder()
    n = 1
    f() = n + 1
end

g = builder()
@test propertynames(g) == (:n,)
@test g.n == 1

end # module test_julia_functions_propertynames

ext: jl


text:
module Transpilers

include("julia/julia_transpiler.jl")

end


text:
module Powerlog

using Distributions, Plots

export powerlog

# !* version 1.0 -- 12/1/06 -- pbe
# !* version .9 -- 9/26/06 -- pbe
# program define powerlog, rclass
# version 9.2

# /*------------------------------------------------------------------*/
# /*    Based on the SAS macro program, POWERLOG.SAS (version 1.1)    */
# /*    Created 1998, revised 2002                                    */
# /*    SAS author:  Michael Friendly   <FRIENDLY@YorkU.CA>           */
# /*------------------------------------------------------------------*/

# /*------------------------------------------------------------------*/
# /*       p1 -- estimated probability of a 1 at the mean             */
# /*             of the predictor                                     */
# /*       p2 -- estimated probability of a 1 at an X-value           */
# /*             equal to the X-mean plus 1 sd                        */
# /*      rsq -- squared multiple correlation of the predictor        */
# /*             with all other predictors.  Use rsq=0 for a          */
# /*             1-predictor model.                                   */
# /*    alpha -- significance level of test of factor A (default=.05) */
# /*------------------------------------------------------------------*/

"""
    powerlog(p1::Float64,p2::Float64; alpha::Float64 = 0.05, rsq::Real = 0, graph = false, help = false)

Perform power analysis for logistic regression designs with a
continuous predictor variable.  Multiple predictor models can be accomodated using
the `rsq` option. This is a port of Stata powerlog.ado, which was originally ported
from the SAS macro program, POWERLOG.SAS.

## Options
- `p1` = estimated probability of an event (1) at the mean of the predictor
- `p2` = estimated probability of a `1` at the mean + 1 sd of the predictor (required)
- `rsq` = squared multiple correlation of the predictor with all other predictors. Use rsq = 0 for a 1-predictor model (default: 0.0).
- `alpha` = the alpha level, probability of type I error (default: 0.05)
- `graph` = set to `true` to generate a graph of the power analysis (default: false)
- `help` = set to `true` to display explanation of terms

## Examples

```jldoctest
julia> powerlog(.25,.35,alpha = 0.01, help = true)

Logistic regression power analysis
One-tailed test: alpha = 0.01, p1 = 0.25, p2 = 0.35, rsq = 0, odds ratio = 1.6153846153846154

 power        n
  0.60      192
  0.65      211
  0.70      232
  0.75      256
  0.80      284
  0.85      319
  0.90      365
  0.95      439

Explanation of terms
p1  -- the probability that the response variable equals 1
         when the predictor is at the mean
p2  -- the probability that the response variable equals 1
         when the predictor is one standard deviation above the mean
rsq -- the squared mulitple correlation between the predictor
         variable and all other variables in the model

julia> powerlog(.25,.35,rsq=0.4)

Logistic regression power analysis
One-tailed test: alpha = 0.05, p1 = 0.25, p2 = 0.35, rsq = 0.4, odds ratio = 1.6153846153846154

 power          n
 0.60         173
 0.65         196
 0.70         223
 0.75         253
 0.80         290
 0.85         335
 0.90         397
 0.95         498

```

## Acknowledgements
Based on (copied from) the SAS macro program, POWERLOG.SAS (version 1.1),
Created 1998, revised 2002, by  Michael Friendly  <FRIENDLY@YorkU.CA>,
and STATA program, powerlog.ado, by Philip B. Ender <ender@ucla.edu>,
UCLA Academic Technology Services.

"""
function powerlog(p1::Float64,p2::Float64; alpha::Float64 = 0.05, rsq::Real = 0, graph = false, help = false)

    pd = p2 - p1
    l1 = p1/(1-p1)
    l2 = p2/(1-p2)
    θ = l2 / l1
    or = θ
    λ = log(θ)
    λ2 = λ^2
    za = quantile(Normal(),1-alpha)

    # start output
    println("\nLogistic regression power analysis")
    println("One-tailed test: alpha = ",alpha,", p1 = ",p1,", p2 = ",p2,", rsq = ",rsq,", odds ratio = ",or)
    println("\n power          n")

    δ = (1 + (1 + λ2)*exp(5 * λ2/4))/(1 + exp(-1*λ2/4))

    pwr = zeros(Float64,8)
    nn = zeros(Int64,8)

    i = 1
    for power = 0.6:.05:.95
        zb = quantile(Normal(),power)

        N = ((za + zb*exp(-1 * λ2/4))^2 * (1 + 2*p1*δ))/(p1*λ2)
        N /= (1 - rsq)

        @printf(" %3.2f  %10.0f\n",power,N)

        # save for a graph
        pwr[i] = power
        nn[i] = ceil(Int64,N)
        i += 1
    end

    if help == true
        println("\nExplanation of terms")
        println("p1  -- the probability that the response variable equals 1")
        println("         when the predictor is at the mean")
        println("p2  -- the probability that the response variable equals 1")
        println("         when the predictor is one standard deviation above the mean")
        println("rsq -- the squared mulitple correlation between the predictor")
        println("         variable and all other variables in the model")
    end

    if graph == true
        gr()
        len = length(nn)
        n_lower = nn[1]
        n_upper = nn[len]
        xlim_lower = nn[len] / 80
        plot(nn,pwr,
            linecolor = :blue,
            marker=(2,.5,:circle,:blue),
            title = "Logistic regression power analysis",
            yaxis = ("Power",(-0.03,1.03),0.0:0.2:1.0),
            xaxis = ("Sample Size",(n_lower - xlim_lower,n_upper + xlim_lower)),
            legend=false)
        vline!([nn[5]],line=(1,:dot,.8,:red))
    end
end

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ 1eb559b3-dbce-4caf-800f-96a2495dee1b
using Revise

# ╔═╡ 77d14c3f-bcaf-432c-9005-f5ec184d84b7
using QuadGK

# ╔═╡ ba2dadc8-d9a2-452e-9f94-406ccfc7404b
using Optim

# ╔═╡ c05cd489-68c8-418e-bbc2-a33bc42781f0
using Plots

# ╔═╡ 1d300e7d-5efb-4ad2-a826-c3f8a86083d8
using PolaronMobility

# ╔═╡ 6a99cd13-68cb-474c-96af-7d78ee4e7095
using DataFrames

# ╔═╡ 49663b33-1604-430b-a815-9eb985802987
using CSV

# ╔═╡ 314af942-8922-40d4-a92b-f75b715c276c
using LaTeXStrings

# ╔═╡ 6ff68ba5-df0c-42a0-8cf1-dedeced77e8d
using ColorSchemes

# ╔═╡ 29d62408-0629-4f39-909e-a6ba8696e955
begin
	# Physical constants
	const ħ = 1.05457162825e-34; # Reduced Planck's constant (kg m^2 s^{-1})
	const eV = 1.602176487e-19; # Electron Volt (kg m^2 s^{-2})
	const m_e = 9.10938188e-31; # Electron Mass (kg)
	const k_B = 1.3806504e-23; # Boltzmann's constant (kg m^2 K^{-1} s^2)
	const ϵ_0 = 8.854e-12 # Dielectric constant (C^2 N^{-1} m^{-2})
	const c = 2.99792458e8 # Speed of light (m s^{-1})
	const amu = 1.66053906660e-27 # Atomic Mass Unit (kg)
end

# ╔═╡ 35390bdf-e260-4566-8a97-7a0b5a3ef07a
MAPI= [
# 96.20813558773261 0.4996300522819191
# 93.13630357703363 1.7139631746083817
# 92.87834578121567 0.60108592692181
# 92.4847918585963 0.0058228799414729
# 92.26701437594754 0.100590086574602
# 89.43972834606603 0.006278895133832249
# 46.89209141511332 0.2460894564364346
# 46.420949316788 0.14174282581124137
# 44.0380222871706 0.1987196948553428
# 42.89702947649343 0.011159939465770681
# 42.67180170168193 0.02557751102757614
# 41.46971205834201 0.012555230726601503
# 37.08982543385215 0.00107488277468418
# 36.53555265689563 0.02126940080871224
# 30.20608114002676 0.009019481779712388
# 27.374810898415028 0.03994453721421388
# 26.363055017011728 0.05011922682554448
# 9.522966890022039 0.00075631870522737
4.016471586720514 0.08168931020200264
3.887605410774121 0.006311654262282101
3.5313112232401513 0.05353548710183397
2.755392921480459 0.021303020776321225
2.4380741812443247 0.23162784335484837
2.2490917637719408 0.2622203718355982
2.079632190634424 0.23382298607799906
2.0336707697261187 0.0623239656843172
1.5673011873879714 0.0367465760261409
1.0188379384951798 0.0126328938653956
1.0022960504442775 0.006817361620021601
0.9970130778462072 0.0103757951973341
0.9201781906386209 0.01095811116040592
0.800604081794174 0.0016830270365341532
0.5738689505255512 0.00646428491253749
# 0.022939578929507105 8.355742795827834e-06   # Acoustic modes!
# 0.04882611767873102 8.309858592685e-06
# 0.07575149723846182 2.778248540373041e-05
]

# ╔═╡ e0ce47a5-368a-4da7-83fc-9f921af0f747
function ϵ_ionic_mode(phonon_mode_freq, ir_activity, volume) # single ionic mode
    ω_j = 2π * phonon_mode_freq * 1e12 # angular phonon freq in Hz
    ϵ_mode = eV^2 * ir_activity / (3 * volume * ω_j^2 * amu) # single ionic mode
    return ϵ_mode / ϵ_0 # normalise with 1 / (4π ϵ_0)
end

# ╔═╡ 053d717a-2b63-40e7-9dc0-487e3892de98
function ϵ_total(freqs_and_ir_activity, volume) # total ionic contribution to dielectric
    phonon_freqs = freqs_and_ir_activity[:, 1]
    ir_activity = freqs_and_ir_activity[:, 2]
    result = 0.0
    for (f, r) in zip(phonon_freqs, ir_activity)
        result += ϵ_ionic_mode(f, r, volume) # sum over all ionic contributions
    end
    return result
end

# ╔═╡ 86fe8a15-4d9e-48fa-8387-3cfaff7aa182
function frohlich_α_j(ϵ_optic, ϵ_ionic, ϵ_total, phonon_mode_freq, m_eff) # Frohlich alpha decomposed into phonon branch contributions
    Ry = eV^4 * m_e / (2 * ħ^2) # Rydberg energy
    ω = 2π * 1e12 * phonon_mode_freq # angular phonon freq (Hz)
    ϵ_static = ϵ_total + ϵ_optic # static dielectric. Calculate here instead of input so that ionic modes properly normalised.
    return (m_eff * Ry / (ħ * ω))^(1 / 2) * ϵ_ionic / (4π * ϵ_0) / (ϵ_optic * ϵ_static) # 1 / (4π ϵ_0) dielectric normalisation
end

# ╔═╡ 16e7d781-bc61-44b3-84a2-e8c8892bf773
function κ_i(i, v, w) # fictitious spring constant, multiple variational params
    κ = v[i]^2 - w[i]^2
    if length(v) > 1
        for j in 1:length(v)
            if j != i
                κ *= (v[j]^2 - w[i]^2) / (w[j]^2 - w[i]^2)
            end
        end
    end
    return κ
end

# ╔═╡ a78017d7-8835-4a8f-8dd3-f3a4a50eb37b
function h_i(i, v, w) # some vector relating to harmonic eigenmodes
    h = v[i]^2 - w[i]^2
    if length(v) > 1
        for j in 1:length(v)
            if j != i
                h *= (w[j]^2 - v[i]^2) / (v[j]^2 - v[i]^2)
            end
        end
    end
    return h
end

# ╔═╡ d7346684-6b68-4dfc-81da-12c435c4acf0
function C_ij(i, j, v, w) # generalised Feynman C variational parameter (inclusive of multiple v and w params)
    C = w[i] * κ_i(i, v, w) * h_i(j, v, w) / (4 * (v[j]^2 - w[i]^2))
    return C
end

# ╔═╡ 798a4f7a-a2a0-48de-82ad-f85b4199496a
function D_j(τ, β, v, w) # log of dynamic structure factor for polaron
    D = τ * (1 - τ / β)
    for i in 1:length(v)
        if v[i] != w[i]
        D += (h_i(i, v, w) / v[i]^2) * (2 * sinh(v[i] * τ / 2) * sinh(v[i] * (β - τ) / 2) / (v[i] * sinh(v[i] * β / 2)) - τ * (1 - τ / β))
        end
    end
    return D
end

# ╔═╡ ec0744c0-2530-11ec-04c3-1b5d83413b78
function multi_free_energy(v_params, w_params, T, ϵ_optic, m_eff, volume, freqs_and_ir_activity)

    setprecision(BigFloat, 32) # Speed up. Stops potential overflows.

    # Extract phonon frequencies and ir activities.
    phonon_freqs = freqs_and_ir_activity[:, 1]
    ir_activity = freqs_and_ir_activity[:, 2]

    num_of_branches = length(phonon_freqs)

    # total dielectric contribution from all phonon branches (used as a normalisation)
    ϵ_tot = ϵ_total(freqs_and_ir_activity, volume)

    # Generalisation of B i.e. Equation 62c in Hellwarth.
    B_j_integrand(τ, β, v, w) = cosh(β / 2 - abs(τ)) / (sinh(β / 2) * sqrt(D_j(abs(τ), β, v, w)))
    B_j(β, α, v, w) = α / √π * quadgk(τ -> B_j_integrand(τ, β, v, w), 0.0, β / 2)[1]

    # Generalisation of C i.e. Equation 62e in Hellwarth.
    function C_j(β, v, w)
        s = 0.0
        for i in 1:length(v)
            for j in 1:length(v)
                s += C_ij(i, j, v, w) / (v[j] * w[i]) * (coth(β * v[j] / 2)  - 2 / (β * v[j]))
            end
        end
        3 * s / num_of_branches
    end

    # Generalisation of A i.e. Equation 62b in Hellwarth.
    function A_j(β, v, w)
        s = -log(2π * β) / 2
        for i in 1:length(v)
            if v[i] != w[i]
                s += log(v[i] / w[i]) - log(sinh(v[i] * β / 2) / sinh(w[i] * β / 2))
            end
        end
        3 / β * s / num_of_branches
    end

	F = 0.0
	for j in 1:num_of_branches
		ω_j = 2π * 1e12 * phonon_freqs[j] # angular phonon freq im 2π Hz
		β_j = BigFloat(ħ * ω_j / (k_B * T)) # reduced thermodynamic beta
		ϵ_ionic_j = ϵ_ionic_mode(phonon_freqs[j], ir_activity[j], volume) # ionic dielectric contribution for current phonon branch
		α_j = frohlich_α_j(ϵ_optic, ϵ_ionic_j, ϵ_tot, phonon_freqs[j], m_eff) # decomposed alpha for current phonon branch

		# F = -(A + B + C) in Hellwarth.
		F += -(B_j(β_j, α_j, v_params, w_params) + C_j(β_j, v_params, w_params) + A_j(β_j, v_params, w_params)) * ω_j # × ħω branch phonon energy
	end

    return F * ħ / eV * 1e3 # change to meV
end

# ╔═╡ f3b91e16-26d3-45c8-b28f-193afac15d7f
function multi_variation(T, ϵ_optic, m_eff, volume, freqs_and_ir_activity; initial_vw = false, N = 1) # N number of v and w params

    setprecision(BigFloat, 32) # Speed up. Stops potential overflows.

    if initial_vw isa Bool
		# Intial guess for v and w.
    	initial = sort(rand(2 * N)) .* 4.0 .+ 1.0 # initial guess around 4 and ≥ 1.

		# Limits of the optimisation.
		lower = repeat([0.1], 2 * N)
		upper = repeat([60.0], 2 * N)
	else
		# Intial guess for v and w.
		initial = sort(vcat(initial_vw...))

		# Limits of the optimisation.
		lower = repeat([0.1], 2 * N)
		upper = repeat([60.0], 2 * N)
	end

	println("Initial guess: ", initial)

	# Osaka Free Energy function to minimise.
	f(x) = multi_free_energy([x[2 * n] for n in 1:Int(N)], [x[2 * n - 1] for n in 1:Int(N)], T, ϵ_optic, m_eff, volume, freqs_and_ir_activity)

	# Use Optim to optimise the free energy function w.r.t v and w.
	solution = Optim.optimize(
		Optim.OnceDifferentiable(f, initial; autodiff = :forward),
		lower,
		upper,
		initial,
		Fminbox(BFGS()),
		Optim.Options(time_limit = 20.0),
	)

	# Get v and w params that minimised free energy.
	var_params = Optim.minimizer(solution)

	# Update matrices for v and w parameters.
	v_params = [var_params[2 * n] for n in 1:N]
	w_params = [var_params[2 * n - 1] for n in 1:N]

	# Show current v and w that minimise jth phonon branch.
	println("Variational parameters: ", var_params)

    # Return variational parameters that minimise the free energy.
    return v_params, w_params
end

# ╔═╡ a941ccc6-90cb-448a-91bb-f80233b630f8
function HellwarthAScheme(phonon_modes; T = 295)

	phonon_mode_freqs = phonon_modes[:, 1]
	ir_activities = phonon_modes[:, 2]

	condition(f) = coth(π * f * 1e12 * ħ / (k_B * T)) / f - sum(ir_activities .* coth.(π .* phonon_mode_freqs .* 1e12 .* ħ ./ (k_B * T)) ./ phonon_mode_freqs) / sum(ir_activities)

	minimum_frequency = minimum(phonon_mode_freqs)
	maximum_frequency = maximum(phonon_mode_freqs)
	middle_frequency = (maximum_frequency + minimum_frequency) / 2
	print("\n")

	while (maximum_frequency - minimum_frequency) / 2 > 1e-6

		if sign(condition(middle_frequency)) == sign(condition(minimum_frequency))
			minimum_frequency = middle_frequency
			middle_frequency = (maximum_frequency + minimum_frequency) / 2
		else
			maximum_frequency = middle_frequency
			middle_frequency = (maximum_frequency + minimum_frequency) / 2
		end
	end

	return middle_frequency
end

# ╔═╡ 1372f247-4622-4388-b458-89aa6e695388
function HellwarthBScheme(LO)
    println("Hellwarth B Scheme... (athermal)")
    H58 = sum( LO[:,2] ./ LO[:,1].^2 )
    println("Hellwarth (58) summation: ",H58)

    H59 = sum( LO[:,2] ) # sum of total ir activity squarred
    println("Hellwarth (59) summation (total ir activity ^2): ", H59)
    println("Hellwarth (59) W_e (total ir activity ): ", sqrt(H59))

    omega = sqrt(H59 / H58)
    println("Hellwarth (61) Omega (freq): ",omega)

	return(omega)
end

# ╔═╡ c28f6d35-c3ed-4984-9c36-b282fb50001b
function ridgeline(x, y, z; shift = 1, jump = 10, ylabel = :none, xlabel = :none, palette = :twilight, linewidth = 1.5, color = :black, size = (595, 842), ymirror = false, fillalpha = 1, step = 1)
	if !ymirror
		p = plot(x, z[:, 1], fillrange = z[:, jump] .- shift, fillalpha = fillalpha, ylabel = ylabel, xlabel = xlabel, xticks = (0:maximum(x)/10:maximum(x), string.(x[1:Int(ceil((length(x)-1)/10)):end])), yticks = (z[1, 1:jump * step:end] .- shift .* range(0, stop = length(y[1:jump:end]) - 1, step = step), string.(y[1:jump * step:end])), legend = false, size = size, grid = false, ymirror = ymirror, palette = palette, tickfontsize = 12, labelfontsize = 12)
	else
		p = plot(x, z[:, 1], fillrange = z[:, jump] .- shift, fillalpha = fillalpha, ylabel = ylabel, xlabel = xlabel, xticks = (0:maximum(x)/10:maximum(x), string.(x[1:Int((length(x))/10):end])), yticks = (z[end, 1:jump * step:end] .- shift .* range(0, stop = length(y[1:jump:end]) - 1, step = step), string.(y[1:jump * step:end])), legend = false, size = size, grid = false, ymirror = ymirror, palette = palette, tickfontsize = 12, labelfontsize = 12)
	end

	plot!(x, z[:, 1], color = color, linewidth = linewidth)
	plot!(x, z[:, jump] .- shift, color = color, linewidth = linewidth)

	for i in 1:Int(floor(length(y)/jump) - 1)
	plot!(x, z[:, i * jump] .- shift * i, fillrange = z[:, (i + 1) * jump] .- shift * (i + 1), fillalpha = fillalpha, palette = palette)
	plot!(x, z[:, i * jump] .- shift * i, color = color, linewidth = linewidth)
	plot!(x, z[:, (i + 1) * jump] .- shift * (i + 1), color = color, linewidth = linewidth)
	end
	return p
end

# ╔═╡ 49480256-2dcc-4d10-a5a1-4202b243cfd2
T_range = 100:100:400

# ╔═╡ 4966ecf1-ef39-4da3-a8a5-9e47a6c6b0ec
phonon_mode_freqs = MAPI[:, 1]

# ╔═╡ 3321246b-76bd-4678-a5ae-332c3543a1ef
ir_activities = MAPI[:, 2]

# ╔═╡ b2944c8a-158e-461b-a83c-842925587179
A_scheme = [HellwarthAScheme(MAPI, T = i) for i in T_range]

# ╔═╡ fe4e5ef5-833a-4b61-9ea9-92fa52d7a21f
B_scheme = HellwarthBScheme(hcat(phonon_mode_freqs, ir_activities))

# ╔═╡ 8350c7a2-e9fd-4312-9561-44e2524734af
ϵ_ionic = [ϵ_ionic_mode(f, r, (6.29e-10)^3) for (f, r) in zip(phonon_mode_freqs, ir_activities)]

# ╔═╡ db2b4169-0737-44c1-9af5-aa2f2403f5db
ϵ_tot = sum(ϵ_ionic)

# ╔═╡ 2b6d00e1-317b-41cc-9489-14a1450012f4
α_j = [frohlich_α_j(4.5, ϵ_i, ϵ_tot, f, 0.12) for (ϵ_i, f) in zip(ϵ_ionic, phonon_mode_freqs)]

# ╔═╡ 3e1418c2-01a1-4651-8029-59f53264193d
begin
	scatter(phonon_mode_freqs, α_j, yaxis = :log, xticks = (phonon_mode_freqs, hcat(["$(round(i, digits = 2))" for i in phonon_mode_freqs]...)), label = L"\textbf{\alpha_j}", size = (500, 500), xlabel = "Phonon Frequencies (THz)", tickfontsize = 9, xrotation = 90)
	scatter!(phonon_mode_freqs, ir_activities, markershape = :diamond, label = "IR")
	plot!(phonon_mode_freqs, α_j ./ ir_activities, linewidth = 2, linestyle = :dash, label = "ratio")
end

# ╔═╡ d4f392ef-968a-439c-ba53-5707464446e7
begin
	multi_data = DataFrames.DataFrame(
			alpha = α_j,
			phonon_freqs = phonon_mode_freqs,
			ir_activities = ir_activities,
			ionic = ϵ_ionic
	    )
	CSV.write("multi_data.csv", multi_data)
end

# ╔═╡ b3c94ea5-d152-4705-9a30-f8d3e4da0874
α_eff = sum(α_j)

# ╔═╡ 3608af5e-7b44-4d50-9770-5e8afd4fca28
α_hellwarth_A = PolaronMobility.frohlichalpha.(4.5, 24.1, A_scheme .* 1e12, 0.12)

# ╔═╡ 1716f263-f8ad-40a1-a805-d8b7ea6f38df
α_hellwarth_B = PolaronMobility.frohlichalpha(4.5, 24.1, B_scheme * 1e12, 0.12)

# ╔═╡ da629bec-b53c-4acf-a0ab-6057c3b973a9
β_j = [ħ * 2π * 1e12 * f / k_B / T for f in phonon_mode_freqs, T in T_range]

# ╔═╡ 0641f475-0f35-4f5c-abf0-85ec56f4e30c
begin
	multi_beta = DataFrames.DataFrame([[0.0, [i for i in T_range]...]'; [phonon_mode_freqs β_j]], :auto)
	CSV.write("multi_beta.csv", multi_beta)
end

# ╔═╡ c49a4dc3-2546-472f-9fde-b59ec3246237
β_j_avg = [sum(β_j[:, T]) / length(phonon_mode_freqs) for T in 1:length(T_range)]

# ╔═╡ 0e3cf449-6e2a-4353-853e-58d15693743d
β_hellwarth_A = [ħ * 2π / k_B / T_range[T] * A_scheme[T] * 1e12 for T in 1:length(T_range)]

# ╔═╡ 16040051-57d5-4c6c-aa13-8280ae964261
β_hellwarth_B = [ħ * 2π / k_B / T * B_scheme * 1e12 for T in T_range]

# ╔═╡ 32f32534-2a6b-4675-8508-890579058e07
begin
	var_params = []
	push!(var_params, multi_variation(T_range[1], 4.5, 0.12, (6.29e-10)^3, MAPI; N = 1))
	for i in 2:length(T_range)
		push!(var_params, multi_variation(T_range[i], 4.5, 0.12, (6.29e-10)^3, MAPI; initial_vw = var_params[i-1], N = 1))
	end
	var_params
end

# ╔═╡ de89f0aa-21fb-42d2-9a75-5b45938b1e65
var_A = PolaronMobility.feynmanvw.(α_hellwarth_A, β_hellwarth_A)

# ╔═╡ a74836cf-f9ab-41d6-9f01-5066f1afac34
var_B = PolaronMobility.feynmanvw.(α_hellwarth_B, β_hellwarth_B)

# ╔═╡ 2ea4fc88-5285-4c8e-baf9-f0b2f75b33e0
v_j = [var_params[T][1] for T in 1:length(T_range)]

# ╔═╡ e5b21223-f7a5-4560-b03d-f015cdc99c41
sum.(v_j) ./ 2

# ╔═╡ ff8beae7-50e2-46cb-a38d-eb65a8cedaa4
w_j = [var_params[T][2] for T in 1:length(T_range)]

# ╔═╡ 757c41b7-e32b-4a67-b02d-c3f25acc266f
v_A = [i[1] for i in var_A]

# ╔═╡ 368ee8d7-84a8-4711-a9f9-f08d3c03c0ec
w_A = [i[2] for i in var_A]

# ╔═╡ b032bda5-4920-4546-a80f-5491b4044f45
v_B = [i[1] for i in var_B]

# ╔═╡ 5dd88a96-75d2-409a-abf8-960c09c74188
w_B = [i[2] for i in var_B]

# ╔═╡ 83ec77a6-fcf8-4025-ac6d-cd7819312427
begin
	vw_temp = plot(T_range, v_A, label = "H-A v", legend = :bottomright, linewidth = 2, tickfontsize = 12, labelfontsize = 12, legendfontsize = 12, xlabel = "T(K)", ylabel = "v & w (THz)", size = (550, 550), minorgrid = true, linestyle = :dashdot, color = theme_palette(:default)[1])
	plot!(T_range, w_A, label = "H-A w", linewidth = 2, linestyle = :dashdot, color = theme_palette(:default)[2])
	plot!(T_range, v_B, label = "H-B v", linewidth = 2, linestyle = :dash, color = theme_palette(:default)[3])
	plot!(T_range, w_B, label = "H-B w", linewidth = 2, linestyle = :dash, color = theme_palette(:default)[4])
	plot!(T_range, sum.(v_j) ./ length(v_j[1]), label = "multi v", linewidth = 2, linestyle = :solid, color = theme_palette(:default)[1])
	plot!(T_range, sum.(w_j) ./ length(w_j[1]), label = "multi w", linewidth = 2, linestyle = :solid, color = theme_palette(:default)[2])
end

# ╔═╡ 3a52ff01-110b-4755-b43a-acc5b3b870aa
savefig(vw_temp, "vw_temp.pdf")

# ╔═╡ 8fafb718-d95a-4ca7-921f-4abe1ab846ac
F_j = [multi_free_energy(v_j[T], w_j[T], T_range[T], 4.5, 0.12, (6.29e-10)^3, MAPI) for T in 1:length(T_range)]

# ╔═╡ 552b01e3-fd72-489a-aa2b-ae00b6d37c05
Hellwarth_energy_A = F.(v_A, w_A, β_hellwarth_A, α_hellwarth_A) .* 1e3 .* ħ .* 2π .* A_scheme .* 1e12 ./ eV

# ╔═╡ a465eeb5-5b8f-4202-bd58-651867c9d423
Hellwarth_energy_B = F.(v_B, w_B, β_hellwarth_B, α_hellwarth_B) .* 1e3 .* ħ .* 2π .* B_scheme .* 1e12 ./ eV

# ╔═╡ 6b60ab62-6706-4c33-ad7a-efb11803eb6b
begin
	free_energy_temp = plot(T_range, Hellwarth_energy_A, linewidth = 2, linestyle = :dashdot, label = "H-A", tickfontsize = 12, labelfontsize = 12, legendfontsize = 12, xlabel = "T(K)", ylabel = "Free energy (meV)", minorgrid = true, size = (550, 550), legend = :topright)
	plot!(T_range, Hellwarth_energy_B, linewidth = 2, linestyle = :dash, label = "H-B")
	plot!(T_range, F_j, linewidth = 2, linestyle = :solid, label = "multi")
end

# ╔═╡ 7ae9767d-d18f-4188-8e59-dea17c9d5ac4
savefig(free_energy_temp, "free_energy_temp.pdf")

# ╔═╡ 523af3a6-7cc9-41cf-9acb-9a1d25086969
begin
	multi_vw = DataFrames.DataFrame(
			temperature = T_range,
			beta = β_j_avg,
			v = v_j,
			w = w_j,
			F = F_j
	    )
	CSV.write("multi_vw.csv", multi_vw)
end

# ╔═╡ 8db8a4fd-fc04-4d45-8d05-8f69eff611df
begin
	A_data = DataFrames.DataFrame(
		alpha = α_hellwarth_A,
		efffreq = A_scheme,
		temp = T_range,
		beta = β_hellwarth_A,
		v = v_A,
		w = w_A,
		F = Hellwarth_energy_A
    )
	CSV.write("A_data.csv", A_data)
end

# ╔═╡ 0e921a7a-a07d-4afc-a991-e06d9be96e3f
begin
	B_data = DataFrames.DataFrame(
			alpha = [α_hellwarth_B for i in 1:length(T_range)],
			efffreq = [B_scheme for i in 1:length(T_range)],
			temp = T_range,
			beta = β_hellwarth_B,
			v = v_B,
			w = w_B,
			F = Hellwarth_energy_B
		)
	CSV.write("B_data.csv", B_data)
end

# ╔═╡ 2d4759ab-18d8-4368-a568-0da25474366c
function multi_conductivity(ν, β, α, v, w, ω, m_eff)
	z_integrand(t, β, ν) = (1 - exp(1im * 2π * ν * t)) * imag(cos(t - 1im * β / 2) / (  sinh(β / 2) * D_j(-1im * t, β, v, w)^(3/2)))
	z = 0.0
	for j in 1:length(ω)
		println("Photon frequency = $ν, Phonon mode frequency = $(ω[j] / 2π)")
		z += -1im * 2π * ν / length(ω) + 1im * 2 * α[j] * ω[j]^2 * quadgk(t -> z_integrand(t, β[j], ν / ω[j]), 0.0, Inf)[1] / (3 * √π * 2π * ν)
	end
	1 / z * eV * 100^2 / (m_eff * m_e * 1e12)
end

# ╔═╡ 2285006a-a115-414a-b788-b05fd6720bcd
ν_range = 0.01:0.01:3.51

# ╔═╡ 0fdcd049-8f06-4651-9864-a76647154680
begin
	σ_j = Array{ComplexF64}(undef, length(ν_range), length(T_range))
	for j in 1:length(T_range), i in 1:length(ν_range)
		println("ν: $(ν_range[i]) THz, T: $(T_range[j]) K")
		σ_j[i, j] = multi_conductivity(ν_range[i], β_j[:, j], α_j, v_j[1], w_j[1], phonon_mode_freqs .* 2π, 0.12)
		CSV.write("multi_conductivity_data.csv", DataFrames.DataFrame([[0.0, [i for i in T_range]...]'; [ν_range σ_j]], :auto))
	end
	σ_j
end

# ╔═╡ f09c4429-9415-49c1-b22d-a8867cbd4249
begin
	plot(ν_range, real.(σ_j), linewidth = 2, label = "Reσ", legend = :outerright, minorgrid = true)
	plot!(ν_range, imag.(σ_j), linewidth = 2, linestyle = :dash, label = "Imσ")
	plot!(ν_range, abs.(σ_j), linewidth = 2, linestyle = :dot, label = "|σ|")
end

# ╔═╡ 0178c326-eeca-4d7a-a111-58b589cba561
multi_contour_real = Plots.contourf(ν_range, T_range, real.(σ_j)', xlabel = "ν (THz)", ylabel = "T (K)", fill = cgrad(:thermal, rev = false, categorical = true), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, right_margin = 7Plots.mm, clims = (0, 600))

# ╔═╡ d2bbb4c9-98e5-4959-a36c-b3e48c37a03c
multi_contour_imag = Plots.contourf(ν_range, T_range, imag.(σ_j)', xlabel = "ν (THz)", ylabel = "T (K)", fill = cgrad(:thermal, rev = false, categorical = true), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, right_margin = 7Plots.mm, clims = (0, 600))

# ╔═╡ 0d78d4ef-6360-43b3-9952-3467410d8908
multi_contour_abs = Plots.contourf(ν_range, T_range, abs.(σ_j)', xlabel = "ν (THz)", ylabel = "T (K)", fill = cgrad(:thermal, rev = false, categorical = true), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, right_margin = 7Plots.mm, clims = (0, 600))

# ╔═╡ 8822d1c3-5bda-4e4e-9b0b-703c75d10816
begin
	savefig(multi_contour_real, "multi_contour_real.pdf")
	savefig(multi_contour_imag, "multi_contour_imag.pdf")
	savefig(multi_contour_abs, "multi_contour_abs.pdf")
end

# ╔═╡ 72d02e9b-d59e-498b-8947-3a86634b7eb4
multi_plot_temp_real = ridgeline(round.(ν_range, digits = 1), Int.(T_range), [i < 600 ? i : 600 for i in real.(σ_j)], jump = 5, shift = 10, xlabel = "ν (THz)", ylabel = "T (K)", size = (550, 550), linewidth = 0.5, palette = :thermal, fillalpha = 0.75, step = 5)

# ╔═╡ 0916b2de-ff08-4c60-9bb1-70efbea6f6f6
multi_plot_temp_imag = ridgeline(round.(ν_range, digits = 1), Int.(T_range), [i < 600 ? i : 600 for i in imag.(σ_j)], jump = 5, shift = 10, xlabel = "ν (THz)", ylabel = "T (K)", size = (550, 550), linewidth = 0.5, palette = :thermal, fillalpha = 0.75, step = 5)

# ╔═╡ dfd78449-6db6-4f6b-9ecf-00a202ad685d
multi_plot_temp_abs = ridgeline(round.(ν_range, digits = 1), Int.(T_range), [i < 600 ? i : 600 for i in abs.(σ_j)], jump = 5, shift = 10, xlabel = "ν (THz)", ylabel = "T (K)", size = (550, 550), linewidth = 0.5, palette = :thermal, fillalpha = 0.75, step = 5)

# ╔═╡ 7d101426-86d9-4925-b4f8-52c2d5415491
begin
	savefig(multi_plot_temp_real, "multi_plot_temp_real.pdf")
	savefig(multi_plot_temp_imag, "multi_plot_temp_imag.pdf")
	savefig(multi_plot_temp_abs, "multi_plot_temp_abs.pdf")
end

# ╔═╡ ad91591f-1ccf-4347-8156-1d0e2e090fdb
begin
	σ_hellwarth_A = Array{ComplexF64}(undef, length(ν_range), length(T_range))
	for j in 1:length(T_range)
		for i in 1:length(ν_range)
			println("ν: $(ν_range[i]) THz, T: $(T_range[j]) K")
			σ_hellwarth_A[i, j] = multi_conductivity(ν_range[i], β_hellwarth_A[j], α_hellwarth_A[j], v_A[1], w_A[1], A_scheme[j] * 2π, 0.12)
		end
		CSV.write("A_conductivity_data.csv", DataFrames.DataFrame([[0.0, [i for i in T_range]...]'; [ν_range σ_hellwarth_A]], :auto))
	end
	σ_hellwarth_A
end

# ╔═╡ 2c4a84f5-6839-4ea4-bc9c-4eb2f8e6281b
begin
	term = 40
	plot(ν_range, real.(σ_hellwarth_A)[:, term], linewidth = 2, label = "Reσ", legend = :outerright, minorgrid = true, ylims = (0, 500))
	plot!(ν_range, imag.(σ_hellwarth_A)[:, term], linewidth = 2, linestyle = :dash, label = "Imσ")
	plot!(ν_range, abs.(σ_hellwarth_A)[:, term], linewidth = 2, linestyle = :dot, label = "|σ|")
end

# ╔═╡ ee677bc3-2e9e-4e35-a66c-4c6052c4e32e
A_contour_real = Plots.contourf(ν_range, T_range, real.(σ_hellwarth_A)', xlabel = "ν (THz)", ylabel = "T (K)", fill = cgrad(:thermal, rev = false, categorical = true), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, right_margin = 7Plots.mm, clims = (0, 600))

# ╔═╡ 87154df5-6996-4b4b-9f63-45606205b945
A_contour_imag = Plots.contourf(ν_range, T_range, imag.(σ_hellwarth_A)', xlabel = "ν (THz)", ylabel = "T (K)", fill = cgrad(:thermal, rev = false, categorical = true), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, right_margin = 7Plots.mm, clims = (-6, 600))

# ╔═╡ f4f84543-d7b8-43b2-93cb-0de8d7fa0a30
A_contour_abs = Plots.contourf(ν_range, T_range, abs.(σ_hellwarth_A)', xlabel = "ν (THz)", ylabel = "T (K)", fill = cgrad(:thermal, rev = false, categorical = true, scale = :log), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, right_margin = 7Plots.mm, clims = (0, 600))

# ╔═╡ 005705cd-a857-431a-8956-8ba0035c7c05
begin
	savefig(A_contour_real, "A_contour_real.pdf")
	savefig(A_contour_imag, "A_contour_imag.pdf")
	savefig(A_contour_abs, "A_contour_abs.pdf")
end

# ╔═╡ 77a6559f-6c40-4673-9a8d-9e0c8fe3476f
A_plot_temp_real = ridgeline(round.(ν_range, digits = 1), Int.(T_range), [i < 600 ? i : 600 for i in abs.(σ_hellwarth_A)], jump = 5, shift = 10, xlabel = "ν (THz)", ylabel = "T (K)", size = (550, 550), linewidth = 0.5, palette = :thermal, fillalpha = 0.75, step = 5)

# ╔═╡ 498c610b-0646-4a84-ac40-ab8c6efc5278
A_plot_temp_imag = ridgeline(round.(ν_range, digits = 1), Int.(T_range), [i < 600 ? i : 600 for i in imag.(σ_hellwarth_A)], jump = 5, shift = 10, xlabel = "ν (THz)", ylabel = "T (K)", size = (550, 550), linewidth = 0.5, palette = :thermal, fillalpha = 0.75, step = 5)

# ╔═╡ bc47027d-21a0-4a02-914d-aba2381d2e01
A_plot_temp_abs = ridgeline(round.(ν_range, digits = 1), Int.(T_range), [i < 600 ? i : 600 for i in abs.(σ_hellwarth_A)], jump = 5, shift = 10, xlabel = "ν (THz)", ylabel = "T (K)", size = (550, 550), linewidth = 0.5, palette = :thermal, fillalpha = 0.75, step = 5)

# ╔═╡ 1e902322-d319-42c5-889c-c3048515674c
begin
	savefig(A_plot_temp_real, "A_plot_temp_real.pdf")
	savefig(A_plot_temp_imag, "A_plot_temp_imag.pdf")
	savefig(A_plot_temp_abs, "A_plot_temp_abs.pdf")
end

# ╔═╡ 905efcbf-b3f5-4da0-a45c-3204cd69c0ff
begin
	σ_hellwarth_B = Array{ComplexF64}(undef, length(ν_range), length(T_range))
	for j in 1:length(T_range)
		for i in 1:length(ν_range)
			println("ν: $(ν_range[i]) THz, T: $(T_range[j]) K")
			σ_hellwarth_B[i, j] = multi_conductivity(ν_range[i], β_hellwarth_B[j], α_hellwarth_B, v_B[1], w_B[1], B_scheme * 2π, 0.12)
		end
		CSV.write("B_conductivity_data.csv", DataFrames.DataFrame([[0.0, [i for i in T_range]...]'; [ν_range σ_hellwarth_B]], :auto))
	end
	σ_hellwarth_B
end

# ╔═╡ 38a23177-633a-4b2d-8083-5e8afd74796a
begin
	term_B = 2
	plot(ν_range, real.(σ_hellwarth_B)[:, term_B], linewidth = 2, label = "Reσ", legend = :outerright, minorgrid = true, ylims = (0, 600))
	plot!(ν_range, imag.(σ_hellwarth_B)[:, term_B], linewidth = 2, linestyle = :dash, label = "Imσ")
	plot!(ν_range, abs.(σ_hellwarth_B)[:, term_B], linewidth = 2, linestyle = :dot, label = "|σ|")
end

# ╔═╡ 64821290-9bbe-4a7b-af4c-79825f6d4015
B_contour_real = Plots.contourf(ν_range, T_range, real.(σ_hellwarth_B)', xlabel = "ν (THz)", ylabel = "T (K)", fill = cgrad(:thermal, rev = false, categorical = true), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, right_margin = 7Plots.mm, clims = (0, 600))

# ╔═╡ d292b963-3fd1-4d05-b72b-d71b66c5a919
B_contour_imag = Plots.contourf(ν_range, T_range, imag.(σ_hellwarth_B)', xlabel = "ν (THz)", ylabel = "T (K)", fill = cgrad(:thermal, rev = false, categorical = true), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, right_margin = 7Plots.mm, clims = (-6, 600))

# ╔═╡ bfdc8a9b-9522-4f48-8b17-da79447994a8
B_contour_abs = Plots.contourf(ν_range, T_range, abs.(σ_hellwarth_B)', xlabel = "ν (THz)", ylabel = "T (K)", fill = cgrad(:thermal, rev = false, categorical = true), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, right_margin = 7Plots.mm, clims = (0, 600))

# ╔═╡ 273c8ea1-0f7c-42b4-8954-67671d157181
begin
	savefig(B_contour_real, "B_contour_real.pdf")
	savefig(B_contour_imag, "B_contour_imag.pdf")
	savefig(B_contour_abs, "B_contour_abs.pdf")
end

# ╔═╡ b96343f4-bb61-4f4b-9ec9-ad2f37a96cfa
B_plot_temp_real = ridgeline(round.(ν_range, digits = 1), Int.(T_range), [i < 600 ? i : 600 for i in real.(σ_hellwarth_A)], jump = 5, shift = 10, xlabel = "ν (THz)", ylabel = "T (K)", size = (550, 550), linewidth = 0.5, palette = :thermal, fillalpha = 0.75, step = 5)

# ╔═╡ b4ee14e9-a89f-4a29-b6ed-d777fe477b9a
B_plot_temp_imag = ridgeline(round.(ν_range, digits = 1), Int.(T_range), [i < 600 ? i : 600 for i in imag.(σ_hellwarth_B)], jump = 5, shift = 10, xlabel = "ν (THz)", ylabel = "T (K)", size = (550, 550), linewidth = 0.5, palette = :thermal, fillalpha = 0.75, step = 5)

# ╔═╡ f711553f-e567-437f-9c38-79e644bd18b9
B_plot_temp_abs = ridgeline(round.(ν_range, digits = 1), Int.(T_range), [i < 600 ? i : 600 for i in abs.(σ_hellwarth_B)], jump = 5, shift = 10, xlabel = "ν (THz)", ylabel = "T (K)", size = (550, 550), linewidth = 0.5, palette = :thermal, fillalpha = 0.75, step = 5)

# ╔═╡ 26e1a5da-6bdb-49b9-973b-768d901a2a60
begin
	savefig(B_plot_temp_real, "B_plot_temp_real.pdf")
	savefig(B_plot_temp_imag, "B_plot_temp_imag.pdf")
	savefig(B_plot_temp_abs, "B_plot_temp_abs.pdf")
end

# ╔═╡ Cell order:
# ╠═29d62408-0629-4f39-909e-a6ba8696e955
# ╠═1eb559b3-dbce-4caf-800f-96a2495dee1b
# ╠═77d14c3f-bcaf-432c-9005-f5ec184d84b7
# ╠═ba2dadc8-d9a2-452e-9f94-406ccfc7404b
# ╠═c05cd489-68c8-418e-bbc2-a33bc42781f0
# ╠═1d300e7d-5efb-4ad2-a826-c3f8a86083d8
# ╠═6a99cd13-68cb-474c-96af-7d78ee4e7095
# ╠═49663b33-1604-430b-a815-9eb985802987
# ╠═314af942-8922-40d4-a92b-f75b715c276c
# ╠═6ff68ba5-df0c-42a0-8cf1-dedeced77e8d
# ╠═35390bdf-e260-4566-8a97-7a0b5a3ef07a
# ╠═e0ce47a5-368a-4da7-83fc-9f921af0f747
# ╠═053d717a-2b63-40e7-9dc0-487e3892de98
# ╠═86fe8a15-4d9e-48fa-8387-3cfaff7aa182
# ╠═16e7d781-bc61-44b3-84a2-e8c8892bf773
# ╠═a78017d7-8835-4a8f-8dd3-f3a4a50eb37b
# ╠═d7346684-6b68-4dfc-81da-12c435c4acf0
# ╠═798a4f7a-a2a0-48de-82ad-f85b4199496a
# ╠═ec0744c0-2530-11ec-04c3-1b5d83413b78
# ╠═f3b91e16-26d3-45c8-b28f-193afac15d7f
# ╠═a941ccc6-90cb-448a-91bb-f80233b630f8
# ╠═1372f247-4622-4388-b458-89aa6e695388
# ╠═c28f6d35-c3ed-4984-9c36-b282fb50001b
# ╠═49480256-2dcc-4d10-a5a1-4202b243cfd2
# ╠═4966ecf1-ef39-4da3-a8a5-9e47a6c6b0ec
# ╠═3321246b-76bd-4678-a5ae-332c3543a1ef
# ╠═b2944c8a-158e-461b-a83c-842925587179
# ╠═fe4e5ef5-833a-4b61-9ea9-92fa52d7a21f
# ╠═8350c7a2-e9fd-4312-9561-44e2524734af
# ╠═db2b4169-0737-44c1-9af5-aa2f2403f5db
# ╠═2b6d00e1-317b-41cc-9489-14a1450012f4
# ╠═3e1418c2-01a1-4651-8029-59f53264193d
# ╠═d4f392ef-968a-439c-ba53-5707464446e7
# ╠═b3c94ea5-d152-4705-9a30-f8d3e4da0874
# ╠═3608af5e-7b44-4d50-9770-5e8afd4fca28
# ╠═1716f263-f8ad-40a1-a805-d8b7ea6f38df
# ╠═da629bec-b53c-4acf-a0ab-6057c3b973a9
# ╠═0641f475-0f35-4f5c-abf0-85ec56f4e30c
# ╠═c49a4dc3-2546-472f-9fde-b59ec3246237
# ╠═0e3cf449-6e2a-4353-853e-58d15693743d
# ╠═16040051-57d5-4c6c-aa13-8280ae964261
# ╠═32f32534-2a6b-4675-8508-890579058e07
# ╠═de89f0aa-21fb-42d2-9a75-5b45938b1e65
# ╠═a74836cf-f9ab-41d6-9f01-5066f1afac34
# ╠═2ea4fc88-5285-4c8e-baf9-f0b2f75b33e0
# ╠═e5b21223-f7a5-4560-b03d-f015cdc99c41
# ╠═ff8beae7-50e2-46cb-a38d-eb65a8cedaa4
# ╠═757c41b7-e32b-4a67-b02d-c3f25acc266f
# ╠═368ee8d7-84a8-4711-a9f9-f08d3c03c0ec
# ╠═b032bda5-4920-4546-a80f-5491b4044f45
# ╠═5dd88a96-75d2-409a-abf8-960c09c74188
# ╠═83ec77a6-fcf8-4025-ac6d-cd7819312427
# ╠═3a52ff01-110b-4755-b43a-acc5b3b870aa
# ╠═8fafb718-d95a-4ca7-921f-4abe1ab846ac
# ╠═552b01e3-fd72-489a-aa2b-ae00b6d37c05
# ╠═a465eeb5-5b8f-4202-bd58-651867c9d423
# ╠═6b60ab62-6706-4c33-ad7a-efb11803eb6b
# ╠═7ae9767d-d18f-4188-8e59-dea17c9d5ac4
# ╠═523af3a6-7cc9-41cf-9acb-9a1d25086969
# ╠═8db8a4fd-fc04-4d45-8d05-8f69eff611df
# ╠═0e921a7a-a07d-4afc-a991-e06d9be96e3f
# ╠═2d4759ab-18d8-4368-a568-0da25474366c
# ╠═2285006a-a115-414a-b788-b05fd6720bcd
# ╠═0fdcd049-8f06-4651-9864-a76647154680
# ╠═f09c4429-9415-49c1-b22d-a8867cbd4249
# ╠═0178c326-eeca-4d7a-a111-58b589cba561
# ╠═d2bbb4c9-98e5-4959-a36c-b3e48c37a03c
# ╠═0d78d4ef-6360-43b3-9952-3467410d8908
# ╠═8822d1c3-5bda-4e4e-9b0b-703c75d10816
# ╠═72d02e9b-d59e-498b-8947-3a86634b7eb4
# ╠═0916b2de-ff08-4c60-9bb1-70efbea6f6f6
# ╠═dfd78449-6db6-4f6b-9ecf-00a202ad685d
# ╠═7d101426-86d9-4925-b4f8-52c2d5415491
# ╠═ad91591f-1ccf-4347-8156-1d0e2e090fdb
# ╠═2c4a84f5-6839-4ea4-bc9c-4eb2f8e6281b
# ╠═ee677bc3-2e9e-4e35-a66c-4c6052c4e32e
# ╠═87154df5-6996-4b4b-9f63-45606205b945
# ╠═f4f84543-d7b8-43b2-93cb-0de8d7fa0a30
# ╠═005705cd-a857-431a-8956-8ba0035c7c05
# ╠═77a6559f-6c40-4673-9a8d-9e0c8fe3476f
# ╠═498c610b-0646-4a84-ac40-ab8c6efc5278
# ╠═bc47027d-21a0-4a02-914d-aba2381d2e01
# ╠═1e902322-d319-42c5-889c-c3048515674c
# ╠═905efcbf-b3f5-4da0-a45c-3204cd69c0ff
# ╠═38a23177-633a-4b2d-8083-5e8afd74796a
# ╠═64821290-9bbe-4a7b-af4c-79825f6d4015
# ╠═d292b963-3fd1-4d05-b72b-d71b66c5a919
# ╠═bfdc8a9b-9522-4f48-8b17-da79447994a8
# ╠═273c8ea1-0f7c-42b4-8954-67671d157181
# ╠═b96343f4-bb61-4f4b-9ec9-ad2f37a96cfa
# ╠═b4ee14e9-a89f-4a29-b6ed-d777fe477b9a
# ╠═f711553f-e567-437f-9c38-79e644bd18b9
# ╠═26e1a5da-6bdb-49b9-973b-768d901a2a60


text:
using ConditionalMoments
using Test

## Observations

const X = [1.,1.,4.,5.,6.,7,1.,2.,4.,5.,6.,7,3.,1.,4.,5.,6.,7]
const dt = 0.1
const ob = Observation(X,dt)

@testset "Observations" begin
    @test ob[3] == 4.0
    @test ob[ob .> 3] == X[X .> 3]
    @test ob[(ob .> 3)'] == X[X .> 3]
    @test ob[:] == X
    @test ob.dt == dt
    @test ob.npoints == 18

    @test_throws MethodError Observation([1.,1.,4.,5.,6.,7+2im],dt)
end

## Moments

const τ_indices = [1,3,4]
const bin_edges = LinRange(1.5,6.5,5)
const histogramSettings = HistogramSettings(τ_indices,bin_edges)

@testset "Moment settings" begin
    @test histogramSettings.n_time_points == 3
    @test histogramSettings.n_evaluation_points == 5

    @test_throws MethodError HistogramSettings([1,2,3.1],bin_edges)
    @test_throws MethodError HistogramSettings([1,2,3+1im],bin_edges)
end

const moment_sol = MomentSolution(ob, histogramSettings)
@testset "Moment calculations" begin
    @test moment_sol.observation.npoints == 18
    @test size(moment_sol.moments.M1) == (3,4)
    @test size(moment_sol.moments.M2) == (3,4)
    #@test size(moment_sol.errors) == (3,4)
    @test moment_sol.moments.M1[2,4] == -4.5
    @test moment_sol.moments.M2[2,4] == 0.25
end

## Drift and noise functions
const model1 = WhiteNoiseProcess(moment_sol)
const model2 = WhiteNoiseProcess(moment_sol,i=2)

@testset "Drift and noise calculations" begin
    @test size(model1.drift) == (4,)
    @test size(model1.noise) == (4,)
    #@test size(model1.errors) == (4,)
    @test size(model2.drift) == (4,)
    @test size(model2.noise) == (4,)
    #@test size(model2.errors) == (4,)
    @test ~isnan(model1.drift[4])
    @test ~isnan(model2.drift[4])
    @test ~isnan(model1.noise[4])
    @test ~isnan(model2.noise[4])
end

## Ensemble observations

const X_ensemble = repeat(X,1,2)

const ensemble_ob = EnsembleObservation(X_ensemble,dt)

@testset "Ensemble observations" begin
    @test ensemble_ob[3] == 4.0
    @test ensemble_ob[:,1] == X
    @test ensemble_ob[:,2] == X
    @test ensemble_ob[3,:] == [4.0,4.0]
    @test ensemble_ob[:] == X_ensemble
    @test ensemble_ob.dt == dt
    @test ensemble_ob.npoints == 18
    @test ensemble_ob.nsample == 2

    @test_throws MethodError Observation([1.,1.,4.,5.,6.,7+2im],dt)
end

## Ensemble moments

const ensemble_sol = MomentSolution(ensemble_ob, histogramSettings)
@testset "Mapping moments" begin
    @test ensemble_sol.observation.npoints == 18
    @test ensemble_sol.observation.nsample == 2
    @test size(ensemble_sol.moments.M1) == (3,4)
    @test size(ensemble_sol.moments.M2) == (3,4)
    #@test size(ensemble_moments.errors) == (3,4)
    @test ensemble_sol.moments.M1[2,4] == -4.5
    @test ensemble_sol.moments.M2[2,4] == 0.25
end

## Larger test

const t = 0:0.001:200
const dt_large = t.step.hi
const X_large = exp.(-t) .*
    cumsum(exp.(t) .* (sqrt(dt_large)*pushfirst!(randn(length(t)-1),0)))
const ob_large = Observation(X_large,dt_large)

@testset "Larger observation" begin
    @test ob_large[ob_large .> 3] == X_large[X_large .> 3]
    @test ob_large[(ob_large .> 3)'] == X_large[X_large .> 3]
    @test ob_large[:] == X_large
    @test ob_large.dt == dt_large
    @test ob_large.npoints == length(t)
end

const τ_indices_large = 1:20
const τ_vec = dt_large*τ_indices_large
const bin_edges_large = LinRange(-1,1,10)
const histogramSettings_large =
        HistogramSettings(τ_indices_large,bin_edges_large)
const moments_large = MomentSolution(ob_large, histogramSettings_large)

const M1_τ = moments_large.moments.M1 ./ τ_vec
const M2_τ = moments_large.moments.M2 ./ (2*τ_vec)

@testset "Larger moments" begin
    @test moments_large.observation.npoints == length(t)
    expected_size = (length(τ_indices_large),length(bin_edges_large)-1)
    @test size(moments_large.moments.M1) == expected_size
    @test size(moments_large.moments.M2) == expected_size
    #@test size(moments_large.errors) == expected_size
end

#NOTE: Not type stable
const model1_large = WhiteNoiseProcess(moments_large)
const model2_large = WhiteNoiseProcess(moments_large,i=2)

@testset "Larger drift and noise functions" begin
    expected_length = length(bin_edges_large)-1
    @test size(model1_large.drift) == (expected_length,)
    @test size(model1_large.noise) == (expected_length,)
    #@test size(model1_large.errors) == (expected_length,)
    @test size(model2_large.drift) == (expected_length,)
    @test size(model2_large.noise) == (expected_length,)
    #@test size(model2_large.errors) == (expected_length,)
    @test ~isnan(model1_large.drift[4])
    @test ~isnan(model2_large.drift[4])
    @test ~isnan(model1_large.noise[4])
    @test ~isnan(model2_large.noise[4])
end


text:
"""
PolaronMobility.jl - https://github.com/jarvist/PolaronMobility.jl
  Codes by Jarvist Moore Frost, 2017-2020
  Calculate Polaron Mobility - by a Osaka/Hellwarth variational solution to the Feynman model
"""
module PolaronMobility

# These codes were developed with Julia 0.5.0 - Julia 0.6.2, and require the Optim and Plots packages.

export Polaron # Type to hold the data
export frohlichalpha, feynmanvw, F, polaronmobility, savepolaron, plotpolaron
export HellwarthBScheme, HellwarthAScheme
export ImX

export frohlichPartial, IRtoDielectric, IRtoalpha, DielectricFromIRmode
export Hellwarth1999mobilityRHS

##### load in library routines... #####
# stdlib
if VERSION > v"0.7.0-DEV.2005"
    using LinearAlgebra
    using Printf
end
# one-dimensional numerical integration in Julia using adaptive Gauss-Kronrod quadrature
import QuadGK.quadgk

# Using the powerful Julia Optim package to optimise the variational parameters
using Optim

# Physical constants
const hbar = const ħ = 1.05457162825e-34;          # kg m2 / s
const eV = const q = const ElectronVolt = 1.602176487e-19;                         # kg m2 / s2
const me=MassElectron = 9.10938188e-31;                          # kg
const Boltzmann = const kB =  1.3806504e-23;                  # kg m2 / K s2
const ϵ_0 = 8.854E-12 #Units: C2N−1m−2, permittivity of free space

const amu = 1.660_539_066_60e-27 # kg

include("types.jl")            # Polaron types
include("FeynmanTheory.jl")    # Actions + variational functions
include("HellwarthTheory.jl")  # multimode -> equivalent mode.
include("MobilityTheories.jl") # Main polaronmobility function
include("Susceptibility.jl")   # ImX calculation
include("OedipusRex.jl")       # Optical Absorption
include("MultipleBranches.jl")  # Oct 2019 extension to multiple phonon branches

end # module


text:
function UnipartiteNetwork(A::M) where {M<:AbstractMatrix{Bool}}
  check_unipartiteness(A)
  S = "s".*string.(1:size(A,1))
  UnipartiteNetwork(A, S)
end

function UnipartiteNetwork(A::M, S::Vector{NT}) where {M<:AbstractMatrix{Bool}, NT}
  is_valid_species(NT) || throw(ArgumentError("Not allowed"))
  check_unipartiteness(A, S)
  UnipartiteNetwork{Bool,NT}(A, S)
end

function BipartiteNetwork(A::M) where {M<:AbstractMatrix{Bool}}
  check_bipartiteness(A)
  T = "t".*string.(1:size(A,1))
  B = "b".*string.(1:size(A,2))
  BipartiteNetwork{Bool,eltype(T)}(A, T, B)
end

function BipartiteNetwork(A::M, T::Vector{NT}, B::Vector{NT}) where {M<:AbstractMatrix{Bool}, NT}
  is_valid_species(NT) || throw(ArgumentError("Not allowed"))
  check_bipartiteness(A, T, B)
  BipartiteNetwork{Bool,eltype(T)}(A, T, B)
end

function BipartiteProbabilisticNetwork(A::Matrix{IT}) where {IT<:AbstractFloat}
  check_bipartiteness(A)
  check_probability_values(A)
  T = "t".*string.(1:size(A,1))
  B = "b".*string.(1:size(A,2))
  BipartiteProbabilisticNetwork{IT,eltype(T)}(A, T, B)
end

function BipartiteProbabilisticNetwork(A::Matrix{IT}, T::Vector{NT}, B::Vector{NT}) where {IT<:AbstractFloat, NT}
  is_valid_species(NT) || throw(ArgumentError("Not allowed"))
  check_bipartiteness(A, T, B)
  check_probability_values(A)
  BipartiteProbabilisticNetwork{IT,NT}(A, T, B)
end

function BipartiteQuantitativeNetwork(A::Matrix{IT}) where {IT <: Number}
  check_bipartiteness(A)
  T = "t".*string.(1:size(A,1))
  B = "b".*string.(1:size(A,2))
  BipartiteQuantitativeNetwork{IT,eltype(T)}(A, T, B)
end

function BipartiteQuantitativeNetwork(A::Matrix{IT}, T::Vector{NT}, B::Vector{NT}) where {IT<:Number,NT}
  is_valid_species(NT) || throw(ArgumentError("Not allowed"))
  check_bipartiteness(A, T, B)
  BipartiteQuantitativeNetwork{IT,NT}(A, T, B)
end

function UnipartiteQuantitativeNetwork(A::Matrix{IT}) where {IT<:Number}
  check_unipartiteness(A)
  S = "s".*string.(1:size(A,1))
  UnipartiteQuantitativeNetwork{IT,eltype(S)}(A, S)
end

function UnipartiteQuantitativeNetwork(A::Matrix{IT}, S::Vector{NT}) where {IT<:Number,NT}
  is_valid_species(NT) || throw(ArgumentError("Not allowed"))
  check_unipartiteness(A, S)
  UnipartiteQuantitativeNetwork{IT,NT}(A, S)
end

function UnipartiteProbabilisticNetwork(A::Matrix{IT}) where {IT<:AbstractFloat}
  check_unipartiteness(A)
  check_probability_values(A)
  S = "s".*string.(1:size(A,1))
  UnipartiteProbabilisticNetwork{IT,eltype(S)}(A, S)
end

function UnipartiteProbabilisticNetwork(A::Matrix{IT}, S::Vector{NT}) where {IT<:AbstractFloat,NT}
  is_valid_species(NT) || throw(ArgumentError("Not allowed"))
  check_unipartiteness(A, S)
  check_probability_values(A)
  UnipartiteProbabilisticNetwork{IT,NT}(A, S)
end


text:
msg = "Hello World"
if length(ARGS) > 0
    msg = msg * " & $(first(ARGS))!"
else
    msg = msg * "!"
end

println(msg)

max_stars_repo_name: MiroK/julia-tools-meetup

max_issues_repo_name: MiroK/julia-tools-meetup

max_forks_repo_name: MiroK/julia-tools-meetup


text:
# define common PETSc constants
# this excludes configurable constants (e.g. PetscScalar) which are set in lib.jl

const PetscErrorCode      = Cint

struct PetscError <: Exception
    code::PetscErrorCode
end

macro chk(expr)
    :((errcode = $(esc(expr))) == 0 || throw(PetscError(errcode)))
end

const PETSC_DEFAULT = -2
const PETSC_DECIDE = PETSC_DETERMINE = -1

@enum DMBoundaryType begin
    DM_BOUNDARY_NONE
    DM_BOUNDARY_GHOSTED
    DM_BOUNDARY_MIRROR
    DM_BOUNDARY_PERIODIC
    DM_BOUNDARY_TWIST
end

@enum DMStagStencilType begin
    DMSTAG_STENCIL_NONE
    DMSTAG_STENCIL_STAR
    DMSTAG_STENCIL_BOX
end

@enum PetscBool PETSC_FALSE PETSC_TRUE

@enum PetscDataType begin
    PETSC_DATATYPE_UNKNOWN  = 0
    PETSC_DOUBLE            = 1
    PETSC_COMPLEX           = 2
    PETSC_LONG              = 3
    PETSC_SHORT             = 4
    PETSC_FLOAT             = 5
    PETSC_CHAR              = 6
    PETSC_BIT_LOGICAL       = 7
    PETSC_ENUM              = 8
    PETSC_BOOL              = 9
    PETSC_FLOAT128          = 10
    PETSC_OBJECT            = 11
    PETSC_FUNCTION          = 12
    PETSC_STRING            = 13
    PETSC___FP16            = 14
    PETSC_STRUCT            = 15
    PETSC_INT               = 16
    PETSC_INT64             = 17
end

const Petsc64bitInt       = Int64
const PetscLogDouble      = Cdouble

@enum InsertMode NOT_SET_VALUES INSERT_VALUES ADD_VALUES MAX_VALUES MIN_VALUES INSERT_ALL_VALUES ADD_ALL_VALUES INSERT_BC_VALUES ADD_BC_VALUES

@enum NormType begin
    NORM_1 = 0
    NORM_2 = 1
    NORM_FROBENIUS = 2
    NORM_INFINITY = 3
    NORM_1_AND_2 = 4
end

@enum MatAssemblyType MAT_FLUSH_ASSEMBLY=1 MAT_FINAL_ASSEMBLY=0

@enum MatFactorType begin
    MAT_FACTOR_NONE     = 0
    MAT_FACTOR_LU       = 1
    MAT_FACTOR_CHOLESKY = 2
    MAT_FACTOR_ILU      = 3
    MAT_FACTOR_ICC      = 4
    MAT_FACTOR_ILUDT    = 5
end

@enum PetscViewerFormat begin
    PETSC_VIEWER_DEFAULT
    PETSC_VIEWER_ASCII_MATLAB
    PETSC_VIEWER_ASCII_MATHEMATICA
    PETSC_VIEWER_ASCII_IMPL
    PETSC_VIEWER_ASCII_INFO
    PETSC_VIEWER_ASCII_INFO_DETAIL
    PETSC_VIEWER_ASCII_COMMON
    PETSC_VIEWER_ASCII_SYMMODU
    PETSC_VIEWER_ASCII_INDEX
    PETSC_VIEWER_ASCII_DENSE
    PETSC_VIEWER_ASCII_MATRIXMARKET
    PETSC_VIEWER_ASCII_VTK
    PETSC_VIEWER_ASCII_VTK_CELL
    PETSC_VIEWER_ASCII_VTK_COORDS
    PETSC_VIEWER_ASCII_PCICE
    PETSC_VIEWER_ASCII_PYTHON
    PETSC_VIEWER_ASCII_FACTOR_INFO
    PETSC_VIEWER_ASCII_LATEX
    PETSC_VIEWER_ASCII_XML
    PETSC_VIEWER_ASCII_GLVIS
    PETSC_VIEWER_ASCII_CSV
    PETSC_VIEWER_DRAW_BASIC
    PETSC_VIEWER_DRAW_LG
    PETSC_VIEWER_DRAW_LG_XRANGE
    PETSC_VIEWER_DRAW_CONTOUR
    PETSC_VIEWER_DRAW_PORTS
    PETSC_VIEWER_VTK_VTS
    PETSC_VIEWER_VTK_VTR
    PETSC_VIEWER_VTK_VTU
    PETSC_VIEWER_BINARY_MATLAB
    PETSC_VIEWER_NATIVE
    PETSC_VIEWER_HDF5_PETSC
    PETSC_VIEWER_HDF5_VIZ
    PETSC_VIEWER_HDF5_XDMF
    PETSC_VIEWER_HDF5_MAT
    PETSC_VIEWER_NOFORMAT
    PETSC_VIEWER_LOAD_BALANCE
end

@enum MatOperation begin
    MATOP_SET_VALUES=0
    MATOP_GET_ROW=1
    MATOP_RESTORE_ROW=2
    MATOP_MULT=3
    MATOP_MULT_ADD=4
    MATOP_MULT_TRANSPOSE=5
    MATOP_MULT_TRANSPOSE_ADD=6
    MATOP_SOLVE=7
    MATOP_SOLVE_ADD=8
    MATOP_SOLVE_TRANSPOSE=9
    MATOP_SOLVE_TRANSPOSE_ADD=10
end

@enum DMStagStencilLocation begin
    DMSTAG_NULL_LOCATION=0
    DMSTAG_BACK_DOWN_LEFT
    DMSTAG_BACK_DOWN
    DMSTAG_BACK_DOWN_RIGHT
    DMSTAG_BACK_LEFT
    DMSTAG_BACK
    DMSTAG_BACK_RIGHT
    DMSTAG_BACK_UP_LEFT
    DMSTAG_BACK_UP
    DMSTAG_BACK_UP_RIGHT
    DMSTAG_DOWN_LEFT
    DMSTAG_DOWN
    DMSTAG_DOWN_RIGHT
    DMSTAG_LEFT
    DMSTAG_ELEMENT
    DMSTAG_RIGHT
    DMSTAG_UP_LEFT
    DMSTAG_UP
    DMSTAG_UP_RIGHT
    DMSTAG_FRONT_DOWN_LEFT
    DMSTAG_FRONT_DOWN
    DMSTAG_FRONT_DOWN_RIGHT
    DMSTAG_FRONT_LEFT
    DMSTAG_FRONT
    DMSTAG_FRONT_RIGHT
    DMSTAG_FRONT_UP_LEFT
    DMSTAG_FRONT_UP
    DMSTAG_FRONT_UP_RIGHT
end

@enum DMBoundaryType begin
  DM_BOUNDARY_NONE=0
  DM_BOUNDARY_GHOSTED=1
  DM_BOUNDARY_MIRROR=2
  DM_BOUNDARY_PERIODIC=3
  DM_BOUNDARY_TWIST=4
end

@enum DMDAStencilType begin
  DMDA_STENCIL_STAR = 0
  DMDA_STENCIL_BOX = 1
end


text:
using ThreeJS
import ThreeJS
using Compat

function ConvertToTuple(p,t)
    points = Array(Tuple{Float64,Float64,Float64},size(p,2))
    for i in 1:size(p,2)
        points[i] = tuple(p[:,i]...)
    end

    faces = Array(Tuple{Int64,Int64,Int64},size(t,2))
    for i in 1:size(t,2)
        faces[i] = tuple(t[:,i]...)
    end

    return points,faces
end

function main(window)

    push!(window.assets,("ThreeJS","threejs"))
    push!(window.assets,"widgets")

    convert(Int,round(3.6))

    vbox(
         vbox(
          outerdiv() <<
          (
           initscene() <<
           [
            # For presentations
             ThreeJS.mesh(0.0, 0.0, 0.0) << [ geometry(ConvertToTuple(points,faces)...), material(Dict(:color=>"white",:kind=>"basic")) ],
             ThreeJS.mesh(0.0, 0.0, 0.0) << [ geometry(ConvertToTuple(points,faces)...), material(Dict(:color=>"black",:wireframe=>true,:kind=>"basic")) ],
            #ThreeJS.mesh(0.0, 0.0, 0.0) << [ geometry(ConvertToTuple(points,faces)...), material(Dict(:color=>"red",:kind=>"normal")) ],
    #            pointlight(0.,10.,0.),
            camera(10.0, 0., 0.),
           ]
          )
          )
    )
end


text:
using Gadfly, Dates

set_default_plot_size(6inch, 3inch)

dates = Array{Date}(undef, 40)
dates[1] = Date(2019,4,19)
for i=2:length(dates)
    dates[i] = today()+Day(i)
end

plot(x=dates,
     y=[0.316371, 0.814995, 2.06458, 3.07804, 3.11306, 3.42541, 4.47924, 7.12984, 7.43285, 8.37889, 8.43873, 9.44133, 10.0552, 10.0665, 11.1433, 11.1758, 11.2476, 11.5081, 12.5083, 13.8607, 16.9887, 17.4843, 18.2277, 18.534, 19.5445, 21.5212, 21.9593, 22.1282, 22.23, 25.6158, 25.7953, 30.0525, 30.9289, 33.9274, 34.3805, 34.9418, 38.0767, 39.0525, 39.4667, 39.9962],
     Geom.bar)


text:
"""
Function to save the wantedlabels as a file in a certain directory
Arguments:
myscratch::Boostscratch -> Container with the data to be written
filepath::ASCIIString -> Path where the data should be saved
"""
function write_wantedlabels(myscratch::Boostscratch, filepath::ASCIIString)
	h5open(filepath, "w") do file
		write(file, "wantedlabels", myscratch.wantedlabels)
		if(!isempty(myscratch.pooledcovarmat.labels))
			write(file, "usedlabels", myscratch.pooledcovarmat.labels)
		else
			write(file, "usedlabels", "")
		end
	end
end
"""
Function to read the values of the covarmat
Arguments:
filepath::ASCIIString -> Path where the data data is saved
"""
function read_covarmat(filepath::ASCIIString)
	covarmat = Covarmat(Array{Float64,2}(), Array{ASCIIString,1}())
	h5open(filepath, "r") do file
		covarmat.covarmat = read(file, "covarmat")
		covarmat.labels = read(file, "labels")
	end
	return covarmat
end
"""
Function to save all for the boostingprocess-relevant data and the results as a file in a certain directory
Arguments:
myscratch::Boostscratch -> Container with the data to be written
filepath::ASCIIString -> Path where the data should be saved
"""
function write_boostscratch(myscratch::Boostscratch, filepath::ASCIIString)
	h5open(filepath, "w") do file
		write(file, "actualstepno", myscratch.actualstepno)
		write(file, "actualbeta", myscratch.actualbeta)
		write(file, "actualscore", myscratch.actualscore)
		write(file, "pooledunibeta", myscratch.pooledunibeta.unibeta)
		write(file, "pooledunibetalabels", myscratch.pooledunibeta.labels)
		write(file, "pooledunibetan", myscratch.pooledunibeta.n)
		write(file, "pooledcovarmat", myscratch.pooledcovarmat.covarmat)
		write(file, "pooledcovarmatlabels", myscratch.pooledcovarmat.labels)
		write(file, "actualupdates", myscratch.actualupdates)
		write(file, "selections", myscratch.selections)
		write(file, "wantedlabels", myscratch.wantedlabels)
		write(file, "nu", myscratch.nu)
		write(file, "stepno", myscratch.stepno)
		write(file, "actualsel", myscratch.actualsel)
		write(file, "actualnom", myscratch.actualnom)
		write(file, "a", myscratch.a)
		write(file, "x", myscratch.x)
	end
end
"""
Function to set the saved data to continue the boosting (since the process could be interrupted to perform a new data call)
Arguments:
myscratch::Boostscratch -> Container with the data which need to be set to perform the boostingsteps
filepath::ASCIIString -> Path where the data data is saved
"""
function set_boostscratch!(myscratch::Boostscratch, filepath::ASCIIString)
	h5open(filepath, "r") do file
		myscratch.actualstepno = read(file, "actualstepno")
		myscratch.actualbeta  = read(file, "actualbeta")
		myscratch.actualscore  = read(file, "actualscore")
		myscratch.pooledunibeta.unibeta  = read(file, "pooledunibeta")
		myscratch.pooledunibeta.labels  = read(file, "pooledunibetalabels")
		myscratch.pooledunibeta.n  = read(file, "pooledunibetan")
		myscratch.actualupdates  = read(file, "actualupdates")
		myscratch.selections  = read(file, "selections")
		myscratch.wantedlabels  = read(file, "wantedlabels")
		myscratch.nu  = read(file, "nu")
		myscratch.stepno  = read(file, "stepno")
		myscratch.actualsel  = read(file, "actualsel")
		myscratch.actualnom  = read(file, "actualnom")
		myscratch.a  = read(file, "a")
		myscratch.x  = read(file, "x")
		if(myscratch.actualsel != 1)
			myscratch.pooledcovarmat.covarmat  = read(file, "pooledcovarmat")
			myscratch.pooledcovarmat.labels  = read(file, "pooledcovarmatlabels")
		end
	end
end
"""
Function to read the univariate estimators from a certain directory
Arguments:
filepath::ASCIIString -> Path where the data data is saved
"""
function read_unibeta(filepath::ASCIIString)
	unibeta = Unibeta(Array{Float64,1}(), Array{Int,1}(), Array{ASCIIString, 1}())
	h5open(filepath, "r") do file
		unibeta.unibeta = read(file, "unibeta")
		unibeta.n = read(file, "n")
		unibeta.labels = read(file, "labels")
	end
	return unibeta
end
"""
Function to save the univarite estimators as a file in a certain directory
Arguments:
unibeta::Unibeta -> Container with the data to be written
filepath::ASCIIString -> Path where the data should be saved
"""
function write_unibeta(unibeta::Unibeta, filepath::ASCIIString)
	h5open(filepath, "w") do file
		write(file, "unibeta", unibeta.unibeta)
		write(file, "n", unibeta.n)
		write(file, "labels", unibeta.labels)
	end
end
"""
Function to save the covariances matrices as a file in a certain directory
Arguments:
covarmat::Covarmat -> Container with the data to be written
filepath::ASCIIString -> Path where the data should be saved
"""
function write_covarmat(covarmat::Covarmat, filepath::ASCIIString)
	h5open(filepath, "w") do file
		write(file, "covarmat", covarmat.covarmat)
		write(file, "labels", covarmat.labels)
	end
end
"""
Function to read the labels which were already used and the labels which are asked for in a new data call.
Arguments:
filepath::ASCIIString -> Path where the data data is saved
"""
function read_labels(filepath::ASCIIString)
	wantedlabels = Array{ASCIIString,1}()
	usedlabels = Array{ASCIIString, 1}()
	h5open(filepath,"r") do file
		wantedlabels = read(file, "wantedlabels")
		usedlabels = read(file, "usedlabels")
	end
	return wantedlabels, usedlabels
end
"""
Function to read the individual data of a cohort from a certain directory
Arguments:
filepath::ASCIIString -> Path where the data data is saved
"""
function read_data(filepath::ASCIIString)
	mycohort = Discohort(Array{Float64,2}(), Array{Float64,1}(), Array{ASCIIString,1}())
	h5open(filepath, "r") do file
		mycohort.x = read(file, "x")
		mycohort.y = read(file, "y")
		mycohort.labels = read(file, "labels")
		#mycohort.id = read(file, "id")
	end
	return mycohort
end

max_stars_repo_name: danielazoeller/DistributedBoosting

max_issues_repo_name: danielazoeller/DistributedBoosting

max_forks_repo_name: danielazoeller/DistributedBoosting


text:
#####
##### setup and utilities for tests
#####

using Random, StaticArrays, FiniteDifferences, ForwardDiff, DiffResults, LogExpFunctions

Random.seed!(0x3b1ac7d1ef4ad7c3eab4e377ca14b76c) # consistent test runs

"""
Random parameters with `N` sectors. `L` is added to positive parameters, to bound away from
`0`. For unit testing.
"""
function random_parameters(::Val{N}, L = 0.1) where N
    σ = rand() * 2.0 .+ L
    Ω̂s = randn(SVector{N})
    ϵs = rand(SVector{N}) .* 2.0 .+ L
    p̂s = randn(SVector{N})
    Ĉ = rand() * 3.0 + L
    Ê = logsumexp((Ĉ .* ϵs .+ p̂s) .* (1-σ) .+ Ω̂s) / (1-σ)
    (; Ê, σ, Ω̂s, ϵs, p̂s, Ĉ)
end

"""
Relative residual from the Newton solver, should be `< tol`. For correctness checks.
"""
function newton_relative_residual(; Ê, Ω̂s, σ, p̂s, ϵs, Ĉ)
    lhs = logsumexp(@. Ω̂s + (1 - σ) * (p̂s + ϵs * Ĉ))
    res = (lhs - Ê * (1 - σ)) / (1 + abs(Ê))
end

"Derivative of (univariate) `f` at `x` (= 0 by default)."
∂(f; x = 0.0, kwargs...) = central_fdm(5, 1; factor = 1e6, kwargs...)(f, x)

"Add `x` at `i`."
function add_at(v::SVector, i::Int, x)
    m = MVector(v)
    m[i] += x
    SVector(m)
end

"Value and derivative using ForwardDiff."
function fwd_d(f::F, x::Real) where F
    r = ForwardDiff.derivative!(DiffResults.DiffResult(x, x), f, x)
    DiffResults.value(r), DiffResults.derivative(r)
end

"Value and gradient using ForwardDiff."
function fwd_∇(f, x::AbstractVector)
    r = DiffResults.DiffResult(x[1], x)
    r = ForwardDiff.gradient!(r, f, x)
    DiffResults.value(r), DiffResults.gradient(r)
end


text:
module LSDYNA

using Printf

export cross_section_extruder, cfile_cross_section_extruder

function cross_section_extruder(xcoords, ycoords, element_info, L, dL, path, scl_filename,lsdyna_filename)

    cd(path)

    io = open(scl_filename, "w")

    println(io,"define:")
    println(io,"void main(void)")
    println(io,"{")

    num_elem = size(element_info)[1]

    for i=1:num_elem
        @printf(io, "ExecuteCommand(\"line param  0 %f %f 0 %f %f\");\n", xcoords[i],ycoords[i],xcoords[i+1],ycoords[i+1])
    end

    println(io, "ExecuteCommand(\"elgenerate shelltype 6\");")
    println(io, "ExecuteCommand(\"genselect target occobject\");")
    println(io, "ExecuteCommand(\"occfilter add Edge\");")
    println(io, "ExecuteCommand(\"genselect whole\");")

    @printf(io, "ExecuteCommand(\"elgenerate shell curvedrag 1 1 %s %s 0 0 0 1 0 0 0",string(L), string(floor(L/dL)))

    for i=1:num_elem
        @printf(io," %s",string(i))
    end

    @printf(io, "\");\n")

    println(io, "ExecuteCommand(\"elgenerate accept\");")

    println(io, "ExecuteCommand(\"genselect whole\");")
    println(io, "ExecuteCommand(\"dupnode merge 0.001000\");")

    @printf(io, "ExecuteCommand(\"geomag del ")

    for i=1:num_elem
        @printf(io," %s",string(i, "e"))
    end

    @printf(io, "\");\n")

    println(io, string("ExecuteCommand(\"save keyword \\","\"",path,"/",lsdyna_filename,"\\","\"\");"))

    println(io,"}")
    println(io,"main();")

    close(io)

end

function cfile_cross_section_extruder(xcoords, ycoords, element_info, L, dL, cfile_filename,lsdyna_filename)

    # cd(path)

    io = open(cfile_filename, "w")

    num_elem = size(element_info)[1]

    for i=1:num_elem
        @printf(io, "line param  0 %f %f 0 %f %f\n", xcoords[i],ycoords[i],xcoords[i+1],ycoords[i+1])
    end

    println(io, "elgenerate shelltype 6")
    println(io, "genselect target occobject")
    println(io, "occfilter add Edge")
    println(io, "genselect whole")

    @printf(io, "elgenerate shell curvedrag 1 1 %s %s 0 0 0 1 0 0 0",string(L), string(floor(L/dL)))

    for i=1:num_elem
        @printf(io," %s",string(i))

        if i==num_elem
            @printf(io," %s\n",string(i))
        end
    end

    println(io, "elgenerate accept")

    # println(io, "genselect whole")
    # println(io, "dupnode merge 0.4000")

    @printf(io, "geomag del ")

    for i=1:num_elem
        @printf(io," %s",string(i, "e"))

        if i==num_elem
            @printf(io," %s\n",string(i, "e"))
        end
    end

    println(io, "dupnode open 1")
    @printf(io, "dupnode showdup %f\n", dL/100)
    println(io, "dupnode select 1")
    println(io, "genselect target node")
    println(io, "genselect node add part 1/0")
    @printf(io, "dupnode merge %f\n", dL/100)  #this dL/100 is a little tricky, might merge nodes not intended to be merged for a very fine mesh
    println(io, "genselect clear")
    println(io, "dupnode accept")

    # kfile_path="/Volumes/GoogleDrive/Shared drives/RunToSolve/Projects/Geometrica/cladding_study_2020/model_generation/all/"

    # lsdyna_filename=kfile_path*"h_purlin_"*string(i)*"_"*string(i)*".k"

    println(io, string("save keyword \"", lsdyna_filename, "\""))

    close(io)

end

end #module


text:
using BinDeps

@BinDeps.setup

libz = library_dependency("libz", aliases=["libz", "libzlib", "zlib1"])

if !Sys.iswindows()
    provides(AptGet, Dict("zlib" => libz, "zlib1g" => libz))
    provides(Yum, "zlib", [libz])
    provides(Pacman, "zlib", [libz])
else
    using WinRPM
    provides(WinRPM.RPM, "zlib1", [libz])
end

@BinDeps.install Dict([:libz => :libz])


text:
module PolynomialTestSystems

    import MultivariatePolynomials
    const MP = MultivariatePolynomials

    import DynamicPolynomials: Polynomial, PolyVar, @polyvar, differentiate
    using LinearAlgebra

    include("testsystem.jl")
    include("systems.jl")

end


text:
using Documenter, LPLib

makedocs()

max_stars_repo_name: abogushov/lplib

max_issues_repo_name: abogushov/lplib

max_forks_repo_name: abogushov/lplib


text:
function Base.read(io::IO, ::Type{SMPSModel})
    raw = read(io, RawSMPS)
    return SMPSModel(raw)
end
"""
    read(io::IO,
         ::Type{SMPSSampler})

Return an `StochasticModel` from the model definition read from `io` in SMPS format.
"""
function Base.read(io::IO, ::Type{StochasticModel})
    smps::SMPSModel{2} = read(io, SMPSModel)
    return stochastic_model(smps)
end
"""
    read(io::IO,
         ::Type{SMPSSampler})

Return an `SMPSSampler` capable of sampling `SMPSScenario` using the model definition read from `io` in SMPS format.
"""
function Base.read(io::IO, ::Type{SMPSSampler})
    smps::SMPSModel{2} = read(io, SMPSModel)
    return SMPSSampler(smps.raw.sto, smps.stages[2])
end
"""
    read(io::IO,
         ::Type{StochasticProgram};
         num_scenarios::Union{Nothing, Integer} = nothing,
         instantiation::StochasticInstantiation = UnspecifiedInstantiation(),
         optimizer = nothing;
         defer::Bool = false,
         kw...)

Instantiate a two-stage stochastic program using the model definition read from `io` in SMPS format, of size `num_scenarios`. If `num_scenarios = nothing`, instantiate using the full support. Optionally, supply an `optimizer`. If no explicit `instantiation` is provided, the structure is induced by the optimizer. The structure is `Deterministic` by default.
"""
function Base.read(io::IO,
                   ::Type{StochasticProgram};
                   num_scenarios::Union{Nothing, Integer} = nothing,
                   instantiation::StochasticInstantiation = StochasticPrograms.UnspecifiedInstantiation(),
                   optimizer = nothing,
                   defer::Bool = false,
                   direct_model::Bool = false,
                   kw...)
    smps::SMPSModel{2} = read(io, SMPSModel)
    sm = stochastic_model(smps)
    sampler = SMPSSampler(smps.raw.sto, smps.stages[2])
    if num_scenarios != nothing
        return instantiate(sm, sampler, num_scenarios; instantiation = instantiation, optimizer = optimizer, defer = defer, direct_model = direct_model, kw...)
    else
        return instantiate(sm, full_support(sampler); instantiation = instantiation, optimizer = optimizer, defer = defer, direct_model = direct_model, kw...)
    end
end
Base.read(filename::AbstractString, ::Type{StochasticProgram}; kw...) = open(io -> read(io, StochasticProgram; kw...), filename)


text:
"""
  vulnerability(N::UnipartiteNetwork)

This function is mesuring the vulnerability of a given network. The vulnerability of a network represents the ink density (L/S) divided by the number of preys. The number of preys: (prop. of basal + prop. of top) x S
"""

function vulnerability(N::UnipartiteNetwork)
    out_degree = collect(values(degree(N; dims=1)))
    in_degree = collect(values(degree(N; dims=2)))
    degrees = DataFrame(out=out_degree, in=in_degree)
    B = nrow(degrees[degrees.out .== 0, :]) / nrow(degrees)
    T = nrow(degrees[degrees.in .== 0, :]) / nrow(degrees)
    I = nrow(degrees[(degrees.out .!= 0) .& (degrees.in .!= 0), :]) / nrow(degrees)
    L = links(N)
    Vulnerability = links_1 ./ (S.^2 .* (B .+ I))
    return Vulnerability
end


text:
"""
    annealing_step(furnace[; cooldown][, cost][, modifier])

Perform one step of annealing procedure. `cooldown` is a function
defining how fast the furnace will lose temperature. `cost`
determines a function we want to minimize. `modifier` determines which
small modifications are made to the system during the step.

**NB:** The updated state of the annealing process is returned by this
function as a new `Furnace` object. Do not discard it.
"""
function annealing_step(furnace  :: Furnace;
                        cooldown :: Function         = exponential_cooldown(),
                        cost     :: Function         = euclid_mean,
                        modifier :: AbstractModifier = InterfaceFlipper())
    # Some statistics
    rejected = false
    accepted = false

    # Compute the cost function which we are trying to minimize
    c1 = cost(furnace.system, furnace.target)

    # Do a small modification to our system
    state = modify!(furnace.system, modifier)

    # Compute the new value for the cost function
    c2 = cost(furnace.system, furnace.target)

    # if c1 < c2 the swap is accepted
    if c2 > c1
        threshold = exp(-(c2 - c1) / furnace.temperature)
        if rand(Float64) > threshold
            reject!(furnace.system, modifier, state)
            @assert c1 ≈ cost(furnace.system, furnace.target)
            rejected = true
        end

        accepted = !rejected
    end

    # Calculate a new temperature
    newtemperature = rejected ? furnace.temperature : cooldown(furnace.temperature, c2)

    # Construct a new furnace
    return Furnace(furnace.system, furnace.target,
                   newtemperature,
                   furnace.steps + 1,
                   furnace.rejected + rejected,
                   furnace.accepted + accepted)
end


text:
using ROS
ROS.@include std_msgs: Float64MultiArray

pub = nothing

function init()
    ROS.init("julia_pub_sub")
    nh = ROS.NodeHandle()

    global pub = ROS.advertise(nh, "test_pub", ROS.std_msgs_Float64MultiArray, 1)
    sub = ROS.subscribe(nh, "test_sub", 1, ROS.std_msgs_Float64MultiArray, callback)

    while ROS.ok()
        ROS.spinOnce()
    end
end

function callback(t)
    t.data = [1,2,3,4]
    ROS.publish(pub,t)
end

init()


text:
function check_opt(ctx, opt)
    mxstp0 = 500
    mord = [0, 12, 5]
    if ctx.state == 0
        ctx.state = 1
    end
    if ctx.state == 1
        opt.h0 = 0.0
        opt.mxordn = mord[2]
        opt.mxords = mord[3]
    end
    if ctx.neq <= 0
        @warn "ERROR(\"[lsoda] neq = %d is less than 1\\n\",ctx->neq)"
        return 0
    end
    rtol = opt.rtol - 1
    atol = opt.atol - 1
    if ctx.state == 1 || ctx.state == 3
        local i
        @cfor i = 1 i <= ctx.neq i += 1 begin
                rtoli = rtol[i + 1]
                atoli = atol[i + 1]
                if rtoli < 0.0
                    @warn "ERROR(\"[lsoda] rtol = %g is less than 0.\\n\",rtoli)"
                end
                if atoli < 0.0
                    @warn "ERROR(\"[lsoda] atol = %g is less than 0.\\n\",atoli)"
                    return 0
                end
            end
    end
    if opt.itask == 0
        opt.itask = 1
    end
    if opt.itask < 1 || opt.itask > 5
        fprintf(__stderrp, "[lsoda] illegal itask = %d\n", opt.itask)
        return 0
    end
    if opt.ixpr < 0 || opt.ixpr > 1
        fprintf(__stderrp, "[lsoda] ixpr = %d is illegal\n", opt.ixpr)
        return 0
    end
    if opt.mxstep < 0
        fprintf(__stderrp, "[lsoda] mxstep < 0\n")
        return 0
    end
    if opt.mxstep == 0
        opt.mxstep = mxstp0
    end
    if opt.mxhnil < 0
        fprintf(__stderrp, "[lsoda] mxhnil < 0\n")
        return 0
    end
    if ctx.state == 1
        if opt.mxordn < 0
            fprintf(__stderrp, "[lsoda] mxordn = %d is less than 0\n", opt.mxordn)
            return 0
        end
        if opt.mxordn == 0
            opt.mxordn = 100
        end
        opt.mxordn = @warn("min(opt->mxordn,mord[1])")
        if opt.mxords < 0
            fprintf(__stderrp, "[lsoda] mxords = %d is less than 0\n", opt.mxords)
            return 0
        end
        if opt.mxords == 0
            opt.mxords = 100
        end
        opt.mxords = @warn("min(opt->mxords,mord[2])")
    end
    if opt.hmax < 0.0
        fprintf(__stderrp, "[lsoda] hmax < 0.\n")
        return 0
    end
    opt.hmxi = 0.0
    if opt.hmax > 0
        opt.hmxi = 1.0 / opt.hmax
    end
    if opt.hmin < 0.0
        fprintf(__stderrp, "[lsoda] hmin < 0.\n")
        return 0
    end
    return 1
end


text:
# createCR3BPODEProb : Creates a DifferentialEquations.jl ODE Problem for CR3BP state and
# costate differential equations. Uses vector continuous callback for switch detection
# and termination if fuel is depleated for impact with body. A copy of parameters is
# created if copyParams = true, which is required if numerically integrating in parallel
function createCR3BPODEProb(y0::AbstractVector, tspan::Tuple, params::AbstractCR3BPIndirectParams,
    integrationTypeFlag::IntegrationFlag, homotopyFlag::MEMF; copyParams = false,
    termCallbacks = false, inPlace = false, ϵ = -1.0, save_positions = (false, false))

    # Copy parameters if desired
    if copyParams
        ps = deepcopy(params)
    else
        ps = params
    end

    # Set ϵ if desired
    if ϵ >= 0
        ps.ϵ = ϵ
    end

    # Check that utype is set appropriately
    cSc = ps.sp.isp*9.81*ps.crp.TU / (ps.crp.LU*1000.0)
    λv = norm(view(y0,11:13))
    S = computeS(y0, λv, cSc)
    if S > ps.ϵ
        ps.utype = 0
    elseif S < -ps.ϵ
        ps.utype = 2
    else
        ps.utype = 1
    end

    # Instantiate callback
    if termCallbacks
        cb = VectorContinuousCallback(
            cr3bpEomsCondition,
            cr3bpEomsAffect!,
            cr3bpEomsAffect!, 4;
            idxs = nothing,
            rootfind = true,
            interp_points = 20,
            abstol = 1e-14,
            reltol = 0.0,
            save_positions = save_positions)
    else
        cb = ContinuousCallback(
            cr3bpEomsConditionNoTerm,
            cr3bpEomsAffectNoTerm!,
            cr3bpEomsAffectNoTerm!;
            idxs = nothing,
            rootfind = true,
            interp_points = 20,
            abstol = 1e-14,
            reltol = 0.0,
            save_positions = save_positions)
    end

    # ODE Problem
    if inPlace
        ff = ODEFunction{true}((du,u,p,t) -> cr3bpEomIndirect!(du,u,p,t,homotopyFlag))
    else
        ff = ODEFunction{false}((u,p,t) -> cr3bpEomIndirect(u,p,t,homotopyFlag))
    end
    return ODEProblem(ff, y0, tspan, ps; callback=cb)
end

function createCR3BPODEProb(y0::AbstractVector, tspan::Tuple, params::AbstractCR3BPIndirectParams,
    integrationTypeFlag::InitializationWithIntegralCost, homotopyFlag::MEMF; copyParams = false,
    termCallbacks = false, inPlace = false, ϵ = -1.0, save_positions = (false, false))

    if length(y0) != 15
        throw(ArgumentError("When computing the integral cost, y0 must be of length 15."))
    end

    # Copy parameters if desired
    if copyParams
        ps = deepcopy(params)
    else
        ps = params
    end

    # Set ϵ if desired
    if ϵ >= 0
        ps.ϵ = ϵ
    end

    # Check that utype is set appropriately
    cSc = ps.sp.isp*9.81*ps.crp.TU / (ps.crp.LU*1000.0)
    λv = norm(view(y0,11:13))
    S = computeS(y0, λv, cSc)
    if S > ps.ϵ
        ps.utype = 0
    elseif S < -ps.ϵ
        ps.utype = 2
    else
        ps.utype = 1
    end

    # Instantiate callback
    if termCallbacks
        cb = VectorContinuousCallback(
            cr3bpEomsCondition,
            cr3bpEomsAffect!,
            cr3bpEomsAffect!, 4;
            idxs = nothing,
            rootfind = true,
            interp_points = 20,
            abstol = 1e-14,
            reltol = 0.0,
            save_positions = save_positions)
    else
        cb = ContinuousCallback(
            cr3bpEomsConditionNoTerm,
            cr3bpEomsAffectNoTerm!,
            cr3bpEomsAffectNoTerm!;
            idxs = nothing,
            rootfind = true,
            interp_points = 20,
            abstol = 1e-14,
            reltol = 0.0,
            save_positions = save_positions)
    end

    # ODE Problem
    if inPlace
        ff = ODEFunction{true}((du,u,p,t) -> cr3bpEomIndirectIntegralCost!(du,u,p,t,homotopyFlag))
    else
        ff = ODEFunction{false}((u,p,t) -> cr3bpEomIndirectIntegralCost(u,p,t,homotopyFlag))
    end
    return ODEProblem(ff, y0, tspan, ps; callback=cb)
end

function createCR3BPODEProb(y0::AbstractVector, tspan::Tuple, scenario::String,
    integrationTypeFlag::IntegrationFlag, homotopyFlag::HomotopyFlag)
    ps = initCR3BPIndirectParams(scenario)
    return createCR3BPODEProb(y0, tspan, ps, integrationTypeFlag, homotopyFlag)
end

function createCR3BPODEWithSTMProb(z0::AbstractVector, tspan::Tuple, params::CR3BPIndirectWithSTMParams,
    homotopyFlag::MEMF; copyParams = false, ϵ = -1.0)

    # Copy parameters if desired
    if copyParams
        ps = deepcopy(params)
    else
        ps = params
    end

    # Set ϵ if desired
    if ϵ >= 0.0
        ps.ϵ = ϵ
    end

    # Check that utype is set appropriately
    cSc = ps.sp.isp*9.81*ps.crp.TU / (ps.crp.LU*1000.0)
    λv = norm(view(z0,11:13))
    S = computeS(z0, λv, cSc)
    if S > ps.ϵ
        ps.utype = 0
    elseif S < -ps.ϵ
        ps.utype = 2
    else
        ps.utype = 1
    end

    # Instantiate callback
    cb = ContinuousCallback(
        cr3bpEomsConditionNoTerm,
        cr3bpEomsAffectNoTermWithSTM!,
        cr3bpEomsAffectNoTermWithSTM!;
        idxs = nothing,
        rootfind = true,
        interp_points = 10,
        abstol = 1e-14,
        reltol = 0.0,
        save_positions = (false, false))

    # ODE Problem
    ff = ODEFunction{true}((du,u,p,t) -> cr3bpEomIndirectWithSTM!(du,u,p,t,homotopyFlag))
    return ODEProblem(ff, z0, tspan, ps; callback=cb)
end

function createCR3BPODEWithSTMProb(z0::AbstractVector, tspan::Tuple, scenario::String, homotopyFlag::HomotopyFlag)
    ps = initCR3BPIndirectWithSTMParams(scenario)
    return createCR3BPODEWithSTMProb(z0, tspan, ps, homotopyFlag)
end

function createCR3BPODENoControlProb(x0::AbstractVector, tspan::Tuple, params::CR3BPParams;
    copyParams = false, inPlace = false)

    # Copy parameters if desired
    if copyParams
        ps = deepcopy(params)
    else
        ps = params
    end

   # ODE Problem
    if inPlace
        ff = ODEFunction{true}(cr3bpEomNoControl!)
    else
        ff = ODEFunction{false}(cr3bpEomNoControl)
    end

    return ODEProblem(ff, x0, tspan, ps)
end

function createCR3BPODENoControlProb(x0::AbstractVector, tspan::Tuple, scenario::String; inPlace = false)
    psFull = initCR3BPIndirectParams(scenario)
    return createCR3BPODENoControlProb(x0, tspan, psFull.crp; copyParams = false, inPlace = inPlace)
end


text:
function NewEst3Ps(g0,m)
#NewEst3Ps    Calculates covariance matrix of sqrt(T)*sample average.
#
#
#  Usage:     Shat = NewEst3Ps(hhat,m)
#
#  Input:     g0          TxK matrix of moment functions
#             m           order of autoregression
#
#  Output:    Shat        covariance matrix of sum( hhat/sqrt(T) ),
#                         that is the covariance matrix of sqrt(T)*mean(hhat)
#
#
#  Note: The CLT typically say that sqrt(T)*sample average ->d N(mu,Shat),
#        where Shat is what this code estimates. Clearly, Var(sample average) =
#        Shat/T.
#
#
#  Reference: Newey and West, 1987, Econometrica
#             Newey, 1985, Journal of Econometrics
#
#
#  Paul Soderlind (Paul.Soderlind@unisg.ch), to Julia Oct 2015
#-----------------------------------------------------------------------

  T = size(g0,1)                     #g is Txq
  m = min(m,T-1)                     #number of lags

  g = g0 .- mean(g0,1)               #Normalizing to Eg=0

  S = g'g/T                          #(qxT)*(Txq)
  for s = 1:m
    Omega_s = g[s+1:T,:]'g[1:T-s,:]/T   #same as Sum[g(t)*g(t-s)',t=s+1,T]
    S       = S + (1 - s/(m+1))*(Omega_s + Omega_s')
  end

  return S

end
#-----------------------------------------------------------------------

max_stars_repo_name: PaulSoderlind/PaulSoderlindCode

max_issues_repo_name: PaulSoderlind/PaulSoderlindCode

max_forks_repo_name: PaulSoderlind/PaulSoderlindCode


text:
# MOI backend AbstractModel
MOIU.@model(ParametronMOIModel, # modelname
    (MOI.ZeroOne, MOI.Integer), # scalarsets
    (MOI.LessThan, MOI.GreaterThan, MOI.EqualTo), # typedscalarsets
    (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives), # vectorsets
    (), # typedvectorsets
    (MOI.SingleVariable,), # scalarfunctions
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction), # typedscalarfunctions
    (), # vectorfunctions
    (MOI.VectorAffineFunction,) # typedvectorfunctions
)

# Variable conversion
Functions.Variable(v::MOI.VariableIndex) = Variable(v.value)
MOI.VariableIndex(v::Variable) = MOI.VariableIndex(v.index)

# Sense conversion
function MOI.OptimizationSense(sense::Sense)
    if sense == Minimize
        MOI.MIN_SENSE
    elseif sense == Maximize
        MOI.MAX_SENSE
    else
        error()
    end
end

# update!
struct IdentityVarMap end
Base.getindex(::IdentityVarMap, i) = MOI.VariableIndex(i)

function update!(moi_f::MOI.ScalarAffineFunction, f::AffineFunction, varmap = IdentityVarMap())
    moi_f.constant = f.constant[]
    resize!(moi_f.terms, length(f.linear))
    @inbounds for i in eachindex(f.linear)
        term = f.linear[i]
        moi_f.terms[i] = MOI.ScalarAffineTerm(term.coeff, varmap[term.var.index])
    end
    moi_f
end

function update!(moi_f::MOI.ScalarQuadraticFunction, f::QuadraticFunction, varmap = IdentityVarMap())
    affine = f.affine
    moi_f.constant = affine.constant[]
    resize!(moi_f.affine_terms, length(affine.linear))
    @inbounds for i in eachindex(affine.linear)
        term = affine.linear[i]
        moi_f.affine_terms[i] = MOI.ScalarAffineTerm(term.coeff, varmap[term.var.index])
    end
    resize!(moi_f.quadratic_terms, length(f.quadratic))
    @inbounds for i in eachindex(f.quadratic)
        term = f.quadratic[i]
        var_index_1 = varmap[term.rowvar.index]
        var_index_2 = varmap[term.colvar.index]
        coeff = ifelse(term.rowvar == term.colvar, 2 * term.coeff, term.coeff)
        moi_f.quadratic_terms[i] = MOI.ScalarQuadraticTerm(coeff, var_index_1, var_index_2)
    end
    moi_f
end

function update!(moi_f::MOI.VectorAffineFunction, fs::Vector{<:AffineFunction}, varmap = IdentityVarMap())
    resize!(moi_f.constants, length(fs))
    num_linear_terms = 0
    for f in fs
        num_linear_terms += length(f.linear)
    end
    resize!(moi_f.terms, num_linear_terms)
    i = 1
    @inbounds for row in eachindex(fs)
        f = fs[row]
        for term in f.linear
            moi_f.terms[i] = MOI.VectorAffineTerm(Int64(row), MOI.ScalarAffineTerm(term.coeff, varmap[term.var.index]))
            i += 1
        end
        moi_f.constants[row] = f.constant[]
    end
    moi_f
end

update!(moi_f::MOI.AbstractFunction, f::Nothing, varmap = IdentityVarMap()) = moi_f

# make_moi_equivalent
make_moi_equivalent(::Type{AffineFunction{T}}) where {T} =
    MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{T}[], zero(T))
make_moi_equivalent(::Type{QuadraticFunction{T}}) where {T} =
    MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{T}[], MOI.ScalarQuadraticTerm{T}[], zero(T))
make_moi_equivalent(::Type{<:AbstractVector{AffineFunction{T}}}) where {T} =
    MOI.VectorAffineFunction(MOI.VectorAffineTerm{T}[], T[])

# canonical_function_type
canonical_function_type(::Type{Variable}, ::Type{T}) where {T} = AffineFunction{T}
canonical_function_type(::Type{<:LinearTerm}, ::Type{T}) where {T} = AffineFunction{T}
canonical_function_type(::Type{<:AffineFunction}, ::Type{T}) where {T} = AffineFunction{T}
canonical_function_type(::Type{<:AbstractVector{<:AffineFunction}}, ::Type{T}) where {T} = Vector{AffineFunction{T}}
canonical_function_type(::Type{<:QuadraticTerm}, ::Type{T}) where {T} = QuadraticFunction{T}
canonical_function_type(::Type{<:QuadraticFunction}, ::Type{T}) where {T} = QuadraticFunction{T}

# moi_to_native_type
moi_to_native_type(::Type{MOI.ScalarAffineFunction}) = AffineFunction{T} where T
moi_to_native_type(::Type{MOI.ScalarQuadraticFunction}) = QuadraticFunction{T} where T
moi_to_native_type(::Type{MOI.VectorAffineFunction}) = Vector{AffineFunction{T}} where T
moi_to_native_type(::Type{MOI.VectorQuadraticFunction}) = Vector{QuadraticFunction{T}} where T
moi_to_native_type(::Type{MOI.SingleVariable}) = Nothing

# Objective
struct Objective{E, F}
    expr::WrappedExpression{E}
    f::F
    isconstant::Bool
end

function Objective(::Type{T}, expr) where T
    val = evalarg(expr)
    E = canonical_function_type(typeof(val), T)
    converted = @expression convert(E, expr)
    isconstant = converted isa E # i.e., it's just a value; not a LazyExpression or a Parameter
    wrapped = convert(WrappedExpression{E}, converted)
    f = make_moi_equivalent(E)
    F = typeof(f)
    update!(f, wrapped())
    Objective{E, F}(wrapped, f, isconstant)
end

function update!(objective::Objective{E, F}, optimizer::MOI.AbstractOptimizer, varmap) where {E, F}
    if !objective.isconstant
        update!(objective.f, objective.expr(), varmap)
        MOI.set(optimizer, MOI.ObjectiveFunction{F}(), objective.f)
    end
    nothing
end

# Constraint
mutable struct Constraint{E, F<:MOI.AbstractFunction, S<:MOI.AbstractSet}
    expr::WrappedExpression{E}
    f::F
    set::S
    isconstant::Bool
    modelindex::MOI.ConstraintIndex{F, S}
    optimizerindex::MOI.ConstraintIndex{F, S}

    function Constraint(::Type{T}, expr, set::S) where {T, S<:MOI.AbstractSet}
        val = evalarg(expr)
        E = canonical_function_type(typeof(val), T)
        converted = @expression convert(E, expr)
        isconstant = converted isa E # i.e., it's just a value; not a LazyExpression or a Parameter
        wrapped = convert(WrappedExpression{E}, converted)
        f = make_moi_equivalent(E)
        F = typeof(f)
        update!(f, wrapped())
        new{E, F, S}(wrapped, f, set, isconstant)
    end

    function Constraint(f::F, set::S) where {F<:MOI.AbstractFunction, S<:MOI.AbstractSet}
        E = Nothing
        wrapped = convert(WrappedExpression{E}, @expression nothing)
        new{E, F, S}(wrapped, f, set, true)
    end
end

function update!(constraint::Constraint, optimizer::MOI.AbstractOptimizer, varmap)
    if !constraint.isconstant
        update!(constraint.f, constraint.expr(), varmap)
        MOI.set(optimizer, MOI.ConstraintFunction(), constraint.optimizerindex, constraint.f)
    end
    nothing
end
update!(constraint::Constraint{Nothing}, optimizer::MOI.AbstractOptimizer, varmap) = nothing

# Constraints
let
    constraint_specs = [(MOI.ScalarAffineFunction{T} where T, MOI.GreaterThan{T} where T),
                        (MOI.ScalarAffineFunction{T} where T, MOI.LessThan{T} where T),
                        (MOI.ScalarAffineFunction{T} where T, MOI.EqualTo{T} where T),
                        (MOI.VectorAffineFunction{T} where T, MOI.Nonnegatives),
                        (MOI.VectorAffineFunction{T} where T, MOI.Nonpositives),
                        (MOI.VectorAffineFunction{T} where T, MOI.Zeros),
                        (MOI.ScalarQuadraticFunction{T} where T, MOI.GreaterThan{T} where T),
                        (MOI.ScalarQuadraticFunction{T} where T, MOI.LessThan{T} where T),
                        (MOI.ScalarQuadraticFunction{T} where T, MOI.EqualTo{T} where T),
                        (MOI.VectorQuadraticFunction{T} where T, MOI.Nonnegatives),
                        (MOI.VectorQuadraticFunction{T} where T, MOI.Nonpositives),
                        (MOI.VectorQuadraticFunction{T} where T, MOI.Zeros),
                        (MOI.SingleVariable, MOI.Integer),
                        (MOI.SingleVariable, MOI.ZeroOne)]
    fieldnames = Symbol[]
    constrainttypes = []
    for (F, S) in constraint_specs
        E = moi_to_native_type(F)
        funcsym = F isa UnionAll ? F.body.name.name : F.name.name
        setsym = S isa UnionAll ? S.body.name.name : S.name.name
        ctype = if F isa UnionAll && S isa UnionAll
            Constraint{E{T}, F{T}, S{T}} where T
        elseif F isa UnionAll
            Constraint{E{T}, F{T}, S} where T
        elseif S isa UnionAll
            Constraint{E, F, S{T}} where T
        else
            Constraint{E, F, S}
        end
        push!(fieldnames, Symbol(lowercase(String(funcsym)) * "_in_" * lowercase(String(setsym))))
        push!(constrainttypes, ctype)
    end

    # Constraints struct
    @eval begin
        struct Constraints{T}
            $([:($name::Vector{$(ctype isa UnionAll ? :($ctype{T}) : ctype)}) for (name, ctype) in zip(fieldnames, constrainttypes)]...)
        end

        function Constraints{T}() where {T}
            Constraints{T}($([:(Vector{$(ctype isa UnionAll ? :($ctype{T}) : ctype)}()) for ctype in constrainttypes]...))
        end
    end

    # Base.push! methods
    for ((F, S), name, ctype) in zip(constraint_specs, fieldnames, constrainttypes)
        if ctype isa UnionAll
            @eval function Base.push!(constraints::Constraints{T}, constraint::$(:($ctype{T}))) where T
                push!(constraints.$name, constraint)
            end
        else
            @eval Base.push!(constraints::Constraints, constraint::$(ctype)) = push!(constraints.$name, constraint)
        end
    end

    # update!
    @eval begin
        function update!(constraints::Constraints, optimizer::MOI.AbstractOptimizer, varmap)
            $(map(fieldnames) do fieldname
                quote
                    for constraint in constraints.$fieldname
                        update!(constraint, optimizer, varmap)
                    end
                end
            end...)
            nothing
        end
    end

    # mapindices!
    @eval begin
        function mapindices!(constraints::Constraints, idxmap)
            $(map(fieldnames) do fieldname
                quote
                    for constraint in constraints.$fieldname
                        constraint.optimizerindex = idxmap[constraint.modelindex]
                    end
                end
            end...)
            nothing
        end
    end
end


text:
"""
    log_value(logger, name, value, step)

Logs a Floating-point variable with name `name` at step `step`
"""
function log_value(logger::Logger, name::String, value::Real; step=nothing)
    summ    = SummaryCollection()
    push!(summ.value, scalar_summary(name, value))
    write_event(logger.file, make_event(logger, summ, step=step))
end

function log_value(logger::Logger, name::String, value::Complex; step=nothing)
    log_value(logger, name*"/re", real(value), step)
    log_value(logger, name*"/im", imag(value), step)
end

function scalar_summary(name::String, value::Real)
    Summary_Value(tag=name, simple_value=value)
end

## Backward compatibility
log_value(logger, name, value, step) =
    log_value(logger, name, value, step=step)


text:
include("../src/Hash.jl")
include("testing.jl")

fname = tempname()
test_strings = ["foo", "foobar", "baz", "this is a test", ""]
function check_hash(algorithm, fn)
    print("checking " * algorithm * " sums ")
    for value in test_strings
        f = open(fname, "w")
        write(f, value)
        close(f)

        command_line = strip(split(readall(`openssl dgst -$algorithm $fname`), " ")[2])
        library_with_passed_vaue = Hash.hex(fn(value))

        @assert_equal(command_line, library_with_passed_vaue)

        f = open(fname, "r")
        library_from_file = Hash.hex(fn(f))
        close(f)

        @assert_equal(command_line, library_from_file)
    end
    println()
end

check_hash("md4", Hash.md4)
check_hash("md5", Hash.md5)
check_hash("mdc2", Hash.mdc2)
check_hash("sha1", Hash.sha1)
check_hash("sha224", Hash.sha224)
check_hash("sha256", Hash.sha256)
check_hash("sha384", Hash.sha384)
check_hash("sha512", Hash.sha512)
println()