text:
using ProgressMeter
using FunctionalCollections
using FluxionArrays
using FluxionFields

immutable WienerNoise <: FieldExpr
    id :: Symbol
    dimensions :: Array{Dimension, 1}

    WienerNoise(dimensions) = new(gensym("WN"), dimensions)
end

immutable IntegrationProblem
    unknown_field :: UnknownField
    propagation_dimension :: UnknownDimension
    transverse_dimensions :: Array{KnownDimension, 1}
    original_dimension_order :: Array{Dimension, 1}
    rhs :: FieldExpr
    noises :: Array{WienerNoise, 1}
    noise_parameters :: Array{UnknownField, 1}
end

function find_noises(fexpr::FieldExpr)
    noises = inspect(PersistentSet(), fexpr) do controls, state, fexpr
        isa(fexpr, WienerNoise) ? conj(state, fexpr) : state
    end
    collect(noises)
end

function normalize_equation(eq::Equation)

    # left-hand side
    lhs = eq.left
    @assert isa(lhs, Diff)
    @assert length(lhs.dimensions) == 1

    ufield = lhs.field
    pdimension = collect(keys(lhs.dimensions))[1]
    tdimensions = known_dimensions(ufield)
    ufield_snapshot = remove_dimensions(ufield, pdimension)

    # right-hand side
    rhs = eq.right
    noises = find_noises(rhs)
    noise_parameters = []
    for noise in noises
        noise_parameter = UnknownField(noise.dimensions)
        noise_parameter = remove_dimensions(noise_parameter, pdimension)
        push!(noise_parameters, noise_parameter)
        rhs = replace_expr(rhs, noise => noise_parameter)
    end

    rhs = replace_expr(rhs, ufield => ufield_snapshot)

    IntegrationProblem(
        ufield_snapshot, pdimension, tdimensions, ufield.dimensions, rhs,
        noises, noise_parameters)
end

immutable SampleNoises
    afunc
end

function get_sample_noises(backend, noise_parameters)
    rs = unknown_random_state()
    step = unknown_array(shape=[], eltype=Float64)

    rs_new = rs
    returns = []
    for noise in noise_parameters
        shape = [length(vertices(FluxionFields.grid(dim))) for dim in noise.dimensions]
        norm = 1 / sqrt(prod(grid_step(FluxionFields.grid(dim)) for dim in noise.dimensions))
        noise_arr, rs_new = FluxionArrays.random_normal(rs_new, 0, norm, shape...)
        push!(returns, noise_arr ./ sqrt(step))
    end

    SampleNoises(compile(backend, array_function([rs, returns...], [rs, step])))
end

function (sn::SampleNoises)(rs, step)
    results = sn.afunc(rs, step)
    rs = results[1]
    noises = results[2:end]
    rs, noises
end

function normalize_initial_field(initial_field, trajectories)

    @assert isa(initial_field, Number) || isa(initial_field, FieldExpr)

    if isa(initial_field, Number)
        initial_field = convert(FieldExpr, initial_field)
    end

    if has_free_random_variables(initial_field)
        if trajectories === nothing
            error(
                "The initial field has random variables, " *
                "but the number of trajectories is not set")
        else
            initial_field = sample(initial_field, trajectories)
        end
    elseif is_stochastic(initial_field)
        given_trajectories = trajectories_number(initial_field)
        if given_trajectories != trajectories
            error(
                "The number of trajectories in the initial field " *
                "is not equal to the given number of trajectories")
        end
    elseif !(trajectories === nothing)
        tdims = []
        initial_field = sample(initial_field, trajectories)
    end
end

function integrate(
        eq::Equation, initial_field, initial_pdim;
        samplers=Dict(), trajectories=nothing, seed=nothing, backend=nothing)

    # In the future this will be passed from the user.
    # Now we only have one stepper.
    stepper = EulerStochasticStepper(0.02)

    if seed === nothing
        rng = MersenneTwister()
    else
        rng = MersenneTwister(seed)
    end

    if backend === nothing
        backend = default_array_backend()
    end

    get_seed(rng) = rand(rng, 1:2^32)

    problem = normalize_equation(eq)
    initial_field = normalize_initial_field(initial_field, trajectories)

    sequences = Dict((key => val[2] for (key, val) in samplers)...)
    samplers = Dict((key => val[1] for (key, val) in samplers)...)

    eseq = event_sequence(sequences)

    results = Dict(
        (key => Dict(
            :pvalue => [],
            :mean => Array{KnownField, 1}(),
            )
        for key in keys(samplers))...)

    rs = random_state(backend; seed=get_seed(rng))
    stepper_state = stepper(backend, problem, initial_field, initial_pdim, get_seed(rng))

    sample_noises = get_sample_noises(backend, problem.noise_parameters)

    # Main integration loop

    p = Progress(100, 1) # minimum update interval: 1 second

    while true
        pdim_val = pdim(stepper_state)

        to_sample, eseq = event_sequence_next_until(pdim_val, eseq)

        # sample
        for (pdim_val, keys) in to_sample

            # Returns a KnownField
            interp_field = interpolate_at(stepper, stepper_state, pdim_val)
            interp_field = known_field(
                interp_field, [problem.transverse_dimensions..., STOCHASTIC_DIMENSION])

            for key in keys
                push!(results[key][:pvalue], pdim_val)
                result = samplers[key](interp_field, pdim_val)

                outcome = Dict(:pdimension => pdim_val)

                if is_stochastic(result)
                    # For the time being we will just use the basic error estimating strategy.
                    # Ensemble splitting and such will be added later.
                    outcome[:mean] = evaluate_field(backend, stochastic_mean(result))
                    outcome[:stderr] = evaluate_field(
                        backend, stochastic_std(result) / sqrt(trajectories(result) - 1))
                else
                    result = evaluate_field(backend, result)
                    outcome[:mean] = result
                end

                results[key] = outcome
            end
        end

        if event_sequence_done(eseq)
            break
        end

        # propagate a step forward
        rs, noises = sample_noises(rs, stepper.propagation_step)
        stepper_state = make_step(stepper, stepper_state, noises)

        update!(p, convert(Integer, round(pdim_val / (eseq.last_event - initial_pdim) * 99)) + 1)
    end

    for (key, result) in results
        generic_dims = find_dimension_order(result[:mean], problem.original_dimension_order)
        pdimension = BasicDimension(
            problem.propagation_dimension.name, ArbitraryGrid(result[:pvalue]))
        snapshots = [as_field(f, generic_dims) for f in result[:mean]]
        results[key][:mean] = join_fields(snapshots, pdimension)
        if haskey(results[key], :stderr)
            results[key][:stderr] = join_fields(snapshots, pdimension)
        end
        results[key][:pdimension] = pdimension
    end

    results
end


text:
# clark_west.jl
# Clark-West
#
# Copyright (C) 2020   Guilherme Bodin, Paul Söderlind
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

export ClarkWestTest

struct ClarkWestTest <: ZTest
    n::Int       # number of observations
    xbar::Real   # estimated mean
    stderr::Real # population standard error
    z::Real      # t-statistic
    μ0::Real     # mean under h_0
end

"""
    ClarkWestTest(e1::AbstractVector{<:Real}, e2::AbstractVector{<:Real};lookahead::Integer=1)

Perform the Clark-West test of equal performance of two nested predition models, in terms of the
out-of-sample mean squared prediction errors.

`e1` is a vector of forecasts from the smaller (nested) model, `e2` is a vector of forecast
errors from the larger model, and `lookahead` is the number of steps ahead of the forecast.
Typically, the null hypothesis is that the two models perform equally well (a two-sided test),
but sometimes we test whether the larger model performs better, which is indicated by a
positive test statistic, for instance, above 1.645 for the 5% significance level (right tail test).

Implements: [`pvalue`](@ref)
# References
 * Clark, T. E., West, K. D. 2006, Using out-of-sample mean squared prediction errors to test
   the martingale difference hypothesis. Journal of Econometrics, 135(1): 155–186.
 * Clark, T. E., West, K. D. 2007, Approximately normal tests for equal predictive accuracy
   in nested models. Journal of Econometrics, 138(1): 291–311.

"""
function ClarkWestTest(e1::AbstractVector{<:Real}, e2::AbstractVector{<:Real};
                       lookahead::Integer=1)
    length(e1) == length(e2) || throw(DimensionMismatch("inputs must have the same length"))
    n            = length(e1)
    d            = 2*e1.*(e1 - e2)
    cw_cov       = HypothesisTests.autocov(d, collect(0:lookahead-1))
    cw_var       = (cw_cov[1] + 2 * sum(cw_cov[2:end]))/n
    statistic_cw = mean(d)/sqrt(cw_var)
    return ClarkWestTest(n,mean(d),sqrt(cw_var),statistic_cw,0.0)
end

testname(::ClarkWestTest) = "Clark West test"
population_param_of_interest(x::ClarkWestTest) = ("Mean",x.μ0, x.xbar)
default_tail(test::ClarkWestTest) = :both

function show_params(io::IO, x::ClarkWestTest, ident)
    println(io, ident, "number of observations:    $(x.n)")
    println(io, ident, "CW statistic:              $(x.z)")
    println(io, ident, "population standard error: $(x.stderr)")
end

ext: jl


text:
module Khepri
using Base.Iterators
using Sockets
using Dates
using ColorTypes
using Colors
using StaticArrays
using LinearAlgebra
using IntervalSets
using Dierckx
#using PyCall
using MeshCat
#using PlotlyJS

import Base.show, Base.zero, Base.iterate
import Base.convert
import Base.getindex, Base.firstindex, Base.lastindex, Base.broadcastable
import Base.union

include("Parameters.jl")
include("Utils.jl")
include("Coords.jl")
include("Paths.jl")
include("Shapes.jl")
include("BIM.jl")
include("Geometry.jl")
include("Encoders.jl")
include("Primitives.jl")
#include("Tikz.jl")
include("AutoCAD.jl")
include("Rhinoceros.jl")
#include("Revit.jl")
#include("ArchiCAD.jl")
include("Radiance.jl")
include("POVRay.jl")
#include("Notebook.jl")
#include("MeshCat.jl")
include("Grasshopper.jl")
include("Unity.jl")
#include("Unreal.jl")
include("Camera.jl")
#include("Frame3DD.jl")
#const com = PyNULL()
#include("Robot.jl")
include("BIMFamilies.jl")
end


text:
using TcrCore, FactCheck

include("test_agents.jl")
include("test_select_agent.jl")
include("test_simulations.jl")
include("test_vote.jl")
include("test_redistribute.jl")

max_stars_repo_name: mlodykaz/tcr-simulation

max_issues_repo_name: mlodykaz/tcr-simulation

max_forks_repo_name: mlodykaz/tcr-simulation


text:
module TestExamples1DEulerMulti

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_1d_dgsem")

@testset "Compressible Euler Multicomponent" begin

  @trixi_testset "elixir_eulermulti_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_ec.jl"),
      l2   = [1.54912508e-01, 4.45233361e-01, 1.70123765e-02, 3.40247530e-02, 6.80495060e-02],
      linf = [2.70251235e-01, 9.94752402e-01, 3.93233778e-02, 7.86467556e-02, 1.57293511e-01])
  end

  @trixi_testset "elixir_eulermulti_es.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_es.jl"),
      l2   = [1.53378791e-01, 4.41560792e-01, 3.93639477e-02, 7.87278954e-02],
      linf = [2.50203604e-01, 7.19643482e-01, 6.36541480e-02, 1.27308296e-01])
  end

  @trixi_testset "elixir_eulermulti_convergence_ec.jl" begin
  @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_ec.jl"),
      l2   = [8.54278172e-05, 1.63112465e-04, 1.87964310e-05, 3.75928621e-05],
      linf = [3.13612234e-04, 5.82916495e-04, 7.05047326e-05, 1.41009465e-04])
  end

  @trixi_testset "elixir_eulermulti_convergence_es.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"),
      l2   = [1.88845048e-05, 5.49106005e-05, 9.42673716e-07, 1.88534743e-06, 3.77069486e-06, 7.54138973e-06],
      linf = [1.16223512e-04, 3.07922197e-04, 3.21774233e-06, 6.43548465e-06, 1.28709693e-05, 2.57419386e-05])
  end

  @trixi_testset "elixir_eulermulti_two_interacting_blast_waves.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_two_interacting_blast_waves.jl"),
      l2   = [1.28874968e+00, 8.28362140e+01, 3.52901633e-03, 1.37049633e-02, 1.91860960e-02],
      linf = [2.96873461e+01, 1.32285301e+03, 9.04336719e-02, 3.10669595e-01, 4.41422167e-01],
      tspan = (0.0, 0.0001))
  end

end

end # module


text:
#!/usr/bin/env julia6
using ReadGlobal

function angle!(w₁::A,w₂::A,w₃::A,f₁::A,f₂::A,f₃::A,out::A) where {A<:Vector{<:Number}}
    Threads.@threads for i in 1:length(w₁)
      @fastmath @inbounds out[i] = (w₁[i]*f₁[i]+w₂[i]*f₂[i]+w₃[i]*f₃[i])/sqrt((w₁[i]^2 + w₂[i]^2 + w₃[i]^2)*(f₁[i]^2 + f₂[i]^2 + f₃[i]^2))
    end
  return 0
end

function angle(u1::String,u2::String,u3::String,w1::String,w2::String,w3::String,out::String)
    doinchunks(angle!,input=(u1,u2,u3,w1,w2,w3),output=(out,))
end

angle(ARGS...)


text:
x, y = generators(MPoly{Float64}, :x, :y)

p = (x+y)^2
@test_approx_eq evaluate(p, 1, 2) 3^2
@test evaluate_basis(p, 1, 2) == [1.0, 4.0, 4.0]

p = (x+2y)^3
@test_approx_eq evaluate(p, 0, 2) 2^3
@test evaluate_basis(p, 0, 2) == [1.0, 4.0, 4.0]


text:
using Protos
using Protos.Parsing
using Test

@testset "Protos.jl" begin
    @testset "parsing" begin
        testfile = joinpath(@__DIR__, "test.proto")
        io = open(testfile)
        parsed = parse_proto(io)
        @show parsed
        @test parsed.comment == "\nThis is a file-wide comment.\n"
        @test parsed.statements[1].package == "some.package.name"
        @test parsed.statements[2].value == "some.other.package.name"

        # TODO: ensure both messages and enum are present
        @test length(parsed.statements) > 5
    end

    @testset "utils" begin
        import Protos: proto_default
        proto_default(UInt32) == zero(UInt32)
        proto_default(String) == ""
        proto_default(Bool) == false
        proto_default(Vector{String}) == String[]
        proto_default(ProtoFile) == missing
    end

    @testset "serialization" begin
        import Protos.Serialization: writefield
        io = IOBuffer()
        idx = 1
        writefield(io, idx, 150, Val(:uint32))
        @test take!(io) == [0x08, 0x96, 0x01]
    end
end


text:
module DtPlot

#using Winston
using PyPlot
using Winston
using PosixCalendar
using Predictors
using DeepThought.DTPredictors
using DeepThought.DTPDFS
using Verification
using DataFrames
using Datetime
import DeepThought.DTPDFS.PDF

export prepareDataFramePdf, prepareDataFrameObs, leaf

function prepareDataFramePdf(F, period)
    fx = asdataframe(F, period)
    fx[:x1] = [x.seconds/3600 for x = fx[:basetime]]
    fx[:x2] = [x.seconds/3600 for x = fx[:prognosis]]
    fx[:x3] = [x.seconds/3600 for x = fx[:basetime] + fx[:prognosis]]
    fx[:x4] = [datetime(PosixCalendar.ymdhms(x)...) for x = fx[:basetime] + fx[:prognosis]]
    fx[:x5] = [string(datetime(PosixCalendar.ymdhms(x)...)) for x = fx[:basetime] + fx[:prognosis]]

    fx[:median] = [median(x) for x = fx[:value]]
    fx[:mean] = [mean(x) for x = fx[:value]]
    fx[:q75] = [quantile(x, 0.75) for x = fx[:value]]
    fx[:q25] = [quantile(x, 0.25) for x = fx[:value]]
    #fx[:y] = dropna(fx[:y])
    sort!(fx,cols=:x2)
end

function prepareDataFrameObs(F, period)
    fx = asdataframe(F, period)
    fx[:x1] = [x.seconds/3600 for x = fx[:basetime]]
    fx[:x2] = [x.seconds/3600 for x = fx[:prognosis]]
    fx[:x3] = [x.seconds/3600 for x = fx[:basetime] + fx[:prognosis]]
    fx[:x4] = [datetime(PosixCalendar.ymdhms(x)...) for x = fx[:basetime] + fx[:prognosis]]
    fx[:x5] = [string(datetime(PosixCalendar.ymdhms(x)...)) for x = fx[:basetime] + fx[:prognosis]]

    fx[:y] = fx[:value]
    sort!(fx,cols=:x2)
end

function leaf(fx, title)

    progstep = 1
    bulge = progstep * 0.3

    #fp = FramedPlot(title = "candidate")
    fp = FramedPlot(title = title)
    for fc in groupby(fx, :basetime)

        #add(fp, Points(fc[:x3], fc[:y]))
        #add(fp, Curve(fc[:x3], fc[:y]))
        for fpp in groupby(fc, :prognosis)
            bt = fpp[:x1][1]
            vt = fpp[:x3][1]
            pdfi = fpp[:value][1]
            p = [[0.01:0.1:0.26]]
            y1 = [quantile(pdfi, p) for p = [.005:.005:.995]]
            y2 = [quantile(pdfi, p) for p = [.25:.005:.75]]
            x1 = [pdf(pdfi,yi)  for yi = y1] /pdf(pdfi, quantile(pdfi, 0.50)) * bulge
            x2 = [pdf(pdfi,yi)  for yi = y2] /pdf(pdfi, quantile(pdfi, 0.50)) * bulge
            #add(fp, Curve(vt + x, y, color="red"))
            #add(fp, Curve(vt - x, y))

            add(fp, FillBetween(vt .- x1, y1,
                                vt .+ x1, y1; color="orange", alpha=0.1))
            add(fp, FillBetween(vt .- x2, y2,
                                vt .+ x2, y2; color="blue", alpha=0.1))
            # add(fp, FillBetween(vt - x1, y1,
            #                     vt + x1, y1; color="#E0E20C", alpha=0.1))
            # add(fp, FillBetween(vt - x2, y2,
            #                     vt + x2, y2; color="#00274C", alpha=0.1))
            ## add(fp, FillBetween(vt - x1, y1,
            ##                     vt + x1, y1; color="#00274C", alpha=0.1))
            ## add(fp, FillBetween(vt - x2, y2,
            ##                     vt + x2, y2; color="#E0E20C", alpha=0.1))

            #plot(vt+x,y,"r-", color = gca().lines[-1].get_color())
            #plot(vt-x,y,"r-", color = gca().lines[-1].get_color())
        end
        b = Points(fc[:x3], fc[:median])
        #setattr(b, label="b points")
        style(b, kind="filled circle", size=0.5, color="white")
        #style(b, size=0.5)
        add(fp, b)
        add(fp, Curve(fc[:x3], fc[:median]))

        epoch_ticklabels = getattr(fp.x1, "ticklabels")
        epoch_ticks = getattr(fp.x1, "ticks")
        println("xxx ", epoch_ticks)
        println("xxx ", epoch_ticklabels)

        nticks = 4
        date_ticklabels = String[]
        [push!(date_ticklabels, PosixCalendar.iso(i)) for i in fc[:basetime] + fc[:prognosis] ]
        date_idxs = map(int, linspace(1, length(date_ticklabels), nticks))
        date_ticklabels2 = String[]
        println([date_ticklabels[i] for i in date_idxs ]  )
        [push!(date_ticklabels2, replace(date_ticklabels[i], "T", "\n  "))
         for i in date_idxs ]
        #date_ticklabels2 = ["aaaa\nbbb","aaaa\nbbb","aaaa\nbbb","aaaa\nbbb"]

        ## println("yyy ", epoch_ticklabels)
        ## println("yyy ", date_ticklabels)
        ## ticklabels = String[]

        ## [push!(ticklabels, i) for i = ["a","b"]]
        setattr(fp.x1, "ticks", length(date_ticklabels2))
        setattr(fp.x1, "ticklabels", date_ticklabels2)
        # setattr(fp.x1, "ticks_style", angle=90)
        #style(fp, angle=90)
    end
    return fp
end

## function plot{T}(p::Predictor{T},period::Request)
##     BT = Time[]
##     PP = Duration[]
##     value = T[]
##     hint(p, period)
##     d = getdata(p, period)
##     for bt in keys(d)
##         for pp in keys(d[bt])
##             push!(BT,bt)
##             push!(PP,pp)
##             push!(value,d[bt][pp])
##         end
##     end

## ## end
## figure()

## #-------------------------------------------------------------------------------
## figure()
## plot(vts, vals, "k-")
## F = PDFForecast(REPOSITORY,STATION,"TTTTT","polishing_TTTTTNZBlend_TTTTT")
## hint(F, period)
## fx  = getdata(F, period)
## vts = None
## pps = None
## vals = None
## bts = sort(collect(keys(fx)))
## for bt = bts
##     pps = sort(collect(keys(fx[bt])))
##     vts = sort([bt + pp for pp = pps])
##     vals = [mean(fx[bt][pp]) for pp = pps]
##     pps = [pp.seconds for pp in pps]
##     vts = [vt.seconds for vt in vts]
##     #plot(vts, vals)
## end

# Need to add obs to dataframe

# figure()
# for fc in groupby(fx, :basetime)
#     plot(fc[:x1], fc[:y],"o")
# end

# figure()
# for fc in groupby(fx, :basetime)
#     plot(fc[:x2], fc[:y],"o")
# end

## #figure()
## title(candidate)
## for fc in groupby(fx, :basetime)
##     plot(fc[:x3], fc[:y],"o-")
##     for fpp in groupby(fc, :prognosis)
##         bt = fpp[:x1][1]
##         vt = fpp[:x3][1]
##         pdfi = fpp[:value][1]
##         p = [[0.01:0.1:0.26]]
##         y = [quantile(pdfi, p) for p = [.005:.005:.995]]
##         x = [pdf(pdfi,yi)  for yi = y] /pdf(pdfi, quantile(pdfi, 0.50)) * 0.5 * 1 #12
##         plot(vt+x,y,"r-")
##         plot(vt-x,y,"r-")
##         #plot(vt+x,y,"r-", color = gca().lines[-1].get_color())
##         #plot(vt-x,y,"r-", color = gca().lines[-1].get_color())
##     end
## end

## title(candidate)
## for fc in groupby(fx, :basetime)
##     #plot(fc[:x3], fc[:y],"o-")
##     #plot(fc[:x3], fc[:q25],"r-")
##     #plot(fc[:x3], fc[:q75],"ro")
##     yerr = zeros((2,length(fc[:q25])))
##     yerr[1,:] = fc[:y] - fc[:q25]
##     yerr[2,:] = fc[:q75] - fc[:y]
##     errorbar(fc[:x3], fc[:y], yerr=yerr,
##              fmt="o", ecolor='g', capthick=2)
## end

## figure()
## for fc in groupby(fx, :basetime)
##     plot(fc[:x3], fc[:y],"o-")
##     ## errorbar(fc[:x3], fc[:y], yerr=[fc[:q25], fc[:q75]],
##     ##          fmt="o", ecolor='g', capthick=2)
## end

# figure()
# for fc in groupby(fx, :basetime)
#     #plot(fc[:x3], fc[:y],"o-")
#     errorbar(fc[:x3], fc[:y], yerr=[fc[:q25], fc[:q75]],
#              fmt="o", ecolor='g', capthick=2)
# end

# p = FramedPlot(
#          title="title!",
#          xlabel="\\Sigma x^2_i",
#          ylabel="\\Theta_i")

# for fc in groupby(fx, :basetime)
#     x = Float64[]
#     y = Float64[]
#     yp = Float64[]
#     ym = Float64[]
#     [push!(x, i)  for i = fc[:x3]]
#     [push!(y, i)  for i = fc[:y]]
#     [push!(yp, i) for i = fc[:q25]]
#     [push!(ym, i) for i = fc[:q75]]
#     # add(p, FillBetween(x, ym, x, yp))
#     # add(p, Curve(x, y))
#     # add(p, Points(x, y))
#     # savefig("fig.png")
#     add(p,Winston.ErrorBarsY(x,ym,yp,color="green"))
# end
# savefig("fig.png")
# # function basicplot(df) # This should be *much* more sophisticated!
# #     p=Winston.FramedPlot(xlabel="$(names(df)[1]) (hours)")
# #     t=map((x)->x.seconds/3600,df[1])
# #     add(p,Winston.ErrorBarsY(t,df[:q05],df[:q95],color="green"))
# #     add(p,Winston.ErrorBarsY(t,df[:q25],df[:q75],color="blue"))
# #     add(p,Winston.ErrorBarsY(t,df[:q50],df[:q50],color="red"))
# #     return p
# # end

# figure()
# for fc in groupby(fx, :basetime)
#     x = Float64[]
#     y = Float64[]
#     [push!(x, xi) for xi = fc[:x3]]
#     [push!(y, yi) for yi = fc[:y]]
#     oplot(x, y,"o-")
# end
# savefig("figure.png")

# figure()
# for fc in groupby(fx, :basetime)
#    plot_date(fc[:x4], fc[:y],"o")
# end

# for fc in groupby(fx, :basetime)
#     p = Gadfly.plot(fc ,x="x5",y="y",Geom.line)
#     #draw(D3("mammals.js", 6inch, 6inch), p)
#     draw(SVG("myplot.svg", 12inch, 6inch), p)
# end

# p = plot(fx, x="x2", y="y",color="basetime", Geom.line)
# draw(SVG("myplot2.svg", 12inch, 6inch), p)

## #-------------------------------------------------------------------------------
## G = PDFForecast(REPOSITORY,STATION,"TTTTT","revising_TTTTTNZBlend_TTTTT")
## hint(G, period)
## gx  = getdata(G, period)
## vts = None
## pps = None
## vals = None
## bts = sort(collect(keys(obs)))
## for bt = bts
##     pps = sort(collect(keys(gx[bt])))
##     vts = sort([bt + pp for pp = pps])
##     vals = [gx[bt][pp] for pp = pps]
##     pps = [pp.seconds for pp in pps]
##     vts = [vt.seconds for vt in vts]
##     plot(vts, vals)
## end

## function fromTimetoSeconds(t::Time)
##     t.seconds
## end

## function fromTimetoSeconds(t::Time[])
##     return [ti.seconds for ti = t]
## end

## function toString(t::Time)
## end

#req = Request((Time(2014,5,1), Time(2014,6,1),), (0,dayinsecs*15))
#hint(OBS,req)
#data=Predictors.Unary{pointtype}(data,pointfunc)
#data=getdata(OBS,req)
end # module


text:
"Describes whether we're serializing (write) or deserializing (read)."
@enum ConverterModes CM_Write CM_Read

println("#TODO: Handle `missing` like `nothing`")
println("#TODO: Handle Sets")
println("#TODO: Handle rationals (as string or as Float64?)")
println("#TODO: Mark specific struct fields as being serialized/deserialized with their literal type, not their declared type, to handle massive unions like `Contiguous{T}`")
println("#TODO: support 'kwargs' thing from StructTypes")

"""
Serializer/deserializer settings.
This type is a functor, converting data into TOML-friendly primitives, or vice-versa.

The functor signature for serializing is `my_converter(x[, field_type_of_x])::TomlType`.
The functor signature for deserializing is `my_converter(x::TomlType, T)::T`.

The settings are:
 * `{ConverterMode}` : The type parameter determines whether to
        serialize (`CM_Write`) or deserialize (`CM_Read`).
 * `null_string` (default is `"null"`) : allows for reading and writing `nothing` values
        using a special string to represent them. Disable by setting the param to `nothing`.
        Note that null values can also be represented in some circumstances
            by simply excluding the field.
 * `write_null_fields` (default is false) : if true, then whenever a struct's field
        has the value `nothing` and thus could have been omitted from serialization,
        this writer will still explicitly write it as the `null_string` value.
        Enabling this feature and disabling `null_string` will throw an error.
 * `default_struct_type` (default is `Struct()`) : if a struct type hasn't been assigned a
        `StructTypes.StructType`, this one is given to it.
        Set to `nothing` to disable this feature, throwing an error when
        a struct doesn't have a `StructType`.
 * `enums_can_be_ints` (default is true) : if true, enums can be read from integer values.
        Otherwise, they have the default behavior for StructTypes, only parsable from strings.
 * `bools_can_be_strings` (default is true) : if true, bools can be read from strings like "T" or "false".
        Note that bools can always be serialized *as* strings, regardless of this setting.
 * `numbers_can_be_strings` (default is true) : if true, numbers can be parsed from strings,
        using `Base.parse(T, x)`. This also enables serialization of `NaN`, `Inf`, and `-Inf`,
        which normally aren't allowed in TOML.
"""
Base.@kwdef struct Converter{ConverterMode}
    null_string::Union{Nothing, AbstractString} = "null"
    write_null_fields::Bool = false
    default_struct_type::Union{Nothing, StructTypes.StructType} = StructTypes.UnorderedStruct()
    enums_can_be_ints::Bool = true
    bools_can_be_strings::Bool = true
    numbers_can_be_strings::Bool = true
end

#####################
#     Interface     #
#####################

"
For serialization of a `StructTypes.ArrayType`, this provides the serialized type of each element.
For deserialization, follow the `StructTypes.ArrayType` documentation.

The first three parameters are the collection, index, and specific element being serialized.
The last parameter is the desired output type; in many cases it's equal to `typeof(data)`,
    but the user might want to, for example, serialize a list of enums as a `Vector{Int}`.
This type should be prioritized over the data's actual type.

By default, this will usually return `eltype(desired_type)`,
    or for tuples, the desired type of the specific element.
If that type is `Any`, then it will try something similar with the collection's actual type.
If that type is also `Any`, then it will return the exact type of the element.
"
function array_el_type(data, index::Int, element, desired_type)
    if (desired_type <: Tuple) && (desired_type.parameters[index] != Any)
        return desired_type.parameters[index]
    elseif eltype(desired_type) != Any
        return eltype(desired_type)
    elseif (data isa Tuple) && (typeof(data).parameters[index] != Any)
        return typeof(data).parameters[index]
    elseif eltype(data) != Any
        return eltype(data)
    else
        return typeof(element)
    end
end

"
For serialization of a `StructTypes.DictType`,
    this provides the type of its elements.

The first three parameters are the dictionary, key, and specific value being serialized.
The last parameter is the desired output type;
    in many cases it's equal to `typeof(data)`, but the user might want to, for example,
    serialize a dictionary of enum values as a dictionary of Int values.
The desired type should be prioritized over the data's actual type.

Default behavior is to check `eltype(desired_type)`, expecting to get some `Pair{K, V}`,
    and return that `V`.
If the element type isn't a `Pair`, or the element type is `Any`,
    then the actual dictionary's value-type is used instead.
If that also fails, then the value's exact type will be returned.
"
function dict_value_type(dict, key, value, desired_type)
    if (eltype(desired_type) <: Pair) && (eltype(desired_type).types[2] != Any)
        return eltype(desired_type).types[2]
    elseif (eltype(dict) <: Pair) && (eltype(dict).types[2] != Any)
        return eltype(dict).types[2]
    else
        return typeof(value)
    end
end

###################
#     Writing     #
###################

# If the type is already TOML-friendly, pass it through.
# Special exceptions are the special float values "NaN", "+Inf", and "-Inf".
(c::Converter{CM_Write})(x::TomlType) = begin
    if x isa AbstractFloat
        if isinf(x)
            if c.numbers_can_be_strings
                return (x > 0) ? "+Inf" : "-Inf"
            else
                error("Can't convert an infinite float to TOML, unless `numbers_can_be_strings` is enabled")
            end
        elseif isnan(x)
            if c.numbers_can_be_strings
                return "NaN"
            else
                error("Cant convert a NaN to TOML, unless `numbers_can_be_strings` is enabled")
            end
        else
            return x
        end
    else
        return x
    end
end

# Get the StructType of the data and dispatch based on that.
# If it has a partcular expected type (e.x. its declared type in a struct),
#    look out for some special cases, like abstract types.
(c::Converter{CM_Write})(x) = c(x, StructTypes.StructType(x), typeof(x))
(c::Converter{CM_Write})(x, expected_type) = begin
    # Get the type of X according to its declaration
    #    (meaning the 'T' in 'x::T', or in 'Vector{T}', or in 'Dict{T, V}', etc).
    # Note that it could be a union of types.
    types::Tuple = union_types(expected_type)
    x_type_idx = findfirst(t -> x isa t, types)

    # If the input data isn't any of the declared types
    #    (e.x. it's an enum and the declared type is Int),
    #    try writing it as each type.
    if isnothing(x_type_idx)
        for T in union_parse_order(expected_type)
            try
                return c(x, StructTypes.StructType(T), expected_type)
            catch e
                # Left here in case debugging is needed.
                #@warn "Parse failed: $T" exception=(e, catch_backtrace())
            end
        end
        error("Couldn't write data of type ", typeof(x), " as a ", expected_type)
    else
        T = types[x_type_idx]
        return c(x, StructTypes.StructType(T), expected_type)
    end
end

(c::Converter{CM_Write})(x, ::StructTypes.NullType, _) =
    if isnothing(c.null_string)
        error("There is no way to write a `nothing` value, as the `null_string` parameter wasn't set")
    else
        c.null_string
    end
(c::Converter{CM_Write})(x, ::StructTypes.BoolType, _) = Bool(x)
(c::Converter{CM_Write})(x, ::StructTypes.NumberType, _) = c(StructTypes.construct(StructTypes.numbertype(typeof(x)), x))
# Special behavioir for enums: they should be written as ints, not Float64.
(c::Converter{CM_Write})(x::Enum{I}, ::StructTypes.NumberType, _) where {I} = c(StructTypes.construct(I, x))
(c::Converter{CM_Write})(x, ::StructTypes.StringType, _) = string(x)
(c::Converter{CM_Write})(x, ::StructTypes.ArrayType, desired_type) = let vec = [ ]
    for (i::Int, el) in enumerate(x)
        push!(vec, c(el, array_el_type(x, i, el, desired_type)))
    end
    vec
end
(c::Converter{CM_Write})(x, ::StructTypes.DictType, desired_type) = Dict(
    (string(c(k)) => c(v, dict_value_type(x, k, v, desired_type)))
      for (k, v) in StructTypes.keyvaluepairs(x)
)
(c::Converter{CM_Write})(x, ::StructTypes.NoStructType, _) =
    if isnothing(c.default_struct_type)
        error("No StructType defined for serializing ", typeof(x),
              " (or maybe its original field type)")
    else
        c(x, c.default_struct_type)
    end
(c::Converter{CM_Write})(x, ::StructTypes.AbstractType, expected_type) = begin
    subtype_data::NamedTuple = StructTypes.subtypes(expected_type)
    subtype_strings::ConstVector{Symbol} = keys(subtype_data)
    subtype_types::ConstVector{Type} = values(subtype_data)

    subtype_idx::Union{Int, Nothing} = findfirst(t -> (x isa t), subtype_types)

    if isnothing(subtype_idx)
        error("Object type ", typeof(x),
              " is not listed as a `StructTypes.subtypes()` of ", expected_type)
    else
        return Dict(
            "type" => subtype_strings[subtype_idx],
            "value" => c(x)
        )
    end
end
(c::Converter{CM_Write})(x, ::StructTypes.CustomStruct, _) = c(StructTypes.lower(x))
(c::Converter{CM_Write})(x, ::Union{StructTypes.Mutable, StructTypes.Struct}, _) = begin
    renamed_fields::NTuple = StructTypes.names(typeof(x))
    renamed_julia_fields::NTuple = map(t -> t[1], renamed_fields)
    renamed_output_fields::NTuple = map(t -> t[2], renamed_fields)

    fields_exclude = StructTypes.excludes(typeof(x))

    # There are a number of reasons why a field's specific value might not be written.
    fields_omitempty = StructTypes.omitempties(typeof(x))
    @inline hide_value(field_name::Symbol, field_value)::Bool = (
        (!c.write_null_fields && isnothing(field_value)) ||  # It's null, and nulls should be omitted
        (
            ((fields_omitempty === true) || (field_name in fields_omitempty)) && # Empties should be omitted, and...
            (hasmethod(iterate, tuple(typeof(field_value))) && isempty(field_value)) #...this value is empty.
        )
    )

    output = Dict{String, Any}()
    for field_name::Symbol in fieldnames(typeof(x))
        # Check if the field should be excluded.
        if !(field_name in fields_exclude)
            field_value = getproperty(x, field_name)
            # Check if the value shouldn't be written.
            if !hide_value(field_name, field_value)
                field_type = fieldtype(typeof(x), field_name)
                # Change the serialized name, if requested.
                new_name_idx = findfirst(x -> x==field_name, renamed_julia_fields)
                if !isnothing(new_name_idx)
                    field_name = renamed_output_fields[new_name_idx]
                end
                # Finally, output the field's value into the dictionary for serialization.
                output[string(field_name)] = c(field_value, field_type)
            end
        end
    end
    return output
end

###################
#     Reading     #
###################

# A catch-all for logic that doesn't fit into Julia's overload resolution order.
(c::Converter{CM_Read})(x::TIn, ::Type{TOut}) where {TIn, TOut} = begin
    # For Union types, try each element of the union individually.
    if TOut isa Union
        for T in union_parse_order(TOut)
            try_val = union_try_parse(c, x, T)
            if try_val isa Some
                @assert(try_val isa Some{<:T},
                        "WOML.union_try_parse() somehow turned a $T into a $(typeof(something(try_val))): $try_val")
                return something(try_val)
            end
        end
        error("Unable to parse a ", typeof(x),
              " into any of the following types: ", union_types(TOut),
              ". The precise value was: ", x)
    # If the output type is non-specific, try to aggressively parse the incoming data.
    elseif TOut == Any
        if x isa AbstractString
            if c.numbers_can_be_strings
                f = tryparse(Float64, x)
                if !isnothing(f)
                    return (isinteger(f) && (f <= typemax(Int)) && (f >= typemin(Int))) ?
                            Int(f) :
                            f
                end
            end

            if c.bools_can_be_strings
                if x == "true"
                    return true
                elseif x == "false"
                    return false
                end
                # Don't try the more expansive options like "f", or "True",
                #    as it feels a bit too eager.
            end

            if x == c.null_string
                return nothing
            end
        elseif x isa AbstractDict
            return c(x, Dict{Any, Any})
        elseif x isa AbstractVector
            return c(x, Vector{Any})
        end
        return x
    elseif TOut <: Union{Dict{Any, Any}, Dict{<:AbstractString, Any}}
        K = eltype(TOut).parameters[1]
        return Dict( (c(k, K) => c(v, Any)) for (k, v) in x )
    elseif TOut == Vector{Any}
        return map(e -> c(e, Any), x)
    # If the data already satisfies the desired output type, pass it through.
    elseif TIn <: TOut
        return x
    # If the desired type is a Symbol, construct it.
    elseif TOut == Symbol
        return StructTypes.construct(Symbol, x)
    end

    # If nothing is special about this case, dispatch based on the desired type's `StructType`.
    return c(x, TOut, StructTypes.StructType(TOut))
end

# Support deserializing a number from a string.
(c::Converter{CM_Read})(x::AbstractString, T::Type{<:Real}) = begin
    if !c.numbers_can_be_strings
        # Was the data supposed to be a number, or was it generic?
        if T <: Real
            error("Can't parse a ", T, " from a ", typeof(x),
                  " unless `numbers_can_be_strings` is enabled")
        else
            return x
        end
    else
        # Make the parsed type concrete.
        if T in (Real, AbstractFloat)
            T = Float64
        elseif T in (Integer, Signed)
            T = Int64
        elseif T == Unsigned
            T = UInt64
        end
        return parse(T, x)
    end
end

# Support deserializing an enum from an int.
(c::Converter{CM_Read})(x::Number, T::Type{<:Enum{I}}) where {I<:Integer} = begin
    if c.enums_can_be_ints && isinteger(x)
        return T(convert(I, x))
    else
        # Pass through to the usual behavior (almost certainly throwing an error).
        return c(x, T, StructTypes.StructType(T))
    end
end

# Support deserializing a Bool from a string.
(c::Converter{CM_Read})(x::AbstractString, ::Type{Bool}) = begin
    x = lowercase(x)
    if c.bools_can_be_strings && (x in ("t", "true"))
        true
    elseif c.bools_can_be_strings && (x in ("f", "false"))
        false
    else
        # Fall through to the usual behavior (almost certainly throwing an error).
        c(x, Bool, StructTypes.StructType(Bool))
    end
end

# Support reading the null-string.
(c::Converter{CM_Read})(x::AbstractString, ::Type{Nothing}) = begin
    if !isnothing(c.null_string) && (x == c.null_string)
        nothing
    else
        # Fall through to the usual behavior (almost certainly throwing an error).
        c(x, Nothing, StructTypes.StructType(Nothing))
    end
end

(c::Converter{CM_Read})(x, T, ::StructTypes.NullType) = error(
    "Cannot convert a ", typeof(x), " to `Nothing`"
)
(c::Converter{CM_Read})(x, T, ::StructTypes.BoolType) = StructTypes.construct(T, x)
(c::Converter{CM_Read})(x, T, ::StructTypes.NumberType) = StructTypes.construct(T, convert(StructTypes.numbertype(T), x))
(c::Converter{CM_Read})(x, T, ::StructTypes.StringType) = StructTypes.construct(T, x)
(c::Converter{CM_Read})(x::AbstractVector, T::Type{<:Tuple}, ::StructTypes.ArrayType) = begin
    if length(x) != length(T.parameters)
        error("Can't parse a ", length(T.parameters), "-tuple from an array of ",
              length(x), " elements")
    else
        return tuple((c(e, T.parameters[i]) for (i, e) in enumerate(x))...)::T
    end
end
(c::Converter{CM_Read})(x::AbstractVector, T, ::StructTypes.ArrayType) = begin
    # Process each element to be as close to the desired type as possible.
    # Even if the element type is Any, we'd like to aggresively parse things
    #    like numbers/bools as the user configured.
    TElement = (Base.IteratorEltype(T) isa Base.HasEltype) ?
                   eltype(T) :
                   Any
    x = map(i -> c(i, TElement), x)

    # Call through to the constructor if necessary.
    return (x isa T) ? x : StructTypes.construct(T, x)
end
(c::Converter{CM_Read})(x::AbstractDict, T::Type, ::StructTypes.DictType) = StructTypes.construct(T, x)
(c::Converter{CM_Read})(x::AbstractDict, T::Type{<:AbstractDict{K, V}}, ::StructTypes.DictType) where {K, V} = begin
    # Convert the incoming key type to the desired type.
    # Do the same for values.
    return T(c(k, K) => c(v, V) for (k, v) in x)
end
(c::Converter{CM_Read})(x, T::Type, ::StructTypes.NoStructType) =
    if isnothing(c.default_struct_type)
        error("No StructType defined for deserializing ", typeof(x))
    else
        c(x, T, c.default_struct_type)
    end
(c::Converter{CM_Read})(x::AbstractDict, T::Type, ::StructTypes.AbstractType) = begin
    if (length(x) != 2) || !haskey(x, "type") || !haskey(x, "value")
        error("For the AbstractType '", T, "', expected a dictionary",
              " with the fields 'type' and 'value'")
    end
    subtype_string = x["type"]
    subtype_data = x["value"]

    subtypes::NamedTuple = StructTypes.subtypes(field_type)
    subtype_strings::ConstVector{Symbol} = keys(subtypes)
    subtype_types::ConstVector{Type} = values(subtypes)

    subtype_idx::Union{Int, Nothing} = findfirst(n -> (n == subtype_string), subtype_strings)
    if isnothing(subtype_idx)
        error("Serialized type '", subtype_string,
              "' is not a known type of ", T, ": ", subtype_strings)
    end

    return c(subtype_data, subtype_types[subtype_idx])
end
(c::Converter{CM_Read})(x, T::Type, ::StructTypes.CustomStruct) = StructTypes.construct(
    T,
    # If asked to "lower" the incoming value, then do so.
    let in_type = StructTypes.lowertype(T)
        if in_type == Any
            if x isa AbstractDict
                x
            else
                error("CustomStruct type '", T, "' didn't implement lowertype(),",
                      " so it expects a dictionary for deserialization, but it got a ",
                      typeof(x))
            end
        else
            c(x, in_type)
        end
    end
)
using InteractiveUtils
(c::Converter{CM_Read})(x::AbstractDict{<:AbstractString, <:Any}, T::Type, ::Union{StructTypes.Mutable, StructTypes.Struct}) = begin
    renamed_fields::NTuple = StructTypes.names(T)
    renamed_julia_fields::NTuple = map(t -> t[1], renamed_fields)
    renamed_toml_fields::NTuple = map(t -> t[2], renamed_fields)
    function get_field_name(toml_name::Symbol)
        i = findfirst(n -> (n == toml_name), renamed_toml_fields)
        return isnothing(i) ? toml_name : renamed_julia_fields[i]
    end

    fields_exclude = StructTypes.excludes(T)

    fields = Dict{Symbol, Any}()
    for toml_name in keys(x)
        field_name = get_field_name(Symbol(toml_name))
        if !in(field_name, fields_exclude)
            fields[field_name] = c(x[toml_name], fieldtype(T, field_name))
        end
    end

    fields_to_set = setdiff(fieldnames(T), fields_exclude)
    field_values = Iterators.map(f -> get(fields, f, nothing), fields_to_set)
    return T(field_values...)
end


text:
E_from_Gν(G,ν) = 2G*(1 + ν)
λ_from_Gν(G,ν) = 2G*ν / (1 - 2ν)

Dᵉ_func(i,j,k,l,G,λ) = λ*(δ(i,j)*δ(k,l)) + G*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k))
Cᵉ_func(i,j,k,l,G,λ) = (λ/(2G*(3λ + 2G)))*(δ(i,j)*δ(k,l)) + (1/2G)*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k))

get_elastic_stiffness_tensor(G,λ) = SymmetricTensor{4, 3}( (i,j,k,l) -> Dᵉ_func(i,j,k,l,G,λ))
get_elastic_stiffness_tensor(r::Rheology) = SymmetricTensor{4, 3}( (i,j,k,l) -> Dᵉ_func(i,j,k,l,r.G,λ_from_Gν(r.G,r.ν)))
get_elastic_compliance_tensor(r::Rheology) = SymmetricTensor{4, 3}( (i,j,k,l) -> Cᵉ_func(i,j,k,l,r.G,λ_from_Gν(r.G,r.ν)))

free_energy_convexity(r::Rheology,D) = 1/compute_Γ(r,compute_A1B1(r,D)...) > 0 ? true : false
# eq 16 Bhat2012 & 2016 & notes (because in Bhat2011 c2 isn't the same form as in Harsha's notes) :
#compute_c1(d::Damage,D) = sqrt(1-cos(d.ψ)^2)/(π*cos(d.ψ)^(3/2)*((D/d.D₀)^(1/3) - 1 + d.β/cos(d.ψ))^(3/2))
function compute_c1(r,D)
  α = cosd(r.ψ)
  #(α < 0) && (α = -α) #### careful with that
  @assert α >= 0

  if r.D₀>0
    isnan(D) && error("D is NaN")
    isnan(r.D₀) && error("r.D₀ is NaN")
    ((D/r.D₀) < 1) && println("D = $D & D₀ = $(r.D₀)")
    @assert (D/r.D₀) >= 1
    return sqrt(1-α^2)/(π*α^(3/2)*((D/r.D₀)^(1/3) - 1 + r.β/α)^(3/2))
  elseif (r.D₀==0) & (D==0)
    return sqrt(1-α^2)/(π*α^(3/2)*(r.β/α)^(3/2)) # lim D->0 D/D0 = 1
  end
end
# Perol&Bhat2016 : 1/α  or  Harsha's notes : 1/α^2 ???
function compute_c2(d::Rheology,D)
  α = cosd(d.ψ)
  (d.D₀==0) && (return 0.0)
  return (sqrt(1 - α^2)/α^2) * (d.D₀^(2/3)/(1 - D^(2/3))) ###
end

function compute_c3(d::Rheology,D)
  α = cosd(d.ψ)
  #(α < 0) && (α = -α) #### careful with that
  @assert α > 0
  (d.D₀ == D) && (return 0.0)
  @assert (D/d.D₀) >= 1
  return (2sqrt(α)/π)*((D/d.D₀)^(1/3) - 1)^(1/2)
end

function compute_c1c2c3(r,D)
  c1 = compute_c1(r,D)
  c2 = compute_c2(r,D)
  c3 = compute_c3(r,D)
  return c1, c2, c3
end
# eq 15 Bhat2012 (A1 : *c2*c3), Perol&Bhat2016 (A1 : ...*c2)*c3):
# Perol&Bhat2016 is the corrected version, and the one implemented
compute_A(r::Rheology,c1,c2,c3) = r.μ*c1 + (1 + r.μ*c2)*c3
compute_B(c1,c2,c3) = c1 + c2*c3

function compute_AB(r::Rheology,c1,c2,c3)
  A = compute_A(r,c1,c2,c3)
  B = compute_B(c1,c2,c3)
  return A, B
end
compute_AB(r::Rheology,D) = compute_AB(r,compute_c1c2c3(r,D)...)

# eq 11 in Harsha's notes :
compute_A1(r::Rheology,A) = A * sqrt((π*r.D₀*(1 - r.ν))/cosd(r.ψ)^3)
compute_B1(r::Rheology,B) = B * sqrt((π*r.D₀*(1 - r.ν))/cosd(r.ψ)^3)

# Bhat 2012
# compute_A1(r::Rheology,A) = A * sqrt((π*r.D₀)/(cosd(r.ψ)^3*(1 - r.ν)))
# compute_B1(r::Rheology,B) = B * sqrt((π*r.D₀)/(cosd(r.ψ)^3*(1 - r.ν)))

function compute_A1B1(r::Rheology,A,B)
  A1 = compute_A1(r,A)
  B1 = compute_B1(r,B)
  return A1, B1
end
compute_A1B1(r::Rheology,D) = compute_A1B1(r,compute_AB(r,D)...)

function compute_a1(B1,Γ)
  return (1/Γ)*(1 + B1^2/2)
end

function compute_b1(A1,B1,Γ)
  return -(1/Γ)*((A1*B1)/2)
end

function compute_b2(r::Rheology,A1,Γ)
  return (1/Γ)*(A1^2/2 + (3*(1-2r.ν))/(2*(1+r.ν)))
end

function compute_Γ(r::Rheology,A₁,B₁)
  return (3*(1-2r.ν))/(2*(1+r.ν)) + (3*(1-2r.ν)*B₁^2)/(4*(1+r.ν)) + A₁^2/2
end

function compute_Γ(r::Rheology,D)
  c1 = compute_c1(r,D)
  c2 = compute_c2(r,D)
  c3 = compute_c3(r,D)
  A = compute_A(r,c1,c2,c3)
  B = compute_B(c1,c2,c3)
  A₁ = compute_A1(r,A)
  B₁ = compute_B1(r,B)
  ν = r.ν
  return (3*(1-2ν))/(2*(1+ν)) + (3*(1-2ν)*B₁^2)/(4*(1+ν)) + A₁^2/2
end

function compute_dc1dD(r::Rheology,D)
  α = cosd(r.ψ)
  term1 = (-sqrt(1-α^2))/(2*π*α^(3/2)*D^(2/3)*r.D₀^(1/3))
  term2 = ((D/r.D₀)^(1/3) - 1 + (r.β/α))^(-5/2)
  return term1 * term2
end

function compute_dc2dD(r::Rheology,D)
  α = cosd(r.ψ)
  term1 = (2*sqrt(1-α^2)*r.D₀^(2/3))/(3*α^2*D^(1/3))
  term2 = (1 - D^(2/3))^(-2)
  return term1 * term2
end

function compute_dc3dD(r::Rheology,D)
  α = cosd(r.ψ)
  term1 = (sqrt(α))/(3*π*D^(2/3)*r.D₀^(1/3))
  term2 = ((D/r.D₀)^(1/3) - 1)^(-1/2)

  # ensure term2 isn't Infinity
  (term2 == Inf) && (term2 = 1e9)
  return term1 * term2
end

function compute_dA1dD(r::Rheology,dc1dD,dc2dD,dc3dD,c2,c3)
  α = cosd(r.ψ)
  term1 = sqrt((π*r.D₀*(1 - r.ν))/(α^3))
  term2 = r.μ*dc1dD + dc3dD + r.μ*c2*dc3dD + r.μ*c3*dc2dD
  return term1 * term2
end

function compute_dB1dD(r::Rheology,dc1dD,dc2dD,dc3dD,c2,c3)
  α = cosd(r.ψ)
  term1 = sqrt((π*r.D₀*(1 - r.ν))/(α^3))
  term2 = dc1dD + c2*dc3dD + c3*dc2dD
  return term1 * term2
end

function compute_dΓdD(r::Rheology,A1,B1,dA1dD,dB1dD)
  return ((3*(1-2r.ν))/(2*(1+r.ν)))*B1*dB1dD + A1*dA1dD
end

function compute_dA1dt(r::Rheology,dc1dt,dc2dt,dc3dt,c2,c3)
  α = cosd(r.ψ)
  term1 = sqrt((π*r.D₀*(1 - r.ν))/(α^3))
  term2 = r.μ*dc1dt + dc3dt + r.μ*c2*dc3dt + r.μ*c3*dc2dt
  return term1 * term2
end

function compute_dB1dt(r::Rheology,dc1dt,dc2dt,dc3dt,c2,c3)
  α = cosd(r.ψ)
  term1 = sqrt((π*r.D₀*(1 - r.ν))/(α^3))
  term2 = dc1dt + c2*dc3dt + c3*dc2dt
  return term1 * term2
end

function compute_dΓdt(r::Rheology,A1,B1,dA1dt,dB1dt)
  return ((3*(1-2r.ν))/(2*(1+r.ν)))*B1*dB1dt + A1*dA1dt
end

# function compute_dΓdt(p::RockParams,A1,B1,dA1dt,dB1dt)
#     return ((3*(1-2p.ν))/(2*(1+p.ν)))*B1*dB1dt + A1*dA1dt
# end

function compute_da1dt(B1,Γ,dB1dt,dΓdt)
  return -(dΓdt/Γ^2)*(1 + B1^2/2) + (B1*dB1dt)/Γ
end

function compute_db1dt(A1,B1,Γ,dA1dt,dB1dt,dΓdt)
  return (dΓdt/Γ^2)*((A1*B1)/2) - (1/2Γ)*(dA1dt*B1 + A1*dB1dt)
end

function compute_db2dt(r::Rheology,A1,Γ,dA1dt,dΓdt)
  return -(dΓdt/Γ^2)*(A1^2/2 + (3*(1-2r.ν))/(2*(1+r.ν))) + (2/Γ)*A1*dA1dt
end

function compute_dσdt(r::Rheology,a1,b1,da1dt,db1dt,ϵ,γ,dϵdt,dγdt)
  return r.G * (da1dt*ϵ + a1*dϵdt + db1dt*γ + b1*dγdt)
end

function compute_dτdt(r::Rheology,b1,b2,db1dt,db2dt,ϵ,γ,dϵdt,dγdt)
  return r.G * (db1dt*ϵ + b1*dϵdt + db2dt*γ + b2*dγdt)
end

function get_τ(tensor::AbstractArray ; input=:s)
  (input == :s) && (return sqrt(0.5 * tensor ⊡ tensor))
  if input == :σ
    s = dev(tensor)
    return sqrt(0.5 * s ⊡ s)
  end
end

function compute_KI(r::Rheology,σ,τ,A,B)
  return (A*σ + B*τ) * sqrt(π*r.a)
end

function compute_KI(r::Rheology,σ,τ,D)
  c1, c2, c3 = compute_c1c2c3(r,D)
  A, B = compute_AB(r,c1,c2,c3)
  # println("c1 : ",c1)
  # println("c2 : ",c2)
  # println("c3 : ",c3)
  # println("σ : ",σ)
  # println("τ : ",τ)
  return (A*σ + B*τ) * sqrt(π*r.a)
end

function compute_KI(r::Rheology,σij::SymmetricTensor,D)
  A, B = compute_AB(r,D)
  p = 1/3 * tr(σij) # trial pressure, negative in compression
  sij = dev(σij) # trial deviatoric stress
  τ = get_τ(sij)
  return (A*p + B*τ) * sqrt(π*r.a)
end

function compute_dDdl(r::Rheology,D)
  return (3*D^(2/3)*r.D₀^(1/3))/(cosd(r.ψ)*r.a)
end

function compute_subcrit_damage_rate(r::Rheology, KI, D ; vmax = :l̇₀)
  ((KI <= 0) || (D >= 1)) && (return 0.0)
  ρ = 2700 ##### TODO better
  Vs = sqrt(r.G/ρ)
  Vr = Vs * (0.862 + 1.14r.ν)/(1 + r.ν)

  dDdl = compute_dDdl(r,D) # damage derivative wrt crack length

  dlmax = Inf
  if vmax == :l̇₀
    dlmax = r.l̇₀
  elseif vmax == :rayleigh
    dlmax = Vr
  else
    @error "vmax keyword argument must be `:l̇₀` or `:rayleigh`"
  end
  dldt = min(r.l̇₀*(KI/r.K₁c)^(r.n),dlmax)  #Vr cracks growth rate
  #@debug dlmax
  @assert dDdl * dldt >= 0
  return dDdl * dldt
end

### TEST ###
function compute_subcrit_damage_rate_dynamic(r::Rheology, KI, D ; vmax = :l̇₀)
  ((KI <= r.K₁c) || (D >= 1)) && (return 0.0)
  ρ = 2700 ##### TODO better
  Vs = sqrt(r.G/ρ)
  Vr = Vs * (0.862 + 1.14r.ν)/(1 + r.ν)

  dDdl = compute_dDdl(r,D) # damage derivative wrt crack length
  dldt = Vr  #Vr cracks growth rate
  #@debug dlmax
  @assert dDdl * dldt >= 0
  return dDdl * dldt
end

function compute_free_energy(r,D,ϵij)
  G = r.G
  ν = r.ν
  ϵ = tr(ϵij)
  e = dev(ϵij)
  γ = sqrt(2.0 * e ⊡ e)

  A1, B1 = compute_A1B1(r,D)
  Γ = compute_Γ(r,A1,B1)
  a1 = compute_a1(B1,Γ)
  b1 = compute_b1(A1,B1,Γ)
  b2 = compute_b2(r,A1,Γ)

  return G * (0.5*a1*ϵ^2 + 0.5*b2*γ^2 + b1*ϵ*γ)
end

function compute_σij(r,A1,B1,Γ,ϵij)
  # TODO make a visco elastic version of this function

  if all(x->x==0,ϵij)
      return Tensor{2,3}(ϵij)
  end

  G = r.G
  ν = r.ν
  Id = SymmetricTensor{2,3}(δ)
  # strain invariants
  ϵ = tr(ϵij)
  e = dev(ϵij)
  γ = sqrt(2.0 * e ⊡ e)

  # stress tensor calculation
  term1 = ( (3*(1-2ν))/(1+ν) + A1^2 - A1*B1*ϵ/γ ) * ϵij
  term2 = (3ν/(1+ν) + B1^2/2 - A1^2/3 + A1*B1*ϵ/(3γ)) * ϵ
  term3 = -A1*B1*γ/2
  return (G/Γ) * (term1 + (term2 + term3)*Id)
end
compute_σij(r::Rheology,D,ϵij) = D==0 ? compute_σij(r,0.0,0.0,compute_Γ(r,0.0,0.0),ϵij) : compute_σij(r,compute_A1B1(r,D)...,compute_Γ(r,D),ϵij)

#TODO :
function compute_ϵ_oop(r,A1,B1,σᵢⱼ)

  σoop = σᵢⱼ[3,3]
  σ = 1/3*tr(σᵢⱼ)
  sᵢⱼ = dev(σᵢⱼ)
  τ = get_τ(sᵢⱼ)
  G = r.G
  ν = r.ν

  # stress tensor calculation
  term1 = ( 1 + A1*B1*σ/(2*τ) + B1^2/2) * σoop
  term2 = (3ν/(1+ν) + A1*B1*σ/(2*τ) - A1^2/3 + B1^2/2) * σ
  term3 = A1*B1*τ/3
  return (1/(2G)) * (term1 - term2 + term3)
end
compute_ϵ_oop(r,D,σᵢⱼ)=compute_ϵ_oop(r,compute_A1B1(r,D)...,σᵢⱼ)

function compute_ϵij(r,A1,B1,σᵢⱼ)

  σ = 1/3*tr(σᵢⱼ)
  sᵢⱼ = dev(σᵢⱼ)
  τ = get_τ(sᵢⱼ)
  G = r.G
  ν = r.ν
  Id = SymmetricTensor{2,3}(δ)
  # strain invariants
  term1 = ( 1 + A1*B1*σ/(2*τ) + B1^2/2) * σᵢⱼ
  term2 = (3ν/(1+ν) + A1*B1*σ/(2*τ) - A1^2/3 + B1^2/2) * σ
  term3 = A1*B1*τ/3
  return (1/(2G)) * (term1 + (term3 - term2)*Id)
end
compute_ϵij(r,D,σᵢⱼ) = D==0 ? compute_ϵij(r,0.0,0.0,σᵢⱼ) : compute_ϵij(r,compute_A1B1(r,D)...,σᵢⱼ)

function compute_damaged_stiffness_tensor(r::Rheology,ϵij,D)

  # unpack
  G = r.G
  ν = r.ν

  # Damage constants
  A₁,B₁ = compute_A1B1(r,D)
  Γ = compute_Γ(r,A₁,B₁)

  # strain invariants
  ϵ = tr(ϵij)
  e = dev(ϵij)
  γ = sqrt(2.0 * e ⊡ e)

  @assert !isnan(G)
  @assert !isnan(A₁)
  @assert !isnan(B₁)
  @assert !isnan(Γ)
  @assert !isnan(ϵ)
  @assert !isnan(γ)
  @assert !isnan((3*(1-2ν))/(2*(1+ν)))
  @assert !isnan(A₁^2/2)
  # println("A₁*B₁*ϵ : ", A₁*B₁*ϵ)
  # println("(2*γ) : ", (2*γ))
  #@assert !isnan(A₁*B₁*ϵ/(2*γ)) # returns NaN, because 0/0

  (γ == 0) && (γ += 1e-15)#zero(typeof(ϵij))) # TODO remove redundancy
  ϵ̂ = ϵij/γ

  # get stiffness factors
  if G/Γ == 0
      Cμ = 0.0
      Cλ = 0.0
      Cσ = 0.0
      Cσσ = 0.0
  else
      Cμ = G/Γ * ( (3*(1-2ν))/(2*(1+ν)) + A₁^2/2 - A₁*B₁*ϵ/(2*γ) )
      Cλ = G/Γ * ( 3*ν/(1+ν) + B₁^2/2 - A₁^2/3 + A₁*B₁*ϵ/γ + 2A₁*B₁*ϵ^3/(9γ^3) )
      Cσ = - G/Γ * ( A₁*B₁ + 2*A₁*B₁*ϵ^2/(3*γ^2) )
      Cσσ = G/Γ * (2A₁*B₁*ϵ/γ)
  end

  if isnan(Cμ)
      println("Cμ = ", Cμ)
      @assert !isnan(G/Γ)
      println("G/Γ = ", G/Γ)
      @assert !isnan( (3*(1-2ν))/(2*(1+ν)) + A₁^2/2 - A₁*B₁*ϵ/(2*γ) )
      println("second term Cμ = ", (3*(1-2ν))/(2*(1+ν)) + A₁^2/2 - A₁*B₁*ϵ/(2*γ) )
      @assert !isnan(G/Γ * ( (3*(1-2ν))/(2*(1+ν)) + A₁^2/2 - A₁*B₁*ϵ/(2*γ) ))
  end
  if isnan(Cλ)
      println("Cλ = ", Cλ)
  end
  @assert !isnan(Cσ)
  @assert !isnan(Cσσ)
  # functional form of the stiffness tensor
  C_func(i,j,k,l) = Cμ * ( δ(k,i)*δ(l,j) + δ(l,i)*δ(k,j) ) +
                    Cλ * ( δ(i,j)*δ(k,l) ) +
                    Cσ * ( ϵ̂[i,j]*δ(k,l) + δ(i,j)*ϵ̂[k,l] ) +
                    Cσσ * (ϵ̂[i,j]*ϵ̂[k,l])

  # assemble the tensor
  return SymmetricTensor{4,3}(C_func)
end

##
function compute_ϵ̇ij_2(r,D,Ḋ,σᵢⱼ,σ̇ᵢⱼ; damaged_allowed=true, R1_allowed=false, keep_2G=true)

  # Convert in case automatic differentiation supplies a Matrix to the function
  if σ̇ᵢⱼ isa Matrix
    sym_test = σ̇ᵢⱼ - σ̇ᵢⱼ'
    if all(.≈(0),sym_test)
      @warn "σ̇ᵢⱼ matrix input is not symmetric"
      σ̇ᵢⱼ = SymmetricTensor{2,3}(σ̇ᵢⱼ)
    else
      println("printed :       ",[elem.value.value for elem in sym_test])
    end
  end

  𝕀 = SymmetricTensor{2,3}(δ) # Second order identity tensor

  # stress at previous timestep
  σ = (1/3)*(tr(σᵢⱼ))
  sᵢⱼ = dev(σᵢⱼ)
  τ = get_τ(sᵢⱼ)

  # stress derivatives
  σ̇ = (1/3)*(tr(σ̇ᵢⱼ))
  τ̇ = sᵢⱼ ⊡ σ̇ᵢⱼ / (2*τ)

  if R1_allowed && (Ḋ == 0)
    KI = compute_KI(r,σ,τ,D)
    (KI <= 0) && return (1/(2*r.G)) * ( σ̇ᵢⱼ - (3r.ν/(1+r.ν))*σ̇*𝕀 )
  end
  #damage constants and derivatives
  if r.D₀ == 0
    A1, B1 = 0.0, 0.0
    dA1dD, dB1dD = 0.0, 0.0
  else
    c1, c2, c3 = compute_c1c2c3(r,D)
    A, B = compute_AB(r,c1,c2,c3)
    A1, B1 = compute_A1B1(r,A,B)
    dc1dD = compute_dc1dD(r,D)
    dc2dD = compute_dc2dD(r,D)
    dc3dD = compute_dc3dD(r,D)
    dA1dD = compute_dA1dD(r,dc1dD,dc2dD,dc3dD,c2,c3)
    dB1dD = compute_dB1dD(r,dc1dD,dc2dD,dc3dD,c2,c3)
  end

  t1 = λ₁(A1,B1,σ,τ)*σ̇ᵢⱼ - ( λ₂(r,A1,B1,σ,τ)*σ̇ - (1/3)*A1*B1*τ̇ )*𝕀
  t2 = ( dλ₁dσ(A1,B1,τ)*σ̇ + dλ₁dτ(A1,B1,σ,τ)*τ̇ )*sᵢⱼ
  if damaged_allowed
    t3 = dλ₁dD(A1,B1,dA1dD,dB1dD,σ,τ)*Ḋ*σᵢⱼ - ( dλ₂dD(A1,B1,dA1dD,dB1dD,σ,τ)*Ḋ*σ - (1/3)*Ḋ*(dA1dD*B1 + A1*dB1dD)*τ )*𝕀
  else
    t3 = zero(SymmetricTensor{2, 3})
  end
  ϵ̇ᵢⱼ = t1 + t2 + t3
  !keep_2G && (ϵ̇ᵢⱼ = (1/(2r.G)) * ϵ̇ᵢⱼ)
  return ϵ̇ᵢⱼ#convert(SymmetricTensor{2,3,eltype(ϵ̇ᵢⱼ)},ϵ̇ᵢⱼ)
end

function compute_ϵ̇ij(r,D,σij,σijnext,Δt ; damaged_allowed=true)

  # Convert in case automatic differentiation supplies a Matrix to the function
  if σijnext isa Matrix
    sym_test = σijnext - σijnext'
    if all(.≈(0),sym_test)
      σijnext = SymmetricTensor{2,3}(σijnext)
    else
      println("printed :       ",[elem.value.value for elem in sym_test])
    end
  end

  𝕀 = SymmetricTensor{2,3}(δ) # Second order identity tensor

  # stress at previous timestep
  σ = (1/3)*(tr(σij))
  sij = σij - σ*𝕀
  τ = get_τ(sij)

  # stress derivatives
  σ̇ij = (σijnext-σij)/Δt
  σ̇ = (1/3)*(tr(σ̇ij))
  #τ̇_2 = get_τ(ṡᵢⱼ)
  τ̇ = sij ⊡ σ̇ij / (2*τ)
  #@show(τ̇,τ̇_i)

  # initialize damage rate
  Ḋ = 0.0

  #damage constants and derivatives
  if r.D₀ == 0
    A1, B1 = 0.0, 0.0
    dA1dD, dB1dD = 0.0, 0.0
  else
    c1, c2, c3 = compute_c1c2c3(r,D)
    A, B = compute_AB(r,c1,c2,c3)
    A1, B1 = compute_A1B1(r,A,B)
    dc1dD = compute_dc1dD(r,D)
    dc2dD = compute_dc2dD(r,D)
    dc3dD = compute_dc3dD(r,D)
    dA1dD = compute_dA1dD(r,dc1dD,dc2dD,dc3dD,c2,c3)
    dB1dD = compute_dB1dD(r,dc1dD,dc2dD,dc3dD,c2,c3)
  end

  t1 = λ₁(A1,B1,σ,τ)*σ̇ij - ( λ₂(r,A1,B1,σ,τ)*σ̇ - (1/3)*A1*B1*τ̇ )*𝕀
  t2 = ( dλ₁dσ(A1,B1,τ)*σ̇ + dλ₁dτ(A1,B1,σ,τ)*τ̇ )*sij
  if damaged_allowed
    KI = compute_KI(r,(σijnext+σij)/2,D) # compute KI at intermediate stress : TO TEST.
    Ḋ = compute_subcrit_damage_rate(r,KI,D)
    t3 = dλ₁dD(A1,B1,dA1dD,dB1dD,σ,τ)*Ḋ*σij - ( dλ₂dD(A1,B1,dA1dD,dB1dD,σ,τ)*Ḋ*σ - (1/3)*Ḋ*(dA1dD*B1 + A1*dB1dD)*τ )*𝕀
  else
    t3 = zero(Tensor{2, 3})
  end
  ϵ̇ij = 1/(2r.G) * (t1 + t2 + t3)

  #ϵ̇ij = insert_into(ϵ̇ij, -1e-5, (1,1)) ######## !!!!!!!!
  return ϵ̇ij, Ḋ
end

λ₁(A1,B1,σ,τ) =  1 + A1*B1*σ/(2*τ) + B1^2/2
λ₂(r,A1,B1,σ,τ) = 3*r.ν/(1+r.ν) + A1*B1*σ/(2*τ) - A1^2/3 + B1^2/2
λ₃(A1,B1) = (A1*B1)/3
dλ₁dσ(A1,B1,τ) = A1*B1/(2*τ)
dλ₁dτ(A1,B1,σ,τ) = -A1*B1*σ/(2*τ^2)
dλ₁dD(A1,B1,dA1dD,dB1dD,σ,τ) = (dA1dD*B1 + A1*dB1dD)*σ/(2*τ) + B1*dB1dD
function dλ₁dD(r,D,σ,τ)
  c1, c2, c3 = compute_c1c2c3(r,D)
  A, B = compute_AB(r,c1,c2,c3)
  A1, B1 = compute_A1B1(r,A,B)
  dc1dD = compute_dc1dD(r,D)
  dc2dD = compute_dc2dD(r,D)
  dc3dD = compute_dc2dD(r,D)
  dA1dD = compute_dA1dD(r,dc1dD,dc2dD,dc3dD,c2,c3)
  dB1dD = compute_dB1dD(r,dc1dD,dc2dD,dc3dD,c2,c3)
  return (dA1dD*B1 + A1*dB1dD)*σ/(2*τ) + B1*dB1dD
end
function dλ₁dt(r,D,σ,τ,Ḋ,σ̇,τ̇)
  c1, c2, c3 = compute_c1c2c3(r,D)
  A, B = compute_AB(r,c1,c2,c3)
  A1, B1 = compute_A1B1(r,A,B)
  dc1dD = compute_dc1dD(r,D)
  dc2dD = compute_dc2dD(r,D)
  dc3dD = compute_dc2dD(r,D)
  dA1dD = compute_dA1dD(r,dc1dD,dc2dD,dc3dD,c2,c3)
  dB1dD = compute_dB1dD(r,dc1dD,dc2dD,dc3dD,c2,c3)
  return dλ₁dD(A1,B1,dA1dD,dB1dD,σ,τ)*Ḋ + dλ₁dσ(A1,B1,τ)*σ̇ + dλ₁dτ(A1,B1,σ,τ)*τ̇
end

dλ₂dD(A1,B1,dA1dD,dB1dD,σ,τ) = dλ₁dD(A1,B1,dA1dD,dB1dD,σ,τ) - (2/3)*A1*dA1dD
dλ₃dD(A1,B1,dA1dD,dB1dD) = (A1*dB1dD + dA1dD*B1)/3


text:
using ColorSchemes
using ColorTypes
using Images
using Optim
using Random

include("functions.jl")
include("../OptMaps.jl")
using .OptMaps

dimensions = [3, 5, 7, 10, 20]
all_funcs = Dict()

func_types = ["sphere", "rastrigin", "rosenbrock", "styblinskitang", "linear"]
funcs = [get_sphere_funcs, get_rastrigin_funcs, get_rosenbrock_funcs, get_styblinskitang_funcs, get_linear_funcs]

for (ftype, fs) in zip(func_types, funcs)
    all_funcs[ftype] = fs(dimensions)
end

solvers = [NelderMead, LBFGS, GradientDescent, SimulatedAnnealing, ConjugateGradient, Newton]
solver_names = ["NelderMead", "LBFGS", "GradientDescent", "SimulatedAnnealing", "ConjugateGradient", "Newton"]

starting_points = ["half", "ones", "two_and_half", "threes", "rand"]
starting_points_elt = [0.5, 1.0, 2.5, 3.0, 0.0]

Random.seed!(787)

for ftype in func_types
    fs = all_funcs[ftype]

    for dim in dimensions
        println("Function: $ftype, dimensions: $dim")

        variable_domains = convert(Array{Tuple{Float64,Float64},1}, [])
        variable_names = convert(Array{String,1}, [])

        for i in 1:dim
            push!(variable_domains, (-5.0,5.0))
            push!(variable_names, string(i))
        end

        O = OptMap(variable_names, variable_domains, fs[dim]; verbose=true)
        create_map!(O)

        for (starting_point, starting_point_elt) in zip(starting_points, starting_points_elt)
            if starting_point == "rand"
                initial_x = (rand(Float64, dim) .* 10.0) .- 5.0
            else
                initial_x = zeros(Float64, dim)
                initial_x = initial_x .+ starting_point_elt
            end

            for (solver, solver_name) in zip(solvers, solver_names)
                result = optimize(fs[dim], initial_x, solver(), Optim.Options(store_trace = true, extended_trace=true); autodiff = :forward)
                result_log = "$(ftype)_$(dim)_$(solver_name)_$(starting_point)_optmap.log"

                open(result_log, "w") do f
                    println(f, "$(ftype), $(dim) dims, $solver_name - starting from $starting_point")
                    println(f, "$(initial_x)")
                    print(f, result)
                end

                println("#######################################################")
                println("$solver_name - $starting_point")

                if solver_name == "NelderMead"
                    tr = Optim.trace(result)
                    points = [state.metadata["centroid"] for state in tr]
                else
                    points = Optim.x_trace(result)
                end

                solution_path = zeros(Float64, length(variable_names), length(points))

                for i in 1:length(points)
                    for j in 1:length(variable_names)
                        solution_path[j, i] = points[i][j]
                    end
                end

                draw_solution!(O, solution_path)
                save_all(O, "$(ftype)_$(dim)_$(solver_name)_$(starting_point)", ".")
            end
        end

        O = nothing

    end
end

max_stars_repo_name: nmaterise/optmap

max_issues_repo_name: nmaterise/optmap

max_forks_repo_name: nmaterise/optmap


text:
"""
    DiscreteTimeTypes

Type for discrete time handling.
"""
DiscreteTimeTypes = Union{Integer, TimeType}

"""
    SupportedTimeTypes

Types supported for the timestamps data in the time series.
"""
SupportedTimeTypes = Union{Real, DiscreteTimeTypes}

"""
    TS{P}

Abstract type to store data. P is the type of the time stamps in the stored time series.

Subtypes are :
    - TimeSeries
    - MarkedTimeSeries
"""
abstract type TS end

"""
    start_time(ts::TS)

Get the start time of a time series.
"""
function start_time(ts::TS)
    return first(ts.times)
end

"""
    end_time(ts::TS)

Get the end time of a time series.
"""
function end_time(ts::TS)
    return last(ts.times)
end

include("TS/TimeSeries.jl")
include("TS/MarkedTimeSeries.jl")


text:
# Raw types from node_api.h
module NapiTypes

import Base: ==

# Types
export NapiPointer, NapiEnv, NapiValue,
    NapiRef, NapiHandleScope, NapiDeferred,
    NapiFinalize, NapiCallback,
    NapiCallbackInfo,
    NapiTypedArrayInfo,
    NapiArrayBufferInfo,
    NapiExtendedErrorInfo,
    NapiPropertyAttributes,
    NapiPropertyDescriptor
# Enums
export NapiValueType, NapiTypedArrayType, NapiStatus

abstract type AbstractNapiPointer end
primitive type NapiPointer <: AbstractNapiPointer 64 end
primitive type NapiEnv <: AbstractNapiPointer 64 end
primitive type NapiValue <: AbstractNapiPointer 64 end
Base.convert(::Type{T}, x::UInt64) where T <: Union{NapiEnv, NapiPointer} = reinterpret(T, x)
Base.convert(::Type{T}, x::AbstractNapiPointer) where T <: AbstractNapiPointer = reinterpret(T, x)
Base.convert(::Type{T}, x::NapiPointer) where T <: Ptr = reinterpret(T, x)
Base.convert(::Type{T}, x::NapiEnv) where T <: Ptr = reinterpret(T, x)
Base.convert(::Type{T}, x::NapiValue) where T <: Ptr = reinterpret(T, x)
Base.Ptr{T}(x::NapiPointer) where T = convert(Ptr{T}, x)
Base.convert(T::Type{NapiPointer}, x::Ptr) = reinterpret(T, UInt64(x))
Base.convert(T::Type{NapiEnv}, x::Ptr) = reinterpret(T, UInt64(x))
Base.convert(T::Type{NapiValue}, x::Ptr) = reinterpret(T, UInt64(x))
Base.convert(::Type{NapiPointer}, x::NapiPointer) = x
Base.convert(::Type{NapiEnv}, x::NapiEnv) = x
Base.convert(::Type{NapiValue}, x::NapiValue) = x
Base.show(io::IO, x::AbstractNapiPointer) = print(io, string(typeof(x), ": ", convert(Ptr{Nothing}, x)))
NapiPointer(x::Ptr) = convert(NapiPointer, x)
==(x::Ptr, y::AbstractNapiPointer) = x == Ptr{Nothing}(y)
const NapiRef = NapiPointer
const NapiHandleScope = NapiPointer
const NapiDeferred = NapiPointer
const NapiFinalize = Ptr{Cvoid}
const NapiCallback = Ptr{Cvoid}

@enum NapiValueType begin
    napi_undefined
    napi_null
    napi_boolean
    napi_number
    napi_string
    napi_symbol
    napi_object
    napi_function
    napi_external
    napi_bigint
end

@enum NapiTypedArrayType begin
    napi_int8_array
    napi_uint8_array
    napi_uint8_clamped_array
    napi_int16_array
    napi_uint16_array
    napi_int32_array
    napi_uint32_array
    napi_float32_array
    napi_float64_array
    napi_bigint64_array
    napi_biguint64_array
end

@enum NapiStatus begin
    napi_ok
    napi_invalid_arg
    napi_object_expected
    napi_string_expected
    napi_name_expected
    napi_function_expected
    napi_number_expected
    napi_boolean_expected
    napi_array_expected
    napi_generic_failure
    napi_pending_exception
    napi_cancelled
    napi_escape_called_twice
    napi_handle_scope_mismatch
    napi_callback_scope_mismatch
    napi_queue_full
    napi_closing
    napi_bigint_expected
    napi_date_expected
    napi_arraybuffer_expected
    napi_detachable_arraybuffer_expected
    napi_would_deadlock
end
Base.getindex(status::NapiStatus) = status

struct NapiTypedArrayInfo
    type::NapiTypedArrayType
    length::Csize_t
    data::Ptr{Cvoid}
    arraybuffer::NapiValue
    byte_offset::Csize_t
end

struct NapiArrayBufferInfo
    data::Ptr{UInt8}
    length::Csize_t
end

struct NapiExtendedErrorInfo
    error_message::Cstring
    engine_reserved::NapiPointer
    engine_error_code::UInt32
    error_code::NapiStatus
end

@enum NapiPropertyAttributes begin
    napi_default = 0
    napi_writable = 1 << 0
    napi_enumerable = 1 << 1
    napi_configurable = 1 << 2

    napi_static = 1 << 10

    # napi_default_method = napi_writable | napi_configurable
    napi_default_method = (1 << 0) | (1 << 2)
  #   napi_default_jsproperty = napi_writable | napi_enumerable | napi_configurable
    napi_default_jsproperty = (1 << 0) | (1 << 1) | (1 << 2)
end

Base.@kwdef struct NapiPropertyDescriptor
    # One of utf8name or name should be NULL.
    utf8name::Cstring = C_NULL
    name::NapiValue = C_NULL

    method::NapiCallback = C_NULL
    getter::NapiCallback = C_NULL
    setter::NapiCallback = C_NULL

    value::NapiValue = C_NULL

    attributes::NapiPropertyAttributes = napi_default_jsproperty
    data::NapiPointer = C_NULL
end

end


text:
module genlasso_test

using Test, ConstrainedLasso, ECOS, Random, LinearAlgebra

@info("Test genlasso")

y = randn(20)
n = p = size(y, 1)
X = Matrix{Float64}(I, n,n)
D = [eye(p-1) zeros(p-1, 1)] - [zeros(p-1, 1) Matrix{Float64}(I, p-1,p-1)]
β̂path, = genlasso(X, y; D = D)
tmp = round.(β̂path, 6)

@testset "fused lasso" begin
for i in 1:(size(tmp, 2) - 1)
  @test length(unique(tmp[:, i])) < n
end
end

end # end of module


text:
"""
Utility module to help define pages (dashboards) as bundles of route, view and model.
"""
module Pages

using Genie.Router
import Genie.Router: Route

using Genie.Renderers
using Genie.Renderers.Html
using Stipple

export Page
export pages

@reactive struct EmptyModel <: ReactiveModel
end

mutable struct Page{M<:ReactiveModel}
  route::Route
  view::Union{Genie.Renderers.FilePath,<:String}
  model::Type{M}
  layout::Union{Genie.Renderers.FilePath,Nothing}
end

const _pages = Page[]
pages() = _pages

function Page(  route::Union{Route,String};
                view::Union{Genie.Renderers.FilePath,<:String,ParsedHTMLString},
                model::Union{M,Function,Nothing} = Stipple.init(EmptyModel),
                layout::Union{Genie.Renderers.FilePath,<:String,Nothing} = nothing,
                context::Module = @__MODULE__
              ) where {M<:ReactiveModel}
  route = isa(route, String) ? Route(; method = GET, path = route) : route
  view = isa(view, String) ? filepath(view) :
          isa(view, ParsedHTMLString) ? string(view) :
            view
  layout = isa(layout, String) ? filepath(layout) :
            isa(layout, ParsedHTMLString) ? string(layout) :
              layout

  route.action = () -> html(view; layout, context, model = (isa(model,Function) ? Base.invokelatest(model) : model))

  Router.route(route)

  push!(_pages, Page(route, view, typeof((isa(model,Function) ? Base.invokelatest(model) : model)), layout))
end

end


text:
"""
    plotsp(X, wl = 1:size(X, 2); color = nothing, kwargs...)
    Plotting spectra.
* `x` : X-data.
* `wl` : Column names of `x`. Must be numeric.
* `colors` : Set a unique color to the spectra.
* `kwargs` : Optional arguments to pass in `Axis`.
Plots lines corresponding to the rows of `x`.
"""
function plotsp(X, wl = 1:size(X, 2); color = nothing, kwargs...)
    n = size(X, 1)
    f = Figure()
    ax = Axis(f; kwargs...)
    @inbounds for i = 1:n
        if isnothing(color)
            lines!(ax, wl, X[i, :])
        else
            lines!(ax, wl, X[i, :]; color = color)
        end
    end
    f[1, 1] = ax
    (f = f, ax = ax)
end


text:
using Makie
using Makie: textslider
using QuadGK

f(x) = exp(-x^3) # anything you want
g(x, λ) = f(x) * cos(λ * x)

function main()
    λ_slider, λ_obs = textslider(range(1, stop = 10000, step=100), "λ", start = 1)
    x = -1:0.001:1
    plotting = lift(λ_obs) do λ
        Point2f0.(x, g.(x, λ))
    end

    integral = lift(λ_obs) do λ
        res=quadgk(x->g(x,λ),-1,1)[1]
        "$res"
    end

    s1 = lines(x, f.(x))
    s2 = lines(plotting)
    Makie.text!(s2,integral,textsize=0.1,position = Point2f0(1.,1.))
    RecordEvents(
        vbox(s1,hbox(s2, λ_slider), parent=Scene()),
        "output"
    )
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
module TestExtras
	using Base.Test
	using DataArrays
	using StatsBase

	##########
	## countmap
	##########

	d = @data [NA,3,3]
	w = weights([1.1,2.2,3.3])
	cm = Dict{Union(Int, NAtype), Int}(@data([NA,3]), [1,2])
	cmw = Dict{Union(Int, NAtype), Real}(@data([NA,3]), [1.1,5.5])
	@assert isequal(countmap(d), cm)
	@assert isequal(countmap(d, w), cmw)

	##########
	## cut
	##########

	@assert isequal(cut([2, 3, 5], [1, 3, 6]), PooledDataArray(["(1,3]", "(1,3]", "(3,6]"]))
	@assert isequal(cut([2, 3, 5], [3, 6]), PooledDataArray(["[2,3]", "[2,3]", "(3,6]"]))
	@assert isequal(cut([2, 3, 5, 6], [3, 6]), PooledDataArray(["[2,3]", "[2,3]", "(3,6]", "(3,6]"]))
	@assert isequal(cut([1, 2, 4], [1, 3, 6]), PooledDataArray(["[1,3]", "[1,3]", "(3,6]"]))
	@assert isequal(cut([1, 2, 4], [3, 6]), PooledDataArray(["[1,3]", "[1,3]", "(3,6]"]))
	@assert isequal(cut([1, 2, 4], [3]), PooledDataArray(["[1,3]", "[1,3]", "(3,4]"]))
	@assert isequal(cut([1, 5, 7], [3, 6]), PooledDataArray(["[1,3]", "(3,6]", "(6,7]"]))

	ages = [20, 22, 25, 27, 21, 23, 37, 31, 61, 45, 41, 32]
	bins = [18, 25, 35, 60, 100]
	cats = cut(ages, bins)
	pdv = PooledDataArray(["(18,25]", "(18,25]", "(18,25]",
	                       "(25,35]", "(18,25]", "(18,25]",
	                       "(35,60]", "(25,35]", "(60,100]",
	                       "(35,60]", "(35,60]", "(25,35]"])
	@assert isequal(cats, pdv)

	##########
	## rep
	##########

	@assert rep(3, 2) == [3,3]
	@assert rep([3,4], 2) == [3,4,3,4]
	@assert rep([3,4], 1, 2) == [3,3,4,4]
	@assert rep([3,4], each = 2) == [3,3,4,4]
	@assert rep([3,4], times = 2) == [3,4,3,4]
	@assert rep([3,4], [2,3]) == [3,3,4,4,4]
	@assert isequal(rep((@data [NA,3,4]), 2), (@data [NA,3,4,NA,3,4]))
	@assert isequal(rep((@data [NA,3,4]), [2,1,2]), (@data [NA,NA,3,4,4]))
	@assert isequal(rep((@data [NA,3,4]), [2,1,0]), (@data [NA,NA,3]))
end


text:
# auxiliary/unique2.jl
#
# Similar to unique, but removes all occurences of elements that appear more than once.

"""unique2(itr)

Returns an array containing only those elements that appear exactly once in itr,
and without any elements that appear more than once.
"""
function unique2(itr)
    out = collect(itr)
    i = 1
    while i < length(out)
        inext = findnext(out,out[i],i+1)
        if inext == 0
            i += 1
            continue
        end
        while inext != 0
            deleteat!(out,inext)
            inext = findnext(out,out[i],i+1)
        end
        deleteat!(out,i)
    end
    out
end


text:
##############################################################################
##
## VcovData stores data you need to compute errors
##
##############################################################################

type VcovData{T, N}
    regressors::Matrix{Float64}       # X
    crossmatrix::T                    # X'X in the simplest case. Can be Matrix but preferably Factorization
    residuals::Array{Float64, N}      # vector or matrix of residuals (matrix in the case of IV, residuals of Xendo on (Z, Xexo))
    df_residual::Int
end

typealias VcovDataVector{T} VcovData{T, 1}
typealias VcovDataMatrix{T} VcovData{T, 2}

##############################################################################
##
## Each method to compute standard errors should define two types
##
## 1. A type that inherits from AbstractVcovMethod (used for the construction)
## It must implement two methods:
## - allvars that specify variables needed
## - VcovMethodData(x, df::AbstractDataFrame) that construct a VcocMethodData object
##   Typicially, it uses a list of symbols and a dataframe to obtain vectors
##
## 2. A type that inherits from AbstractVcovMethodData (used for the computation)
## It must implement two methods:
## - shat! returns a S hat matrix. It may change regressors in place
## - vcov! returns a covariance matrix. It may change regressors in matrix
##
##############################################################################

##############################################################################
##
## VcovMethod stores data you need to compute errors
##
##############################################################################

abstract AbstractVcovMethod
allvars(::AbstractVcovMethod) = Symbol[]

##############################################################################
##
## VcovDataMethod stores data you need to compute errors
##
##############################################################################

abstract AbstractVcovMethodData

# utilities

sandwich(H, S::Matrix{Float64}) = H \ S * inv(H)
df_FStat(::AbstractVcovMethodData, x::VcovData, hasintercept::Bool) = x.df_residual - hasintercept


text:
include("products/executable_generators.jl")
include("products/file_generators.jl")
include("products/library_generators.jl")

macro generate_init_header(dependencies...)
    deps_path_add = Expr[]
    if !isempty(dependencies)
        for dep in dependencies
            push!(deps_path_add, quote
                append!(PATH_list, $(dep).PATH_list)
                append!(LIBPATH_list, $(dep).LIBPATH_list)
            end)
        end
    end

    return excat(
        # This either calls `@artifact_str()`, or returns a constant string if we're overridden.
        :(global artifact_dir = find_artifact_dir()),

        # Initialize PATH_list and LIBPATH_list
        deps_path_add...,
    )
end

macro generate_init_footer()
    return esc(quote
        # Filter out duplicate and empty entries in our PATH and LIBPATH entries
        unique!(PATH_list)
        unique!(LIBPATH_list)
        global PATH = join(PATH_list, $(pathsep))
        global LIBPATH = join(vcat(LIBPATH_list, Base.invokelatest(JLLWrappers.get_julia_libpaths)), $(pathsep))
    end)
end


text:
function besselj_ELZOUKA(n::Int, x::NN) where NN <: Number
end

function bessely_ELZOUKA(n::Int, x::NN) where NN <: Number
end

"""
Calculate Bessel functions using recurrence relation
"""
function recurrence_relation_bessel(function_here::Function, n::Int, x::NN) where NN <: Number
    return (2*(n-1)/x) * function_here(n-1, x) - function_here(n-2, x)
end

"""
    Recurrence formula for ALL spherical Bessel functions.
I mean by spherical Bessel functions the following:
jₙ, spherical Bessel of the first kind
yₙ, spherical Bessel of the second kind
hₙ⁽¹⁾
hₙ⁽²⁾
http://dlmf.nist.gov/10.51.E1
"""
function recurrence_relation_spherical_bessel(function_here::Function, n::Int, x::NN) where NN <: Number
    return ((2*(n-1)+1)/x) * function_here(n-1, x) - function_here(n-2, x)
end

function recurrence_relation_spherical_bessel(n::Int, x::NN, bessel_n_minus_1::NN, bessel_n_minus_2::NN) where NN <: Number
    return (2*(n-1)+1)/x * bessel_n_minus_1 - bessel_n_minus_2
end

"""
    Calculate spherical Bessel function, given the Bessel functions at n=0 and n=1
Works for ALL spherical Bessel functions
"""
function spherical_Bessel_n_ELZOUKA(n::Int, x::NN, bessel_0::NN, bessel_1::NN) where NN <: Number
    if n == 0
        return bessel_0
    elseif n == 1
        return bessel_1
    else
        # use the recurrence
        n_now = 2
        bessel_n_minus_2 = bessel_0
        bessel_n_minus_1 = bessel_1
        while n_now <= n
            bessel_n_now = recurrence_relation_spherical_bessel(n_now, x, bessel_n_minus_1, bessel_n_minus_2)
            if n_now == n
                return bessel_n_now
            else
                bessel_n_minus_2 = bessel_n_minus_1
                bessel_n_minus_1 = bessel_n_now
                n_now += 1
            end
        end
    end
end

# jₙ, spherical Bessel of the first kind =======================================
"""
    Spherical Bessel function of the first kind, at n=0
from http://dlmf.nist.gov/10.49.E3
"""
function spherical_Bessel_j_0_ELZOUKA(x::NN) where NN <: Number
    return sin(x)/x
end

"""
    Spherical Bessel function of the first kind, at n=1
from http://dlmf.nist.gov/10.49.E3
"""
function spherical_Bessel_j_1_ELZOUKA(x::NN) where NN <: Number
    return sin(x)/(x^2) - cos(x)/x
end

"""
    Spherical Bessel function of the first kind, using recurrence relation
"""
function spherical_Bessel_j_n_ELZOUKA(n::Int, x::NN) where NN <: Number
    if n == 0
        return spherical_Bessel_j_0_ELZOUKA(x)
    elseif n == 1
        return spherical_Bessel_j_1_ELZOUKA(x)
    else
        bessel_0 = spherical_Bessel_j_0_ELZOUKA(x)
        bessel_1 = spherical_Bessel_j_1_ELZOUKA(x)
        return spherical_Bessel_n_ELZOUKA(n, x, bessel_0, bessel_1)
    end
end

ChainRulesCore.@scalar_rule(
    spherical_Bessel_j_n_ELZOUKA(n::Integer, x::Real),
    (
        ChainRulesCore.ZeroTangent(),
        (spherical_Bessel_j_n_ELZOUKA(n - 1, x) - spherical_Bessel_j_n_ELZOUKA(n + 1, x)) / 2,
    )
)

# yₙ, spherical Bessel of the second kind =======================================
"""
    Spherical Bessel function of the second kind, at n=0
from http://dlmf.nist.gov/10.49.E5
"""
function spherical_Bessel_y_0_ELZOUKA(x::NN) where NN <: Number
    return -cos(x)/x
end

"""
    Spherical Bessel function of the second kind, at n=1
from http://dlmf.nist.gov/10.49.E5
"""
function spherical_Bessel_y_1_ELZOUKA(x::NN) where NN <: Number
    return -cos(x)/(x^2) - sin(x)/x
end

"""
    Spherical Bessel function of the second kind, using recurrence relation
"""
function spherical_Bessel_y_n_ELZOUKA(n::Int, x::NN) where NN <: Number
    if n == 0
        return spherical_Bessel_y_0_ELZOUKA(x)
    elseif n == 1
        return spherical_Bessel_y_1_ELZOUKA(x)
    else
        bessel_0 = spherical_Bessel_y_0_ELZOUKA(x)
        bessel_1 = spherical_Bessel_y_1_ELZOUKA(x)
        return spherical_Bessel_n_ELZOUKA(n, x, bessel_0, bessel_1)
    end
end

ChainRulesCore.@scalar_rule(
    spherical_Bessel_y_n_ELZOUKA(n::Integer, x::Real),
    (
        ChainRulesCore.ZeroTangent(),
        (spherical_Bessel_y_n_ELZOUKA(n - 1, x) - spherical_Bessel_y_n_ELZOUKA(n + 1, x)) / 2,
    )
)


text:
module TriangulationsTests

using Gridap.Geometry
using Gridap.Arrays
using Gridap.ReferenceFEs
using Test
using FillArrays

domain = (-1,1,-1,1)
cells = (10,10)
model = CartesianDiscreteModel(domain,cells)
Ω = Triangulation(model)
test_triangulation(Ω)
Ω_ = Triangulation(Ω)
@test Ω_ === Ω

@test model === get_background_model(Ω)
@test model === get_active_model(Ω)
@test get_grid(Ω) === get_grid(model)
glue = get_glue(Ω,Val(2))
@test isa(glue.tface_to_mface,IdentityVector)
@test isa(glue.mface_to_tface,IdentityVector)
glue.mface_to_tface === glue.tface_to_mface
@test isa(glue.tface_to_mface_map,Fill)

glue_0 = get_glue(Ω,Val(0))
glue_1 = get_glue(Ω,Val(1))
glue_2 = get_glue(Ω,Val(2))

Θ = GenericTriangulation(
  get_grid(Ω),
  get_background_model(Ω),
  (glue_0,glue_1,glue_2))
test_triangulation(Θ)
@test get_glue(Θ,Val(0)) === glue_0
@test get_glue(Θ,Val(1)) === glue_1
@test get_glue(Θ,Val(2)) === glue_2

Θ = GenericTriangulation(
  get_grid(Ω),
  get_background_model(Ω))
test_triangulation(Θ)

Θ = GenericTriangulation(get_grid(Ω))
@test isa(Θ,Triangulation)

Γ = Triangulation(ReferenceFE{1},model)
@test model === get_background_model(Γ)
glue = get_glue(Γ,Val(1))
@test isa(glue.tface_to_mface,IdentityVector)
@test isa(glue.mface_to_tface,IdentityVector)
glue.mface_to_tface === glue.tface_to_mface
@test isa(glue.tface_to_mface_map,Fill)

cell_xs = get_cell_coordinates(Ω)
cell_mask = lazy_map(cell_xs) do xs
  R = 0.7
  n = length(xs)
  x = (1/n)*sum(xs)
  d = x[1]^2 + x[2]^2 - R^2
  d < 0
end

Ω1 = Triangulation(model,cell_mask)
@test model === get_background_model(Ω1)
@test model !== get_active_model(Ω1)
glue = get_glue(Ω1,Val(2))
@test glue.tface_to_mface == findall(collect1d(cell_mask))
@test isa(glue.tface_to_mface_map,Fill)
@test model === get_background_model(Ω)
@test isa(glue.mface_to_tface,PosNegPartition)
@test glue.mface_to_tface[glue.tface_to_mface] == 1:length(glue.tface_to_mface)

tface_to_val = [ rand(3,4) for i in 1:num_cells(Ω1) ]
mface_to_val = extend(tface_to_val,glue.mface_to_tface)
@test mface_to_val[glue.tface_to_mface] == tface_to_val

cell_mcell = findall(collect1d(cell_mask))
Ω1 = Triangulation(model,cell_mcell)
glue = get_glue(Ω1,Val(2))
@test glue.tface_to_mface == cell_mcell
@test isa(glue.tface_to_mface_map,Fill)
@test model === get_background_model(Ω)
@test isa(glue.mface_to_tface,PosNegPartition)
@test glue.mface_to_tface[glue.tface_to_mface] == 1:length(glue.tface_to_mface)

Ω1 = view(Ω,cell_mcell)
glue = get_glue(Ω1,Val(2))
@test glue.tface_to_mface == cell_mcell
@test isa(glue.tface_to_mface_map,Fill)
@test model === get_background_model(Ω)
@test isa(glue.mface_to_tface,PosNegPartition)
@test glue.mface_to_tface[glue.tface_to_mface] == 1:length(glue.tface_to_mface)

labels = get_face_labeling(model)
entity = num_entities(labels)+1
labels.d_to_dface_to_entity[end][cell_mcell] .= entity
add_tag!(labels,"Ω1",[entity])
Ω1 = Triangulation(model,tags="Ω1")
glue = get_glue(Ω1,Val(2))
@test glue.tface_to_mface == cell_mcell
@test isa(glue.tface_to_mface_map,Fill)
@test model === get_background_model(Ω)
@test isa(glue.mface_to_tface,PosNegPartition)
@test glue.mface_to_tface[glue.tface_to_mface] == 1:length(glue.tface_to_mface)

glue = get_glue(Ω,Val(2))
tface_to_data = rand(num_cells(Ω))
mface_to_data = extend(tface_to_data,glue.mface_to_tface)
@test tface_to_data === mface_to_data

glue = get_glue(Ω1,Val(2))
tface_to_data = rand(num_cells(Ω1))
mface_to_data = extend(tface_to_data,glue.mface_to_tface)
@test mface_to_data[glue.tface_to_mface] == tface_to_data
tface_to_data = get_cell_shapefuns(Ω1)
mface_to_data = extend(tface_to_data,glue.mface_to_tface)
@test mface_to_data[glue.tface_to_mface] == tface_to_data

@test is_change_possible(Ω,Ω1)
@test is_change_possible(Ω1,Ω)

Ω2 = best_target(Ω1,Ω)
glue = get_glue(Ω2,Val(2))
@test isa(glue.tface_to_mface,IdentityVector)
@test isa(glue.mface_to_tface,IdentityVector)

end # module


text:
function plot_results(;paperplots=false)
    if paperplots
        close(3)
        figure(3, figsize=(4,3))
    end

    figure(3)
    clf()

    #N = 1./hgrid
    N = Ngrid_list

    region = (l2_rel_err .> 1e-13) .& (l2_rel_err .< 1e-5)
    p = -10
    c = match_convergence(N[region], l2_rel_err[region], p)
    loglog(N, c*N.^(p*1.0), "-k", label=L"\mathcal{O}(N^{-10})", zorder=-1)

    loglog(N, grad_max_rel_err, "*-", label=L"||\nabla \tilde u-\nabla u||_{\infty}")
    loglog(N, grad_l2_rel_err, "+-", label=L"||\nabla \tilde u-\nabla u||_{2}")
    loglog(N, max_rel_err, ".-", label=L"||\tilde u-u||_{\infty}")
    loglog(N, l2_rel_err, "x-", label=L"||\tilde u-u||_{2}")

    region = (grad_l2_rel_err .> 1e-10) .& (grad_l2_rel_err .< 1e-4)
    p = -10
    c = match_convergence(N[region], grad_l2_rel_err[region], p)
    loglog(N, c*N.^(p*1.0), "-k", zorder=-1)

    legend()
    xlabel("N")
    grid("on")
    minorticks_off()
    ticks = [50,100,200,400,800]
    xticks(ticks, round.(Int, ticks))
    ylim(1e-15, 1)
    ylabel("Relative error")

    tight_layout(0)
end


text:
function point_attribs(mh, points_robj, line_robj)
    isoverpoints, isoverlines = mh[1] == points_robj.id, mh[1] == line_robj.id
    points_robj[:glow_color] = isoverpoints ? RGBA{Float32}(0.9,.1,0.2,0.9) : RGBA{Float32}(0.,0.,0.,0.)
    points_robj[:visible] = isoverpoints || isoverlines
end
function point_edit(past, mousediff_index, point_gpu)
    mousediff_past, index_past = past
    mousediff, index = mousediff_index
    if checkbounds(Bool, point_gpu, index)
        point_gpu[index] = point_gpu[index] - eltype(point_gpu)(mousediff-mousediff_past)
    end
    mousediff_index
end
function edit_line(
        line, direction_restriction::Vec2f0, clampto, window;
        color=default(RGBA{Float32}), knob_scale=9f0, kw_args...
    )
    mouse_hover = mouse2id(window)
    inds = reinterpret(Cuint, collect(partition(Cuint(0):Cuint(length(line)-1),2,1)))
    line_robj = visualize(
        line, :linesegment;
        color=color, thickness=1f0,
        indices=inds,
        kw_args...
    ).children[]
    point_gpu = line_robj[:vertex]
    points = visualize(
        (Circle{Float32}(Point2f0(0), knob_scale), point_gpu);
        color=RGBA{Float32}(0.7, 0.7, 0.7, 1.0),
        kw_args...
    )
    point_robj = points.children[]
    gpu_position = point_robj[:position]
    gpu_color = point_robj[:color]
    gpu_scale = point_robj[:scale]

    m2id = mouse2id(window)
    ids = (point_robj.id, line_robj.id)
    isoverpoint = droprepeats(const_lift(is_same_id, m2id, ids))

    @materialize mouse_buttons_pressed, mouseposition = window.inputs

    key_pressed = const_lift(
        GLAbstraction.singlepressed,
        mouse_buttons_pressed, GLFW.MOUSE_BUTTON_LEFT
    )
    mousedragg = GLAbstraction.dragged(
        mouseposition, key_pressed, isoverpoint
    )
    T = Point2f0
    startvalue = (0, 0, Point2f0(0), Point2f0(0))
    sig = foldp(startvalue, mousedragg) do v0, dragg
        id, index, p0, np = v0
        if dragg == Vec2f0(0) # if drag just started. Not the best way, maybe dragged should return a tuple of (draggvalue, started)
            id, index = value(m2id)
            if id==point_robj.id && length(gpu_position) >= index
                p0 = gpu_position[index]
            end
        else
            if id==point_robj.id && length(gpu_position) >= index
                np = p0 + T(dragg).*T(direction_restriction)
                np = T(np[1], clamp(np[2], clampto...))
                gpu_position[index] = np
            end
        end
        return id, index, p0, np
    end

    Context(line_robj, point_robj), map(x->(x[2], x[4]), sig)
end

function widget{T}(
        points::Vector{Point{2,T}}, window;
        color=default(RGBA), kw_args...
    )
    line = visualize(points, :lines; thickness=2f0, kw_args...)
    line_robj = line.children[]
    point_gpu = line_robj[:vertex]
    points = visualize(
        (Circle(Point2f0(0), 8f0), point_gpu);
        glow_width=2f0, kw_args...
    )
    points_robj = points.children[]
    preserve(const_lift(point_attribs, window.inputs[:mouse_hover], points_robj, line_robj))
    mousediff_index = dragged_on(points_robj, MOUSE_LEFT, window)
    preserve(foldp(point_edit, mousediff_index.value, mousediff_index, Signal(point_gpu)))
    Context(line, points)
end

function mouse_drag_diff(past, drag_index)
    index_past, _, drag_past = past
    drag, index = drag_index
    index, drag-drag_past, drag
end

function dragg_gpu(v0, dragg)
    if dragg == Vec2f0(0) # if drag just started. Not the best way, maybe dragged should return a tuple of (draggvalue, started)
        id, index = value(m2id)
        if id==point_robj.id && length(gpu_position) >= index
            p0 = gpu_position[index]
        else
            p0 = v0[3]
        end
    else
        id, index, p0 = v0
        if id==point_robj.id && length(gpu_position) >= index
            gpu_position[index] = Point2f0(p0) + Point2f0(dragg)
        end
    end
    return id, index, p0
end

immutable ClampFunctor{T}
    a::T
    b::T
end
@compat (c::ClampFunctor)(elem) = clamp(elem, c.a, c.b)
Base.clamp(x::FixedVector, a, b) = map(ClampFunctor(a,b) , x)

clampU8{T}(x::RGBA{T}) = RGBA{T}(ntuple(i->clamp(getfield(x, i), 0.,1.), Val{4})...)
channel_color(channel, value) = RGBA{Float32}(ntuple(i->i==channel ? value : 0.0f0, Val{3})..., 1f0)
function c_setindex{T}(color::RGBA{T}, val, channel)
    v = clamp(val, 0, 1)
    RGBA{T}(
        1==channel ? v : comp1(color),
        2==channel ? v : comp2(color),
        3==channel ? v : comp3(color),
        4==channel ? v : alpha(color)
    )
end
function edit_color(tex, buff, index_value, channel, maxval)
    index, value = index_value
    if checkbounds(Bool, tex, index)
        color = c_setindex(buff[index], value[2]/maxval, channel) # we need buff, since getindex is very slow for textures
        buff[index] = color
        tex[index] = color
    end
    nothing
end

function widget{T<:Colorant}(colormap::VecTypes{T}, window;
        area = (300, 30),
        slider_colors = (
            RGBA{Float32}(0.78125,0.1796875,0.41796875),
            RGBA{Float32}(0.41796875,0.78125,0.1796875),
            RGBA{Float32}(0.1796875,0.41796875,0.78125),
            RGBA{Float32}(0.9,0.9,0.9)
        ), knob_scale=9f0,
        kw_args...
    )
    colors = to_cpu_mem(value(colormap))
    N = length(colors)
    color_tex = GLAbstraction.gl_convert(Texture, colormap)
    @assert colors == to_cpu_mem(color_tex)
    scale = Point2f0(area) .* Point2f0(1, 6)
    dir_restrict = Vec2f0(0,1)
    vis = ntuple(Val{4}) do i
        c_channel = Point2f0[
            Point2f0(x, getfield(c, i)) .* scale for (x, c) in zip(linspace(0,1,N), colors)
        ]
        c_i, diff = edit_line(
            c_channel, dir_restrict,
            (0, scale[2]), window, color=slider_colors[i],
            knob_scale=knob_scale
        )
        preserve(const_lift(edit_color, color_tex, colors, diff, i, scale[2]))
        c_i
    end
    tex = visualize(
        color_tex;
        is_fully_opaque=false,
        primitive=SimpleRectangle{Float32}(0, scale[2]+6, scale[1], 10),
        kw_args...
    )
    color_tex, Context(tex, vis...)
end