text:
using DataFrames, CSV, Statistics, StatsBase

df = CSV.read("./preprocessing/input/land-sustainable.csv")
vars = CSV.read("./preprocessing/input/land-sustainable-attributes.csv")

# remove last row
df = df[.!ismissing.(df.Scenario),:]
# remove region col
select!(df, Not(:Region))
df = stack(df, 5:size(df)[2])
rename!(df, [:year, :value, :model, :scenario, :variable, :unit])
df[occursin.("N2O", df.variable), :value] /= 1000
df[occursin.("CO2", df.variable), :value] /= 1000
df[occursin.("Sequestration", df.variable), :value] /= 1000
df[occursin.("N2O", df.variable), :unit] = "Mt N₂O/yr"
df[occursin.("CO2", df.variable), :unit] = "Gt CO₂/yr"
df[occursin.("Sequestration", df.variable), :unit] = "Gt CO₂/yr"
df[occursin.("CH4", df.variable), :unit] = "Mt CH₄/yr"
df = unstack(df, :year, :value)

df[occursin.("Def", df.scenario), :scenario] = "default"
df[occursin.("Sust", df.scenario), :scenario] = "sustainable"

# sort!(df, [order(:variable, rev=true), order(:scenario)])
df = join(vars,df,on = :variable, kind= :outer)
CSV.write("./src/assets/data/land-sustainable.csv", df)

max_stars_repo_name: SensesProject/land-use

max_issues_repo_name: SensesProject/land-use

max_forks_repo_name: SensesProject/land-use


text:
using Ska
using Test

@testset "Ska.jl" begin
    # Write your tests here.
end

ext: jl


text:
function second_moment(q::Normal)
    return abs2(mean(q)) + var(q)
end

function second_moment(q::Normal, y::Real)
    return abs2(mean(q) - y) + var(q)
end

# This is not exactly an approximation but corresponds anyway
# to the expectation of the function σ(f)
function approx_expected_logistic(μ::Real, c::Real)
    lower, upper = LogExpFunctions._logistic_bounds(μ)
    return μ < lower ? zero(μ) : (μ > upper ? one(μ) : exp(μ / 2) * sech(c / 2) / 2)
end

# Same thing as above for softmax
function approx_expected_logisticsoftmax(
    μ::AbstractVector, c::AbstractVector, θ::AbstractVector
)
    σs = θ[1:(end - 1)] .* approx_expected_logistic.(μ, c)
    return σs / (θ[end] * logistic(0) + sum(σs))
end

transpose_nested(x::ArrayOfSimilarArrays) = nestedview(flatview(x)')


text:
module GLU
#Autogenerated, hopefully stays that way!

import GetC.@getCFun

typealias GLenum Uint32
typealias GLboolean Uint8
typealias GLbitfield Uint32
typealias GLvoid Void
typealias GLbyte Uint8
typealias GLshort Int16
typealias GLint Int32
typealias GLubyte Uint8
typealias GLushort Uint16
typealias GLuint Uint32
typealias GLsizei Int32
typealias GLfloat Float32
typealias GLclampf Float32
typealias GLdouble Float64
typealias GLclampd Float64
typealias ptrdiff_t Int64
typealias size_t Int64
typealias wchar_t Int32
typealias GLchar Uint8
typealias GLintptr ptrdiff_t
typealias GLsizeiptr ptrdiff_t
typealias GLintptrARB ptrdiff_t
typealias GLsizeiptrARB ptrdiff_t
typealias GLcharARB Uint8
typealias GLhandleARB Uint32
typealias GLhalfARB Uint16
typealias GLhalfNV Uint16
typealias int8_t Uint8
typealias int16_t Int16
typealias int32_t Int32
typealias int64_t Int64
typealias uint8_t Uint8
typealias uint16_t Uint16
typealias uint32_t Uint32
typealias uint64_t Uint64
typealias int_least8_t Uint8
typealias int_least16_t Int16
typealias int_least32_t Int32
typealias int_least64_t Int64
typealias uint_least8_t Uint8
typealias uint_least16_t Uint16
typealias uint_least32_t Uint32
typealias uint_least64_t Uint64
typealias int_fast8_t Uint8
typealias int_fast16_t Int64
typealias int_fast32_t Int64
typealias int_fast64_t Int64
typealias uint_fast8_t Uint8
typealias uint_fast16_t Uint64
typealias uint_fast32_t Uint64
typealias uint_fast64_t Uint64
typealias intptr_t Int64
typealias uintptr_t Uint64
typealias intmax_t Int64
typealias uintmax_t Uint64
typealias __gwchar_t Int32
export GLenum
export GLboolean
export GLbitfield
export GLvoid
export GLbyte
export GLshort
export GLint
export GLubyte
export GLushort
export GLuint
export GLsizei
export GLfloat
export GLclampf
export GLdouble
export GLclampd
export ptrdiff_t
export size_t
export wchar_t
export GLchar
export GLintptr
export GLsizeiptr
export GLintptrARB
export GLsizeiptrARB
export GLcharARB
export GLhandleARBt32
export GLhalfARB
export GLhalfNV
export int8_t
export int16_t
export int32_t
export int64_t
export uint8_t
export uint16_t
export uint32_t
export uint64_t
export int_least8_t
export int_least16_t
export int_least32_t
export int_least64_t
export uint_least8_t
export uint_least16_t
export uint_least32_t
export uint_least64_t
export int_fast8_t
export int_fast16_t
export int_fast32_t
export int_fast64_t
export uint_fast8_t
export uint_fast16_t
export uint_fast32_t
export uint_fast64_t
export intptr_t
export uintptr_t
export intmax_t
export uintmax_t
export __gwchar_t
@getCFun "libGLU" imaxabs imaxabs(__n::intmax_t)::intmax_t
export imaxabs
@getCFun "libGLU" imaxdiv imaxdiv(__numer::intmax_t,__denom::intmax_t)::imaxdiv_t
export imaxdiv
@getCFun "libGLU" strtoimax strtoimax(__nptr::Ptr,__endptr::Ptr,__base::Int32)::intmax_t
export strtoimax
@getCFun "libGLU" strtoumax strtoumax(__nptr::Ptr,__endptr::Ptr,__base::Int32)::uintmax_t
export strtoumax
@getCFun "libGLU" wcstoimax wcstoimax(__nptr::Ptr,__endptr::Ptr,__base::Int32)::intmax_t
export wcstoimax
@getCFun "libGLU" wcstoumax wcstoumax(__nptr::Ptr,__endptr::Ptr,__base::Int32)::uintmax_t
export wcstoumax
typealias GLint64EXT int64_t
typealias GLuint64EXT uint64_t
typealias GLint64 int64_t
typealias GLuint64 uint64_t
typealias GLvdpauSurfaceNV GLintptr
typealias GLeglImageOES Ptr
export GLint64EXT
export GLuint64EXT
export GLint64
export GLuint64
export GLvdpauSurfaceNV
export GLeglImageOES
@getCFun "libGLU" gluPerspective gluPerspective(fovy::GLdouble,aspect::GLdouble,znear::GLdouble,zfar::GLdouble)::Void
export gluPerspective
@getCFun "libGLU" gluBuild2DMipmaps gluBuild2DMipmaps(target::GLenum,internalFormat::GLint,width::GLsizei,height::GLsizei,format::GLenum,thetype::GLenum,data::Ptr{Uint8})::GLint
export gluBuild2DMipmaps
@getCFun "libGLU" gluNewQuadric gluNewQuadric()::Ptr{Void}
export gluNewQuadric
@getCFun "libGLU" gluQuadricNormals gluQuadricNormals(quad::Ptr{Void},normal::GLenum)::Void
export gluQuadricNormals
@getCFun "libGLU" gluQuadricTexture gluQuadricTexture(quad::Ptr{Void},texture::GLboolean)::Void
export gluQuadricTexture
@getCFun "libGLU" gluCylinder gluCylinder(quad::Ptr{Void},base::GLdouble,top::GLdouble,height::GLdouble,slices::GLint,stacks::GLint)::Void
export gluCylinder
@getCFun "libGLU" gluDisk gluDisk(quad::Ptr{Void},inner::GLdouble,outer::GLdouble,slices::GLint,stacks::GLint)::Void
export gluDisk
@getCFun "libGLU" gluSphere gluSphere(quad::Ptr{Void},radius::GLdouble,slices::GLint,stacks::GLint)::Void
export gluSphere
@getCFun "libGLU" gluPartialDisk gluPartialDisk(quad::Ptr{Void},inner::GLdouble,outer::GLdouble,slices::GLint,loops::GLint,start::GLdouble,sweep::GLdouble)::Void
export gluPartialDisk
const GLU_EXT_object_space_tess       = 1
const GLU_EXT_nurbs_tessellator       = 1
const GLU_FALSE                       = 0
const GLU_TRUE                        = 1
const GLU_VERSION_1_1                 = 1
const GLU_VERSION_1_2                 = 1
const GLU_VERSION_1_3                 = 1
const GLU_VERSION                     = 100800
const GLU_EXTENSIONS                  = 100801
const GLU_INVALID_ENUM                = 100900
const GLU_INVALID_VALUE               = 100901
const GLU_OUT_OF_MEMORY               = 100902
const GLU_INCOMPATIBLE_GL_VERSION     = 100903
const GLU_INVALID_OPERATION           = 100904
const GLU_OUTLINE_POLYGON             = 100240
const GLU_OUTLINE_PATCH               = 100241
const GLU_NURBS_ERROR                 = 100103
const GLU_ERROR                       = 100103
const GLU_NURBS_BEGIN                 = 100164
const GLU_NURBS_BEGIN_EXT             = 100164
const GLU_NURBS_VERTEX                = 100165
const GLU_NURBS_VERTEX_EXT            = 100165
const GLU_NURBS_NORMAL                = 100166
const GLU_NURBS_NORMAL_EXT            = 100166
const GLU_NURBS_COLOR                 = 100167
const GLU_NURBS_COLOR_EXT             = 100167
const GLU_NURBS_TEXTURE_COORD         = 100168
const GLU_NURBS_TEX_COORD_EXT         = 100168
const GLU_NURBS_END                   = 100169
const GLU_NURBS_END_EXT               = 100169
const GLU_NURBS_BEGIN_DATA            = 100170
const GLU_NURBS_BEGIN_DATA_EXT        = 100170
const GLU_NURBS_VERTEX_DATA           = 100171
const GLU_NURBS_VERTEX_DATA_EXT       = 100171
const GLU_NURBS_NORMAL_DATA           = 100172
const GLU_NURBS_NORMAL_DATA_EXT       = 100172
const GLU_NURBS_COLOR_DATA            = 100173
const GLU_NURBS_COLOR_DATA_EXT        = 100173
const GLU_NURBS_TEXTURE_COORD_DATA    = 100174
const GLU_NURBS_TEX_COORD_DATA_EXT    = 100174
const GLU_NURBS_END_DATA              = 100175
const GLU_NURBS_END_DATA_EXT          = 100175
const GLU_NURBS_ERROR1                = 100251
const GLU_NURBS_ERROR2                = 100252
const GLU_NURBS_ERROR3                = 100253
const GLU_NURBS_ERROR4                = 100254
const GLU_NURBS_ERROR5                = 100255
const GLU_NURBS_ERROR6                = 100256
const GLU_NURBS_ERROR7                = 100257
const GLU_NURBS_ERROR8                = 100258
const GLU_NURBS_ERROR9                = 100259
const GLU_NURBS_ERROR10               = 100260
const GLU_NURBS_ERROR11               = 100261
const GLU_NURBS_ERROR12               = 100262
const GLU_NURBS_ERROR13               = 100263
const GLU_NURBS_ERROR14               = 100264
const GLU_NURBS_ERROR15               = 100265
const GLU_NURBS_ERROR16               = 100266
const GLU_NURBS_ERROR17               = 100267
const GLU_NURBS_ERROR18               = 100268
const GLU_NURBS_ERROR19               = 100269
const GLU_NURBS_ERROR20               = 100270
const GLU_NURBS_ERROR21               = 100271
const GLU_NURBS_ERROR22               = 100272
const GLU_NURBS_ERROR23               = 100273
const GLU_NURBS_ERROR24               = 100274
const GLU_NURBS_ERROR25               = 100275
const GLU_NURBS_ERROR26               = 100276
const GLU_NURBS_ERROR27               = 100277
const GLU_NURBS_ERROR28               = 100278
const GLU_NURBS_ERROR29               = 100279
const GLU_NURBS_ERROR30               = 100280
const GLU_NURBS_ERROR31               = 100281
const GLU_NURBS_ERROR32               = 100282
const GLU_NURBS_ERROR33               = 100283
const GLU_NURBS_ERROR34               = 100284
const GLU_NURBS_ERROR35               = 100285
const GLU_NURBS_ERROR36               = 100286
const GLU_NURBS_ERROR37               = 100287
const GLU_AUTO_LOAD_MATRIX            = 100200
const GLU_CULLING                     = 100201
const GLU_SAMPLING_TOLERANCE          = 100203
const GLU_DISPLAY_MODE                = 100204
const GLU_PARAMETRIC_TOLERANCE        = 100202
const GLU_SAMPLING_METHOD             = 100205
const GLU_U_STEP                      = 100206
const GLU_V_STEP                      = 100207
const GLU_NURBS_MODE                  = 100160
const GLU_NURBS_MODE_EXT              = 100160
const GLU_NURBS_TESSELLATOR           = 100161
const GLU_NURBS_TESSELLATOR_EXT       = 100161
const GLU_NURBS_RENDERER              = 100162
const GLU_NURBS_RENDERER_EXT          = 100162
const GLU_OBJECT_PARAMETRIC_ERROR     = 100208
const GLU_OBJECT_PARAMETRIC_ERROR_EXT = 100208
const GLU_OBJECT_PATH_LENGTH          = 100209
const GLU_OBJECT_PATH_LENGTH_EXT      = 100209
const GLU_PATH_LENGTH                 = 100215
const GLU_PARAMETRIC_ERROR            = 100216
const GLU_DOMAIN_DISTANCE             = 100217
const GLU_MAP1_TRIM_2                 = 100210
const GLU_MAP1_TRIM_3                 = 100211
const GLU_POINT                       = 100010
const GLU_LINE                        = 100011
const GLU_FILL                        = 100012
const GLU_SILHOUETTE                  = 100013
const GLU_SMOOTH                      = 100000
const GLU_FLAT                        = 100001
const GLU_NONE                        = 100002
const GLU_OUTSIDE                     = 100020
const GLU_INSIDE                      = 100021
const GLU_TESS_BEGIN                  = 100100
const GLU_BEGIN                       = 100100
const GLU_TESS_VERTEX                 = 100101
const GLU_VERTEX                      = 100101
const GLU_TESS_END                    = 100102
const GLU_END                         = 100102
const GLU_TESS_ERROR                  = 100103
const GLU_TESS_EDGE_FLAG              = 100104
const GLU_EDGE_FLAG                   = 100104
const GLU_TESS_COMBINE                = 100105
const GLU_TESS_BEGIN_DATA             = 100106
const GLU_TESS_VERTEX_DATA            = 100107
const GLU_TESS_END_DATA               = 100108
const GLU_TESS_ERROR_DATA             = 100109
const GLU_TESS_EDGE_FLAG_DATA         = 100110
const GLU_TESS_COMBINE_DATA           = 100111
const GLU_CW                          = 100120
const GLU_CCW                         = 100121
const GLU_INTERIOR                    = 100122
const GLU_EXTERIOR                    = 100123
const GLU_UNKNOWN                     = 100124
const GLU_TESS_WINDING_RULE           = 100140
const GLU_TESS_BOUNDARY_ONLY          = 100141
const GLU_TESS_TOLERANCE              = 100142
const GLU_TESS_ERROR1                 = 100151
const GLU_TESS_ERROR2                 = 100152
const GLU_TESS_ERROR3                 = 100153
const GLU_TESS_ERROR4                 = 100154
const GLU_TESS_ERROR5                 = 100155
const GLU_TESS_ERROR6                 = 100156
const GLU_TESS_ERROR7                 = 100157
const GLU_TESS_ERROR8                 = 100158
const GLU_TESS_MISSING_BEGIN_POLYGON  = 100151
const GLU_TESS_MISSING_BEGIN_CONTOUR  = 100152
const GLU_TESS_MISSING_END_POLYGON    = 100153
const GLU_TESS_MISSING_END_CONTOUR    = 100154
const GLU_TESS_COORD_TOO_LARGE        = 100155
const GLU_TESS_NEED_COMBINE_CALLBACK  = 100156
const GLU_TESS_WINDING_ODD            = 100130
const GLU_TESS_WINDING_NONZERO        = 100131
const GLU_TESS_WINDING_POSITIVE       = 100132
const GLU_TESS_WINDING_NEGATIVE       = 100133
const GLU_TESS_WINDING_ABS_GEQ_TWO    = 100134
const GLU_TESS_MAX_COORD              = 1.0e150
export GLU_EXT_object_space_tess
export GLU_EXT_nurbs_tessellator
export GLU_FALSE
export GLU_TRUE
export GLU_VERSION_1_1
export GLU_VERSION_1_2
export GLU_VERSION_1_3
export GLU_VERSION
export GLU_EXTENSIONS
export GLU_INVALID_ENUM
export GLU_INVALID_VALUE
export GLU_OUT_OF_MEMORY
export GLU_INCOMPATIBLE_GL_VERSION
export GLU_INVALID_OPERATION
export GLU_OUTLINE_POLYGON
export GLU_OUTLINE_PATCH
export GLU_NURBS_ERROR
export GLU_ERROR
export GLU_NURBS_BEGIN
export GLU_NURBS_BEGIN_EXT
export GLU_NURBS_VERTEX
export GLU_NURBS_VERTEX_EXT
export GLU_NURBS_NORMAL
export GLU_NURBS_NORMAL_EXT
export GLU_NURBS_COLOR
export GLU_NURBS_COLOR_EXT
export GLU_NURBS_TEXTURE_COORD
export GLU_NURBS_TEX_COORD_EXT
export GLU_NURBS_END
export GLU_NURBS_END_EXT
export GLU_NURBS_BEGIN_DATA
export GLU_NURBS_BEGIN_DATA_EXT
export GLU_NURBS_VERTEX_DATA
export GLU_NURBS_VERTEX_DATA_EXT
export GLU_NURBS_NORMAL_DATA
export GLU_NURBS_NORMAL_DATA_EXT
export GLU_NURBS_COLOR_DATA
export GLU_NURBS_COLOR_DATA_EXT
export GLU_NURBS_TEXTURE_COORD_DATA
export GLU_NURBS_TEX_COORD_DATA_EXT
export GLU_NURBS_END_DATA
export GLU_NURBS_END_DATA_EXT
export GLU_NURBS_ERROR1
export GLU_NURBS_ERROR2
export GLU_NURBS_ERROR3
export GLU_NURBS_ERROR4
export GLU_NURBS_ERROR5
export GLU_NURBS_ERROR6
export GLU_NURBS_ERROR7
export GLU_NURBS_ERROR8
export GLU_NURBS_ERROR9
export GLU_NURBS_ERROR10
export GLU_NURBS_ERROR11
export GLU_NURBS_ERROR12
export GLU_NURBS_ERROR13
export GLU_NURBS_ERROR14
export GLU_NURBS_ERROR15
export GLU_NURBS_ERROR16
export GLU_NURBS_ERROR17
export GLU_NURBS_ERROR18
export GLU_NURBS_ERROR19
export GLU_NURBS_ERROR20
export GLU_NURBS_ERROR21
export GLU_NURBS_ERROR22
export GLU_NURBS_ERROR23
export GLU_NURBS_ERROR24
export GLU_NURBS_ERROR25
export GLU_NURBS_ERROR26
export GLU_NURBS_ERROR27
export GLU_NURBS_ERROR28
export GLU_NURBS_ERROR29
export GLU_NURBS_ERROR30
export GLU_NURBS_ERROR31
export GLU_NURBS_ERROR32
export GLU_NURBS_ERROR33
export GLU_NURBS_ERROR34
export GLU_NURBS_ERROR35
export GLU_NURBS_ERROR36
export GLU_NURBS_ERROR37
export GLU_AUTO_LOAD_MATRIX
export GLU_CULLING
export GLU_SAMPLING_TOLERANCE
export GLU_DISPLAY_MODE
export GLU_PARAMETRIC_TOLERANCE
export GLU_SAMPLING_METHOD
export GLU_U_STEP
export GLU_V_STEP
export GLU_NURBS_MODE
export GLU_NURBS_MODE_EXT
export GLU_NURBS_TESSELLATOR
export GLU_NURBS_TESSELLATOR_EXT
export GLU_NURBS_RENDERER
export GLU_NURBS_RENDERER_EXT
export GLU_OBJECT_PARAMETRIC_ERROR
export GLU_OBJECT_PARAMETRIC_ERROR_EXT
export GLU_OBJECT_PATH_LENGTH
export GLU_OBJECT_PATH_LENGTH_EXT
export GLU_PATH_LENGTH
export GLU_PARAMETRIC_ERROR
export GLU_DOMAIN_DISTANCE
export GLU_MAP1_TRIM_2
export GLU_MAP1_TRIM_3
export GLU_POINT
export GLU_LINE
export GLU_FILL
export GLU_SILHOUETTE
export GLU_SMOOTH
export GLU_FLAT
export GLU_NONE
export GLU_OUTSIDE
export GLU_INSIDE
export GLU_TESS_BEGIN
export GLU_BEGIN
export GLU_TESS_VERTEX
export GLU_VERTEX
export GLU_TESS_END
export GLU_END
export GLU_TESS_ERROR
export GLU_TESS_EDGE_FLAG
export GLU_EDGE_FLAG
export GLU_TESS_COMBINE
export GLU_TESS_BEGIN_DATA
export GLU_TESS_VERTEX_DATA
export GLU_TESS_END_DATA
export GLU_TESS_ERROR_DATA
export GLU_TESS_EDGE_FLAG_DATA
export GLU_TESS_COMBINE_DATA
export GLU_CW
export GLU_CCW
export GLU_INTERIOR
export GLU_EXTERIOR
export GLU_UNKNOWN
export GLU_TESS_WINDING_RULE
export GLU_TESS_BOUNDARY_ONLY
export GLU_TESS_TOLERANCE
export GLU_TESS_ERROR1
export GLU_TESS_ERROR2
export GLU_TESS_ERROR3
export GLU_TESS_ERROR4
export GLU_TESS_ERROR5
export GLU_TESS_ERROR6
export GLU_TESS_ERROR7
export GLU_TESS_ERROR8
export GLU_TESS_MISSING_BEGIN_POLYGON
export GLU_TESS_MISSING_BEGIN_CONTOUR
export GLU_TESS_MISSING_END_POLYGON
export GLU_TESS_MISSING_END_CONTOUR
export GLU_TESS_COORD_TOO_LARGE
export GLU_TESS_NEED_COMBINE_CALLBACK
export GLU_TESS_WINDING_ODD
export GLU_TESS_WINDING_NONZERO
export GLU_TESS_WINDING_POSITIVE
export GLU_TESS_WINDING_NEGATIVE
export GLU_TESS_WINDING_ABS_GEQ_TWO
export GLU_TESS_MAX_COORD

end


text:
# variation.jl

"""
variation(α::Float64; v = 7.0, w = 6.0)

    Calculate v and w variational Polaron parameters that minimise the free energy, for the supplied α Frohlich coupling.
    This version uses the original athermal action (Feynman 1955).
    Returns v, w.
"""
function variation(α; v = 0.0, w = 0.0)

    # Intial guess for v and w.
    if v == 0.0 || w == 0.0 # Default values to start with. Generates a random float between 1.0 and 11.0
        initial = sort(rand(2), rev=true) .* 10.0 .+ 1.0
    else
        initial = [v, w]
    end

    # Limits of the optimisation.
    lower = [0.0, 0.0]
    upper = [Inf, Inf]

    # Osaka Free Energy function to minimise.
    f(x) = free_energy(x[1], x[2], α)

    # 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()),
    )

    # Get v and w values that minimise the free energy.
    v, w = Optim.minimizer(solution)

    # If optimisation does not converge or if v ≤ w, pick new random starting guesses for v and w between 1.0 and 11.0. Repeat until the optimisation converges with v > w.
    while Optim.converged(solution) == false || v <= w
        initial = sort(rand(2), rev=true) .* 10.0 .+ 1.0
        solution = Optim.optimize(
            Optim.OnceDifferentiable(f, initial; autodiff = :forward),
            lower,
            upper,
            initial,
            Fminbox(BFGS()),
        )
        v, w = Optim.minimizer(solution)
    end

    # Return variational parameters that minimise the free energy.
    return v, w
end

"""
variation(α::Float64, β::Float64; v = 3.0, w = 2.0)

    Calculate v and w variational Polaron parameters that minimise the free energy, for the supplied α Frohlich coupling and β reduced thermodynamic beta.
    This version uses temperature dependent action (Osaka 1959).
    Returns v, w.
"""
function variation(α, β; v = 0.0, w = 0.0)

    # Intial guess for v and w.
    if v == 0.0 || w == 0.0 # Default values to start with. Generates a random float between 1.0 and 11.0
        initial = sort(rand(2), rev=true) .* 10.0 .+ 1.0
    else
        initial = [v, w]
    end

    # Limits of the optimisation.
    lower = [0.0, 0.0]
    upper = [Inf, Inf]

    # Osaka Free Energy function to minimise.
    f(x) = free_energy(x[1], x[2], α, β)

    # 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()),
    )

    # Get v and w values that minimise the free energy.
    v, w = Optim.minimizer(solution)

    # If optimisation does not converge or if v ≤ w, pick new random starting guesses for v and w between 1.0 and 11.0. Repeat until the optimisation converges with v > w.
    while Optim.converged(solution) == false || v <= w
        initial = sort(rand(2), rev=true) .* 10.0 .+ 1.0
        solution = Optim.optimize(
            Optim.OnceDifferentiable(f, initial; autodiff = :forward),
            lower,
            upper,
            initial,
            Fminbox(BFGS()),
        )
        v, w = Optim.minimizer(solution)
    end

    # Return variational parameters that minimise the free energy.
    return v, w
end


text:
# FunctionAlgebra.jl

# [[file:~/Documents/Julia/scrap.org::*FunctionAlgebra.jl][FunctionAlgebra.jl:1]]
function Base.:+(f1::Union{Function,Operator}, f2::Union{Function,Operator})
    if f1 == f2
	t -> 2*f1(t)
    else
	t -> f1(t) + f2(t)
    end
end

Base.:+(a::Number, f::Union{Function,Operator}) = t -> a + f(t)
Base.:+(f::Union{Function,Operator}, a::Number) = t -> f(t) + a

#_____________________________________________
# Subtraction
function Base.:-(f1::Union{Function,Operator}, f2::Union{Function,Operator})
    if f1 == f2
	function (t)
	    0
	end
    else
	function (t)
	    f1(t) - f2(t)
	end
    end
end

Base.:-(a::Number, f::Union{Function,Operator}) = t -> a - f(t)
Base.:-(f::Union{Function,Operator}, a::Number) = t -> f(t) - a

#_____________________________________________
# Multiplication
function Base.:*(f1::Function, f2::Function)
    if f1 == f2
	t -> f1(t)*f1(t)
    else
	t -> f1(t)*f2(t)
    end
end

function Base.:*(f1::Operator, f2::Operator)
    t -> (f1(f2))(t)
end

Base.:*(a::Number, f::Union{Function,Operator}) = t -> a*f(t)
Base.:*(f::Union{Function,Operator}, a::Number) = *(a, f)

#_____________________________________________
# Division

function Base.:/(f1::Function, f2::Function)
    if f1 == f2
	t -> 1
    else
	t -> f1(t)/f2(t)
    end
end

Base.:/(a::Number, f::Function) = t -> a/f(t)
Base.:/(f::Function, a::Number) = t -> f(t)/a

#_____________________________________________
# Exponents

function Base.:^(a::Operator, b::Integer)
    function (t)
	foldl((x,y)->a(x),1:b,init=t)
    end
end
# FunctionAlgebra.jl:1 ends here


text:
using Test
using KnetMetrics

@testset "classification-metrics" begin
    random_true = rand(2:6,1000)
    random_pred = rand(2:6,1000)
    random_labels = [2,3,4,5,6]
    random_confusion_matrix = confusion_matrix(random_true, random_pred, labels = random_labels)

    @testset "binary-confusion-matrix" begin
        y_true = [1,1,2,2]
        y_pred = [1,1,2,2]
        x = confusion_matrix(y_true, y_pred)
        @test x.matrix == [2 0; 0 2]
    end

    y_true = [3, 2, 4, 3, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5, 3, 4, 4, 2]
    y_pred = [3, 2, 5, 3, 2, 4, 4, 3, 2, 3, 4, 2, 2, 3, 4, 2, 3, 3, 4, 4]
    labels = [2, 3, 4, 5]

    x = confusion_matrix(y_true, y_pred, labels = labels)

    @testset "multiclass-confusion-matrix" begin
        @test x.matrix == [4 0 1 0; 0 5 1 0; 0 1 4 1; 2 1 0 0]
    end

    @testset "class-confusion" begin
        @test class_confusion(x,ith_class =2) == [5 2; 1 12]
        @test class_confusion(random_confusion_matrix,ith_class =2) == class_confusion(random_confusion_matrix, class_name = 3)
    end

    @testset "condition-positive" begin
        @test condition_positive(x) == [5, 6, 6, 3]
        @test condition_positive(x, ith_class = 2) == 6
        @test condition_positive(random_confusion_matrix, class_name = 3) == condition_positive(random_confusion_matrix, ith_class = 2)
    end

    @testset "condition-negative" begin
        @test condition_negative(x) == [15, 14, 14, 17]
        @test condition_negative(x,ith_class = 3) == 14
        @test condition_negative(random_confusion_matrix, class_name = 3) == condition_negative(random_confusion_matrix, ith_class = 2)
    end

    @testset "sensitivity-recall-score" begin
        @test isapprox(sensitivity_score(x), [0.8, 0.8333333333333334, 0.6666666666666666, 0.0])
        @test sensitivity_score(random_confusion_matrix) == recall_score(random_confusion_matrix)
    end

    @testset "specificity-score" begin
        @test isapprox(specificity_score(x),[0.8666666666666667, 0.8571428571428571, 0.8571428571428571, 0.9411764705882353])
    end

    @testset "precision-score" begin
        @test isapprox(precision_score(x), [0.6666666666666666, 0.7142857142857143, 0.6666666666666666, 0.0])
    end

    @testset "accuracy-score" begin
        @test isapprox(accuracy_score(x), 0.16499999999999998)
        @test accuracy_score(x, normalize = false) ==  3.3
    end

    @testset "balanced-accuracy-score" begin
        @test isapprox(balanced_accuracy_score(x),[0.8333333333333334, 0.8452380952380952, 0.7619047619047619, 0.47058823529411764])
    end

    @testset "negative-predictive-value" begin
        @test isapprox(negative_predictive_value(x), [0.9285714285714286, 0.9230769230769231, 0.8571428571428571, 0.8421052631578947])
    end

    @testset "false-positive-rate" begin
        @test isapprox(false_positive_rate(x), [0.13333333333333333, 0.14285714285714285, 0.14285714285714285, 0.058823529411764705])
    end

    @testset "false-negative-rate" begin
        @test isapprox(false_negative_rate(x),  [0.2, 0.16666666666666666, 0.3333333333333333, 1.0])
    end

    @testset "false-discovery-rate" begin
        @test isapprox(false_discovery_rate(x) , [0.13333333333333333, 0.14285714285714285, 0.14285714285714285, 0.058823529411764705])
    end

    @testset "false-omission-rate" begin
        @test isapprox(false_omission_rate(x),[0.0714285714285714, 0.07692307692307687, 0.1428571428571429, 0.1578947368421053])
    end

    @testset "f1-score" begin
        @test isapprox(f1_score(x),[0.7272727272727273, 0.7692307692307693, 0.6666666666666666, 0.0])
    end

    @testset "prevalence-threshold" begin
        @test isapprox(prevalence_threshold(x),[0.24494897427831755, 0.22347381718647807, 0.4165977904505312, -4.12310562561766])
    end

    @testset "threat-score" begin
        @test isapprox(threat_score(x),[0.5714285714285714, 0.625, 0.5, 0.0])
    end

    @testset "matthews-correlation-coefficient" begin
        @test isapprox(matthews_correlation_coeff(x),[0.629940788348712, 0.6633880657639324, 0.5238095238095238, -0.09637388493048533])
    end

    @testset "fowlkes-mallows-index" begin
        @test isapprox( fowlkes_mallows_index(x),[0.7302967433402214, 0.7715167498104596, 0.6666666666666666, 0.0])
    end

    @testset "informedness" begin
        @test isapprox(informedness(x) ,[0.6666666666666667, 0.6904761904761905, 0.5238095238095237, -0.05882352941176472])
    end

    @testset "markedness" begin
        @test isapprox(markedness(x), [-0.38095238095238093, -0.34065934065934056, -0.4285714285714286, -1.0])
    end

    @testset "hamming-loss" begin
        @test hamming_loss(x) == 0.35
    end

    @testset "jaccard-score" begin
        @test jaccard_score(x) == [1.25, 1.5, 1.5, 0.75]
        @test isapprox(0.48148148148148145, jaccard_score(x, average = "micro"))
        @test isapprox(0.42410714285714285,  jaccard_score(x, average = "macro"))
    end

    @testset "cohen-kappa-score" begin
        @test isapprox(0.5155709342560555, cohen_kappa_score(x))
    end

end


text:
#using Plots; plotly()
using SpikingNeuralNetworks
const SNN = SpikingNeuralNetworks
#using Unitful
#using Unitful.DefaultSymbols
include(joinpath(@__DIR__, "..", "src", "units.jl")) # FIXME

include("chain.jl")
include("hh_net.jl")
include("hh_neuron.jl")
include("if_net.jl")
include("if_neuron.jl")
include("iz_net.jl")
include("iz_neuron.jl")
include("oja.jl")
include("rate_net.jl")
include("stdp_demo.jl")
# FIXME: include("stdp_song.jl")


text:
using Distributions

abstract type Bandit end
pull(b::Bandit, action_idx::Integer) = 0
regret(b::Bandit, action_idx::Integer) = 0

struct StaticBandit <: Bandit
  action_distributions::AbstractVector{Distribution}
  num_arms::Integer
  optimal_expected_return::Float64 #here the arms are static so the best arm is always the same
end

# TODO: weird typing issue here - have a sense it's a Julia bug or I'm misunderstanding something in v0.6
staticbandit(action_distributions) = StaticBandit(action_distributions, length(action_distributions), maximum([mean(distr) for distr in action_distributions]))

pull(sb::StaticBandit, action_idx::Integer) = rand(sb.action_distributions[action_idx])

regret(sb::StaticBandit, action_idx::Integer) = sb.optimal_expected_return - mean(sb.action_distributions[action_idx])


text:
# packaged version of nested_optim_d_only.jl

import NMRDataSetup
import NMRSpectraSimulator
import NMRSpecifyRegions

include("../src/NMRCalibrate.jl")
import .NMRCalibrate

using LinearAlgebra
using FFTW
import PyPlot
import BSON
import JSON

import OffsetArrays
import Interpolations

import Statistics
import Random

import NLopt
import MultistartOptimization
import FiniteDiff

PyPlot.close("all")
fig_num = 1

Random.seed!(25)
PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])

# for convinence.
T = Float64
combinevectors = NMRSpectraSimulator.combinevectors

### user inputs.
region_min_dist = 0.1 # in ppm.

fit_config_path = "/home/roy/Documents/repo/NMRData/input/fit_configs/calibrate_700MHz_type1_select_compounds.json"

project_name = "NRC-glucose-2018"

project_base_folder = "/home/roy/MEGAsync/outputs/NMR/calibrate/NRC"
project_folder = joinpath(project_base_folder, project_name)
#mkpath(project_folder)
### end inputs.

### load block.
load_path = joinpath(project_folder, "proxy.bson")
dict = BSON.load(load_path)

u_min = dict[:u_min]
u_max = dict[:u_max]
y = convert(Vector{Complex{Float64}}, dict[:y])
U_y = convert(Vector{Float64}, dict[:U_y])
SW = dict[:SW]
fs = dict[:fs]
ν_0ppm  = dict[:ν_0ppm]
λ0 = dict[:λ0]
molecule_names = dict[:molecule_names]
#w = dict[:w]
w = ones(length(dict[:As]))

As = collect( dict[:As][i] for i = 1:length(dict[:As]) )
Bs = collect( dict[:Bs][i] for i = 1:length(dict[:Bs]) )

### end load.

### alternatively, load from experiment.

### prepare positions.
hz2ppmfunc = uu->(uu - ν_0ppm)*SW/fs
ppm2hzfunc = pp->(ν_0ppm + pp*fs/SW)
P_y = hz2ppmfunc.(U_y)

ΩS0 = NMRSpectraSimulator.getΩS(As)
ΩS0_ppm = NMRSpectraSimulator.getPs(ΩS0, hz2ppmfunc)

Δsys_cs, y_cost_all, U_cost_all, P_cost_all, exp_info, cost_inds, cost_inds_set,
    λ_lbs, λ_ubs, κs_β_orderings,
    κs_β_DOFs = NMRCalibrate.prepareoptim(fit_config_path, molecule_names,
    hz2ppmfunc, U_y, y, As; region_min_dist = region_min_dist)

# visualize cost regions.project_folder
# sort(P_cost_set[1]) # etc..
P_cost_set = collect( P_y[cost_inds_set[r]] for r = 1:length(cost_inds_set) )

a_setp, b_setp, minxs,
    rets = NMRCalibrate.setupitpab(0.1, 10, 0.7; optim_algorithm = :LN_BOBYQA)
#

#include("align.jl")
#include("pkg_align.jl")


text:
export Region

typealias Region hsa_region_t

import Base.convert

convert(::Type{UInt64}, r::Region) = r.handle
convert(::Type{Region}, r::UInt64) = Region(r)

import Base.==

(==)(a::Region, b::Region) = a.handle == b.handle

import Base.hash

hash(r::Region) = r.handle

const iterate_regions_cb_ptr = cfunction(iterate_cb, hsa_status_t, (hsa_region_t, Ptr{Void}))

function iterate_regions(a :: Agent, callback :: Function)
    state = IterState(
        function(x::hsa_region_t)
            r = Region(x)
            callback(r)
        end)

    state_ptr = pointer_from_objref(state)

    err = ccall((:hsa_agent_iterate_regions, libhsa), hsa_status_t, (hsa_agent_t, Ptr{Void}, Ptr{Void}),
                a, iterate_regions_cb_ptr, state_ptr)

    if !isnull(state.err)
        rethrow(state.err.value)
    end

    test_status(err)
end

function regions(a :: Agent)
    reg = Array(Region, 0)

    iterate_regions(a,
        function(r)
            push!(reg, r)
        end)

    return reg
end

function memory_register(
    ptr,
    size)

    err = ccall((:hsa_memory_register, libhsa), hsa_status_t, (Ptr{Void}, Csize_t),
    ptr, size)

    test_status(err)
end

function memory_deregister(
    ptr,
    size)

    err = ccall((:hsa_memory_deregister, libhsa), hsa_status_t, (Ptr{Void}, Csize_t),
    ptr, size)

    test_status(err)
end

type Allocation
    ptr::Ptr{Void}

    function Allocation(ptr)
        a = new(ptr)
        finalizer(a, deallocate)

        return a
    end
end

function deallocate(a :: Allocation)
    if a.ptr != C_NULL
        memory_free(a.ptr)
        a.ptr = C_NULL
    end
end

convert(::Type{Ptr{Void}}, a::Allocation) = a.ptr

function memory_allocate(
    region,
    size)

    res = Ref{Ptr{Void}}(0)

    err = ccall((:hsa_memory_allocate, libhsa), hsa_status_t, (
    hsa_region_t,
    Csize_t,
    Ptr{Ptr{Void}}),
    region, size, res)

    test_status(err)

    return Allocation(res.x)
end

function memory_free(ptr)
    err = ccall((:hsa_memory_free, libhsa), hsa_status_t, (
    Ptr{Void},),
    ptr)

    test_status(err)
end


text:
using SimJulia
using Distributions

function visit(customer::Process,
 			time_in_bank::Float64,
 			counter::Resource,
 			max_in_queue::Int)
  arrive = now(customer);
  @printf("%8.4f %s: Here I am\n", arrive, customer);
  if length(counter.wait_queue) < max_in_queue
      request(customer, counter);
      wait = now(customer) – arrive;
      @printf("%8.4f %s: Waited %6.3f\n",
 			now(customer), customer, wait);
      hold(customer, time_in_bank)
      release(customer, counter);
      @printf("%8.4f %s: Finished\n", now(customer), customer)
  else
      @printf("%8.4f %s: Waiting\n", now(customer), customer)
  end
end

function generate(source::Process, number::Int,
 			mean_time_between_arrivals::Float64,
 			mean_time_in_bank::Float64,
 			counter::Resource,
 			max_in_queue::Int)
  d_tba = Exponential(mean_time_of_arrivals);
  d_tib = Exponential(mean_time_in_bank);
  for i = 1:number
      c = Process(simulation(source),
 			@sprintf("Customer%02d", i));
    tib = rand(d_tib);
    activate(c, now(source), visit, tib,
    			counter, max_in_queue);
    tba = rand(d_tba);
    hold(source, tba);
  end
end

function number_in_system(counter::Resource)
  length(counter.active_set) + length(counter.wait_queue);
end

# Set simulation data
max_number = 8;
max_time = 3900.0;
max_queue = 1;
mean_arrival_time = 10.0;
mean_service_time = 12.0;
number_of_tellers = 1;
seed = randomize();	# Use function defined early or pick
# a value to make run reproducible

# Model/Experiment
srand(seed);
sim = Simulation(uint(16));
k = Resource(sim, "Counter", uint(1), false);
s = Process(sim, "Source");
activate(s, 0.0,  generate,
max_number, mean_arrival_time,
 	mean_service_time, k, max_queue);
run(sim, max_time);

max_stars_repo_path: Module 3/Chapter06/banky.jl

max_issues_repo_path: Module 3/Chapter06/banky.jl

max_forks_repo_path: Module 3/Chapter06/banky.jl


text:
using LinearAlgebra, Plots, JSON
include("../src/SFFM.jl")
include("METools.jl")

# define SFM
T = [-2.0 2.0; 1.0 -1.0]
C = [1.0; -2.0]
fn(x) = [ones(size(x)) ones(size(x))]
model = SFFM.Model(;T=T,C=C,r=(r=fn,R=fn),Bounds=[0 10;-Inf Inf])
N₋ = sum(C.<=0)
N₊ = sum(C.>=0)
NPhases = length(C)

Δ = 1
NBases = 11
order = NBases

Erlang = MakeErlang(order,mean=Δ)

Coxian = MakeSomeCoxian(order,mean=Δ)

Coxianbkwd = reversal(Coxian)

t = 0:0.01:4
orbitCoxian = orbit(t,Coxian)
orbitCoxianbkwd = orbit(t,Coxianbkwd)
plot(t,orbitCoxian,label = "fwd")
plot!(t,orbitCoxianbkwd,label = "bkwd")

plot(orbitCoxian[:,2],orbitCoxian[:,3],label = "fwd",legend=:outerright)
plot!(orbitCoxianbkwd[:,2],orbitCoxianbkwd[:,3],label = "bkwd")

αME, QME, ~ = MakeME(CMEParams[NBases],mean=Δ)#MakeSomeCoxian(order,mean=Δ)#
CME = (α = αME, Q = QME, q = sum(QME,dims=2), π = αME*QME^(-1)./sum(αME*QME^(-1)))
CMEr = reversal(CME)
πME = (α = CME.π, Q = QME, π = CME.π)
πMEr = (α = CMEr.π, Q = CMEr.Q, π = CME.π)
idx = [1; [3:2:NBases 2:2:NBases]'[:]]
CMEt = (α = αME[idx]', Q = QME')
D = jumpMatrixD(CME)
MEj = (α = αME*D, Q=CME.Q)

orbitME = orbit(t,CME)
orbitMEr = orbit(t,CMEr)
orbitMED = orbit(t,CME)*D^-1
orbitMEj = orbit(t,MEj,norm=2)
# orbitMEQ = (orbit(t,CME)*-CME.Q^-1)./sum(orbit(t,πME)*-CME.Q^-1,dims=2)
orbitπME = orbit(t,πME)
orbitπMEQ = (orbit(t,πME)*-CME.Q)./sum(orbit(t,πME)*-CME.Q,dims=2)
orbitπMED = orbit(t,πME)*D^-1
orbitπMEr = orbit(t,πMEr)
orbitCMEt = orbit(t,CMEt)
densityME = density(t,CME)
densityMEj = density(t,MEj)
densityπME = density(t,πME)
densityCMEt = density(t,CMEt)
intensityME = intensity(t,CME)
intensityπME = intensity(t,πME)
intensityCMEt = intensity(t,CMEt)
plot(t,orbitME,label = "fwd")
# plot(t,orbitMED,label = "D")
# plot(t,orbitπME,label = "π")
plot(t,orbitCMEt,label = "tr")
plot(orbitME[:,2],orbitME[:,3],label = "ME fwd orbit",legend=:outerright)
# plot!(orbitMED[:,4],orbitMED[:,5],label = "MED fwd",legend=:outerright)
# plot!(orbitπMEQ[:,4],orbitπMEQ[:,5],label = "MEQ fwd",legend=:outerright)
# plot!(orbitMEQ[:,4],orbitMEQ[:,5],label = "MEQ fwd",legend=:outerright)
plot!(orbitMEr[:,4],orbitMEr[:,5],label = "MEr",legend=:outerright)
# plot!(orbitMEj[:,4],orbitMEj[:,5],label = "MEj",legend=:outerright)
# plot!(orbitπME[:,4],orbitπME[:,5],label = "π orbit",legend=:outerright)
plot!(orbitπMEr[:,4],orbitπMEr[:,5],label = "πr orbit",legend=:outerright)
# plot!(orbitπMED[:,4],orbitπMED[:,5],label = "πD orbit",legend=:outerright)
plot!(orbitCMEt[:,4],orbitCMEt[:,5],label = "tr orbit",legend=:outerbottomright)
scatter!([CME.α[4]],[CME.α[5]],label = :false,series_annotations = ["     α"])
shift = CME.α*exp(CME.Q)./sum(exp(CME.Q[1,1]))#orbit(-1,CME,norm=2)
scatter!([shift[4]],[shift[5]],label = :false,series_annotations = ["     α+"])

scatter!([(CME.α*D^-1)[4]],[(CME.α*D^-1)[5]],label = :false,series_annotations = ["     α"])
scatter!([CMEr.α[4]],[CMEr.α[5]],label = :false,series_annotations = ["     α̃"])
scatter!([orbit(1,CME)[4]],[orbit(1,CME)[5]],label = :false, series_annotations = ["       t=1"])
scatter!([(orbit(0.1,CME))[4]],[(orbit(0.1,CME))[5]],label = :false, series_annotations = ["           t=0.1"])
scatter!([(orbit(1,CME)*D^-1)[4]],[(orbit(1,CME)*D^-1)[5]],label = :false, series_annotations = ["                A(t)D, t=1"])
scatter!([orbit(1,CMEt)[4]],[orbit(1,CMEt)[5]],label = :false, series_annotations = ["       t̃=1"])
# τ = 0.5
# scatter!([orbit(τ,CME)[2]],[orbit(τ,CME)[3]],label = :false, series_annotations = ["         t="*string(τ)])
# # scatter!([orbit(1,πME)[2]],[orbit(1,πME)[3]],label = :false, series_annotations = ["        t=1"])
# # scatter!([orbit(2.2885,πME)[2]],[orbit(2.2885,πME)[3]],label = :false, series_annotations = ["                                πexp(Qt)Q, t=2.2885"])
# tstar = orbit_zero(CME,orbit(1,CMEt))
# val = orbit(τ,CME)*exp(CME.Q*tstar)./sum(exp(CME.Q*tstar)[1,1])
# scatter!([val[2]],[val[3]],label = :false,
#         series_annotations = ["                         A(t)exp(Q) t=0.1"])
# Atplus = val./sum(val)
# scatter!([Atplus[2]],[Atplus[3]],label = :false,
#         series_annotations = ["                         A(t)exp(Q) t=0.1"])
#
# scatter!([orbit(2,CME)[2]],[orbit(2,CME)[3]],label = :false, series_annotations = ["       t=2"])
# # scatter!([orbit(0.1,CMEr)[2]],[orbit(0.1,CMEr)[3]],label = :false, series_annotations = ["          t̃=0.1"])
# # scatter!([orbit(2.05,CME)[2]],[orbit(2.05,CME)[3]], label = "t=2")
# # scatter!([CME.π[2]],[CME.π[3]],label = :false, series_annotations = ["      π"])

CME = (α = αME, Q = QME, q = sum(QME,dims=2), π = αME*QME^(-1)./sum(αME*QME^(-1)))
CMEminus = (α = αME, Q = QME'-12*I, q = sum(QME,dims=2), π = αME*QME^(-1)./sum(αME*QME^(-1)))

orbitME = orbit(t,CME)
orbitMEminus = orbit(t,CMEminus)
plot(orbitME[:,4],orbitME[:,5],label = "ME fwd orbit",legend=:outerright)
plot!(orbitMEminus[:,4],orbitMEminus[:,5],label = "ME fwd orbit",legend=:outerright)

display(CMEminus.Q)

densityCME = density(t,CME)
intensityCME = intensity(t,CME)
densityCMEminus = density(t,CMEminus)
intensityCMEminus = intensity(t,CMEminus)

plot(t,intensityCME,label = "ME")
plot!(t,intensityCMEminus,label = "ME-")

plot(t,densityCME,label = "ME")
plot!(t,densityCMEminus,label = "ME-")

#
plot(t,densityME,label = "ME")
plot!(t,densityπME,label = "πME")
plot!(t,densityCMEt,label = "CMEt")
plot!(t,densityMEj,label = "MEj")

# CCDF(t) = sum(CME.α*exp(CME.Q*t))
# plot!(t,CCDF.(t),label = "πME")
#
# S = Coxian.Q
# s = -sum(Coxian.Q,dims=2)
# α = Coxian.α
# m₁ = α*(-S)^-1
# π₁ = m₁./sum(m₁)
# m₂ = -π₁*S^-1
# ρ₁ = sum(m₁)
# α̂₁ = [s'*diagm(m₂[:]) zeros(1,order)]
# Ŝ₁ = [diagm(m₂[:])^-1*S'*diagm(m₂[:]) ρ₁^-1*diagm(m₂[:])^-1*diagm(m₁[:]);
#         zeros(order,order) diagm(m₁[:])^-1*S'*diagm(m₁[:])]
# Spreadbkwd = (α = α̂₁, Q = Ŝ₁)
#
# vecρ₁ = sum(-S^-1,dims=2)
# α₁ = [ρ₁^-1*α*diagm(vecρ₁[:]) zeros(1,order)]
# S₁ = [diagm(vecρ₁[:])^-1*S*diagm(vecρ₁[:]) diagm(vecρ₁[:])^-1;
#         zeros(order,order) S]
# Spread = (α = α₁, Q = S₁)
#
#
# plot(t,density(t,Coxian))
# plot!(t,density(t,Spreadbkwd))
# plot!(t,density(t,Spread))
#
# MEα₁ = [CME.α*S^-1/(sum(CME.α*S^-1)) 0*CME.α]
# MES₁ = [CME.Q -CME.Q; 0*CME.Q CME.Q]
# MEs₁ = sum(-MES₁,dims=1)

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
using Mocking

cancel_key_deletion_1_patch = @patch function cancel_key_deletion(a...; b...)
    return Dict{String, Any}("KeyId" => "arn:aws:kms:us-east-2:111122223333:key/1234abcd-12ab-34cd-56ef-1234567890ab")
end

create_grant_1_patch = @patch function create_grant(a...; b...)
    return Dict{String, Any}("GrantToken" => "AQpAM2RhZTk1MGMyNTk2ZmZmMzEyYWVhOWViN2I1MWM4Mzc0MWFiYjc0ZDE1ODkyNGFlNTIzODZhMzgyZjBlNGY3NiKIAgEBAgB4Pa6VDCWW__MSrqnre1HIN0Grt00ViSSuUjhqOC8OT3YAAADfMIHcBgkqhkiG9w0BBwaggc4wgcsCAQAwgcUGCSqGSIb3DQEHATAeBglghkgBZQMEAS4wEQQMmqLyBTAegIn9XlK5AgEQgIGXZQjkBcl1dykDdqZBUQ6L1OfUivQy7JVYO2-ZJP7m6f1g8GzV47HX5phdtONAP7K_HQIflcgpkoCqd_fUnE114mSmiagWkbQ5sqAVV3ov-VeqgrvMe5ZFEWLMSluvBAqdjHEdMIkHMlhlj4ENZbzBfo9Wxk8b8SnwP4kc4gGivedzFXo-dwN8fxjjq_ZZ9JFOj2ijIbj5FyogDCN0drOfi8RORSEuCEmPvjFRMFAwcmwFkN2NPp89amA", "GrantId" => "0c237476b39f8bc44e45212e08498fbe3151305030726c0590dd8d3e9f3d6a60")
end

create_key_1_patch = @patch function create_key(a...; b...)
    return Dict{String, Any}("KeyMetadata" => Dict{String, Any}("Arn" => "arn:aws:kms:us-east-2:111122223333:key/1234abcd-12ab-34cd-56ef-1234567890ab", "Description" => "", "KeyState" => "Enabled", "KeyUsage" => "ENCRYPT_DECRYPT", "KeyManager" => "CUSTOMER", "Origin" => "AWS_KMS", "CreationDate" => "2017-07-05T14:04:55-07:00", "KeyId" => "1234abcd-12ab-34cd-56ef-1234567890ab", "AWSAccountId" => "111122223333", "Enabled" => true))
end

decrypt_1_patch = @patch function decrypt(a...; b...)
    return Dict{String, Any}("Plaintext" => "<binary data>", "KeyId" => "arn:aws:kms:us-west-2:111122223333:key/1234abcd-12ab-34cd-56ef-1234567890ab")
end

describe_key_1_patch = @patch function describe_key(a...; b...)
    return Dict{String, Any}("KeyMetadata" => Dict{String, Any}("Arn" => "arn:aws:kms:us-east-2:111122223333:key/1234abcd-12ab-34cd-56ef-1234567890ab", "Description" => "", "KeyState" => "Enabled", "KeyUsage" => "ENCRYPT_DECRYPT", "KeyManager" => "CUSTOMER", "Origin" => "AWS_KMS", "CreationDate" => "2017-07-05T14:04:55-07:00", "KeyId" => "1234abcd-12ab-34cd-56ef-1234567890ab", "AWSAccountId" => "111122223333", "Enabled" => true))
end

encrypt_1_patch = @patch function encrypt(a...; b...)
    return Dict{String, Any}("CiphertextBlob" => "<binary data>", "KeyId" => "arn:aws:kms:us-west-2:111122223333:key/1234abcd-12ab-34cd-56ef-1234567890ab")
end

generate_data_key_1_patch = @patch function generate_data_key(a...; b...)
    return Dict{String, Any}("CiphertextBlob" => "<binary data>", "Plaintext" => "<binary data>", "KeyId" => "arn:aws:kms:us-east-2:111122223333:key/1234abcd-12ab-34cd-56ef-1234567890ab")
end

generate_data_key_without_plaintext_1_patch = @patch function generate_data_key_without_plaintext(a...; b...)
    return Dict{String, Any}("CiphertextBlob" => "<binary data>", "KeyId" => "arn:aws:kms:us-east-2:111122223333:key/1234abcd-12ab-34cd-56ef-1234567890ab")
end

generate_random_1_patch = @patch function generate_random(a...; b...)
    return Dict{String, Any}("Plaintext" => "<binary data>")
end

get_key_policy_1_patch = @patch function get_key_policy(a...; b...)
    return Dict{String, Any}("Policy" => "{\n  \"Version\" : \"2012-10-17\",\n  \"Id\" : \"key-default-1\",\n  \"Statement\" : [ {\n    \"Sid\" : \"Enable IAM User Permissions\",\n    \"Effect\" : \"Allow\",\n    \"Principal\" : {\n      \"AWS\" : \"arn:aws:iam::111122223333:root\"\n    },\n    \"Action\" : \"kms:*\",\n    \"Resource\" : \"*\"\n  } ]\n}")
end

get_key_rotation_status_1_patch = @patch function get_key_rotation_status(a...; b...)
    return Dict{String, Any}("KeyRotationEnabled" => true)
end

get_parameters_for_import_1_patch = @patch function get_parameters_for_import(a...; b...)
    return Dict{String, Any}("ParametersValidTo" => "2016-12-01T14:52:17-08:00", "ImportToken" => "<binary data>", "KeyId" => "arn:aws:kms:us-east-2:111122223333:key/1234abcd-12ab-34cd-56ef-1234567890ab", "PublicKey" => "<binary data>")
end

list_aliases_1_patch = @patch function list_aliases(a...; b...)
    return Dict{String, Any}("Truncated" => false, "Aliases" => Any[Dict{String, Any}("AliasArn" => "arn:aws:kms:us-east-2:111122223333:alias/aws/acm", "AliasName" => "alias/aws/acm", "TargetKeyId" => "da03f6f7-d279-427a-9cae-de48d07e5b66"), Dict{String, Any}("AliasArn" => "arn:aws:kms:us-east-2:111122223333:alias/aws/ebs", "AliasName" => "alias/aws/ebs", "TargetKeyId" => "25a217e7-7170-4b8c-8bf6-045ea5f70e5b"), Dict{String, Any}("AliasArn" => "arn:aws:kms:us-east-2:111122223333:alias/aws/rds", "AliasName" => "alias/aws/rds", "TargetKeyId" => "7ec3104e-c3f2-4b5c-bf42-bfc4772c6685"), Dict{String, Any}("AliasArn" => "arn:aws:kms:us-east-2:111122223333:alias/aws/redshift", "AliasName" => "alias/aws/redshift", "TargetKeyId" => "08f7a25a-69e2-4fb5-8f10-393db27326fa"), Dict{String, Any}("AliasArn" => "arn:aws:kms:us-east-2:111122223333:alias/aws/s3", "AliasName" => "alias/aws/s3", "TargetKeyId" => "d2b0f1a3-580d-4f79-b836-bc983be8cfa5"), Dict{String, Any}("AliasArn" => "arn:aws:kms:us-east-2:111122223333:alias/example1", "AliasName" => "alias/example1", "TargetKeyId" => "4da1e216-62d0-46c5-a7c0-5f3a3d2f8046"), Dict{String, Any}("AliasArn" => "arn:aws:kms:us-east-2:111122223333:alias/example2", "AliasName" => "alias/example2", "TargetKeyId" => "f32fef59-2cc2-445b-8573-2d73328acbee"), Dict{String, Any}("AliasArn" => "arn:aws:kms:us-east-2:111122223333:alias/example3", "AliasName" => "alias/example3", "TargetKeyId" => "1374ef38-d34e-4d5f-b2c9-4e0daee38855")])
end

list_grants_1_patch = @patch function list_grants(a...; b...)
    return Dict{String, Any}("Truncated" => true, "Grants" => Any[Dict{String, Any}("Operations" => Any["Encrypt", "ReEncryptFrom", "ReEncryptTo"], "RetiringPrincipal" => "acm.us-east-2.amazonaws.com", "CreationDate" => "2016-10-25T14:37:41-07:00", "GrantId" => "91ad875e49b04a9d1f3bdeb84d821f9db6ea95e1098813f6d47f0c65fbe2a172", "IssuingAccount" => "arn:aws:iam::111122223333:root", "GranteePrincipal" => "acm.us-east-2.amazonaws.com", "KeyId" => "arn:aws:kms:us-east-2:111122223333:key/1234abcd-12ab-34cd-56ef-1234567890ab"), Dict{String, Any}("Operations" => Any["ReEncryptFrom", "ReEncryptTo"], "RetiringPrincipal" => "acm.us-east-2.amazonaws.com", "CreationDate" => "2016-10-25T14:37:41-07:00", "GrantId" => "a5d67d3e207a8fc1f4928749ee3e52eb0440493a8b9cf05bbfad91655b056200", "IssuingAccount" => "arn:aws:iam::111122223333:root", "GranteePrincipal" => "acm.us-east-2.amazonaws.com", "KeyId" => "arn:aws:kms:us-east-2:111122223333:key/1234abcd-12ab-34cd-56ef-1234567890ab"), Dict{String, Any}("Operations" => Any["Encrypt", "ReEncryptFrom", "ReEncryptTo"], "RetiringPrincipal" => "acm.us-east-2.amazonaws.com", "CreationDate" => "2016-10-25T14:37:41-07:00", "GrantId" => "c541aaf05d90cb78846a73b346fc43e65be28b7163129488c738e0c9e0628f4f", "IssuingAccount" => "arn:aws:iam::111122223333:root", "GranteePrincipal" => "acm.us-east-2.amazonaws.com", "KeyId" => "arn:aws:kms:us-east-2:111122223333:key/1234abcd-12ab-34cd-56ef-1234567890ab"), Dict{String, Any}("Operations" => Any["Encrypt", "ReEncryptFrom", "ReEncryptTo"], "RetiringPrincipal" => "acm.us-east-2.amazonaws.com", "CreationDate" => "2016-10-25T14:37:41-07:00", "GrantId" => "dd2052c67b4c76ee45caf1dc6a1e2d24e8dc744a51b36ae2f067dc540ce0105c", "IssuingAccount" => "arn:aws:iam::111122223333:root", "GranteePrincipal" => "acm.us-east-2.amazonaws.com", "KeyId" => "arn:aws:kms:us-east-2:111122223333:key/1234abcd-12ab-34cd-56ef-1234567890ab")])
end

list_key_policies_1_patch = @patch function list_key_policies(a...; b...)
    return Dict{String, Any}("PolicyNames" => Any["default"], "Truncated" => false)
end

list_keys_1_patch = @patch function list_keys(a...; b...)
    return Dict{String, Any}("Truncated" => false, "Keys" => Any[Dict{String, Any}("KeyId" => "0d990263-018e-4e65-a703-eff731de951e", "KeyArn" => "arn:aws:kms:us-east-2:111122223333:key/0d990263-018e-4e65-a703-eff731de951e"), Dict{String, Any}("KeyId" => "144be297-0ae1-44ac-9c8f-93cd8c82f841", "KeyArn" => "arn:aws:kms:us-east-2:111122223333:key/144be297-0ae1-44ac-9c8f-93cd8c82f841"), Dict{String, Any}("KeyId" => "21184251-b765-428e-b852-2c7353e72571", "KeyArn" => "arn:aws:kms:us-east-2:111122223333:key/21184251-b765-428e-b852-2c7353e72571"), Dict{String, Any}("KeyId" => "214fe92f-5b03-4ae1-b350-db2a45dbe10c", "KeyArn" => "arn:aws:kms:us-east-2:111122223333:key/214fe92f-5b03-4ae1-b350-db2a45dbe10c"), Dict{String, Any}("KeyId" => "339963f2-e523-49d3-af24-a0fe752aa458", "KeyArn" => "arn:aws:kms:us-east-2:111122223333:key/339963f2-e523-49d3-af24-a0fe752aa458"), Dict{String, Any}("KeyId" => "b776a44b-df37-4438-9be4-a27494e4271a", "KeyArn" => "arn:aws:kms:us-east-2:111122223333:key/b776a44b-df37-4438-9be4-a27494e4271a"), Dict{String, Any}("KeyId" => "deaf6c9e-cf2c-46a6-bf6d-0b6d487cffbb", "KeyArn" => "arn:aws:kms:us-east-2:111122223333:key/deaf6c9e-cf2c-46a6-bf6d-0b6d487cffbb")])
end

list_resource_tags_1_patch = @patch function list_resource_tags(a...; b...)
    return Dict{String, Any}("Truncated" => false, "Tags" => Any[Dict{String, Any}("TagValue" => "87654", "TagKey" => "CostCenter"), Dict{String, Any}("TagValue" => "ExampleUser", "TagKey" => "CreatedBy"), Dict{String, Any}("TagValue" => "Test", "TagKey" => "Purpose")])
end

list_retirable_grants_1_patch = @patch function list_retirable_grants(a...; b...)
    return Dict{String, Any}("Truncated" => false, "Grants" => Any[Dict{String, Any}("Operations" => Any["Decrypt", "Encrypt"], "RetiringPrincipal" => "arn:aws:iam::111122223333:role/ExampleRole", "CreationDate" => "2016-12-07T11:09:35-08:00", "GrantId" => "0c237476b39f8bc44e45212e08498fbe3151305030726c0590dd8d3e9f3d6a60", "IssuingAccount" => "arn:aws:iam::444455556666:root", "GranteePrincipal" => "arn:aws:iam::111122223333:role/ExampleRole", "KeyId" => "arn:aws:kms:us-east-2:444455556666:key/1234abcd-12ab-34cd-56ef-1234567890ab")])
end

re_encrypt_1_patch = @patch function re_encrypt(a...; b...)
    return Dict{String, Any}("CiphertextBlob" => "<binary data>", "SourceKeyId" => "arn:aws:kms:us-east-2:111122223333:key/1234abcd-12ab-34cd-56ef-1234567890ab", "KeyId" => "arn:aws:kms:us-east-2:111122223333:key/0987dcba-09fe-87dc-65ba-ab0987654321")
end

schedule_key_deletion_1_patch = @patch function schedule_key_deletion(a...; b...)
    return Dict{String, Any}("KeyId" => "arn:aws:kms:us-east-2:111122223333:key/1234abcd-12ab-34cd-56ef-1234567890ab", "DeletionDate" => "2016-12-17T16:00:00-08:00")
end


text:
using ImperoHooks
using ClimateMachine
using ClimateMachine.Mesh.Grids
using Impero, MPI
using Makie, GLMakie, AbstractPlotting
using ImageTransformations, Colors
using AbstractPlotting.MakieLayout
using GaussQuadrature, Random, Test
using Statistics, Printf
using ClimateMachine.MPIStateArrays
using Base.Threads, LinearAlgebra
import Impero: Circle
import ImperoHooks: ScalarField
# include(pwd() * "/examples/gridhelper.jl")
##
ClimateMachine.init()
const ArrayType = ClimateMachine.array_type()
const mpicomm = MPI.COMM_WORLD
const FT = Float64
# Domain
Ω = Circle(-1,1) × Circle(-1,1) × Circle(-1,1)
# Discretization
ex, ey, ez = (3,5,7)
(npx, npy, npz) = (1,1,1)
ClimateMachine.gpu_allowscalar(true)
grid = DiscontinuousSpectralElementGrid(Ω, elements = (ex, ey, ez), polynomialorder = (npx, npy, npz), array = ArrayType)
gridhelper = GridHelper(grid)
x, y, z = coordinates(grid)
xC, yC, zC = cellcenters(grid)
ϕ = ScalarField(copy(x), gridhelper)
##
f = MPIStateArray{FT}(mpicomm, ArrayType, size(x)...,1)
g(x,y,z) = 1 - x^2 - y^2 - z^2
@. f =  g(x,y,z)
@. ϕ = f[:, :, 1]
newsize = (20, 20, 20)
nx, ny, nz = newsize
newx = range(-1,1, length = nx)
newy = range(-1,1, length = ny)
newz = range(-1,1, length = nz)
# can use
ϕf = ϕ(newx, newy, newz) # interpolate to the new grid

newf = zeros((nx,ny,nz))
checkf = zeros((nx,ny,nz))
tic = time()
p = gridhelper.element.permutation
lin = gridhelper.element.cartesianindex
@threads for i in 1:nx
    for j in 1:ny
        for k in 1:nz
            location = (newx[i], newy[j], newz[k])
            #findelement(xC, yC, zC, location, p, lin)
            newf[i,j,k] = getvalue(f, location, gridhelper)
            checkf[i,j,k] = g(location[1], location[2], location[3])
        end
    end
end
toc = time()

println("Time to interpolate is ", toc - tic)
println("The error is ", norm(newf - checkf)/ norm(checkf))
println("The maximum difference is ", maximum(abs.(newf - checkf)))
println("The largest error is at ", argmax(abs.(newf - checkf)))
errorloc = argmax(abs.(newf - checkf))
i, j, k = errorloc[1], errorloc[2], errorloc[3]
println("The largest error is at the point ",(x[i], y[j], z[k]))

## Note always put LScenes before LText for now
resolution = (3000, 860)
width = round(Int, resolution[1] / 6)
scene, layout = layoutscene(resolution = resolution)
# LScenes first! (bad plots otherwise)
lscene = layout[2:4, 2:4] = LScene(scene)
lscene3 = layout[2:4, 8:10] = LScene(scene)
lscene2 = layout[2:4, 5:7] = LScene(scene)
# LTexts next! (bad plots otherwise)
layout[1, 8:10] = LText(scene, "|Difference|", textsize = 50)
layout[1, 2:4] = LText(scene, "MPIStateArray", textsize = 50)
layout[1, 5:7] = LText(scene, "Exact", textsize = 50)

layout[1,1] = LText(scene, "Menu", width = width, textsize = 50)
layout[1,11] = LText(scene, "Menu", width = width, textsize = 50)

upperclim_slider = LSlider(scene, range = range(0, 1, length = 101), startvalue = 0.99)
upperclim_node = upperclim_slider.value
lowerclim_slider = LSlider(scene, range = range(0, 1, length = 101), startvalue = 0.01)
lowerclim_node = lowerclim_slider.value

upperclim_slider2 = LSlider(scene, range = range(0, 1, length = 101), startvalue = 0.99)
upperclim_node2 = upperclim_slider2.value
lowerclim_slider2 = LSlider(scene, range = range(0, 1, length = 101), startvalue = 0.01)
lowerclim_node2 = lowerclim_slider2.value

clims = @lift((quantile(f[:], $lowerclim_node) , quantile(f[:], $upperclim_node)))
cmap_rgb = to_colormap(:balance);

vs1 = volume!(lscene, -1..1, -1..1, -1..1, newf, colormap = cmap_rgb, colorrange = clims, camera = cam3d!, width = width)
# axis = scene.children[1].plots[1]
# axis[:names"][:axisnames] = ("up", "down", "left")
# axis[:names][:textsize] = (10.0, 10.0, 20.0)

volume!(lscene2, -1..1, -1..1, -1..1, checkf, camera = cam3d!, colormap = cmap_rgb, colorrange = clims)

Δf = abs.(checkf - newf)
clims2 = @lift((quantile(Δf[:], $lowerclim_node2) , quantile(Δf[:], $upperclim_node2)))
volume!(lscene3, -1..1, -1..1, -1..1, Δf, colormap = cmap_rgb, colorrange = clims2, camera = cam3d!)

lowerclim_string = @lift("lower clim quantile = " *  @sprintf("%0.2f", $lowerclim_node) * ", value = " * @sprintf("%0.1e", $clims[1]))
upperclim_string = @lift("upper clim quantile = " *  @sprintf("%0.2f", $upperclim_node) * ", value = " * @sprintf("%0.1e", $clims[2]))
layout[2, 1] = vgrid!(
    LText(scene, lowerclim_string, width = nothing),
    lowerclim_slider,
    LText(scene, upperclim_string, width = nothing),
    upperclim_slider,
)

lowerclim_string2 = @lift("lower clim quantile = " *  @sprintf("%0.2f", $lowerclim_node2) * ", value = " * @sprintf("%0.1e", $clims2[1]))
upperclim_string2 = @lift("upper clim quantile = " *  @sprintf("%0.2f", $upperclim_node2) * ", value = " * @sprintf("%0.1e", $clims2[2]))
layout[2, 11] = vgrid!(
    LText(scene, lowerclim_string2, width = nothing),
    lowerclim_slider2,
    LText(scene, upperclim_string2, width = nothing),
    upperclim_slider2,
)

display(scene)

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks


text:
module DLMReader

using InMemoryDatasets
using Reexport
@reexport using Dates, InlineStrings
export filereader,
    filewriter,
    # informats
    COMMA!,
    COMMAX!,
    NUM_NA!,
    CHAR_NA!,
    BOOL!

include("linebuffer.jl")
include("util.jl")
include("informats.jl")
include("parsers.jl")
include("reader.jl")
include("writer.jl")

# needs work
# include("precompile.jl")

end


text:
function extract_coef_info(coefs, ix)
    # 1 = basisname, 2 = coefname, 3 = colname
    return [c[ix] for c in split.(coefs, " : ")]
end

function get_coefnames(uf::Union{UnfoldModel,DesignMatrix})
    coefnames = Unfold.coefnames(formula(uf))
    coefnames = vcat(coefnames...) # gets rid of an empty Any() XXX not sure where it comes from, only in MixedModel Timexpanded case
end

modelfit(uf::UnfoldModel) = uf.modelfit
StatsModels.coef(uf::UnfoldModel) = coef(modelfit(uf))
StatsModels.coef(mf::LinearModelFit) = mf.estimate

function StatsModels.coeftable(
    uf::Union{UnfoldLinearModelContinuousTime,UnfoldLinearMixedModelContinuousTime},
)
    coefsRaw = get_coefnames(uf)
    coefs = extract_coef_info(coefsRaw, 2)
    #colnames_basis_raw = get_colnames_basis(formula(uf))# this is unconverted basisfunction basis,
    colnames_basis = extract_coef_info(coefsRaw, 3) # this is converted to strings!
    basisnames = extract_coef_info(coefsRaw, 1)
    @debug coefs
    @debug colnames_basis

    nchan = size(coef(uf), 1)

    coefs_rep = permutedims(repeat(coefs, 1, nchan), [2, 1])

    if length(colnames_basis) == 1
        colnames_basis = [colnames_basis]
    end

    #colnames_basis_rep = permutedims(repeat(colnames_basis_raw,Int(length(colnames_basis)/length(colnames_basis_raw)),nchan),[2,1])
    # XXX HOTFIX for above lines
    colnames_basis_rep = permutedims(repeat(colnames_basis, 1, nchan), [2, 1])
    try
        colnames_basis_rep = parse.(Float64, colnames_basis_rep)
    catch
        #do nothing
    end
    chan_rep = repeat(1:nchan, 1, size(colnames_basis_rep, 2))

    basisnames_rep = permutedims(repeat(basisnames, 1, nchan), [2, 1])
    return make_long_df(
        uf,
        coefs_rep,
        chan_rep,
        colnames_basis_rep,
        basisnames_rep,
        collabel(uf),
    )
end

function StatsModels.coeftable(uf::Union{UnfoldLinearModel,UnfoldLinearMixedModel})
    # Mass Univariate Case
    coefnames = get_coefnames(uf)

    colnames_basis = collect(values(design(uf)))[1][2]#times#get_colnames_basis(m.X.formulas,times)

    @debug "coefs: $(size(coefnames)),colnames_basis:$(size(colnames_basis)))"
    @debug "coefs: $coefnames,colnames_basis:$colnames_basis"
    nchan = size(coef(uf), 1)
    ncols = length(colnames_basis)
    ncoefs = length(coefnames)

    # replicate
    coefs_rep = permutedims(repeat(coefnames, outer = [1, nchan, ncols]), [2, 3, 1])
    colnames_basis_rep = permutedims(repeat(colnames_basis, 1, nchan, ncoefs), [2 1 3])
    chan_rep = repeat(1:nchan, 1, ncols, ncoefs)
    basisnames_rep = repeat(["mass-univariate"], nchan, ncols, ncoefs)
    #
    results =
        make_long_df(uf, coefs_rep, chan_rep, colnames_basis_rep, basisnames_rep, :time)

    return results
end

#---
# Returns a long df given the already matched
function make_long_df(m, coefs, chans, colnames, basisnames, collabel)
    @assert all(size(coefs) .== size(chans)) "coefs, chans and colnames need to have the same size at this point, $(size(coefs)),$(size(chans)),$(size(colnames)), should be $(size(coef(m))))"
    @assert all(size(coefs) .== size(colnames)) "coefs, chans and colnames need to have the same size at this point"
    estimate, stderror, group = make_estimate(m)

    return DataFrame(
        Dict(
            :coefname => String.(linearize(coefs)),
            :channel => linearize(chans),
            :basisname => linearize(basisnames),
            collabel => linearize(colnames),
            :estimate => linearize(estimate),
            :stderror => linearize(stderror),
            :group => linearize(group),
        ),
    )
end
#---------

# extracts betas (and sigma's for mixed models) with string grouping indicator
# returns as a ch x beta, or ch x time x beta (for mass univariate)
function make_estimate(
    m::Union{UnfoldLinearMixedModel,UnfoldLinearMixedModelContinuousTime},
)
    estimate = cat(coef(m), ranef(m), dims = ndims(coef(m)))
    if ndims(coef(m)) == 3
        group_f = repeat(
            [nothing],
            size(coef(m), 1),
            size(coef(m), 2),
            size(coef(m), ndims(coef(m))),
        )

        ranef_group = [x.group for x in MixedModels.tidyσs(modelfit(m))]

        # reshape to pred x time x chan and then invert to chan x time x pred
        ranef_group = permutedims(
            reshape(ranef_group, :, size(coef(m), 2), size(coef(m), 1)),
            [3 2 1],
        )

        group_s = ranef_group
        stderror_fixef = Unfold.stderror(m)
        stderror_ranef = fill(nothing, size(ranef(m)))
        stderror = cat(stderror_fixef, stderror_ranef, dims = 3)
    else
        group_f = repeat([nothing], size(coef(m), 1), size(coef(m), ndims(coef(m))))
        group_s = repeat(["ranef"], size(coef(m), 1), size(ranef(m), ndims(coef(m))))
        stderror = fill(nothing, size(estimate))
    end
    group = cat(group_f, group_s, dims = ndims(coef(m)))
    return Float64.(estimate), stderror, group
end

stderror(m::UnfoldModel) = stderror(modelfit(m))
function stderror(m::Union{UnfoldLinearMixedModel,UnfoldLinearMixedModelContinuousTime})
    return permutedims(
        reshape(vcat([[b.se...] for b in modelfit(m).fits]...), reverse(size(coef(m)))),
        [3, 2, 1],
    )
end
function stderror(m::LinearModelFit)
    if isempty(m.standarderror)
        stderror = fill(nothing, size(coef(m)))
    else
        stderror = Float64.(m.standarderror)
    end
end

function make_estimate(uf::UnfoldModel)
    return Float64.(coef(uf)), stderror(uf), fill(nothing, size(coef(uf)))
end

# Return the column names of the basis functions.
function get_colnames_basis(formula::FormulaTerm)
    return get_colnames_basis(formula.rhs)
end

function get_colnames_basis(rhs::Tuple)
    return colnames(rhs[1].basisfunction)
end

function get_colnames_basis(rhs::AbstractTerm)
    return colnames(rhs.basisfunction)
end

function get_basis_name(m::UnfoldModel)
    return extract_coef_info(Unfold.get_coefnames(m), 1)
end
function get_basis_name(rhs::AbstractTerm)
    return name(rhs.basisfunction)
end

function get_colnames_basis(formulas::AbstractArray{<:FormulaTerm})
    # In case of multiple basisfunction we can have an array of formulas.
    # in that case we have to add an unique identifier
    colnames_all = []
    for formula in formulas
        append!(colnames_all, get_colnames_basis(formula.rhs))
    end
    #get_colnames_basis_name(formula.rhs) .*
    return colnames_all
end


text:
path = "docs/src/assets/"

using ISA, Plots
ψ₀ = AMFMcomp(t->exp(-t^2),t->200.0,0.0)
ψ₁ = AMFMcomp(t->1.0,t->100*t,0.1)
ψ₂ = AMFMcomp(t->1+0.8*cos(2t),t->100 + 70.5*sin(t),π)
fs = 16_000
t = 0:1/fs:1
𝚿₀ = numComp( ψ₀(t), fs )
𝚿₁ = numComp( ψ₁(t), fs )
𝚿₂ = numComp( ψ₂(t), fs )
𝐂₀ = AMFMdemod(𝚿₀)
𝐂₁ = AMFMdemod(𝚿₁)
𝐂₂ = AMFMdemod(𝚿₂)
𝐒 = numSet([𝐂₀,𝐂₁,𝐂₂])
plot(𝐒)
png(path*"IS_num_spectra1.png")

using ISA, Plots
ψ₀ = AMFMcomp(t->exp(-t^2),t->20.0,0.0)
ψ₁ = AMFMcomp(t->1.0,t->100*t,0.1)
ψ₂ = AMFMcomp(t->t,t->10 + 7.5*sin(t),π)
fs = 16_000
t = 0:1/fs:2
𝚿₀ = numComp( ψ₀(t), fs )
𝚿₁ = numComp( ψ₁(t), fs )
𝚿₂ = numComp( ψ₂(t), fs )
𝐂₀ = AMFMdemod(𝚿₀)
𝐂₁ = AMFMdemod(𝚿₁)
𝐂₂ = AMFMdemod(𝚿₂)
𝐒 = numSet([𝐂₀,𝐂₁,𝐂₂])
plot(𝐒,view="TF")
png(path*"IS_num_spectra2.png")


text:
using JuliaWeBWorK
using  Roots

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

meta=(AuthorText="JuliaWeBWorK",
      Keywords="numeric derivatives",
      Question="1"
      )

qs = JuliaWeBWorK.QUESTIONS()
letters = JuliaWeBWorK.letters()

intro = raw"""

**Derivatives**

"""

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

numericq(jmt"""

### Problem $(letters())

Calculate the slope of the secant line of \( f(x)= {{:a1}} x^2 + {{:a2}} \)
between \( ({{:a3}},f({{:a3}})) \) and \( ({{:a4}},f({{:a4}})) \).
""",
         (a₂,a₀,a,b) -> begin
         fn = x  ->  a₂*x^2 + a₀
         m  = (fn(a)-fn(b))/(a-b)
         m
         end,
         (2:4, 6:8, 2:4, 5:9)
         ) |> qs

## ----

numericq(raw"""

###  Problem $(letters())

For the function \( f(x)= {{:a1}} x^2 + {{:a2}} \) between
\( ({{:a3}},f({{:a3}})) \) and \( ({{:a4}},f({{:a4}})) \),
plot the function and the secant line.

Estimate from the graph the largest distance between the two
functions between {{:a3}} and {{:a4}}.
""",
         (a₂,a₀,a,b) ->
         fn =  x  ->  a₂*x^2 + a₀
         ## could work harde rhere
         m  = (fn(a)-fn(b))/(a-b)
         M = maximum(abs, [fn(x) - (a +  m*(x-a)) for  x ∈ range(min(a,b),  max(a,b), length=100)])
         M
         end,
         (2:4, 6:8, 2:4, 5:9),
         tolerance=5.0
         ) |> qs

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

p  =  Page(intro, qs; meta...)