text:
module eqws_qa

    using FFTW
    using Plots
    using Measures
    theme(:dark)

    export plot_filter_bank_QA
    export plot_filter_bank_QAxy
    export plot_filter_bank_QAxz

    function plot_filter_bank_QA(filt, info; fname="filter_bank_QA.png",p=2)

        # -------- define plot1 to append plot to a list
        function plot1(ps, image; clim=nothing, bin=1.0, fsz=12, label=nothing)
            # ps    - list of plots to append to
            # image - image to heatmap
            # clim  - tuple of color limits (min, max)
            # bin   - bin down by this factor, centered on nx/2+1
            # fsz   - font size
            marg   = -1mm
            nx, ny = size(image)
            nxb    = nx/round(Integer, 2*bin)

            # -------- center on nx/2+1
            i0 = max(1,round(Integer, (nx/2+2)-nxb-1))
            i1 = min(nx,round(Integer, (nx/2)+nxb+1))
            lims = [i0,i1]
            subim = image[i0:i1,i0:i1]
            push!(ps, heatmap(image,aspect_ratio=:equal,clim=clim,
                xlims=lims, ylims=lims, size=(400,400),
                legend=false,xtickfontsize=fsz,ytickfontsize=fsz,#tick_direction=:out,
                rightmargin=marg,leftmargin=marg,topmargin=marg,bottommargin=marg))
            if label != nothing
                annotate!(lims'*[.96,.04],lims'*[.09,.91],text(label,:left,:white,16))
            end
            return
        end

        # -------- initialize array of plots
        ps   = []
        ind  = info["psi_index"]
        jval = info["j_value"]
        J    = length(jval)
        L    = length(info["theta_value"])
        # plot the first 3 and last 3 j values
        jplt = [1:3;J-2:J]

        # -------- loop over j
        for j_ind in jplt
            # first 3 filters in Fourier space
            j = jval[j_ind]
            bfac = max(1,2.0^(j-1))
            for l = 1:3
                label=string("j=",j_ind," ℓ=",l-1)
                #label=latexstring("\\tilde{\\psi}_{$j_ind,$(l-1)}")
                plot1(ps, fftshift(filt[:,:,ind[j_ind,l]]), clim=(0,1),
                      bin=bfac, label=label)
            end
            fac = max(1,2.0^(5-j))
            # real part of config space of 3rd filter
            plot1(ps, fftshift(real(ifft(filt[:,:,ind[j_ind,3]]))),
                  bin=fac, label=string("j=",jval[j_ind]," ℓ=",3-1))

            # sum all filters to form (partial) ring
            ring = dropdims(sum(filt[:,:,ind[j_ind,:]].^2, dims=3), dims=3)
            plot1(ps, fftshift(ring[:,:,1]), clim=(0,1),
                  bin=bfac, label=string("j=",jval[j_ind]," ℓ=",0,":",L-1))
        end

        wavepow = fftshift(dropdims(sum(filt[:,:,ind].^p, dims=(3,4)), dims=(3,4)))
        plot1(ps, wavepow, label=string("j=1:",jval[end]," ℓ=0:",L-1)) #clim=(-0.1,1)

        if info["t"]==1
            wavepow ./= 2
            wavepow += circshift(wavepow[end:-1:1,end:-1:1],(1,1))
        end
        plot1(ps, wavepow, bin=32, label="missing power") #clim=(-0.1,1)
        phi = filt[:,:,info["phi_index"]]
        phi_shift = fftshift(phi)
        plot1(ps, phi_shift.^p, bin=32, label="ϕ") #clim=(-0.1,1)
        disc = wavepow+(phi_shift.^p)
        plot1(ps, fftshift(real(ifft(phi))), label="ϕ")
        plot1(ps, disc, bin=32, label="all") #clim=(-0.1,1)

        myplot = plot(ps..., layout=(7,5), size=(1400,2000))
        savefig(myplot, fname)
        return disc
    end

    #AKS added 3d plotting functions 2021_02_22; They are not fully complete
    #but do give an ok picture of what is going on. Enough to proceed with some applications
    function plot_filter_bank_QAxy(hash; fname="filter_bank_QAxy.png")

        # -------- define plot1 to append plot to a list
        function plot1(ps, image; clim=nothing, bin=1.0, fsz=12, label=nothing,color=:white)
            # ps    - list of plots to append to
            # image - image to heatmap
            # clim  - tuple of color limits (min, max)
            # bin   - bin down by this factor, centered on nx/2+1
            # fsz   - font size
            marg   = -1mm
            nx, ny = size(image)
            nxb    = nx/round(Integer, 2*bin)

            # -------- center on nx/2+1
            i0 = max(1,round(Integer, (nx/2+2)-nxb-1))
            i1 = min(nx,round(Integer, (nx/2)+nxb+1))
            lims = [i0,i1]
            subim = image[i0:i1,i0:i1]
            push!(ps, heatmap(image,aspect_ratio=:equal,clim=clim,
                xlims=lims, ylims=lims, size=(400,400),
                legend=false,xtickfontsize=fsz,ytickfontsize=fsz,#tick_direction=:out,
                rightmargin=marg,leftmargin=marg,topmargin=marg,bottommargin=marg))
            if label != nothing
                annotate!(lims'*[.96,.04],lims'*[.09,.91],text(label,:left,color,16))
            end
            return
        end

        # -------- initialize array of plots
        ps   = []
        jval = hash["j_value"]
        kval = hash["k_value"]
        J    = length(hash["j_value"])
        L    = length(hash["theta_value"])
        K    = length(hash["k_value"])
        nx   = hash["npix"]
        nz   = hash["nz"]
        lup  = hash["psi_index"]
        index= hash["filt_index"]
        value= hash["filt_value"]
        # plot the first 3 and last 3 j values
        jplt = [1:3;J-2:J]
        lplt = 1:3
        kplt = [1,2,K]

        # -------- loop over k
        for k_ind in kplt
            # -------- loop over j
            k = kval[k_ind]
            krs = convert(Int64,(max(1,2.0^(7-k))+nz/2))
            kfs = convert(Int64,(max(1,2.0^(7-k))+1))
            for j_ind in jplt
                # first 3 filters in Fourier space
                j = jval[j_ind]
                bfac = max(1,2.0^(j-1))
                # -------- loop over l
                for l_ind in lplt
                    label=string("j=",j_ind," ℓ=",l_ind-1," k=",k_ind)
                    filt = zeros(nx,nx,nz)
                    ind = lup[j_ind,l_ind,k_ind]
                    filt[index[ind]] = value[ind]
                    #temp = dropdims(sum(filt,dims=3),dims=3)
                    temp = filt[:,:,kfs]
                    plot1(ps, fftshift(temp[:,:]),
                          bin=bfac, label=label)
                end
                fac = max(1,2.0^(5-j))
                # real part of config space of 3rd filter
                l_ind = 3
                filt = zeros(nx,nx,nz)
                ind = lup[j_ind,l_ind,k_ind]
                filt[index[ind]] = value[ind]
                rs = fftshift(real(ifft(filt)))
                rs_slice = rs[:,:,krs]
                #rs_sum = dropdims(sum(rs,dims=3),dims=3)
                plot1(ps, rs_slice,
                      bin=fac, label=string("j=",j_ind," ℓ=",l_ind-1," k=",k_ind))

                # sum all filters to form (partial) ring
                ring = zeros(nx,nx)
                filt = zeros(nx,nx,nz)
                for indx in lup[j_ind,:,k_ind]
                    filt[index[indx]] .+= value[indx].^2
                end
                ring = filt[:,:,kfs]
                plot1(ps, fftshift(ring[:,:]),
                      bin=bfac, label=string("j=",j_ind," ℓ=",0,":",L-1," k=",k_ind))
            end
        end

        filt = zeros(nx,nx,nz)
        for indx in lup[1:J,:,2]
            filt[index[indx]] .+= value[indx].^2
        end

        wavepow = fftshift(dropdims(sum(filt, dims=3), dims=3))
        plot1(ps, wavepow, label=string("j=1:",J," ℓ=0:",L-1," K=1:",K),color=:white)

        if hash["2d_t"]==1
            wavepow ./= 2
            wavepow += circshift(wavepow[end:-1:1,end:-1:1],(1,1))
        end
        plot1(ps, wavepow, #clim=(0,K),
            bin=32, label="missing power",color=:black)

        phi_z = zeros(nx,nx,nz)
        for indx in lup[J+1,1,2]
            phi_z[index[indx]] .+= value[indx].^2
        end

        phi_rs = zeros(nx,nx,nz)
        for indx in lup[J+1,1,2]
            phi_rs[index[indx]] .+= value[indx]
        end
        phi_rst = fftshift(real(ifft(phi_rs)))[:,:,krs]

        k = kval[K]
        krs = convert(Int64,(max(1,2.0^(7-k))+nz/2))
        kfs = convert(Int64,(max(1,2.0^(7-k))+1))

        phi_shift = fftshift(dropdims(sum(phi_z, dims=3), dims=3))
        plot1(ps, phi_shift, #clim=(0,K),
            bin=32, label="ϕ")
        disc = wavepow.+phi_shift
        plot1(ps, phi_rst, label=string("ϕ k=",convert(Int64,k)))
        plot1(ps, disc, #clim=(0,K),
            bin=32, label="all", color=:black)

        myplot = plot(ps..., layout=(19,5), size=(3000,8000))
        savefig(myplot, fname)
        return phi_z
    end

    function plot_filter_bank_QAxz(hash; fname="filter_bank_QAxz.png")

        # -------- define plot1 to append plot to a list
        function plot1(ps, image; clim=nothing, binx=1.0, biny=1.0, fsz=12, label=nothing,color=:white)
            # ps    - list of plots to append to
            # image - image to heatmap
            # clim  - tuple of color limits (min, max)
            # bin   - bin down by this factor, centered on nx/2+1
            # fsz   - font size
            marg   = -1mm
            nx, ny = size(image)
            nxb    = nx/round(Integer, 2*binx)
            nyb    = ny/round(Integer, 2*biny)

            # -------- center on nx/2+1
            xi0 = max(1,round(Integer, (nx/2+2)-nxb-1))
            xi1 = min(nx,round(Integer, (nx/2)+nxb+1))
            xlims = [xi0,xi1]

            yi0 = max(1,round(Integer, (ny/2+2)-nyb-1))
            yi1 = min(ny,round(Integer, (ny/2)+nyb+1))
            ylims = [yi0,yi1]

            subim = image[xi0:xi1,yi0:yi1]
            push!(ps, heatmap(image,aspect_ratio=biny/binx,clim=clim,
                xlims=xlims, ylims=ylims, size=(400,400),
                legend=false,xtickfontsize=fsz,ytickfontsize=fsz,#tick_direction=:out,
                rightmargin=marg,leftmargin=marg,topmargin=marg,bottommargin=marg))
            if label != nothing
                annotate!(xlims'*[.96,.04],ylims'*[.09,.91],text(label,:left,color,16))
            end
            return
        end

        # -------- initialize array of plots
        ps   = []
        jval = hash["j_value"]
        kval = hash["k_value"]
        J    = length(hash["j_value"])
        L    = length(hash["theta_value"])
        K    = length(hash["k_value"])
        nx   = hash["npix"]
        nz   = hash["nz"]
        lup  = hash["psi_index"]
        index= hash["filt_index"]
        value= hash["filt_value"]
        # plot the first 3 and last 3 j values
        jplt = [1:3;J-2:J]
        lplt = [1,2,L]
        kplt = [1,2,K]

        # -------- loop over k
        for k_ind in kplt
            # -------- loop over j
            k = kval[k_ind]
            krs = convert(Int64,(max(1,2.0^(7-k))+nz/2))
            kfs = convert(Int64,(max(1,2.0^(7-k))+1))
            kfac = max(1,2.0^(k-1))
            for j_ind in jplt
                # first 3 filters in Fourier space
                j = jval[j_ind]
                jrs = convert(Int64,(max(1,2.0^(7-j))+nx/2))
                jfs = convert(Int64,(max(1,2.0^(7-j))+1))

                bfac = max(1,2.0^(j-1))
                # -------- loop over l
                for l_ind in lplt
                    label=string("j=",j_ind," ℓ=",l_ind-1," k=",k_ind)
                    filt = zeros(nx,nx,nz)
                    ind = lup[j_ind,l_ind,k_ind]
                    filt[index[ind]] = value[ind]
                    #temp = dropdims(sum(filt,dims=2),dims=2)
                    #temp = filt[:,nx÷2,:]
                    plot1(ps, fftshift(filt)[nx÷2,:,:], #clim=(0,1),
                          binx=kfac, biny=bfac, label=label)
                end
                j_fac = max(1,2.0^(5-j))
                k_fac = max(1,2.0^(5-k))
                # real part of config space of 1st filter
                l_ind = 1

                filt = zeros(nx,nx,nz)
                ind = lup[j_ind,l_ind,k_ind]
                filt[index[ind]] = value[ind]
                rs = fftshift(real(ifft(filt)))
                rs_slice = rs[nx÷2,:,:]
                plot1(ps, rs_slice,
                      binx=k_fac, biny=j_fac, label=string("j=",j_ind," ℓ=",l_ind-1," k=",k_ind))

                filt = zeros(nx,nx,nz)
                ind = lup[j_ind,l_ind,k_ind]
                filt[index[ind]] = value[ind]
                rs = fftshift(real(ifft(filt)))
                rs_slice = rs[:,nx÷2,:]
                plot1(ps, rs_slice,
                    binx=k_fac, biny=j_fac, label=string("j=",j_ind," ℓ=",l_ind-1," k=",k_ind))

                # sum all filters to form (partial) ring
                ring = zeros(nx,nx)
                filt = zeros(nx,nx,nz)
                for indx in lup[j_ind,:,k_ind]
                    filt[index[indx]] .+= value[indx].^2
                end
                ring = filt
                plot1(ps, fftshift(ring)[nx÷2,:,:], #clim=(0,1),
                      binx=kfac, biny=bfac, label=string("j=",j_ind," ℓ=",0,":",L-1," k=",k_ind))
            end
        end

        filt = zeros(nx,nx,nz)
        for indx in lup[1:J,:,:]
            filt[index[indx]] .+= value[indx].^2
        end

        wavepow = fftshift(filt)
        plot1(ps, wavepow[:,nx÷2,:], #clim=(-0.1,1),
        label=string("j=1:",J," ℓ=0:",L-1," K=1:",K))

        wavepow2 = wavepow[:,nx÷2,:]
        if hash["2d_t"]==1
            wavepow2 += reverse(wavepow[:,nx÷2,:], dims = 2)
            wavepow2 += circshift(wavepow2[end:-1:1,end:-1:1],(1,1))
        end
        plot1(ps, wavepow2, #clim=(0,K),
            binx=16, label="missing power",color=:black)

        phi_z = zeros(nx,nx,nz)
        for indx in lup[J+1,1,:]
            phi_z[index[indx]] .+= value[indx].^2
        end

        phi_rs = zeros(nx,nx,nz)
        for indx in lup[J+1,1,:]
            phi_rs[index[indx]] .+= value[indx]
        end
        phi_rst = fftshift(real(ifft(phi_rs)))[:,:,krs]

        k = kval[K]
        krs = convert(Int64,(max(1,2.0^(7-k))+nz/2))
        kfs = convert(Int64,(max(1,2.0^(7-k))+1))

        phi_shift = fftshift(dropdims(sum(phi_z, dims=1), dims=1))
        plot1(ps, phi_shift, #clim=(0,K),
            binx=32, label="ϕ")
        disc = wavepow[:,nx÷2,:].+phi_shift
        plot1(ps, phi_rst, label=string("ϕ k=",convert(Int64,k)))
        plot1(ps, disc, #clim=(0,K),
            binx=32, label="all", color=:black)

        plot1(ps, disc, #clim=(0,K),
            binx=32, label="all", color=:black)

        myplot = plot(ps..., layout=(19,6), size=(3000,8000))
        savefig(myplot, fname)

        filt = zeros(nx,nx,nz)
        ind = lup[1,1,1]
        filt[index[ind]] = value[ind]

        return wavepow
    end
end # of module

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
mutable struct BoltzmannModel{T <: SSDFloat} <: AbstractTemperatureModel{T}

    #temperatures needed for the drift velocity rescaling
    temperature::T      # actual temperature of the crystal
    reftemperature::T   # temperature at which the reference values were measured

    #fit parameters for electrons (e) and holes (h) for <100> and <111> axes
    p1e100::T
    p2e100::T
    p3e100::T

    p1e111::T
    p2e111::T
    p3e111::T

    p1h100::T
    p2h100::T
    p3h100::T

    p1h111::T
    p2h111::T
    p3h111::T
end

function BoltzmannModel(config_file::Dict; T::Type{<:AbstractFloat} = Float32)::BoltzmannModel
    return BoltzmannModel{T}(config_file::Dict)
end

function BoltzmannModel{T}(config_file::Dict; temperature::Union{Missing, T} = missing)::BoltzmannModel where T <: SSDFloat
    config_file["temperature_dependence"]["model"] != "Boltzmann" ? error() : nothing
    ismissing(temperature) ? temperature = config_file["temperature_dependence"]["temperature"] : nothing
    m = BoltzmannModel{T}(
        temperature, #temperature = config_file["temperature_dependence"]["temperature"]
        config_file["drift"]["velocity"]["temperature"],

        config_file["temperature_dependence"]["parameters"]["e100"]["p1"],
        config_file["temperature_dependence"]["parameters"]["e100"]["p2"],
        config_file["temperature_dependence"]["parameters"]["e100"]["p3"],

        config_file["temperature_dependence"]["parameters"]["e111"]["p1"],
        config_file["temperature_dependence"]["parameters"]["e111"]["p2"],
        config_file["temperature_dependence"]["parameters"]["e111"]["p3"],

        config_file["temperature_dependence"]["parameters"]["h100"]["p1"],
        config_file["temperature_dependence"]["parameters"]["h100"]["p2"],
        config_file["temperature_dependence"]["parameters"]["h100"]["p3"],

        config_file["temperature_dependence"]["parameters"]["h111"]["p1"],
        config_file["temperature_dependence"]["parameters"]["h111"]["p2"],
        config_file["temperature_dependence"]["parameters"]["h111"]["p3"]
    )
end

function scale_to_given_temperature(m::BoltzmannModel{T})::NTuple{4,T} where T <: SSDFloat

    function f(p1::T, p2::T, p3::T)::T where T <: SSDFloat
        return (p1+ p2*exp(-p3/m.reftemperature))/(p1 + p2*exp(-p3/m.temperature))
    end

    scale_e100 = f(m.p1e100, m.p2e100, m.p3e100)
    scale_e111 = f(m.p1e111, m.p2e111, m.p3e111)
    scale_h100 = f(m.p1h100, m.p2h100, m.p3h100)
    scale_h111 = f(m.p1h111, m.p2h111, m.p3h111)

    return scale_e100, scale_e111, scale_h100, scale_h111
end

ext: jl


text:
# Copyright (c) 2021, Eric Sabo
# All rights reserved.
#
# This source code is licensed under the BSD-style license found in the
# LICENSE file in the root directory of this source tree.

function ord(n::Integer, q::Integer)
    if q <= 0 || n <= 0
        error("q and n both need to be positive. Passed: q = $q, n = $n")
    end

    # finite stop instead of while
    for i = 1:200
        if mod(BigInt(q)^i, n) == 1
            return i
        end
    end

    error("Unable to compute order($q, $n).")
end

function cyclotomiccoset(x::Integer, q::Integer, n::Integer, verbose::Bool=false)
    # keep small order assumption for now
    temp = [mod(x, n)]
    for i = 0:(n - 1)
        y = mod(temp[end] * q, n)
        if y ∉ temp
            append!(temp, y)
        else
            break
        end
    end
    sort!(temp)
    len = length(temp)

    if verbose
        print("C_$x = {")
        for (i, y) in enumerate(temp)
            if i != len
                print("$y, ")
            else
                println("$y}")
            end
        end
    end

    return temp
end

function allcyclotomiccosets(q::Integer, n::Integer, verbose::Bool=false)
    if n % q == 0
        error("Cyclotomic coset requires gcd(n, q) = 1")
    end

    arr = [[0]]
    for x in 1:(n - 1)
        found = false
        for a in arr
            if x ∈ a
                found = true
                break
            end
        end

        if !found
            Cx = cyclotomiccoset(x, q, n, false)
            push!(arr, Cx)
        end
    end
    #
    # sort!(arr, by=x->x[1])

    if verbose
        for Cx in arr
            len = length(Cx)
            print("C_$(Cx[1]) = {")
            for (i, y) in enumerate(Cx)
                if i != len
                    print("$y, ")
                else
                    println("$y}")
                end
            end
        end
    end

    if sort!(vcat(arr...)) != [i for i in 0:(n - 1)]
        error("Missed some")
    end

    return arr
end

function complementqcosets(q::Integer, n::Integer, qcosets::Vector{Vector{Int64}})
    all = allcyclotomiccosets(q, n)
    compcosets = Vector{Vector{Int64}}()
    for a in all
        # if a != [0]
            found = false
            for b in qcosets
                if a[1] == b[1]
                    found = true
                    break
                end
            end

            if !found
                push!(compcosets, a)
            end
        # end
    end

    return compcosets
end

function qcosetpairings(arr::Vector{Vector{Int64}}, n::Integer)
    cosetreplist = Vector{Tuple{Int64, Int64}}()
    cosetpairlist = Vector{Tuple{Vector{Int64}, Vector{Int64}}}()
    for a in arr
        found = false
        for pair in cosetreplist
            if a[1] == pair[1] || a[1] == pair[2]
                found = true
            end
        end

        if !found
            neg = sort!([mod(n - i, n) for i in a])
            if neg == a
                push!(cosetreplist, (a[1], a[1]))
                push!(cosetpairlist, (a, a))
            else
                for b in arr
                    if neg == b
                        push!(cosetreplist, (a[1], b[1]))
                        push!(cosetpairlist, (a, b))
                        break
                    end
                end
            end
        end
    end

    return cosetpairlist, cosetreplist
end

function qcosetpairings(q::Integer, n::Integer)
    arr = allcyclotomiccosets(q, n, false)
    return qcosetpairings(arr, n)
end

function qcosettable(a::Integer, b::Integer, q::Integer)
    for n in a:b
        if n % q != 0
            println("n = $n")
            allcyclotomiccosets(q, n, true)
            println(" ")
        end
    end
end

function dualqcosets(q::Integer, n::Integer, qcosets::Vector{Vector{Int64}})
    compcosets = complementqcosets(q, n, qcosets)
    for a in compcosets
        for (i, x) in enumerate(a)
            a[i] = mod(n - a[i], n)
        end
        sort!(a)
    end

    return compcosets
end

max_stars_repo_name: esabo/CodingTheory

max_issues_repo_name: esabo/CodingTheory

max_forks_repo_name: esabo/CodingTheory


text:
@testset "stereo" begin

@testset "stereohydrogens" begin
    LAla1 = parse(SMILES, "N[C@@H](C)C(=O)O")
    @test LAla1.nodeattrs[2].stereo == :clockwise
    LAla1 = removestereohydrogens(LAla1)
    @test nodecount(LAla1) == 6
    @test LAla1.nodeattrs[2].stereo == :clockwise
    @test LAla1.nodeattrs[3].symbol == :C
    LAla1 = addstereohydrogens(LAla1)
    @test nodecount(LAla1) == 7
    @test LAla1.nodeattrs[2].stereo == :clockwise
    @test LAla1.nodeattrs[3].symbol == :H

    LAla2 = parse(SMILES, "N[C@@](C)([H])C(=O)O")
    LAla2 = removestereohydrogens(LAla2)
    @test nodecount(LAla2) == 6
    @test LAla2.nodeattrs[2].stereo == :anticlockwise

    LAla3 = parse(SMILES, "[H][C@@](C)(N)C(=O)O")
    LAla3 = removestereohydrogens(LAla3)
    @test nodecount(LAla3) == 6
    @test LAla3.nodeattrs[1].stereo == :clockwise

    LAla4 = parse(SMILES, "N[C@@](C)(C(=O)O)[H]")
    LAla4 = removestereohydrogens(LAla4)
    @test nodecount(LAla4) == 6
    @test LAla4.nodeattrs[2].stereo == :clockwise
end

@testset "setdiastereosmiles" begin
    mol1 = parse(SMILES, "C/C=C\\C")
    setdiastereo!(mol1)
    @test mol1.edgeattrs[2].stereo === :cis

    mol2 = parse(SMILES, "C/C=C/C")
    setdiastereo!(mol2)
    @test mol2.edgeattrs[2].stereo === :trans

    cis = parse(SMILES, "C(/C)=C/C")
    setdiastereo!(cis)
    @test cis.edgeattrs[2].stereo === :cis

    mol3 = parse(SMILES, "C\\C([H])=C([H])/C")
    setdiastereo!(mol3)
    @test mol3.edgeattrs[3].stereo === :trans

    mol4 = parse(SMILES, "C/C=C(\\C=C)/C=C(/C)C")
    setdiastereo!(mol4)
    @test mol4.edgeattrs[2].stereo === :cis
    @test mol4.edgeattrs[4].stereo === :unspecified
    @test mol4.edgeattrs[6].stereo === :trans
end

@testset "setdiastereosdfile" begin
    bonds = [
        SDFileBond(1), SDFileBond(2), SDFileBond(3),
    ]

    atoms = [
        SDFileAtom(:C, 0, 1, nothing, [-0.5, 1.41]),
        SDFileAtom(:C, 0, 1, nothing, [0.0, 0.0]),
        SDFileAtom(:C, 0, 1, nothing, [1.0, 0.0]),
        SDFileAtom(:C, 0, 1, nothing, [1.5, 1.41])
    ]
    mol1 = graphmol([(1,2), (2,3), (3,4)], atoms, bonds)
    setdiastereo!(mol1)
    @test mol1.edgeattrs[2].stereo === :cis

    atoms[4] = SDFileAtom(:C, 0, 1, nothing, [1.5, -1.41])
    mol2 = graphmol([(1,2), (2,3), (3,4)], atoms, bonds)
    setdiastereo!(mol2)
    @test mol2.edgeattrs[2].stereo === :trans

    atoms[4] = SDFileAtom(:C, 0, 1, nothing, [2.0, 0.0])
    mol3 = graphmol([(1,2), (2,3), (3,4)], atoms, bonds)
    setdiastereo!(mol3)
    @test mol3.edgeattrs[2].stereo === :unspecified
end

@testset "setstereocenter" begin
    # degree=4
    atoms = [
        SDFileAtom(:C, 0, 1, nothing, [0.0, 0.0]),
        SDFileAtom(:C, 0, 1, nothing, [-1.41, 0.5]),
        SDFileAtom(:C, 0, 1, nothing, [1.41, 0.5]),
        SDFileAtom(:C, 0, 1, nothing, [-0.5, -1.41]),
        SDFileAtom(:C, 0, 1, nothing, [0.5, -1.41])
    ]

    bonds = [
        SDFileBond(1, 0),
        SDFileBond(1, 0),
        SDFileBond(1, 0),
        SDFileBond(1, 0)
    ]
    uns1 = graphmol([(1,2), (1,3), (1,4), (1,5)], atoms, bonds)
    setstereocenter!(uns1)
    @test uns1.nodeattrs[1].stereo === :unspecified

    bonds = [
        SDFileBond(1, 0),
        SDFileBond(1, 0),
        SDFileBond(1, 1),
        SDFileBond(1, 6)
    ]
    mol1 = graphmol([(1,2), (1,3), (1,4), (1,5)], atoms, bonds)
    setstereocenter!(mol1)
    @test mol1.nodeattrs[1].stereo === :clockwise

    bonds = [
        SDFileBond(1, 1),
        SDFileBond(1, 0),
        SDFileBond(1, 6),
        SDFileBond(1, 0)
    ]
    mol2 = graphmol([(1,2), (1,3), (1,4), (1,5)], atoms, bonds)
    setstereocenter!(mol2)
    @test mol2.nodeattrs[1].stereo === :anticlockwise

    bonds = [
        SDFileBond(1, 1),
        SDFileBond(1, 0),
        SDFileBond(1, 0),
        SDFileBond(1, 1)
    ]
    mol3 = graphmol([(1,2), (1,3), (1,4), (1,5)], atoms, bonds)
    setstereocenter!(mol3)
    @test mol3.nodeattrs[1].stereo === :anticlockwise

    bonds = [
        SDFileBond(1, 6),
        SDFileBond(1, 0),
        SDFileBond(1, 0),
        SDFileBond(1, 6)
    ]
    mol4 = graphmol([(1,2), (1,3), (1,4), (1,5)], atoms, bonds)
    setstereocenter!(mol4)
    @test mol4.nodeattrs[1].stereo === :clockwise

    bonds = [
        SDFileBond(1, 0),
        SDFileBond(1, 1),
        SDFileBond(1, 0),
        SDFileBond(1, 0)
    ]
    mol5 = graphmol([(1,2), (1,3), (1,4), (1,5)], atoms, bonds)
    setstereocenter!(mol5)
    @test mol5.nodeattrs[1].stereo === :clockwise

    bonds = [
        SDFileBond(1, 0),
        SDFileBond(1, 0),
        SDFileBond(1, 6),
        SDFileBond(1, 0)
    ]
    mol6 = graphmol([(1,2), (1,3), (1,4), (1,5)], atoms, bonds)
    setstereocenter!(mol6)
    @test mol6.nodeattrs[1].stereo === :anticlockwise

    bonds = [
        SDFileBond(1, 6),
        SDFileBond(1, 6),
        SDFileBond(1, 6),
        SDFileBond(1, 0)
    ]
    wrong1 = graphmol([(1,2), (1,3), (1,4), (1,5)], atoms, bonds)
    setstereocenter!(wrong1)
    @test wrong1.nodeattrs[1].stereo === :atypical

    bonds = [
        SDFileBond(1, 1),
        SDFileBond(1, 0),
        SDFileBond(1, 1),
        SDFileBond(1, 1)
    ]
    wrong2 = graphmol([(1,2), (1,3), (1,4), (1,5)], atoms, bonds)
    setstereocenter!(wrong2)
    @test wrong2.nodeattrs[1].stereo === :atypical

    bonds = [
        SDFileBond(1, 1),
        SDFileBond(1, 0),
        SDFileBond(1, 0),
        SDFileBond(1, 6)
    ]
    wrong3 = graphmol([(1,2), (1,3), (1,4), (1,5)], atoms, bonds)
    setstereocenter!(wrong3)
    @test wrong3.nodeattrs[1].stereo === :atypical

    bonds = [
        SDFileBond(1, 1),
        SDFileBond(1, 1),
        SDFileBond(1, 0),
        SDFileBond(1, 0)
    ]
    wrong4 = graphmol([(1,2), (1,3), (1,4), (1,5)], atoms, bonds)
    setstereocenter!(wrong4)
    @test wrong4.nodeattrs[1].stereo === :atypical

    # degree=3
    atoms = [
        SDFileAtom(:C, 0, 1, nothing, [0.0, 0.0]),
        SDFileAtom(:C, 0, 1, nothing, [0.0, 1.0]),
        SDFileAtom(:C, 0, 1, nothing, [1.41, -0.5]),
        SDFileAtom(:C, 0, 1, nothing, [-1.41, -0.5])
    ]

    bonds = [
        SDFileBond(1, 0),
        SDFileBond(1, 0),
        SDFileBond(1, 0)
    ]
    uns2 = graphmol([(1,2), (1,3), (1,4)], atoms, bonds)
    setstereocenter!(uns2)
    @test uns2.nodeattrs[1].stereo === :unspecified

    bonds = [
        SDFileBond(1, 1), SDFileBond(1, 0), SDFileBond(1, 0)
    ]
    implh1 = graphmol([(1,2), (1,3), (1,4)], atoms, bonds)
    setstereocenter!(implh1)
    @test implh1.nodeattrs[1].stereo === :anticlockwise

    bonds = [
        SDFileBond(1, 0), SDFileBond(1, 6), SDFileBond(1, 0)
    ]
    implh2 = graphmol([(1,2), (1,3), (1,4)], atoms, bonds)
    setstereocenter!(implh2)
    @test implh2.nodeattrs[1].stereo === :clockwise

    bonds = [
        SDFileBond(1, 0), SDFileBond(1, 6), SDFileBond(1, 1)
    ]
    wrong5 = graphmol([(1,2), (1,3), (1,4)], atoms, bonds)
    setstereocenter!(wrong5)
    @test wrong5.nodeattrs[1].stereo === :atypical
end

@testset "sdfhydrogens" begin
    atoms = [
        SDFileAtom(:C, 0, 1, nothing, [0.0, 0.0]),
        SDFileAtom(:H, 0, 1, nothing, [0.5, 1.41]),
        SDFileAtom(:C, 0, 1, nothing, [-0.5, 1.41]),
        SDFileAtom(:N, 0, 1, nothing, [1.41, -0.5]),
        SDFileAtom(:O, 0, 1, nothing, [-1.41, -0.5])
    ]
    bonds = [
        SDFileBond(1, 1),
        SDFileBond(1, 6),
        SDFileBond(1, 0),
        SDFileBond(1, 0)
    ]
    mol1 = graphmol([(1,3), (1,2), (1,4), (1,5)], atoms, bonds)
    setstereocenter!(mol1)
    mol1 = removestereohydrogens(mol1)
    @test nodecount(mol1) == 4
    @test mol1.nodeattrs[1].stereo === :anticlockwise
    mol1 = addstereohydrogens(mol1)
    @test nodecount(mol1) == 5
    @test mol1.nodeattrs[1].stereo === :anticlockwise
end

end # stereo


text:
using DistributedFactorGraphs
using IncrementalInference
using Test

@testset "FileDFG Tests" begin
    for filename in ["/tmp/fileDFG", "/tmp/FileDFGExtension.tar.gz"]
        global dfg
        dfg = DistributedFactorGraphs._getDuplicatedEmptyDFG(dfg)

        if typeof(dfg) <: CloudGraphsDFG
            @warn "TEST: Nuking all data for user '$(dfg.userId)', robot '$(dfg.robotId)'!"
            clearRobot!!(dfg)
        end

        # Same graph as iifInterfaceTests.jl
        numNodes = 5

        #change ready and solvable for x7,x8 for improved tests on x7x8f1
        verts = map(n -> addVariable!(dfg, Symbol("x$n"), ContinuousScalar, labels = [:POSE]), 1:numNodes)
        #TODO fix this to use accessors
        verts[3].solvable = 1
        verts[4].solvable = 0
        getSolverData(verts[4]).solveInProgress = 1
        #call update to set it on cloud
        updateVariable!(dfg, verts[3])
        updateVariable!(dfg, verts[4])

        setSolvedCount!(verts[1], 5)
        # Add some bigData to x1, x2
        addBigDataEntry!(verts[1], GeneralBigDataEntry(:testing, :testing; mimeType="application/nuthin!"))
        addBigDataEntry!(verts[2], FileBigDataEntry(:testing2, "/dev/null"))
        #call update to set it on cloud
        updateVariable!(dfg, verts[1])
        updateVariable!(dfg, verts[2])

        facts = map(n -> addFactor!(dfg, [verts[n], verts[n+1]], LinearConditional(Normal(50.0,2.0))), 1:(numNodes-1))

        #FIXME a lot of these makes test fail
        # facts[4].solverData.solveInProgress = 1
        # facts[4].solverData.multihypo = [1, 0.1, 0.9]
        # facts[4].solverData.eliminated = true
        # facts[4].solverData.potentialused = true
        # updateFactor!(dfg, facts[4])
        # Save and load the graph to test.
        saveDFG(dfg, filename)

        copyDfg = DistributedFactorGraphs._getDuplicatedEmptyDFG(dfg)
        @info "Going to load $filename"

        #FIXME loadDFG!(copyDfg, filename), but use this to test deprecation
        retDFG = loadDFG(filename, Main, copyDfg) #, loaddir="/tmp")

        @test issetequal(ls(dfg), ls(retDFG))
        @test issetequal(lsf(dfg), lsf(retDFG))
        for var in ls(dfg)
            @test getVariable(dfg, var) == getVariable(retDFG, var)
        end
        for fact in lsf(dfg)
            @test getFactor(dfg, fact) == getFactor(retDFG, fact)
        end

        @test length(getBigDataEntries(getVariable(retDFG, :x1))) == 1
        @test typeof(getBigDataEntry(getVariable(retDFG, :x1),:testing)) == GeneralBigDataEntry
        @test length(getBigDataEntries(getVariable(retDFG, :x2))) == 1
        @test typeof(getBigDataEntry(getVariable(retDFG, :x2),:testing2)) == FileBigDataEntry
    end
end


text:
# Experiment to see how the motion-compensated diffusion sequences
# affects microstructural signals.

using MRIsim, JLD2, ProgressMeter
# using PlotlyJS: savefig

Ls = (1:4:120)*1e-6
bs = 0:.25:10
E = zeros(ComplexF64,length(Ls),length(bs))
Emc = zeros(ComplexF64,length(Ls),length(bs))
@load "M0MC.jld2"
DIF1 = DIF
@load "M2MXMC.jld2"
DIF2 = DIF
##Timings
@showprogress for (j,b) in enumerate(bs)
    for (i,L) in enumerate(Ls)
        μ = MRIsim.Planes(L)
        b1 = MRIsim.get_bvalue(DIF1)*1e-6
        seq = DIF1 * b
        E[i,j] = MRIsim.SignalE(μ, seq)
        b2 = MRIsim.get_bvalue(DIF2)*1e-6
        seq = DIF2 * sqrt(b1/b2) * b
        Emc[i,j] = MRIsim.SignalE(μ, seq)
    end
end
@save "ResultsMC.jld2" E Emc
##
# p = MRIsim.plot_grads_moments(DIF1,title="M0, b=$(round(MRIsim.get_bvalue(DIF1)*1e-6, digits=2)) s/mm2",mode="hd")
# PlotlyJS.savefig(p,"M0.pdf";width=900,height=250)
# ##
# b1 = MRIsim.get_bvalue(DIF1)*1e-6
# b2 = MRIsim.get_bvalue(DIF2)*1e-6
# DIF = DIF2 * sqrt(b1/b2)
# p = MRIsim.plot_grads_moments(DIF,title="M2 M1 MX MC, b=$(round(MRIsim.get_bvalue(DIF)*1e-6, digits=2)) s/mm2",mode="hd")
# PlotlyJS.savefig(p,"M1M2MC.pdf";width=900,height=250)
##
@load "ResultsMC.jld2" E Emc

using Plots#, Interact, Blink
b1 = MRIsim.get_bvalue(DIF1)*1e-6
# w = Window()
anim = @animate for l = 1:9
    D = 2e-9
    b = b1 * 1e6 * bs[l]
    p = plot(Ls*1e6,abs.(E[:,l]),label="M0",legend=:topright,ylim=(0,1),xlim=(0,30))
    plot!(Ls*1e6,abs.(Emc[:,l]),label="M2-Maxwell")
    plot!(Ls*1e6,abs.(ones(size(Emc[:,l])) .* exp(-b*D)),label="Free diffusion")
    xlabel!("L [μm]")
    ylabel!("E(q)")
    title!("Comparison (b=$(round(b*1e-6)) s/mm2)")
end
# body!(w,p)
gif(anim, "MC_L.gif", fps = 1)
# Plots.savefig(p,"Eq.pdf")

## Plot
anim = @animate for l = 1:10
    D = 2e-9
    p = plot(b1*bs, abs.(E[l,:]),label="M0",legend=:topright,xlim=(0,2000))# ,yaxis=:log)
    plot!(b1*bs, abs.(Emc[l,:]),label="M2-Maxwell")
    plot!(b1*bs, exp.(-b1*bs*1e6*D),label="Free diffusion")
    xlabel!("b [s/mm2]")
    ylabel!("E(q)")
    title!("Comparison (L=$(round(Ls[l]*1e6)) um)")
end
gif(anim, "MC_B.gif", fps = 2)

## Effect of rotationa angle
θs = -π/2:.1:π/2
E = zeros(ComplexF64,length(θs))
Emc = zeros(ComplexF64,length(θs))
f = 2
L = 20e-6
@showprogress for (i,θ) in enumerate(θs)
        μ = MRIsim.Planes(L)
        R = MRIsim.rotz(θ)
        b1 = MRIsim.get_bvalue(DIF1)*1e-6
        seq = DIF1 * f * R
        E[i] = MRIsim.SignalE(μ, seq)
        b2 = MRIsim.get_bvalue(DIF2)*1e-6
        seq = DIF2 * sqrt(b1/b2) * f * R
        Emc[i] = MRIsim.SignalE(μ, seq)
end
##
using Plots
f = 2
D = 2e-9
b = MRIsim.get_bvalue(DIF1)
p = plot(θs,abs.(E[:]),label="M0",legend=:topright,ylim=(0,1))
plot!(θs,abs.(Emc[:]),label="M2-Maxwell")
plot!(θs,abs.(ones(size(Emc[:])) .* exp(-b*f*D)),label="Free diffusion")
xlabel!("θ [rad]")
ylabel!("E(q)")
title!("Comparison (L=20 um, b=$(round(b*1e-6*f)) s/mm2)")


text:
# Good source of mappings from MiniZinc:
# https://github.com/radsz/jacop/blob/develop/src/main/java/org/jacop/fz/constraints/GlobalConstraints.java

# CP.AllDifferent
function MOI.supports_constraint(
    ::Optimizer,
    ::Type{F},
    ::Type{S},
) where {
    F <: MOI.VectorOfVariables,
    S <: CP.AllDifferent,
}
    return true
end

function _build_constraint(
    model::Optimizer,
    f::MOI.VectorOfVariables,
    ::CP.AllDifferent,
)
    # This could also be Alldiff or Alldistinct, with various algorithms
    # for filtering behind. TODO: add a parameter to choose?
    # http://jacopguide.osolpro.com/guideJaCoP.pdf#subsection.3.3.1
    # fzn-jacop uses Alldiff:
    # https://github.com/radsz/jacop/blob/develop/src/main/java/org/jacop/fz/constraints/GlobalConstraints.java#L279
    return Alldifferent(_parse_to_vars(model, f))
end

# CP.Domain
# Cannot be reified, hence no _build_constraint!
# Only for VariableIndex-in-Domain.
# Not implemented using In, because IntDomain is too limited (only a few
# values are possible in SmallDenseDomain or intervals).

# # CP.AntiDomain
# function MOI.supports_constraint(
#     ::Optimizer,
#     ::Type{F},
#     ::Type{S},
# ) where {
#     T <: Int,
#     F <: Union{MOI.VariableIndex, MOI.ScalarAffineFunction{T}},
#     S <: CP.AntiDomain{T},
# }
#     return true
# end

# function _build_constraint(
#     model::Optimizer,
#     f::Union{MOI.VariableIndex, MOI.ScalarAffineFunction{T}},
#     s::CP.AntiDomain{T},
# ) where {T <: Int}
#     return cpo_java_forbiddenassignments(
#         model.inner,
#         _parse(model, f),
#         collect(Int32(v) for v in s.values),
#     )
# end

# # CP.DifferentFrom
# function MOI.supports_constraint(
#     ::Optimizer,
#     ::Type{F},
#     ::Type{S},
# ) where {
#     T <: Int,
#     F <: Union{MOI.VariableIndex, MOI.ScalarAffineFunction{T}},
#     S <: CP.DifferentFrom{T},
# }
#     return true
# end

# function _build_constraint(
#     model::Optimizer,
#     f::Union{MOI.VariableIndex, MOI.ScalarAffineFunction{T}},
#     s::CP.DifferentFrom{T},
# ) where {T <: Int}
#     return cpo_java_neq(model.inner, _parse(model, f), Int32(s.value))
# end

# # CP.Count
# function MOI.supports_constraint(
#     ::Optimizer,
#     ::Type{F},
#     ::Type{S},
# ) where {
#     T <: Int,
#     F <: Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
#     S <: CP.Count{MOI.EqualTo{T}},
# }
#     return true
# end

# function _build_constraint(
#     model::Optimizer,
#     f::Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
#     s::CP.Count{MOI.EqualTo{T}},
# ) where {T <: Int}
#     @assert MOI.output_dimension(f) >= 2

#     f_parsed = _parse(model, f)
#     count_assign = f_parsed[1]
#     count_values = f_parsed[2:end]

#     expr = cpo_java_count(model.inner, count_values, Int32(s.set.value))
#     return cpo_java_eq(model.inner, count_assign, expr)
# end

# # CP.CountDistinct
# function MOI.supports_constraint(
#     ::Optimizer,
#     ::Type{F},
#     ::Type{S},
# ) where {
#     T <: Int,
#     F <: Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
#     S <: CP.CountDistinct,
# }
#     return true
# end

# function _build_constraint(
#     model::Optimizer,
#     f::Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
#     ::CP.CountDistinct,
# ) where {T <: Int}
#     @assert MOI.output_dimension(f) >= 2

#     f_parsed = _parse(model, f)
#     count_assign = f_parsed[1]
#     count_values = f_parsed[2:end]

#     expr = cpo_java_countdifferent(model.inner, count_values)
#     return cpo_java_eq(model.inner, count_assign, expr)
# end

# # CP.Strictly
# function MOI.supports_constraint(
#     ::Optimizer,
#     ::Type{F},
#     ::Type{S},
# ) where {
#     T <: Int,
#     Sense <: Union{MOI.LessThan{T}, MOI.GreaterThan{T}},
#     F <: Union{MOI.VariableIndex, MOI.ScalarAffineFunction{T}},
#     S <: CP.Strictly{Sense, T},
# }
#     return true
# end

# function _build_constraint(
#     model::Optimizer,
#     f::Union{MOI.VariableIndex, MOI.ScalarAffineFunction{T}},
#     s::CP.Strictly{MOI.LessThan{T}, T},
# ) where {T <: Int}
#     return cpo_java_lt(
#         model.inner,
#         _parse(model, f),
#         cpo_java_constant(model.inner, s.set.upper),
#     )
# end

# function _build_constraint(
#     model::Optimizer,
#     f::Union{MOI.VariableIndex, MOI.ScalarAffineFunction{T}},
#     s::CP.Strictly{MOI.GreaterThan{T}, T},
# ) where {T <: Int}
#     return cpo_java_gt(
#         model.inner,
#         _parse(model, f),
#         cpo_java_constant(model.inner, s.set.lower),
#     )
# end

# # CP.Element
# function MOI.supports_constraint(
#     ::Optimizer,
#     ::Type{F},
#     ::Type{S},
# ) where {
#     T <: Int,
#     F <: Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
#     S <: CP.Element{T},
# }
#     return true
# end

# function _build_constraint(
#     model::Optimizer,
#     f::Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
#     s::CP.Element{T},
# ) where {T <: Int}
#     # Split the dimensions in the right parts.
#     f_parsed = _parse(model, f)
#     element_assign = f_parsed[2]
#     element_index = f_parsed[1]

#     # Build the constraint.
#     expr = cpo_java_element(
#         model.inner,
#         convert(Vector{Int32}, s.values),
#         element_index,
#     )
#     return cpo_java_eq(model.inner, element_assign, expr)
# end

# # CP.BinPacking
# # Unlike other CP solvers, the item weights are fixed (i.e. constant expressions, not variables).
# function MOI.supports_constraint(
#     ::Optimizer,
#     ::Type{F},
#     ::Type{S},
# ) where {
#     T <: Int,
#     F <: Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
#     S <: CP.BinPacking{CP.NO_CAPACITY_BINPACKING, T},
# }
#     return true
# end

# function _build_constraint(
#     model::Optimizer,
#     f::Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
#     s::CP.BinPacking{CP.NO_CAPACITY_BINPACKING, T},
# ) where {T <: Int}
#     f = MOI.Utilities.canonical(f)

#     @assert MOI.output_dimension(f) == s.n_bins + s.n_items

#     # Split the dimensions in the right parts.
#     f_parsed = _parse(model, f)
#     pack_load = f_parsed[1:(s.n_bins)]
#     pack_assigned = f_parsed[(s.n_bins + 1):(s.n_bins + s.n_items)]

#     # Build the constraint.
#     return cpo_java_pack(
#         model.inner,
#         pack_load,
#         pack_assigned,
#         convert(Vector{Int32}, s.weights),
#     )
# end

# # CP.MinimumDistance
# function MOI.supports_constraint(
#     ::Optimizer,
#     ::Type{F},
#     ::Type{S},
# ) where {
#     T <: Int,
#     F <: Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
#     S <: CP.MinimumDistance{T},
# }
#     return true
# end

# function _build_constraint(
#     model::Optimizer,
#     f::MOI.VectorOfVariables,
#     s::CP.MinimumDistance{T},
# ) where {T <: Int}
#     # The Java API does not support single variables as argument, only expressions.
#     return _build_constraint(model, MOI.VectorAffineFunction{T}(f), s)
# end

# function _build_constraint(
#     model::Optimizer,
#     f::MOI.VectorAffineFunction{T},
#     s::CP.MinimumDistance{T},
# ) where {T <: Int}
#     f_parsed = _parse(model, f)
#     return cpo_java_allmindistance(model.inner, f_parsed, Int32(s.k))
# end

# # CP.Inverse
# function MOI.supports_constraint(
#     ::Optimizer,
#     ::Type{F},
#     ::Type{S},
# ) where {
#     T <: Int,
#     F <: Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
#     S <: CP.Inverse,
# }
#     return true
# end

# function _build_constraint(
#     model::Optimizer,
#     f::Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
#     s::CP.Inverse,
# ) where {T <: Int}
#     f_parsed = _parse(model, f)

#     # Map indices between Julia (starts at 1) and CPLEX (starts at 0).
#     f_parsed = [cpo_java_diff(model.inner, f_p, Int32(1)) for f_p in f_parsed]

#     inverse_first = f_parsed[1:(s.dimension)]
#     inverse_second = f_parsed[(1 + s.dimension):(2 * s.dimension)]

#     return cpo_java_inverse(model.inner, inverse_first, inverse_second)
# end

# # CP.LexicographicallyLessThan
# function MOI.supports_constraint(
#     ::Optimizer,
#     ::Type{F},
#     ::Type{S},
# ) where {
#     T <: Int,
#     F <: Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}},
#     S <: Union{
#         CP.LexicographicallyLessThan,
#         CP.Strictly{CP.LexicographicallyLessThan},
#     },
# }
#     return true
# end

# function _build_constraint(
#     model::Optimizer,
#     f::MOI.VectorOfVariables,
#     s::CP.LexicographicallyLessThan,
# )
#     # The Java API does not support single variables as argument, only expressions.
#     return _build_constraint(model, MOI.VectorAffineFunction{Int}(f), s)
# end

# function _build_constraint(
#     model::Optimizer,
#     f::MOI.VectorAffineFunction{T},
#     s::CP.LexicographicallyLessThan,
# ) where {T <: Int}
#     # _build_constraint can only return one constraint at a time, CPLEX only
#     # has lexicographic constraints between two vectors.
#     @assert s.row_dim == 2

#     f_parsed = _parse(model, f)
#     lex_first = f_parsed[1:(s.column_dim)]
#     lex_second = f_parsed[(s.column_dim + 1):(2 * s.column_dim)]

#     return cpo_java_lexicographic(model.inner, lex_first, lex_second)
# end

# function _build_constraint(
#     model::Optimizer,
#     f::MOI.VectorOfVariables,
#     s::CP.Strictly{CP.LexicographicallyLessThan},
# )
#     # The Java API does not support single variables as argument, only expressions.
#     return _build_constraint(model, MOI.VectorAffineFunction{Int}(f), s)
# end

# function _build_constraint(
#     model::Optimizer,
#     f::MOI.VectorAffineFunction{T},
#     s::CP.Strictly{CP.LexicographicallyLessThan},
# ) where {T <: Int}
#     # _build_constraint can only return one constraint at a time, CPLEX only
#     # has strict lexicographic constraints between two vectors.
#     @assert s.set.row_dim == 2

#     f_parsed = _parse(model, f)
#     lex_first = f_parsed[1:(s.set.column_dim)]
#     lex_second = f_parsed[(s.set.column_dim + 1):(2 * s.set.column_dim)]

#     return cpo_java_strictlexicographic(model.inner, lex_first, lex_second)
# end


text:
systemid=parse(Int8,ARGS[1])
tend=parse(Float64,ARGS[2])
atol=parse(Float64,ARGS[3])
rtol=parse(Float64,ARGS[4])

using LinearAlgebra: BLAS
BLAS.set_num_threads(1) # Recording CPU time means we should
                        # stick to one thread

# Load data
load_data = :(begin
    using SQLite
    db = SQLite.DB(ENV["KURABENCH_DB"])

    # First read the system info
    results = DBInterface.execute(db,string(
            "select ",
                "name, nodes, edges, coupling_constant " ,
            "from ",
                "systems ",
            "where " ,
                "id=$systemid;"))
    name = ""
    nodes = 0
    edges = 0
    coupling_constant = 0.0
    for row in results
        global name = row.name
        global nodes = row.nodes
        global edges = row.edges
        global coupling_constant = row.coupling_constant
    end
    # Now read connectivity
    # Allocate
    sources = Array{Int32}(undef,edges)
    dests = Array{Int32}(undef,edges)
    data = ones(Int32,edges)

    results = DBInterface.execute(db,string(
            "select ",
                "id,source,dest ",
             "from " ,
                "edges " ,
             "where " ,
                "system=$systemid;"))

    for (i,row) in enumerate(results)
        sources[i] = row.source
        dests[i]   = row.dest
    end
    # Store the node omegas
    omegas = Array{Float64}(undef,nodes)
    x0     = Array{Float64}(undef,nodes)
    results = DBInterface.execute(db,string(
            "select ",
                "id,omega,initialCondition ",
            "from " ,
                "nodes " ,
            "where " ,
                "system=$systemid;"))
    for (i,row) in enumerate(results)
        omegas[i] = row.omega
        x0[i] = row.initialCondition
    end
    my_matrix = sparse(sources,dests,data)
end)

prepare_ode = :(begin
    g = SimpleGraph(my_matrix')
    integrator = radau
    ### Network dynamics vertex and edge functions

    @inline function kuramoto_vertex!(dv, v, edges, p, t)
        dv[1] = p # parameter for vertices are eigenfrequencies omega
        @inbounds for e in edges
            dv[1] += e[1]
        end
        nothing
    end

    @inline function kuramoto_edge!(e,v_s,v_d,p,t)
        e[1] = p * sin(v_s[1] - v_d[1])
        nothing
    end

    ### Constructing the network dynamics

    odevertex = ODEVertex(f! = kuramoto_vertex!, dim = 1)
    staticedge = StaticEdge(f! = kuramoto_edge!, dim = 1, coupling = :antisymmetric)

    # generating random values between -0.5 and 0.5 for the parameter value ω_0 of the vertices
    e_pars = [coupling_constant for i in 1:g.ne]
    parameters = (omegas, e_pars)

    # setting up the  network_dynamics
    kuramoto_network! = network_dynamics(odevertex, staticedge, g)

    ### Simulation
    # constructing random initial conditions for nodes (variable θ)

    prob = ODEProblem(kuramoto_network!, x0, (0.,tend), parameters, abstol=atol, reltol=rtol)
end)

# Precompile eval
eval(:nothing)

using Pkg
Pkg.activate(@__DIR__)
println("Enviroment Activated, starting program...")
using CPUTime

# Times
tstartup = @CPUelapsed begin
	using NetworkDynamics: ODEVertex, StaticEdge, network_dynamics
	using LightGraphs: SimpleGraph
	using OrdinaryDiffEq: ODEProblem, solve, DP5
    using SparseArrays
    using ODEInterfaceDiffEq
end
times = 0.:2.:tend
println("time to startup is $tstartup")
tload    = @CPUelapsed(eval(load_data))
println("time to load is $tload")
tjit     = @CPUelapsed(eval(prepare_ode)) +
           @CPUelapsed(solve(prob, integrator(), saveat=times))
           # Note: We don't need ot run a whole simulation to compile the
           # specalized functions. Just one small step will suffice.
           # This makes the comparison to the JIT step in jitcode more direct.

println("time to jit is $tjit")
tintegrate = @CPUelapsed(sol = solve(prob, integrator(), saveat=times))

println("time to solve is $tintegrate")

ttotal = tstartup + tload + tjit + tintegrate
# Now write the config to the experiments table and the states to the states table
#

## Config
tidy = cmd -> chomp(String(read(cmd)))

githash = tidy(`git rev-parse HEAD`)
hostname = tidy(`hostname`)
datetime = tidy(`date`)
gitchanges = String(read(`git status -s -uno`))
argsstring = join(ARGS," ")
intstring = string(integrator)
DBInterface.execute(db,string(
    "insert into experiments ",
    "(runtime,datetime,githash,gitchanges,commandline,host,",
    "system,solver,atol,rtol,tstart,tend,",
    "tstartup,tload,tjit,tintegration,ttotal)",
    "values ",
    "(\"networkdynamics.jl\",",
    "\"$datetime\",\"$githash\",\"$gitchanges\",\"$argsstring\",\"$hostname\",",
    "$systemid,\"$intstring\",$atol,$rtol,0,$tend,",
    "$tstartup,$tload,$tjit,$tintegrate,$ttotal);"
                              )
                   )

results = DBInterface.execute(db,"select last_insert_rowid()")
experimentid = 0
for row in results
    # TODO: How to get a single element from julia iterator?
    global experimentid = row[1]
end

for (i,t) in enumerate(sol.t)
    for (k,state) in enumerate(sol[i])
        state = mod(state,2*pi)
        DBInterface.execute(db,string(
            "insert into states",
            "(experimentid,time,idx,node,state)",
            "values",
            "($experimentid,$t,$k,$k,$state);"
                                     )
                           )
    end
end

max_stars_repo_name: PIK-ICoNe/NetworkDynamicsBenchmarks

max_issues_repo_name: PIK-ICoNe/NetworkDynamicsBenchmarks

max_forks_repo_name: PIK-ICoNe/NetworkDynamicsBenchmarks


text:
using DelaunayMeshes
using Base.Test

# write your own tests here
include("TestHelpers.jl")
#include("QuadEdgeTests.jl")
#include("TriangulationTests.jl")
#include("DelaunayMeshTests.jl")
include("PostProcessTests.jl")


text:
import DCELs._pseudoangle, DCELs._pseudoangle_unsigned, DCELs._is_sorted_ccw

@testset "Pseudoangle" begin
    @test _pseudoangle(1, 0) == 0
    @test _pseudoangle(1, 1) == 0.5
    @test _pseudoangle(0, 1) == 1
    @test _pseudoangle(-1, 1) == 1.5
    @test _pseudoangle(-1, 0) == 2
    @test _pseudoangle(-1, -1) == -1.5
    @test _pseudoangle(0, -1) == -1
    @test _pseudoangle(1, -1) == -0.5
end

@testset "Pseudoangle unsigned" begin
    @test _pseudoangle_unsigned(1, 0) == 0
    @test _pseudoangle_unsigned(1, 1) == 0.5
    @test _pseudoangle_unsigned(0, 1) == 1
    @test _pseudoangle_unsigned(-1, 1) == 1.5
    @test _pseudoangle_unsigned(-1, 0) == 2
    @test _pseudoangle_unsigned(-1, -1) == 2.5
    @test _pseudoangle_unsigned(0, -1) == 3
    @test _pseudoangle_unsigned(1, -1) == 3.5
end

@testset "Sorting" begin
    u = Vertex(1, 0)
    v = Vertex(0, 1)
    w = Vertex(-1, 0)
    center = Vertex(0, 0)

    @test _is_sorted_ccw(u, v, w, center)
    @test !_is_sorted_ccw(u, w, v, center)
end


text:
abstract type Kroki <: ShortCode end

function Base.show(io::IO, ::MIME"text/plain", kroki::Kroki)
    print(io, kroki.text * "\n" * kroki(kroki.text, lowercase(String(nameof(typeof(kroki)))), "svg"))
end

function Base.show(io::IO, ::MIME"image/svg+xml", kroki::Kroki)
    write(io, fetch_kroki(kroki.text, lowercase(String(nameof(typeof(kroki)))), "svg"))
end

struct GraphViz <: Kroki
    text::String
end

struct Mermaid <: Kroki
    text::String
end

struct PlantUML <: Kroki
    text::String
end

struct BlockDiag <: Kroki
    text::String
end

function kroki(krokitext, generator="graphviz", format="svg")
    encoded = base64urlencode(compress(krokitext))
    url = "https://kroki.io/$generator/$format/$encoded"
end

@memoize function fetch_kroki(krokitext, generator="graphviz", format="svg")
    response = HTTP.get(kroki(krokitext, generator, format))
    String(response.body)
end

function base64urlencode(text)
    str = base64encode(text)
    return replace(replace(str, "/"=>"_"),"+"=>"-")
end

function compress(text)
    return transcode(ZlibCompressor, text)
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ e9ad66b0-0d6b-11eb-26c0-27413c19dd32
begin
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.add([
		"Plots",
		"PlutoUI",
		"LaTeXStrings",
	])
	using LaTeXStrings
	using Plots
	using PlutoUI

	md"#### Package dependencies"
end

# ╔═╡ 05031b60-1df4-11eb-2b61-956e526b3d4a
md"## Lecture 21: Snowball Earth, the ice-albedo feedback, and multiple equilibria"

# ╔═╡ a05bd3ce-1e06-11eb-1af2-65886cab38ef
md"![](https://static01.nyt.com/images/2019/12/17/science/02TB-SNOWBALLEARTH1/02TB-SNOWBALLEARTH1-superJumbo-v2.jpg?quality=90&auto=webp)
[Source (New York Times)](https://static01.nyt.com/images/2019/12/17/science/02TB-SNOWBALLEARTH1/02TB-SNOWBALLEARTH1-superJumbo-v2.jpg?quality=90&auto=webp)
"

# ╔═╡ a6b6ea46-2519-11eb-1134-8994582c5546
md"""
### 0) Review of Lecture 20

Recall from [Lecture 20 (Part I)](https://www.youtube.com/watch?v=Gi4ZZVS2GLA&t=15s) that the the **zero-dimensional energy balance equation** is

\begin{gather}
\color{brown}{C \frac{dT}{dt}}
\; \color{black}{=} \; \color{orange}{\frac{(1 - α)S}{4}}
\; \color{black}{-} \; \color{blue}{(A - BT)}
\; \color{black}{+} \; \color{grey}{a \ln \left( \frac{[\text{CO}₂]}{[\text{CO}₂]_{\text{PI}}} \right)},
\end{gather}
"""

# ╔═╡ 9d3a6312-2519-11eb-3550-89ea034bf119
html""""""

# ╔═╡ c50b6b24-2506-11eb-00e1-0f183526ed4e
md"""
Today, we will ignore changes in CO₂, so that

$\ln \left( \frac{ [\text{CO}₂]_{\text{PI}} }{[\text{CO}₂]_{\text{PI}}} \right) = \ln(1) = 0$

and the model simplifies to

\begin{gather}
\color{brown}{C \frac{dT}{dt}}
\; \color{black}{=} \; \color{orange}{\frac{(1 - α)S}{4}}
\; \color{black}{-} \; \color{blue}{(A - BT)}.
\end{gather}

The dynamics of this **Ordinary Differential Equation (ODE)** are quite simple because it is *linear*: we can rewrite it in the form

$\dot{T} = f(T(t))$ where $f(x) = \alpha x + \beta$ is a *linear* function of x. A linear ODE permits only one stable solution, $\dot{T} = 0$, which in Lecture 20 we found was Earth's pre-industrial temperature $T_{0} = 14$°C.

In this lecture, we show how a small modification that makes one term in our simple climate model non-linear completely changes its dynamics, allowing us to explain the existence of both "Snowball Earth" and the relatively warm pre-industrial climate that allowed humans to thrive.
"""

# ╔═╡ 4f5f3038-1e06-11eb-16a2-b11035701fb8
md"""
### 1) Background: Snowball Earth

Geological evidence shows that the Neoproterozoic Era (550 to 1000 million years ago) is marked by two global glaciation events, in which Earth's surface was covered in ice and snow from the Equator to the poles (see review by [Pierrehumbert et al. 2011](https://www.annualreviews.org/doi/full/10.1146/annurev-earth-040809-152447)).
"""

# ╔═╡ 68364804-2517-11eb-00de-d1655dd7754b
html"""

"""

# ╔═╡ 9c118f9a-1df0-11eb-22dd-b14428994076
html"""

"""

# ╔═╡ 38346e6a-0d98-11eb-280b-f79787a3c788
md"""
We can represent the ice-albedo feedback crudely in our energy balance model by allowing the albedo to depend on temperature:

$\alpha(T) = \begin{cases}
\alpha_{i} & \mbox{if }\;\; T \leq -10\text{°C} &\text{(completely frozen)}\\
\alpha_{i} + (\alpha_{0}-\alpha_{i})\frac{T + 10}{20} & \mbox{if }\;\; -10\text{°C} \leq T \leq 10\text{°C} &\text{(partially frozen)}\\
\alpha_{0} &\mbox{if }\;\; T \geq 10\text{°C} &\text{(no ice)}
\end{cases}$
"""

# ╔═╡ 816f1d96-2508-11eb-0873-c3b564a31dea
md"""
##### 1.2) Adding the ice-albedo feedback to our simple climate model

First, we program albedo as a function of temperature.
"""

# ╔═╡ a8dcc0fc-1df8-11eb-21fd-1fdebe5dabfc
md"""
To add this function into our energy balance model from [Lecture 20 (Part I)](https://www.youtube.com/watch?v=Gi4ZZVS2GLA&t=15s) (which we've copied into the cell below), all we have to do is overwrite the definition of the `timestep!` method to specify that the temperature-dependent albedo should be updated based on the current state:
"""

# ╔═╡ 96ed2f9a-1e29-11eb-09f4-23df52152b2f
module Model

const S = 1368; # solar insolation [W/m^2]  (energy per unit time per unit area)
const α = 0.3; # albedo, or planetary reflectivity [unitless]
const B = -1.3; # climate feedback parameter [W/m^2/°C],
const T0 = 14.; # preindustrial temperature [°C]

absorbed_solar_radiation(; α=α, S=S) = S*(1 - α)/4; # [W/m^2]
outgoing_thermal_radiation(T; A=A, B=B) = A - B*T;

const A = S*(1. - α)/4 + B*T0; # [W/m^2].

greenhouse_effect(CO2; a=a, CO2_PI=CO2_PI) = a*log(CO2/CO2_PI);

const a = 5.0; # CO2 forcing coefficient [W/m^2]
const CO2_PI = 280.; # preindustrial CO2 concentration [parts per million; ppm];
CO2_const(t) = CO2_PI; # constant CO2 concentrations

const C = 51.; # atmosphere and upper-ocean heat capacity [J/m^2/°C]

function timestep!(ebm)
	append!(ebm.T, ebm.T[end] + ebm.Δt*tendency(ebm));
	append!(ebm.t, ebm.t[end] + ebm.Δt);
end;

tendency(ebm) = (1. /ebm.C) * (
	+ absorbed_solar_radiation(α=ebm.α, S=ebm.S)
	- outgoing_thermal_radiation(ebm.T[end], A=ebm.A, B=ebm.B)
	+ greenhouse_effect(ebm.CO2(ebm.t[end]), a=ebm.a, CO2_PI=ebm.CO2_PI)
);

begin
	mutable struct EBM
		T::Array{Float64, 1}

		t::Array{Float64, 1}
		Δt::Float64

		CO2::Function

		C::Float64
		a::Float64
		A::Float64
		B::Float64
		CO2_PI::Float64

		α::Float64
		S::Float64
	end;

	# Make constant parameters optional kwargs
	EBM(T::Array{Float64, 1}, t::Array{Float64, 1}, Δt::Float64, CO2::Function;
		C=C, a=a, A=A, B=B, CO2_PI=CO2_PI, α=α, S=S) = (
		EBM(T, t, Δt, CO2, C, a, A, B, CO2_PI, α, S)
	);

	# Construct from float inputs for convenience
	EBM(T0::Float64, t0::Float64, Δt::Float64, CO2::Function;
		C=C, a=a, A=A, B=B, CO2_PI=CO2_PI, α=α, S=S) = (
		EBM([T0], [t0], Δt, CO2;
			C=C, a=a, A=A, B=B, CO2_PI=CO2_PI, α=α, S=S);
	);
end;

begin
	function run!(ebm::EBM, end_year::Real)
		while ebm.t[end] < end_year
			timestep!(ebm)
		end
	end;

	run!(ebm) = run!(ebm, 200.) # run for 200 years by default
end

CO2_hist(t) = CO2_PI * (1 .+ fractional_increase(t));
fractional_increase(t) = ((t .- 1850.)/220).^3;

begin
	hist = EBM(T0, 1850., 1., CO2_hist,
		C=C, B=B, A=A
	)
	run!(hist, 2020.)
end

begin
	CO2_RCP26(t) = CO2_PI * (1 .+ fractional_increase(t) .* min.(1., exp.(-((t .-1850.).-170)/100))) ;
	RCP26 = EBM(T0, 1850., 1., CO2_RCP26)
	run!(RCP26, 2100.)

	CO2_RCP85(t) = CO2_PI * (1 .+ fractional_increase(t) .* max.(1., exp.(((t .-1850.).-170)/100)));
	RCP85 = EBM(T0, 1850., 1., CO2_RCP85)
	run!(RCP85, 2100.)
end

end

# ╔═╡ 016c1074-1df4-11eb-2da8-578e25d9456b
md"""##### 1.1) The ice-albedo feedback

In Lecture 20, we used a **constant** value $α =$ $(Model.hist.α) for Earth's planetary albedo, which is a reasonable thing to do for small climate variations relative to the present (such as the difference between the present-day and preindustrial climates). In the case of large variations, however, this approximation is not very reliable.

While oceans are dark and absorbant, $α_{ocean} \approx 0.05$, ice and snow are bright and reflective: $\alpha_{ice,\,snow} \approx 0.5$ to $0.9$. Thus, if much of the ocean's surface freezes over, we expect Earth's albedo to rise dramatically, causing more sunlight to be reflected to space, which in turn causes even more cooling and more of the ocean to freeze, etc. This *non-linear positive feedback effect* is referred to as the **ice-albedo feedback** (see illustration below).
"""

# ╔═╡ 262fc3c6-1df2-11eb-332d-c1c9561b3710
function α(T; α0=Model.α, αi=0.5, ΔT=10.)
	if T < -ΔT
		return αi
	elseif -ΔT <= T < ΔT
		return αi + (α0-αi)*(T+ΔT)/(2ΔT)
	elseif T >= ΔT
		return α0
	end
end

# ╔═╡ f7761e40-1e23-11eb-2741-cfebfaf434ec
begin
	T_example = -20.:1.:20.
	plot(size=(500, 230), ylims=(0.2, 0.6))
	plot!([-20, -10], [0.2, 0.2], fillrange=[0.6, 0.6], color=:lightblue, alpha=0.2, label=nothing)
	plot!([10, 20], [0.2, 0.2], fillrange=[0.6, 0.6], color=:red, alpha=0.12, label=nothing)
	plot!(T_example, α.(T_example), lw=3., label="α(T)", color=:black)
	plot!(ylabel="albedo α\n(planetary reflectivity)", xlabel="Temperature [°C]")
	annotate!(-15.5, 0.252, text("completely\nfrozen", 10, :darkblue))
	annotate!(15.5, 0.252, text("no ice", 10, :darkred))
	annotate!(-0.3, 0.252, text("partially frozen", 10, :darkgrey))
end

# ╔═╡ 872c8f2a-1df1-11eb-3cfc-3dd568926442
function Model.timestep!(ebm)
	ebm.α = α(ebm.T[end]) # Added this line
	append!(ebm.T, ebm.T[end] + ebm.Δt*Model.tendency(ebm));
	append!(ebm.t, ebm.t[end] + ebm.Δt);
end

# ╔═╡ 13f42334-1e27-11eb-11a0-f51af4574a6b
md"""### 2) Multiple Equilibria
**OR: the existence of "alternate Earths"**

Human civilization flourished over the last several thousand years in part because Earth's global climate has been remarkably stable and forgiving. The preindustrial combination of natural greenhouse effect and incoming solar radiation yielded temperatures between the freezing and boiling points of water across most of the planet, allowing ecoystems based on liquid water to thrive.

The climate system, however, is rife with non-linear effects like the **ice-albedo effect**, which reveal just how fragile our habitable planet is and just how unique our stable pre-industrial climate was.

We learned in Lecture 20 that in response to temperature fluctuations, *net-negative feedbacks* act to restore Earth's temperature back towards a single equilibrium state in which absorbed solar radiation is balanced by outgoing thermal radiation. Here, we explore how *non-linear positive feedbacks* can temporarily result in a *net-positive feedback* and modify Earth's state space.
"""

# ╔═╡ 03292b48-2503-11eb-1514-5d1923d1d9a2
md"""
##### 2.1) Exploring the non-linear ice-albedo feedback

In [Lecture 20 (Part II)](https://www.youtube.com/watch?v=D3jpfeQCISU), we learned how introducing non-linear terms in *ordinary differential equations* can lead to complex state spaces that allow for multiple fixed points (e.g. for $\dot{x} = \mu + x^{2}$).

Let's explore how this plays out with the non-linear ice-albedo feedback by varying the initial condition $T_{0} \equiv T(t=0)$ and allowing the system to evolve for 200 years.
"""

# ╔═╡ 7a1b3138-2503-11eb-1165-c399836b66a7
begin
	T0_slider = @bind T0_interact Slider(-60.:0.25:30., default=24., show_value=true)
	md""" T₀ = $(T0_slider) °C"""
end

# ╔═╡ 1be1cce0-251b-11eb-35f7-e15741b0a712
begin
	ebm_interact = Model.EBM(Float64(T0_interact), 0., 1., Model.CO2_const)
	Model.run!(ebm_interact, 200)
end

# ╔═╡ a954ce70-2510-11eb-0d8c-c7b3f0fb3a06
md"We can get an overview of the behavior by drawing a set of these curves all on the same graph:"

# ╔═╡ 94852946-1e25-11eb-3425-210be17c23cd
begin
	p_equil = plot(xlabel="year", ylabel="temperature [°C]", legend=:bottomright, xlims=(0,205), ylims=(-60, 30.))

	plot!([0, 200], [-60, -60], fillrange=[-10., -10.], fillalpha=0.3, c=:lightblue, label=nothing)
	annotate!(120, -20, text("completely frozen", 10, :darkblue))

	plot!([0, 200], [10, 10], fillrange=[30., 30.], fillalpha=0.09, c=:red, lw=0., label=nothing)
	annotate!(120, 25, text("no ice", 10, :darkred))
	for T0_sample in (-60.:5.:30.)
		ebm = Model.EBM(T0_sample, 0., 1., Model.CO2_const)
		Model.run!(ebm, 200)

		plot!(p_equil, ebm.t, ebm.T, label=nothing)
	end

	T_un = 7.5472
	for δT in 1.e-2*[-2, -1., 0., 1., 2.]
		ebm_un = Model.EBM(T_un+δT, 0., 1., Model.CO2_const)
		Model.run!(ebm_un, 200)

		plot!(p_equil, ebm_un.t, ebm_un.T, label=nothing, linestyle=:dash)
	end

	plot!(p_equil, [200], [Model.T0], markershape=:circle, label="Our pre-industrial climate (stable ''warm'' branch)", color=:orange, markersize=8)
	plot!(p_equil, [200], [-38.3], markershape=:circle, label="Alternate universe pre-industrial climate (stable ''cold'' branch)", color=:aqua, markersize=8)
	plot!(p_equil, [200], [T_un], markershape=:diamond, label="Impossible alternate climate (unstable branch)", color=:lightgrey, markersize=8, markerstrokecolor=:white, alpha=1., markerstrokestyle=:dash)
	p_equil
end

# ╔═╡ c81b4956-2510-11eb-388c-87e93923e45a
md"We see that for T₀ ⪆  $(round(T_un, digits=2)) °C, all of the curves seem to converge on the T = 14°C *equilibrium* (or *fixed point*) that we saw in Lecture 20. Curves that start below this value warm up and while curves that start above this value will cool down. For T₀ ⪅ $(round(T_un, digits=2)) °C, however, the temperatures converge on a much colder equilibrium around T = -40°C. This is the **Snowball Earth** equilibrium. These two states are referred to as **stable equilibria** because even if the state gets temporarily pushed slightly away from its equilibrium, it will eventually converge right back to its equilibrium.

So what happens is T₀ ≈ $(round(T_un, digits=2)) °C? For some exact temperature near there, there is indeed an equilibrim state: if you start with that temperature you will stay there forever. However, if the temperature starts off even *one one-hundredth of a degree* above or below this exact value, we see that temperatures eventually converge to one of the other two equilibria. Thus, we call this intermediate equilibrium an **unstable equilibrium**, because any *infinitesimal* push away will cause it to careen away towards another state.
"

# ╔═╡ 67f43076-1fa5-11eb-093f-75e406d054c6
md"""
##### 2.2) Radiative stability analysis

We can understand why our model has two stable equilibria and one unstable equilibrium by applying concepts from dynamical systems theory.

Recall that, with fixed CO₂ concentrations, our energy balance model differential equation can be expressed as:

$C\frac{d\,T}{d\,t} = \text{ASR}(T) - \text{OTR}(T),$

where now the Absorbed Solar Radiation (ASR) is also temperature dependent because the albedo $α(T)$ is.

In particular, by plotting the right-hand-side tendency terms as a function of the state variable $T$, we can plot a stability diagram for our system that tells us whether the planet will warm ($C \frac{d\,T}{d\,t} > 0$) or cool ($C \frac{d\,T}{d\,t} < 0$).
"""

# ╔═╡ 83819644-1fa5-11eb-0152-a9fab2f1730c
begin
	T_samples = -60.:1.:30.
	OTR = Model.outgoing_thermal_radiation.(T_samples)
	ASR = [Model.absorbed_solar_radiation.(α=α(T_sample)) for T_sample in T_samples]
	imbalance = ASR .- OTR
end;

# ╔═╡ e7b5ea16-1fa5-11eb-192e-8d50120c805b
begin
	p1_stability = plot(legend=:topleft, ylabel="energy flux [W/m²]", xlabel="temperature [°C]")
	plot!(p1_stability, T_samples, OTR, label="Outgoing Thermal Radiation", color=:blue, lw=2.)
	plot!(p1_stability, T_samples, ASR, label="Absorbed Solar Radiation", color=:orange, lw=2.)

	p2_stability = plot(ylims=(-50, 45), ylabel="energy flux [W/m²]", xlabel="temperature [°C]")
	plot!([-60., 30.], [-100, -100], fillrange=[0, 0], color=:blue, alpha=0.1, label=nothing)
	plot!([-60., 30.], [100, 100], fillrange=[0, 0], color=:red, alpha=0.1, label=nothing)
	annotate!(-58, -40, text("cooling", :left, :darkblue))
	annotate!(-58, 38, text("warming", :left, :darkred))
	plot!(p2_stability, T_samples, imbalance, label="Radiative Imbalance\n(ASR - OTR)", color=:black, lw=2.)
	plot!([7.542], [0], markershape=:diamond, markersize=6, color=:lightgrey, markerstrokecolor=:darkgrey, label=nothing)
	plot!([Model.T0], [0], markershape=:circle, markersize=6, color=:orange, markerstrokecolor=:black, label=nothing)
	plot!([-38.3], [0], markershape=:circle, markersize=6, color=:aqua, markerstrokecolor=:black, label=nothing)

	p_stability = plot(p1_stability, p2_stability, layout=(1,2), size=(680, 300))
end

# ╔═╡ 58e9b802-1e11-11eb-3479-0f7eb69b2c3a
md"
### 3) Transitioning to and from Snowball Earth

##### 3.1) Turning up the Sun

Over the entire history of the Earth, the Sun is thought to have brightened by about 40%."

# ╔═╡ 846618ea-1f7e-11eb-15e8-6176acb9278b
html"""

"""

# ╔═╡ 873befc2-1f7e-11eb-16af-fba7f1e52eeb
md"In the Neoproterozoic (~700 million years ago), the Sun was 93% as bright as it is today, such that the incoming solar radiation was $S =$ 1272 W/m², Earth's average temperature plunged to $T = -50$°C, and Earth's ice-covered surface had a high albedo (reflectivity) of $α_{i} = 0.5$.

##### 3.2) Did the increasing brightness of the Sun melt the Snowball?
If we start out in the Neoproterozoic climate and all we do is increase solar insolation to today's value of $S =$ 1368 W/m², can we warm the planet up to the pre-industrial temperature of $T=14$°C?
"

# ╔═╡ d423b466-1e2a-11eb-3bb0-77f8151cdeea
md"""*Extend upper-limit of insolation* $(@bind extend_S CheckBox(default=false))"""

# ╔═╡ e2e08386-1e22-11eb-1e02-059ce290e80f
if extend_S
	md"""
	*"Cold" branch* $(@bind show_cold CheckBox(default=false))    |
	*"Warm" branch* $(@bind show_warm CheckBox(default=false))    |
	*Unstable branch* $(@bind show_unstable CheckBox(default=false))
	"""
else
	show_cold = true;
	nothing
end

# ╔═╡ 9474f8e4-1f80-11eb-3c9e-c54e662cc29c
if extend_S
	md"""
	##### Abrupt climate transitions

	In this model, temperature variations are fairly smooth unless temperatures rise above -10°C or fall below 10°C, in which case the *ice-albedo positive feedback* kicks in and causes an **abrupt climate transition**. While this is just a simple hypothetical model, these kinds of abrupt climate transitions show up all the time in the paleoclimate record and in more realistic climate model simulations.

	![](https://www.pnas.org/content/pnas/115/52/13288/F1.large.jpg)

	This simulation teaches us that we should not take the stability of our climate for granted and that pushing our present climate outside of its historical regime of stability could trigger catastrophic abrupt climate transitions.
	"""
end

# ╔═╡ 6eca000c-1f81-11eb-068e-01d06c1beeb9
md"""
##### 3.3) If not the Sun, how did Snowball Earth melt?

The leading theory is that a slow but steady outgassing of CO₂ from volcanoes eventually caused a strong enough greenhouse gas effect to offset the cooling effect of the frozen surface's high albedo and raise temperatures above the melting point $-10$°C.
"""

# ╔═╡ c0c26cf0-1fa4-11eb-2a1b-77a9d03b9c25
html"""

"""

# ╔═╡ bda79af4-1fa4-11eb-17da-21b452f43d92
md"""
In **Homework 9**, you will extend the above model to include the effect of CO₂ and determine how much CO2 would need to be added to the snowball for it to melt.
"""

# ╔═╡ da4df78a-1e2c-11eb-1d4c-69b86e196526
md"""### 4) Towards realistic climate modelling

In this simple model, the preindustrial climate of $T=14$°C is so warm that there is no ice anywhere on the planet. Indeed, the only two valid stable climates are one with *no ice* or one with *ice everywhere*.

So how did Earth's preindustrial climate, which was relatively stable for thousands of years, have substantial ice caps at the poles?

"""

# ╔═╡ 908b50a2-1f76-11eb-0eb5-11318be6896b
md"""
**The "Aquaplanet", a simple three-dimensional ocean-world climate model**

An "Aquaplanet" is a three-dimensional global climate simulation of a hypothetical planet covered entirely by a single global Ocean. While this is of course very different from Earth, where 27% of the planet is covered in land, the "Aquaplanet" exhibits many of the same characteristics as Earth and is much more realistic than our zero-dimensional climate model above.

The video below shows that the Aquaplanet simulation exhibits a third equilibrium state, with a *mostly-liquid ocean but ice caps at the poles*, in addition to the two we found in our zero-dimensional model.

In **Homework 10**, you will build a simple two-dimensional version of the aqua-planet and explore its stability.
"""

# ╔═╡ 507342e0-1f76-11eb-098a-0973155652e2
html"""

"""

# ╔═╡ d6a7e652-1f7d-11eb-2a5f-13ff13b9850c
md"""
## Appendix
"""

# ╔═╡ d2f4d882-1e2c-11eb-206b-b9ef5031a1bb
md"#### Pluto Magic"

# ╔═╡ 1dc709e2-1de8-11eb-28da-fd7469ec50c2
function ingredients(path::String)
	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol(basename(path))
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
end

# ╔═╡ a212680c-1fa2-11eb-1fb9-f910081c562a
as_svg(x) = PlutoUI.Show(MIME"image/svg+xml"(), repr(MIME"image/svg+xml"(), x))

# ╔═╡ 8e28ef88-1f88-11eb-2d5d-bff2f59ea998
md"""
#### Bifurcation diagram helper functions
"""

# ╔═╡ 61fc91ec-1df8-11eb-13c1-098c113b46ec
function restart_ebm!(ebm)
	ebm.T = [ebm.T[end]]
	ebm.t = [ebm.t[1]]
end

# ╔═╡ 9b39df12-1df9-11eb-2eb0-f138980be597
function plot_trajectory!(p, x, y; lw=8)
	n = size(x,1)
	plot!(x, y, color="black", alpha=collect(1/n:1/n:1.), linewidth=collect(0.:lw/n:lw), label=nothing)
	plot!((x[end], y[end]), color="black", marker=:c, markersize=lw/2*1.2, label=nothing, markerstrokecolor=nothing, markerstrokewidth=0.)
	return p
end;

# ╔═╡ b5849000-0d68-11eb-3beb-c575e8d0ce8e
begin
	Smin = 1200
	Smax = 1800
	Smax_limited = 1650
	Svec = Smin:1.:Smax
	Svec = vcat(Svec, reverse(Svec))
	Tvec = zeros(size(Svec))

	local T_restart = -100.
	for (i, S) = enumerate(Svec)
		ebm = Model.EBM(T_restart, 0., 5., Model.CO2_const);
		ebm.S = S
		Model.run!(ebm, 400.)
		T_restart = ebm.T[end]
		Tvec[i] = deepcopy(T_restart)
	end

	md"**Data structures for storing warm & cool branch climates**"
end

# ╔═╡ 364e1d2a-1f78-11eb-177b-9dcec4b6da38
md"""
For insolations $S$ between $(Smin) W/m² and $(Smax_limited) W/m², temperatures always remain below -10°C and the planet remains completely frozen. What if we extend the upper limit on insolation so that the Sun becomes bright enough to start melting ice?
"""

# ╔═╡ b4357da4-1f9e-11eb-2925-55cd6f37be22
begin
	T_unstable_branch(S, A, B, αi, α0) = (
		(A/B-S/(4B)*(1-0.5(αi+α0))) /
		(1+S*(αi-α0)/(80B))
	)
	S_unstable = Smin:2:Smax
	T_unstable = T_unstable_branch.(S_unstable, Model.A, Model.B, 0.5, 0.3)
	T_unstable[.~(0.3 .< α.(T_unstable) .< 0.5)] .= NaN
	md"**Unstable branch solution**"
end

# ╔═╡ 90ae18dc-0db8-11eb-0d73-c3b7efaef9b0
begin
	Sneo = Model.S*0.93
	Tneo = -48.
	md"**Initial conditions**"
end

# ╔═╡ 5ca7ac14-1df9-11eb-13f3-e5c86333ef83
begin
	if extend_S
		solarSlider = @bind S Slider(Smin:2.:Smax, default=Sneo);
		md""" $(Smin) W/m² $(solarSlider) $(Smax) W/m²"""
	else
		solarSlider = @bind S Slider(Smin:2.:Smax_limited, default=Sneo);
		md""" $(Smin) W/m² $(solarSlider) $(Smax_limited) W/m²"""
	end
end

# ╔═╡ 0bbcdf5a-0dba-11eb-3e81-2b075d4f67ea
begin
	md"""
	*Move the slider below to change the brightness of the Sun (solar insolation):* S = $(S) [W/m²]
	"""
end

# ╔═╡ 7765f834-0db0-11eb-2c46-ef65f4a1fd1d
begin
	ebm = Model.EBM(Tneo, 0., 5., Model.CO2_const)
	ebm.S = Sneo

	ntraj = 10
	Ttraj = repeat([NaN], ntraj)
	Straj = repeat([NaN], ntraj)

	md"**Data structures for storing trajectories of recent climates**"
end

# ╔═╡ 0f222836-0d6c-11eb-2ee8-45545da73cfd
begin
	S
	warming_mask = (1:size(Svec)[1]) .< size(Svec)./2
	p = plot(xlims=(Smin, Smax_limited), ylims=(-55, 75), title="Earth's solar insolation bifurcation diagram")
	plot!([Model.S, Model.S], [-55, 75], color=:yellow, alpha=0.3, lw=8, label="Pre-industrial / present insolation")
	if extend_S
		plot!(p, xlims=(Smin, Smax))
		if show_cold
			plot!(Svec[warming_mask], Tvec[warming_mask], color=:blue,lw=3., alpha=0.5, label="cool branch")
		end
		if show_warm
			plot!(Svec[.!warming_mask], Tvec[.!warming_mask], color="red", lw=3., alpha=0.5, label="warm branch")
		end
		if show_unstable
			plot!(S_unstable, T_unstable, color=:darkgray, lw=3., alpha=0.4, ls=:dash, label="unstable branch")
		end
	end
	plot!(legend=:topleft)
	plot!(xlabel="solar insolation S [W/m²]", ylabel="Global temperature T [°C]")
	plot!([Model.S], [Model.T0], markershape=:circle, label="Our preindustrial climate", color=:orange, markersize=8)
	plot!([Model.S], [-38.3], markershape=:circle, label="Alternate preindustrial climate", color=:aqua, markersize=8)
	plot!([Sneo], [Tneo], markershape=:circle, label="neoproterozoic (700 Mya)", color=:lightblue, markersize=8)
	plot_trajectory!(p, reverse(Straj), reverse(Ttraj), lw=9)

	plot!([Smin, Smax], [-60, -60], fillrange=[-10., -10.], fillalpha=0.3, c=:lightblue, label=nothing)
	annotate!(Smin+12, -19, text("completely\nfrozen", 10, :darkblue, :left))

	plot!([Smin, Smax], [10, 10], fillrange=[80., 80.], fillalpha=0.09, c=:red, lw=0., label=nothing)
	annotate!(Smin+12, 15, text("no ice", 10, :darkred, :left))
end |> as_svg

# ╔═╡ 477732c4-0daf-11eb-1422-cf0f55cd835b
begin
	S
	restart_ebm!(ebm)
	ebm.S = S
	Model.run!(ebm, 500)

	insert!(Straj, 1, copy(S))
	pop!(Straj)

	insert!(Ttraj, 1, copy(ebm.T[end]))
	pop!(Ttraj)
end;

# ╔═╡ f4b24d20-251a-11eb-04cf-fd69abdcfd54
function add_labels!(p)
	plot!(p, xlabel="year", ylabel="temperature [°C]", legend=:bottomright, xlims=(-5,205), ylims=(-60, 30.))

	plot!(p, [-5, 200], [-60, -60], fillrange=[-10., -10.], fillalpha=0.3, c=:lightblue, label=nothing)
	annotate!(120, -20, text("completely frozen", 10, :darkblue))

	plot!(p, [-5, 200], [10, 10], fillrange=[30., 30.], fillalpha=0.09, c=:red, lw=0., label=nothing)
	annotate!(p, 120, 25, text("no ice", 10, :darkred))
end

# ╔═╡ af197a52-2503-11eb-1da1-d36bab6ef2d3
begin
	p_interact = plot(ebm_interact.t, ebm_interact.T, label=nothing, lw=3)
	plot!([0.], [T0_interact], label=nothing, markersize=4, markershape=:circle)

	add_labels!(p_interact)
end |> as_svg

# ╔═╡ Cell order:
# ╟─05031b60-1df4-11eb-2b61-956e526b3d4a
# ╟─a05bd3ce-1e06-11eb-1af2-65886cab38ef
# ╟─a6b6ea46-2519-11eb-1134-8994582c5546
# ╟─9d3a6312-2519-11eb-3550-89ea034bf119
# ╟─c50b6b24-2506-11eb-00e1-0f183526ed4e
# ╟─4f5f3038-1e06-11eb-16a2-b11035701fb8
# ╟─68364804-2517-11eb-00de-d1655dd7754b
# ╟─016c1074-1df4-11eb-2da8-578e25d9456b
# ╟─9c118f9a-1df0-11eb-22dd-b14428994076
# ╟─38346e6a-0d98-11eb-280b-f79787a3c788
# ╟─f7761e40-1e23-11eb-2741-cfebfaf434ec
# ╟─816f1d96-2508-11eb-0873-c3b564a31dea
# ╠═262fc3c6-1df2-11eb-332d-c1c9561b3710
# ╟─a8dcc0fc-1df8-11eb-21fd-1fdebe5dabfc
# ╟─96ed2f9a-1e29-11eb-09f4-23df52152b2f
# ╠═872c8f2a-1df1-11eb-3cfc-3dd568926442
# ╟─13f42334-1e27-11eb-11a0-f51af4574a6b
# ╟─03292b48-2503-11eb-1514-5d1923d1d9a2
# ╟─7a1b3138-2503-11eb-1165-c399836b66a7
# ╠═1be1cce0-251b-11eb-35f7-e15741b0a712
# ╟─af197a52-2503-11eb-1da1-d36bab6ef2d3
# ╟─a954ce70-2510-11eb-0d8c-c7b3f0fb3a06
# ╟─94852946-1e25-11eb-3425-210be17c23cd
# ╟─c81b4956-2510-11eb-388c-87e93923e45a
# ╟─67f43076-1fa5-11eb-093f-75e406d054c6
# ╠═83819644-1fa5-11eb-0152-a9fab2f1730c
# ╟─e7b5ea16-1fa5-11eb-192e-8d50120c805b
# ╟─58e9b802-1e11-11eb-3479-0f7eb69b2c3a
# ╟─846618ea-1f7e-11eb-15e8-6176acb9278b
# ╟─873befc2-1f7e-11eb-16af-fba7f1e52eeb
# ╟─0bbcdf5a-0dba-11eb-3e81-2b075d4f67ea
# ╟─5ca7ac14-1df9-11eb-13f3-e5c86333ef83
# ╟─0f222836-0d6c-11eb-2ee8-45545da73cfd
# ╟─364e1d2a-1f78-11eb-177b-9dcec4b6da38
# ╟─d423b466-1e2a-11eb-3bb0-77f8151cdeea
# ╟─e2e08386-1e22-11eb-1e02-059ce290e80f
# ╟─9474f8e4-1f80-11eb-3c9e-c54e662cc29c
# ╟─6eca000c-1f81-11eb-068e-01d06c1beeb9
# ╟─c0c26cf0-1fa4-11eb-2a1b-77a9d03b9c25
# ╟─bda79af4-1fa4-11eb-17da-21b452f43d92
# ╟─da4df78a-1e2c-11eb-1d4c-69b86e196526
# ╟─908b50a2-1f76-11eb-0eb5-11318be6896b
# ╟─507342e0-1f76-11eb-098a-0973155652e2
# ╟─d6a7e652-1f7d-11eb-2a5f-13ff13b9850c
# ╟─d2f4d882-1e2c-11eb-206b-b9ef5031a1bb
# ╟─1dc709e2-1de8-11eb-28da-fd7469ec50c2
# ╟─a212680c-1fa2-11eb-1fb9-f910081c562a
# ╟─8e28ef88-1f88-11eb-2d5d-bff2f59ea998
# ╠═61fc91ec-1df8-11eb-13c1-098c113b46ec
# ╠═7765f834-0db0-11eb-2c46-ef65f4a1fd1d
# ╠═477732c4-0daf-11eb-1422-cf0f55cd835b
# ╠═9b39df12-1df9-11eb-2eb0-f138980be597
# ╠═b5849000-0d68-11eb-3beb-c575e8d0ce8e
# ╠═b4357da4-1f9e-11eb-2925-55cd6f37be22
# ╠═90ae18dc-0db8-11eb-0d73-c3b7efaef9b0
# ╟─f4b24d20-251a-11eb-04cf-fd69abdcfd54
# ╠═e9ad66b0-0d6b-11eb-26c0-27413c19dd32

max_stars_repo_name: mathijsvdv/ComputationalThinking

max_issues_repo_name: mathijsvdv/ComputationalThinking

max_forks_repo_name: mathijsvdv/ComputationalThinking


text:
function get_right_endpoint(
    theta_init::Vector{Float64}, # initial point of parameters
    theta_num::Int, # number of parameter to scan
    loss_func::Function, # lambda(theta) - labmbda_min - delta_lambda
    method::Val{:LIN_EXTRAPOL}; # function works only for method LIN_INTER

    theta_bounds::Vector{Tuple{Float64,Float64}} = fill(
        (-Inf, Inf), length(theta_init)
        ),
    scan_bound::Float64 = 9.0,
    scan_tol::Float64 = 1e-3,
    loss_tol::Float64 = 0., # 1e-3,
    # method args
    scan_hini = 1.,
    scan_hmax = Inf,
    # local alg args
    local_alg::Symbol = :LN_NELDERMEAD,
    max_iter::Int = 10^5,
    kwargs... # options for local fitter
    )
    # dim of the theta vector
    n_theta = length(theta_init)

    # checking arguments
    # methods which are not supported
    if local_alg in [:LN_BOBYQA, :LN_SBPLX, :LN_NEWUOA]
        @warn "Using local_alg = :"*String(local_alg)*" may result in wrong output."
    end
    # when using :LN_NELDERMEAD initial parameters should not be zero
    if local_alg == :LN_NELDERMEAD
        zeroParameter = [ isapprox(theta_init[i], 0., atol=1e-2) for i in 1:n_theta]
        if any(zeroParameter)
            @warn "Close-to-zero parameters found when using :LN_NELDERMEAD."
            show(findall(zeroParameter))
        end
    end

    # to count loss function calls inside this function, accumulation
    accum_counter::Int = 0
    # empty container
    pps = ProfilePoint[]

    prof = _profile(
        theta_init,
        theta_num,
        loss_func;
        theta_bounds = theta_bounds,
        local_alg = local_alg,
        ftol_abs = loss_tol
    )

    # first iteration
    x_2 = theta_init[theta_num]
    theta_init_2 = theta_init
    iteration_count = 0

    # other iterations
    while true
        # preparation
        global x_1 # not initialized for the first iteration
        global point_1 # not initialized for the first iteration
        iteration_count += 1 # to understand if this is a first iteration

        # get profile point
        point_2 = prof(
            x_2;
            theta_init_i = theta_init_2, # hypothetically this makes optimization more effective
            maxeval = max_iter - accum_counter # how many calls left
            )
        push!(pps, point_2)
        accum_counter += point_2.counter # update counter

        if point_2.ret == :MAXEVAL_REACHED
            return (nothing, pps, :MAX_ITER_STOP)
        elseif point_2.ret == :FORCED_STOP
            return (nothing, pps, :LOSS_ERROR_STOP)
        elseif x_2 >= scan_bound && point_2.loss < 0. # successfull result
            return (nothing, pps, :SCAN_BOUND_REACHED)
        # no checking for the first iteration
        elseif iteration_count>1 && isapprox((x_2 - x_1) * point_2.loss / (point_2.loss - point_1.loss), 0., atol = scan_tol) # successfull result
            return (x_2, pps, :BORDER_FOUND_BY_SCAN_TOL)
        elseif isapprox(point_2.loss, 0., atol = loss_tol) # successfull result
            return (x_2, pps, :BORDER_FOUND_BY_LOSS_TOL)
        end

        # next step
        if iteration_count==1
            x_3_extrapol = x_2 + scan_hini
        else
            if (point_2.loss - point_1.loss) / (x_2 - x_1) <= 0.
                x_3_extrapol = x_2 + scan_hini
            else
                x_3_extrapol = x_2 + (x_2 - x_1) * (0. - point_2.loss) / (point_2.loss - point_1.loss)
            end
        end
        x_3 = minimum([x_2+scan_hmax, x_3_extrapol, theta_bounds[theta_num][2]])
        # preparation for the next iteration
        (point_1, x_1, x_2, theta_init_2) = (point_2, x_2, x_3, point_2.params)
    end
end