text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
"""
	add_outage_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict)

Adds the Outages results to the dictionary passed back from `run_reopt` using the solved model `m` and the `REoptInputs`.
Only added to results when multiple outages are modeled via the `ElectricUtility.outage_durations` input.

!!! note
	When modeling PV the name of the PV system is used for the output keys to allow for modeling multiple PV systems. The default PV name is `PV`.

!!! warn
	The Outage results can be very large when many outages are modeled and can take a long time to generate.

Outages results:
- `expected_outage_cost` The expected outage cost over the random outages modeled.
- `max_outage_cost_per_outage_duration_series` The maximum outage cost in every outage duration modeled.
- `unserved_load_series` The amount of unserved load in each outage and each time step.
- `unserved_load_per_outage_series` The total unserved load in each outage.
- `mg_storage_upgrade_cost` The cost to include the storage system in the microgrid.
- `storage_upgraded` Boolean that is true if it is cost optimal to include the storage system in the microgrid.
= `discharge_from_storage_series` Array of storage power discharged in every outage modeled.
- `PVmg_kw` Optimal microgrid PV capacity. Note that the name `PV` can change based on user provided `PV.name`.
- `mg_PV_upgrade_cost` The cost to include the PV system in the microgrid.
- `mgPV_to_storage_series` Array of PV power sent to the battery in every outage modeled.
- `mgPV_curtailed_series` Array of PV curtailed in every outage modeled.
- `mgPV_to_load_series` Array of PV power used to meet load in every outage modeled.
- `Generatormg_kw` Optimal microgrid Generator capacity. Note that the name `Generator` can change based on user provided `Generator.name`.
- `mg_Generator_upgrade_cost` The cost to include the Generator system in the microgrid.
- `mgGenerator_to_storage_series` Array of Generator power sent to the battery in every outage modeled.
- `mgGenerator_curtailed_series` Array of Generator curtailed in every outage modeled.
- `mgGenerator_to_load_series` Array of Generator power used to meet load in every outage modeled.
- `mg_Generator_fuel_used` Array of Generator fuel used in every outage modeled.

!!! warn
	The output keys for "Outages" are subject to change.
"""
function add_outage_results(m, p, d::Dict)
	# TODO with many outages the dispatch arrays are so large that it can take hours to create them
	# (eg. 8760 * 12 hour outages with PV, storage and diesel makes 7*12*8760 = 735,840 values)
	# For now the outage dispatch outputs are not created (commented out below). Perhaps make a new
	# function to optionally get the outage dispatch values so that we don't slow down returning the
	# other results.
	r = Dict{String, Any}()
	r["expected_outage_cost"] = value(m[:ExpectedOutageCost])
	r["max_outage_cost_per_outage_duration_series"] = value.(m[:dvMaxOutageCost]).data
	r["unserved_load_series"] = value.(m[:dvUnservedLoad]).data
	S = length(p.s.electric_utility.scenarios)
	T = length(p.s.electric_utility.outage_start_time_steps)
	unserved_load_per_outage = Array{Float64}(undef, S, T)
	for s in 1:S, t in 1:T
		unserved_load_per_outage[s, t] = sum(r["unserved_load_series"][s, t, ts] for
											 ts in 1:p.s.electric_utility.outage_durations[s])
	# need the ts in 1:p.s.electric_utility.outage_durations[s] b/c dvUnservedLoad has unused values in third dimension
	end
	r["unserved_load_per_outage_series"] = round.(unserved_load_per_outage, digits=2)
	r["mg_storage_upgrade_cost"] = value(m[:dvMGStorageUpgradeCost])
	r["discharge_from_storage_series"] = value.(m[:dvMGDischargeFromStorage]).data

	for t in p.techs.all
		r[t * "_upgraded"] = value(m[:binMGTechUsed][t])
	end
	r["storage_upgraded"] = value(m[:binMGStorageUsed])

	if !isempty(p.techs.pv)
		for t in p.techs.pv

			# need the following logic b/c can have non-zero mg capacity when not using the capacity
			# due to the constraint for setting the mg capacities equal to the grid connected capacities
			if Bool(round(r[t * "_upgraded"], digits=1))
				r[string(t, "mg_kw")] = round(value(m[:dvMGsize][t]), digits=4)
			else
				r[string(t, "mg_kw")] = 0
			end
			r[string("mg_", t, "_upgrade_cost")] = round(value(m[:dvMGTechUpgradeCost][t]), digits=2)

			if !isempty(p.s.storage.types.elec)
				PVtoBatt = (m[:dvMGProductionToStorage][t, s, tz, ts] for
					s in p.s.electric_utility.scenarios,
					tz in p.s.electric_utility.outage_start_time_steps,
					ts in p.s.electric_utility.outage_time_steps)
			else
				PVtoBatt = []
			end
			r[string("mg", t, "_to_storage_series")] = round.(value.(PVtoBatt), digits=3)

			PVtoCUR = (m[:dvMGCurtail][t, s, tz, ts] for
				s in p.s.electric_utility.scenarios,
				tz in p.s.electric_utility.outage_start_time_steps,
				ts in p.s.electric_utility.outage_time_steps)
			r[string("mg", t, "_curtailed_series")] = round.(value.(PVtoCUR), digits=3)

			PVtoLoad = (
				m[:dvMGRatedProduction][t, s, tz, ts] * p.production_factor[t, tz+ts]
						* p.levelization_factor[t]
				- m[:dvMGCurtail][t, s, tz, ts]
				- m[:dvMGProductionToStorage][t, s, tz, ts] for
					s in p.s.electric_utility.scenarios,
					tz in p.s.electric_utility.outage_start_time_steps,
					ts in p.s.electric_utility.outage_time_steps
			)
			r[string("mg", t, "_to_load_series")] = round.(value.(PVtoLoad), digits=3)
		end
	end

	if !isempty(p.techs.gen)
		for t in p.techs.gen

			# need the following logic b/c can have non-zero mg capacity when not using the capacity
			# due to the constraint for setting the mg capacities equal to the grid connected capacities
			if Bool(round(r[t * "_upgraded"], digits=1))
				r[string(t, "_mg_kw")] = round(value(m[:dvMGsize][t]), digits=4)
			else
				r[string(t, "mg_kw")] = 0
			end

			r[string("mg_", t, "_fuel_used_per_outage_series")] = value.(m[:dvMGFuelUsed][t, :, :]).data
			r[string("mg_", t, "_upgrade_cost")] = round(value(m[:dvMGTechUpgradeCost][t]), digits=2)

			if !isempty(p.s.storage.types.elec)
				GenToBatt = (m[:dvMGProductionToStorage][t, s, tz, ts] for
					s in p.s.electric_utility.scenarios,
					tz in p.s.electric_utility.outage_start_time_steps,
					ts in p.s.electric_utility.outage_time_steps)
			else
				GenToBatt = []
			end
			r[string("mg", t, "_to_storage_series")] = round.(value.(GenToBatt), digits=3)

			GENtoCUR = (m[:dvMGCurtail][t, s, tz, ts] for
				s in p.s.electric_utility.scenarios,
				tz in p.s.electric_utility.outage_start_time_steps,
				ts in p.s.electric_utility.outage_time_steps)
			r[string("mg", t, "_curtailed_series")] = round.(value.(GENtoCUR), digits=3)

			GENtoLoad = (
				m[:dvMGRatedProduction][t, s, tz, ts] * p.production_factor[t, tz+ts]
						* p.levelization_factor[t]
				- m[:dvMGCurtail][t, s, tz, ts]
				- m[:dvMGProductionToStorage][t, s, tz, ts] for
					s in p.s.electric_utility.scenarios,
					tz in p.s.electric_utility.outage_start_time_steps,
					ts in p.s.electric_utility.outage_time_steps
			)
			r[string("mg", t, "_to_load_series")] = round.(value.(GENtoLoad), digits=3)
		end
	end
	d["Outages"] = r
end

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
using GeneralizedGenerated
using JuliaVariables
using Test
using BenchmarkTools
using DataStructures

rmlines = NGG.rmlines

@testset "no kwargs" begin

@gg function f1(a)
    quote
        x -> a + x
    end |> rmlines
end

@test f1(1)(2) == 3

@gg function f2(a)
    quote
        x -> begin
            a += 2
            x + a
        end
    end |> rmlines
end

@test f2(1)(2) == 5

@gg function f3(a)
    quote
        k = 20
        x -> begin
            a += 2
            x + a + k
        end
    end
end

@test f3(1)(2) == 25

end

@testset "kwargs" begin

@gg function f4(a)
    quote
        function z(x, k=1)
            x + 20 + a + k
        end
    end
end

@test f4(10)(2) == 33
end

@testset "namedtuple" begin

@gg function f5(a)
    quote
        function z(x, k=1)
            (x=x, k=k, q=a)
        end
    end
end

@test f5(10)(2) == (x=2, k=1, q=10)
end

@testset "mk funcs" begin

f_ = mk_function(:((x, y) -> x + y))
@test f_(1, 2) == 3

f_ = mk_function(:(function (x, y) x + y end))
@test f_(1, 2) == 3

end

@testset "type encoding more datatypes" begin

@gg function f6(a)
    tp = (1, 2, 3)
    quote
        function z(x, k=$tp)
            (x=x, k=k, q=a)
        end
    end
end

@test f6(10)(2) == (x=2, k=(1, 2, 3), q=10)

@gg function f7(a)
    tp = (a1=1, a2=2, a3=3)
    quote
        function z(x, k=$tp)
            (x=x, k=k, q=a)
        end
    end
end

@test f7(10)(2) == (x=2, k=(a1=1, a2=2, a3=3), q=10)

@gg function f8(a)
    tp = "233"
    quote
        function z(x, k=$tp)
            (x=x, k=k, q=a)
        end
    end
end

@test f8(10)(2) == (x=2, k="233", q=10)

@gg function f9(a)
    tp = list(1, 2, 3)
    quote
        function z(x; k=$tp)
            (x=x, k=k, q=a)
        end
    end
end

@test f9(10)(2) == (x=2, k=list(1, 2, 3), q=10)
@test f9(10)(2; k=10) == (x=2, k=10, q=10)

end

@testset "runtime eval" begin

a = to_type(:(1 + 2))
@test :(1 + 2) == from_type(a)
@test string(from_type(a)) == string(:(1 + 2))

@test runtime_eval(1) == 1
@test mk_function(:(
    x -> x + 1
))(2) == 3

@test_throws Any mk_function(quote
    x -> x + 1
end)

@test runtime_eval(quote
    x -> x + 1
end)(1) == 2

end

@testset "self recursive" begin
    to_test = quote
        g(x, r=0) = x === 0 ? r : begin
            g = g # required for self recur
            g(x-1, r + x)
        end
        g(10)
    end |> runtime_eval

    g(x, r=0) = x === 0 ? r : g(x-1, r + x)
    expected = g(10)
    @test expected == to_test
end

@testset "self recursive" begin
    to_test = quote
        g(x, r=0) = x === 0 ? r : begin
            g = g # required for self recur
            g(x-1, r + x)
        end
        g(10)
    end |> runtime_eval

    g(x, r=0) = x === 0 ? r : g(x-1, r + x)
    expected = g(10)
    @test expected == to_test
end

@testset "support where clauses and return type annotations for @gg" begin
    @gg function foo(x::T) where T
        :(x, T)
    end
    @test foo(1) == (1, Int)
    @gg function bar(x::T) where T
        quote
            g = x + 20
            x = 10
            () -> begin
                x = g
                x
            end
        end
    end
    @test bar(2)() == 2 + 20

    @gg function foobar(x::T, y::A) where {T <: Number, A <: AbstractArray{T}}
        quote
            g = x + 20
            x = 10
            () -> begin
                x = g
                (A, x + y[1])
            end
        end
    end
    @test foobar(2, [3])() == (Vector{Int}, 2 + 20 + 3)
end

@testset "support default arguments" begin
    @gg function h(x, c)
        quote
            d = x + 10
            function g(x, y=c)
                x + y + d
            end
        end
    end
    @test h(1, 2)(3) == 16
end

module S
    run(y) = y + 1
end

struct K
    f1::Function
    f2::Function
end
@testset "specifying evaluation modules" begin
    @gg function g(m::Module, y)
        @under_global :m :(run(y))
    end
    @test g(S, 1) == 2

    @gg function h(m, y)
        @under_global :m quote
            c = m.f1(y)
            () -> begin c = m.f2(c) end
        end
    end
    k = K(x -> x + 1, x -> x * 9)
    next = h(k, 1)
    @test next() == 18
    @test next() == 18 * 9
end

@testset "test free variables of let bindings" begin
    @gg function test_free_of_let()
        quote
            let x = 1
                f = () -> begin
                    x * 3
                end
                x = 2
                f
            end
        end
    end
    @test test_free_of_let()() == 6
end

@testset "show something" begin
    f1 = mk_function(:(x -> x + 1))
    f2 = mk_function(:((x :: Int = 2, ) -> x + 1))
    @test f1(1) == 2
    @test f2() == 3
    println(f1)
    println(f2)
end

@testset "omit func argname: #34" begin
   f1 = mk_function(:( (:: Int) -> 0 ))
   @test f1(1) == 0
   @test_throws MethodError f1("")
end

@testset "support for 1-d arrays' type encoding: #42" begin
    from_type(to_type([1, 2, 3])) == [1, 2, 3]
    from_type(to_type([1, 2, "3"])) == [1, 2, "3"]
end

@testset "#54: fastmath/typeable Val" begin
    f = mk_function([:x],[],Base.FastMath.make_fastmath(:(0.5+1.0*x^2)))
    @test f(1) == 1.5
    @test f(10) == 100.5
end

struct PseudoModule{X, F}
    x :: X
    (+) :: F
end

@testset "struct as module" begin
    @gg f(y) = @under_global PseudoModule(10, *) :(x + y)
    @test f(2) == 20
end

@testset "#57, #59: big expression" begin

using GeneralizedGenerated: mk_function

function rand_ex()
    u = rand()
    c = randn()
    if u < 0.5
        return :(m1 + m2 + r² + $c)
    else
        return :(m1 * m2 / r² * $c)
    end
end

function gen_f_mkfunc(ex)
    return mk_function(@__MODULE__, rmlines(:((m1, m2, r²) -> $ex)))
end

function main(N::Int)
    for i in 1:N
        m1, m2, r² = rand(3)
        ex = rand_ex()
        f_mkfunc = gen_f_mkfunc(ex)
        f_mkfunc(m1, m2, r²)
    end
end

main(5_000)

@test true

end

ext: jl


text:
"""
    AbstractImageSmoothAlgorithm <: AbstractImageFilter

The root type for `ImageSmooth` package.

Any image smoothing algorithm shall subtype it to support
[`smooth`](@ref) and [`smooth!`](@ref) APIs.

# Examples

All algorithm in `ImageSmooth` are called in the following pattern,
take `L0Smooth <: AbstractImageSmoothAlgorithm` as an example:

```julia
# First generate an algorithm instance
fₛ = L0Smooth()

# Then pass the algorithm to `smooth`
imgₛ = smooth(img, fₛ)

# Or use in-place version `smooth!`
imgₛ = similar(img)

smooth!(imgₛ, img, fₛ)
```

Some algorithms also receive parameters to control the smoothing process and
to get an expected smooth degree.

```julia
# You could explicit specify it
f = L0Smooth(λ=0.04, κ=1.5, βmax=1e6)

# Or infer the default value
f = L0Smooth()
```

For more examples, please check [`smooth`](@ref) and [`smooth!`](@ref) and concret
algorithms.
"""
abstract type AbstractImageSmoothAlgorithm <: AbstractImageFilter end

smooth!(out::GenericImage,
        img::GenericImage,
        f::AbstractImageSmoothAlgorithm,
        args...; kwargs...) =
    f(out, img, args...; kwargs...)

function smooth(img::GenericGrayImage,
                f::AbstractImageSmoothAlgorithm,
                args...; kwargs...)
    out = similar(img)
    smooth!(out, img, f, args...; kwargs...)
    return out
end

function smooth(img::GenericImage{<:AbstractRGB},
                f::AbstractImageSmoothAlgorithm,
                args...; kwargs...)
    out = similar(img)
    smooth!(out, img, f, args...; kwargs...)
    return out
end

### Docstrings

"""
    smooth!(out::GenericImage, img::GenericImage, fₛ::AbstractImageSmoothAlgorithm, args...; kwargs...)

Smooth `img::GenericImage` using algorithm `fₛ`

# Output

`out` will be changed in place.

# Examples

Just simply pass an algorithm to `smooth!`:

```julia
# First generate an algorithm instance
fₛ = L0Smooth()

## For Gray or RGB images
imgₛ = similar(img)

smooth!(imgₛ, img, fₛ)
```

See also: [`smooth`](@ref)
"""
smooth!

"""
    smooth(img, fₛ::AbstractImageSmoothAlgorithm, args...; kwargs...)

Smooth `img` using algorithm `fₛ`

# Output

The return image `imgₛ` is an `Array{Gray{N0f8}} for `Gray` image, and `Array{RGB{N0f8}}` for `RGB` image.

# Examples

Just simply pass the input image and algorithm to `smooth`

```julia
fₛ = L0Smooth()

imgₛ = smooth(img, fₛ)
```

This reads as "`smooth` image `img` using smoothing algorithm `fₛ`".

See also [`smooth!`](@ref) for in-place smoothing.
"""
smooth


text:
module TestSDPA

import MathOptInterface
using Test

const MOI = MathOptInterface
const MOIU = MOI.Utilities
const SDPA = MOI.FileFormats.SDPA
const SDPA_TEST_FILE = "test.sdpa"
const SDPA_MODELS_DIR = joinpath(@__DIR__, "models")

function _set_var_and_con_names(model::MOI.ModelLike)
    variable_names = String[]
    for j in MOI.get(model, MOI.ListOfVariableIndices())
        var_name_j = "v" * string(j.value)
        push!(variable_names, var_name_j)
        MOI.set(model, MOI.VariableName(), j, var_name_j)
    end

    idx = 0
    constraint_names = String[]
    single_variable_constraints = Tuple[]
    for i in MOI.get(
        model,
        MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.Integer}(),
    )
        idx += 1
        x = MOI.get(model, MOI.VariableName(), MOI.VariableIndex(i.value))
        push!(single_variable_constraints, (x, MOI.Integer()))
    end
    for i in Iterators.flatten((
        MOI.get(
            model,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.Nonnegatives,
            }(),
        ),
        MOI.get(
            model,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.PositiveSemidefiniteConeTriangle,
            }(),
        ),
    ))
        idx += 1
        con_name_i = "c" * string(idx)
        push!(constraint_names, con_name_i)
        MOI.set(model, MOI.ConstraintName(), i, con_name_i)
    end

    return variable_names, constraint_names, single_variable_constraints
end

function _test_write_then_read(model_string::String)
    model1 = SDPA.Model()
    MOIU.loadfromstring!(model1, model_string)
    args = _set_var_and_con_names(model1)

    MOI.write_to_file(model1, SDPA_TEST_FILE)
    model2 = SDPA.Model()
    MOI.read_from_file(model2, SDPA_TEST_FILE)
    _set_var_and_con_names(model2)
    if MOI.get(model1, MOI.ObjectiveSense()) == MOI.MAX_SENSE
        MOI.set(model2, MOI.ObjectiveSense(), MOI.MAX_SENSE)
        attr = MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}()
        obj = MOI.get(model2, attr)
        MOI.set(model2, attr, MOIU.operate(-, Float64, obj))
    end

    return MOI.Test.util_test_models_equal(model1, model2, args...)
end

function _test_read(filename::String, model_string::String)
    model1 = MOI.FileFormats.Model(filename = filename)
    MOIU.loadfromstring!(model1, model_string)
    args = _set_var_and_con_names(model1)
    model2 = SDPA.Model()
    MOI.read_from_file(model2, filename)
    _set_var_and_con_names(model2)
    return MOI.Test.util_test_models_equal(model1, model2, args...)
end

function test_show()
    @test sprint(show, SDPA.Model()) ==
          "A SemiDefinite Programming Algorithm Format (SDPA) model"
end

function test_support()
    for set in [
        MOI.EqualTo(1.0),
        MOI.LessThan(1.0),
        MOI.GreaterThan(1.0),
        MOI.Interval(1.0, 2.0),
        MOI.Semiinteger(1.0, 2.0),
        MOI.Semicontinuous(1.0, 2.0),
        MOI.ZeroOne(),
    ]
        model_string = """
        variables: x
        minobjective: 1x
        x in $set
        """
        model = SDPA.Model()
        @test !MOI.supports_constraint(model, MOI.VariableIndex, typeof(set))
        err = MOI.UnsupportedConstraint{MOI.VariableIndex,typeof(set)}
        @test_throws err MOIU.loadfromstring!(model, model_string)
    end
end

function test_delete()
    for T in [Int, Float64]
        model = SDPA.Model(; number_type = T)
        x = MOI.add_variable(model)
        MOI.delete(model, x)
        y = MOI.add_variable(model)
        MOI.add_constraint(
            model,
            MOIU.vectorize([one(T) * y]),
            MOI.Nonnegatives(1),
        )
        err = ErrorException(
            "Non-contiguous variable indices not supported. This might be due to deleted variables.",
        )
        @test_throws err MOI.write_to_file(model, SDPA_TEST_FILE)
    end
end

function test_objective()
    for T in [Int, Float64]
        model = SDPA.Model(; number_type = T)
        @test !MOI.supports(model, MOI.ObjectiveFunction{MOI.VariableIndex}())
        @test !MOI.supports(
            model,
            MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{T}}(),
        )
    end
end

function test_nonempty()
    model = SDPA.Model()
    MOI.add_variable(model)
    err = ErrorException("Cannot read in file because model is not empty.")
    @test_throws err MOI.read_from_file(
        model,
        joinpath(SDPA_MODELS_DIR, "example_A.dat-s"),
    )
end

function test_bad_blocks()
    model = SDPA.Model()
    err = ErrorException(
        "The number of blocks (3) does not match the length of the list of blocks dimensions (2).",
    )
    @test_throws err MOI.read_from_file(
        model,
        joinpath(SDPA_MODELS_DIR, "bad_blocks.sdpa"),
    )
end

function test_bad_number_variables()
    model = SDPA.Model()
    err = ErrorException(
        "The number of variables (3) does not match the length of the list of coefficients for the objective function vector of coefficients (2).",
    )
    @test_throws err MOI.read_from_file(
        model,
        joinpath(SDPA_MODELS_DIR, "bad_vars.sdpa"),
    )
end

function test_wrong_number_of_values_in_entry()
    model = SDPA.Model()
    err = ErrorException(
        "Invalid line specifying entry: 0 1 2 2. There are 4 values instead of 5.",
    )
    @test_throws err MOI.read_from_file(
        model,
        joinpath(SDPA_MODELS_DIR, "bad_entry.sdpa"),
    )
end

function test_nondiagonal_entry()
    model = SDPA.Model()
    err = ErrorException(
        "Invalid line specifying entry: 0 1 1 2 1.0. `1 != 2` while block 1 has dimension 2 so it is a diagonal block.",
    )
    @test_throws err MOI.read_from_file(
        model,
        joinpath(SDPA_MODELS_DIR, "bad_diag.sdpa"),
    )
end

function test_nonzero_in_objective()
    model = SDPA.Model()
    MOIU.loadfromstring!(
        model,
        """
variables: x
minobjective: x + 1
""",
    )
    err = ErrorException(
        "Nonzero constant in objective function not supported. Note that " *
        "the constant may be added by the substitution of a bridged variable.",
    )
    @test_throws err MOI.write_to_file(model, SDPA_TEST_FILE)
end

function test_model_name()
    model = SDPA.Model()
    MOI.set(model, MOI.Name(), "FooBar")
    MOI.write_to_file(model, SDPA_TEST_FILE)
    @test readlines(SDPA_TEST_FILE) == ["\"FooBar", "0", "0", "", ""]
end

const _WRITE_READ_MODELS = [
    (
        "min ScalarAffine",
        """
    variables: x, y
    minobjective: 1.2x + -1y
""",
    ),
    (
        "max ScalarAffine",
        """
    variables: x, y
    maxobjective: 1.2x + -1y
""",
    ),
    (
        "VectorAffineFunction in Nonnegatives",
        """
    variables: x, y
    minobjective: 1.2x
    c1: [1.1 * x, y + 1] in Nonnegatives(2)
""",
    ),
    (
        "VectorAffineFunction in PositiveSemidefiniteConeTriangle",
        """
    variables: x, y, z
    minobjective: 1.2x
    c1: [1.1x, y + 1, 2x + z] in PositiveSemidefiniteConeTriangle(2)
""",
    ),
]

function test_write_read_models()
    for (model_name, model_string) in _WRITE_READ_MODELS
        _test_write_then_read(model_string)
    end
end

const _EXAMPLE_MODELS = [
    (
        "example_A.dat-s",
        """
    variables: x, y
    minobjective: 10x + 20y
    c1: [x + -1, 0, x + y + -2] in PositiveSemidefiniteConeTriangle(2)
    c2: [5y + -3, 2y, 6y + -4] in PositiveSemidefiniteConeTriangle(2)
""",
    ),
    (
        "example_B.sdpa",
        """
    variables: x
    minobjective: 1x
    c1: [0, 1x + -1, 0] in PositiveSemidefiniteConeTriangle(2)
""",
    ),
    (
        "example_integer.sdpa",
        """
    variables: x, y, z
    minobjective: 1x + -2y + -1z
    c1: [1x, 1y, 1z] in PositiveSemidefiniteConeTriangle(2)
    c2: [1z, 1x, 2.1] in PositiveSemidefiniteConeTriangle(2)
    c3: [1x + 1y + 1z + -1, -1x + -1y + -1z + 8] in Nonnegatives(2)
    x in Integer()
    y in Integer()
    z in Integer()
""",
    ),
]

function test_examples()
    for (model_name, model_string) in _EXAMPLE_MODELS
        _test_read(joinpath(SDPA_MODELS_DIR, model_name), model_string)
        _test_write_then_read(model_string)
    end
end

function runtests()
    for name in names(@__MODULE__, all = true)
        if startswith("$(name)", "test_")
            @testset "$name" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    sleep(1.0)  # Allow time for unlink to happen.
    rm(SDPA_TEST_FILE, force = true)
    return
end

end

TestSDPA.runtests()


text:
@testset "Utils" begin
    Random.seed!(0)
    n, r = 100, 10
    A = randn(n, r)
    A = A*A'

    # Check PSD
    @test RP.check_psd(A)
    @test_throws ArgumentError !RP.check_psd(A - 1e-1I)

    μ = 1e-2
    @test rank(A) - RP.deff(A, μ) < 2μ

end


text:
# Autogenerated wrapper script for HHsuite_jll for x86_64-linux-gnu-cxx03
export a3m_database_extract, a3m_database_filter, a3m_database_reduce, a3m_extract, a3m_reduce, cstranslate, ffindex_apply, ffindex_build, ffindex_from_fasta, ffindex_from_fasta_with_split, ffindex_get, ffindex_modify, ffindex_order, ffindex_reduce, ffindex_unpack, hhalign, hhalign_omp, hhblits, hhblits_ca3m, hhblits_omp, hhconsensus, hhfilter, hhmake, hhsearch, hhsearch_omp

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("HHsuite")
JLLWrappers.@declare_executable_product(a3m_database_extract)
JLLWrappers.@declare_executable_product(a3m_database_filter)
JLLWrappers.@declare_executable_product(a3m_database_reduce)
JLLWrappers.@declare_executable_product(a3m_extract)
JLLWrappers.@declare_executable_product(a3m_reduce)
JLLWrappers.@declare_executable_product(cstranslate)
JLLWrappers.@declare_executable_product(ffindex_apply)
JLLWrappers.@declare_executable_product(ffindex_build)
JLLWrappers.@declare_executable_product(ffindex_from_fasta)
JLLWrappers.@declare_executable_product(ffindex_from_fasta_with_split)
JLLWrappers.@declare_executable_product(ffindex_get)
JLLWrappers.@declare_executable_product(ffindex_modify)
JLLWrappers.@declare_executable_product(ffindex_order)
JLLWrappers.@declare_executable_product(ffindex_reduce)
JLLWrappers.@declare_executable_product(ffindex_unpack)
JLLWrappers.@declare_executable_product(hhalign)
JLLWrappers.@declare_executable_product(hhalign_omp)
JLLWrappers.@declare_executable_product(hhblits)
JLLWrappers.@declare_executable_product(hhblits_ca3m)
JLLWrappers.@declare_executable_product(hhblits_omp)
JLLWrappers.@declare_executable_product(hhconsensus)
JLLWrappers.@declare_executable_product(hhfilter)
JLLWrappers.@declare_executable_product(hhmake)
JLLWrappers.@declare_executable_product(hhsearch)
JLLWrappers.@declare_executable_product(hhsearch_omp)
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_executable_product(
        a3m_database_extract,
        "bin/a3m_database_extract",
    )

    JLLWrappers.@init_executable_product(
        a3m_database_filter,
        "bin/a3m_database_filter",
    )

    JLLWrappers.@init_executable_product(
        a3m_database_reduce,
        "bin/a3m_database_reduce",
    )

    JLLWrappers.@init_executable_product(
        a3m_extract,
        "bin/a3m_extract",
    )

    JLLWrappers.@init_executable_product(
        a3m_reduce,
        "bin/a3m_reduce",
    )

    JLLWrappers.@init_executable_product(
        cstranslate,
        "bin/cstranslate",
    )

    JLLWrappers.@init_executable_product(
        ffindex_apply,
        "bin/ffindex_apply",
    )

    JLLWrappers.@init_executable_product(
        ffindex_build,
        "bin/ffindex_build",
    )

    JLLWrappers.@init_executable_product(
        ffindex_from_fasta,
        "bin/ffindex_from_fasta",
    )

    JLLWrappers.@init_executable_product(
        ffindex_from_fasta_with_split,
        "bin/ffindex_from_fasta_with_split",
    )

    JLLWrappers.@init_executable_product(
        ffindex_get,
        "bin/ffindex_get",
    )

    JLLWrappers.@init_executable_product(
        ffindex_modify,
        "bin/ffindex_modify",
    )

    JLLWrappers.@init_executable_product(
        ffindex_order,
        "bin/ffindex_order",
    )

    JLLWrappers.@init_executable_product(
        ffindex_reduce,
        "bin/ffindex_reduce",
    )

    JLLWrappers.@init_executable_product(
        ffindex_unpack,
        "bin/ffindex_unpack",
    )

    JLLWrappers.@init_executable_product(
        hhalign,
        "bin/hhalign",
    )

    JLLWrappers.@init_executable_product(
        hhalign_omp,
        "bin/hhalign_omp",
    )

    JLLWrappers.@init_executable_product(
        hhblits,
        "bin/hhblits",
    )

    JLLWrappers.@init_executable_product(
        hhblits_ca3m,
        "bin/hhblits_ca3m",
    )

    JLLWrappers.@init_executable_product(
        hhblits_omp,
        "bin/hhblits_omp",
    )

    JLLWrappers.@init_executable_product(
        hhconsensus,
        "bin/hhconsensus",
    )

    JLLWrappers.@init_executable_product(
        hhfilter,
        "bin/hhfilter",
    )

    JLLWrappers.@init_executable_product(
        hhmake,
        "bin/hhmake",
    )

    JLLWrappers.@init_executable_product(
        hhsearch,
        "bin/hhsearch",
    )

    JLLWrappers.@init_executable_product(
        hhsearch_omp,
        "bin/hhsearch_omp",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""The integration dimension determines the in- and outputs of the transition function."""
abstract type IntegrationDimension end
const IntDim = IntegrationDimension
abstract type All <: IntDim end
abstract type Separable <: IntDim end
abstract type Separable_States <: IntDim end
abstract type Separable_ExogStates <: IntDim end
const Separable_Union = Union{Separable, Separable_States, Separable_ExogStates}

"""Transition function depending on parameters for dispatch."""
struct Transition{ID<:IntDim,nChoices,nShocks,nExogStates,F}
	f::F
end
(p::Transition)(args...) = p.f(args...)

"""Modify vStatesPrime in place."""
function get_transition!(transfunc::Transition{All}, vStatesPrime,
	vStates, vChoices, vShocks)
	vStatesPrime .= transfunc(vStates, vChoices, vShocks)
end

function get_transition!(transfunc::Transition{Separable,1}, vStatesPrime,
		vStates, vChoices, vShocks)
	vStatesPrime[1] = vChoices
	vStatesPrime[2:end] .= transfunc(vStates, vChoices, vShocks)
end
function get_transition!(transfunc::Transition{Separable,NC}, vStatesPrime,
		vStates, vChoices, vShocks) where NC
	vStatesPrime[1:NC] .= vChoices
	vStatesPrime[1+NC:end] .= transfunc(vStates, vChoices, vShocks)
end

function get_transition!(transfunc::Transition{Separable_States,1}, vStatesPrime,
		vStates, vChoices, vShocks)
	vStatesPrime[1] = vChoices
	vStatesPrime[2:end] .= transfunc(vStates, vShocks)
end
function get_transition!(transfunc::Transition{Separable_States,NC}, vStatesPrime,
		vStates, vChoices, vShocks) where NC
	vStatesPrime[1:NC] .= vChoices
	vStatesPrime[1+NC:end] .= transfunc(vStates, vShocks)
end

function get_transition!(transfunc::Transition{Separable_ExogStates,1,NSh,NSE}, vStatesPrime,
		vStates, vChoices, vShocks) where {NSh,NSE}
	vStatesPrime[1] = vChoices
	if NSE == 1
		vStatesPrime[2] = transfunc(vStates[2], vShocks)
	else
		vStatesPrime[2:end] .= transfunc(vStates[2:end], vShocks)
	end
end
function get_transition!(transfunc::Transition{Separable_ExogStates,NC,NSh,NSE}, vStatesPrime,
		vStates, vChoices, vShocks) where {NC,NSh,NSE}
	vStatesPrime[1:NC] .= vChoices
	if NSE == 1
		vStatesPrime[1+NC:end] .= transfunc(vStates[1+NC], vShocks)
	else
		vSim_it1[1+NC:end] .= transfunc(vSim_it[1+NC:end], mShocks_it)
	end
end


text:
using Base.Test

function solve(chars, offset::Int)
  sum = 0
  len = length(chars)
  for i in 1:len
    # 1-base indexing is a bit weird for this
    if chars[i] == chars[(i - 1 + offset)%len + 1]
      sum += parse(Int, chars[i])
    end
  end
  sum
end

@test solve("1122", 1) == 3
@test solve("1234", 1) == 0
@test solve("91212129", 1) == 9
@test solve("1212", 2) == 6
@test solve("1221", 2) == 0
@test solve("123425", 3) == 4
@test solve("123123", 3) == 12
@test solve("12131415", 4) == 4

open("/Users/mariosangiorgio/Downloads/input") do input
  chars = strip(readstring(input))
  println(solve(chars, 1))
  println(solve(chars, length(chars)/2))
end


text:
##### Discrete operators

export ∇!
export ∇²!

#####
##### cent <-> face
#####

function ∇!(face::CellFace{∇_dir}, cent::CellCenter, grid::Grid) where {∇_dir}
    iR = over_points(grid.c[∇_dir], Boundary(), Vertex1D())
    for (Ipre, i, Ipost) in sweep(cent, ∇_dir, iR)
        D = grid.c[∇_dir].deriv.c2v.diag[i-1]
        U = grid.c[∇_dir].deriv.c2v.upper[i-1]
        face[Ipre, i, Ipost] = D*cent[Ipre, i-1, Ipost] + U*cent[Ipre, i, Ipost]
    end
    nothing
end

function ∇!(cent::CellCenter, face::CellFace{∇_dir}, grid::Grid) where {∇_dir}
    for (Ipre, i, Ipost) in sweep(cent, ∇_dir)
        D = grid.c[∇_dir].deriv.v2c.diag[i]
        U = grid.c[∇_dir].deriv.v2c.upper[i]
        cent[Ipre, i, Ipost] = D*face[Ipre, i, Ipost] + U*face[Ipre, i+1, Ipost]
    end
    nothing
end

#####
##### corn <-> edge
#####

function ∇!(edge::CellEdge{∇_dir}, corn::CellCorner, grid::Grid) where {∇_dir}
    for (Ipre, i, Ipost) in sweep(edge, ∇_dir)
        D = grid.c[∇_dir].deriv.v2c.diag[i]
        U = grid.c[∇_dir].deriv.v2c.upper[i]
        edge[Ipre, i, Ipost] = D*corn[Ipre, i, Ipost] + U*corn[Ipre, i+1, Ipost]
    end
    nothing
end

function ∇!(corn::CellCorner, edge::CellEdge{∇_dir}, grid::Grid) where {∇_dir}
    iR = over_points(grid.c[∇_dir], Boundary(), Vertex1D())
    for (Ipre, i, Ipost) in sweep(edge, ∇_dir, iR)
        D = grid.c[∇_dir].deriv.c2v.diag[i-1]
        U = grid.c[∇_dir].deriv.c2v.upper[i-1]
        corn[Ipre, i, Ipost] = D*edge[Ipre, i-1, Ipost] + U*edge[Ipre, i, Ipost]
    end
    nothing
end

#####
##### face <-> edge
#####

∇!(edge::CellEdge, face::CellFace, grid::Grid, ∇_dir::Int) =
    ∇!(edge, face, grid, Val(∇_dir))

∇!(edge::CellEdge{dir}, face::CellFace{dir}, grid::Grid, ::Val{dir}) where {dir} =
    MethodError("∇(CellFace{dir}) does not live on CellEdge{dir}!")

function ∇!(edge::CellEdge{edgedir}, face::CellFace{facedir}, grid::Grid, ::Val{∇_dir}) where {edgedir,facedir,∇_dir}
    iR = over_points(grid.c[∇_dir], Boundary(), Vertex1D())
    for (Ipre, i, Ipost) in sweep(edge, ∇_dir, iR)
        D = grid.c[∇_dir].deriv.c2v.diag[i-1]
        U = grid.c[∇_dir].deriv.c2v.upper[i-1]
        edge[Ipre, i, Ipost] = D*face[Ipre, i-1, Ipost] + U*face[Ipre, i, Ipost]
    end
    nothing
end

#####
##### Generic forwarding
#####

VecTup(T) = Tuple{<:T,<:T,<:T}

function ∇!(∇f::VecTup(T∇F), f::TF, grid::Grid) where {T∇F,TF}
    ∇!(∇f[1], f, grid)
    ∇!(∇f[2], f, grid)
    ∇!(∇f[3], f, grid)
end

#####
##### ∇²!
#####

real_domain(::Center1D) = Interior()
real_domain(::Vertex1D) = Boundary()
is_cent(::Vertex1D) = 0
is_cent(::Center1D) = 1

# TODO: this can be much better optimized.
function ∇²!(∇²f::T, f::T, grid::Grid) where {T<:AbstractField}
    p_along = primary_along(f)
    d_along = dual_along(f)
    i_cent = is_cent.(d_along)

    @inbounds for k in over_points(grid.c[3], real_domain(p_along[3]), p_along[3])
    @inbounds for j in over_points(grid.c[2], real_domain(p_along[2]), p_along[2])
    @inbounds for i in over_points(grid.c[1], real_domain(p_along[1]), p_along[1])
        Δx⁺ = coords(grid.c[1], p_along[1]).Δh[i]
        Δy⁺ = coords(grid.c[2], p_along[2]).Δh[j]
        Δz⁺ = coords(grid.c[3], p_along[3]).Δh[k]
        Δx⁻ = coords(grid.c[1], p_along[1]).Δh[i-1]
        Δy⁻ = coords(grid.c[2], p_along[2]).Δh[j-1]
        Δz⁻ = coords(grid.c[3], p_along[3]).Δh[k-1]
        Δx = coords(grid.c[1], d_along[1]).Δh[i-i_cent[1]]
        Δy = coords(grid.c[2], d_along[2]).Δh[j-i_cent[2]]
        Δz = coords(grid.c[3], d_along[3]).Δh[k-i_cent[3]]
        f_ijk = f[i,j,k]
        ∇fx⁺ = (f[i+1,j,k]-f_ijk)/Δx⁺
        ∇fy⁺ = (f[i,j+1,k]-f_ijk)/Δy⁺
        ∇fz⁺ = (f[i,j,k+1]-f_ijk)/Δz⁺
        ∇fx⁻ = (f_ijk-f[i-1,j,k])/Δx⁻
        ∇fy⁻ = (f_ijk-f[i,j-1,k])/Δy⁻
        ∇fz⁻ = (f_ijk-f[i,j,k-1])/Δz⁻
        ∇²f_x = (∇fx⁺-∇fx⁻)/Δx
        ∇²f_y = (∇fy⁺-∇fy⁻)/Δy
        ∇²f_z = (∇fz⁺-∇fz⁻)/Δz
        ∇²f[i,j,k] = ∇²f_x+∇²f_y+∇²f_z
    end
    end
    end
    nothing
end


text:
#=
papers

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4788749/
Systems Biophysics: Multiscale Biophysical Modeling of Organ Systems

https://onlinelibrary.wiley.com/doi/full/10.1002/gamm.201900013
The dynamics of the skeletal muscle: A systems biophysics perspective on muscle modeling with the focus on Hill‐type muscle models

https://dev.biologists.org/content/develop/144/12/2104.full.pdf
Making muscle: skeletal myogenesis in vivo and in vitro

=#


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

using ISA, Plots
𝐶₀ = AMFMtriplet(t->exp(-t^2/5),t->200.0,0.0)
𝐶₁ = AMFMtriplet(t->1.0,t->100*t,0.1)
𝐶₂ = AMFMtriplet(t->0.8*cos(11t),t->100 + 70.5*sin(5t),π)
𝑆 = compSet([𝐶₀,𝐶₁,𝐶₂])
z = AMFMmodel(𝑆)
plot(z; timeaxis=-1.0:0.001:1.0)
png(path*"IS_exAMFMmodel1.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->exp(-t^2/5),t->200.0,0.0)
𝐶₁ = AMFMtriplet(t->1.0,t->100*t,0.1)
𝐶₂ = AMFMtriplet(t->0.8*cos(11t),t->100 + 70.5*sin(5t),π)
𝑆 = compSet([𝐶₀,𝐶₁,𝐶₂])
z = AMFMmodel(𝑆)
plot(𝑆; timeaxis=-1.0:0.001:1.0)
png(path*"IS_exAMFMmodel2.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->exp(-t^2/5),t->200.0,0.0)
𝐶₁ = AMFMtriplet(t->1.0,t->100*t,0.1)
𝐶₂ = AMFMtriplet(t->0.8*cos(11t),t->100 + 70.5*sin(5t),π)
𝑆 = compSet([𝐶₀,𝐶₁,𝐶₂])
z = AMFMmodel(𝑆)
plot(𝑆, timeaxis=-1.0:0.001:1.0, view="TF",
     left_margin=15Plots.mm, margin=5Plots.mm)
png(path*"IS_exAMFMmodel3.png")

#--------------------------------
using ISA, Plots
𝐶₀ = AMFMtriplet(t->exp(-t^2),t->100,0.0)
𝐶₁ = AMFMtriplet(t->1,t->25+20*sin(t),0.0)
𝐶₂ = AMFMtriplet(t->sin(11t),t->250,0.0)
𝑆 = compSet([𝐶₀,𝐶₁,𝐶₂])
z = AMFMmodel(𝑆)
plot(z; timeaxis=-1.0:0.001:1.0)
png(path*"IS_exAMFMmodel4.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->exp(-t^2),t->100,0.0)
𝐶₁ = AMFMtriplet(t->1,t->25+20*sin(t),0.0)
𝐶₂ = AMFMtriplet(t->sin(11t),t->250,0.0)
𝑆 = compSet([𝐶₀,𝐶₁,𝐶₂])
z = AMFMmodel(𝑆)
plot(𝑆; timeaxis=-1.0:0.001:1.0)
png(path*"IS_exAMFMmodel5.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->exp(-t^2),t->100,0.0)
𝐶₁ = AMFMtriplet(t->1,t->25+20*sin(t),0.0)
𝐶₂ = AMFMtriplet(t->sin(11t),t->250,0.0)
𝑆 = compSet([𝐶₀,𝐶₁,𝐶₂])
z = AMFMmodel(𝑆)
plot(𝑆, timeaxis=-1.0:0.001:1.0, view="TF",
     left_margin=15Plots.mm, margin=5Plots.mm)
png(path*"IS_exAMFMmodel6.png")

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

using ISA, Plots
𝐶₀ = AMFMtriplet(t->10,t->10t^2+2t+100,1.0)
𝐶₁ = AMFMtriplet(t->5exp(0.1t),t->250,0.0)
𝐶₂ = AMFMtriplet(t->10t,t->5cos(50t)+50,0.0)
𝑆 = compSet([𝐶₀,𝐶₁,𝐶₂])
z = AMFMmodel(𝑆)
plot(z; timeaxis=-1.0:0.001:1.0)
png(path*"IS_exAMFMmodel7.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->10,t->10t^2+2t+100,1.0)
𝐶₁ = AMFMtriplet(t->5exp(0.1t),t->250,0.0)
𝐶₂ = AMFMtriplet(t->10t,t->5cos(50t)+50,0.0)
𝑆 = compSet([𝐶₀,𝐶₁,𝐶₂])
z = AMFMmodel(𝑆)
plot(𝑆; timeaxis=-1.0:0.001:1.0)
png(path*"IS_exAMFMmodel8.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->10,t->10t^2+2t+100,1.0)
𝐶₁ = AMFMtriplet(t->5exp(0.1t),t->250,0.0)
𝐶₂ = AMFMtriplet(t->10t,t->5cos(50t)+50,0.0)
𝑆 = compSet([𝐶₀,𝐶₁,𝐶₂])
z = AMFMmodel(𝑆)
plot(𝑆, timeaxis=-1.0:0.001:1.0, view="TF",
     left_margin=15Plots.mm, margin=5Plots.mm)
png(path*"IS_exAMFMmodel9.png")


text:
"""
    localmaxmin!(x, maxes, mins)

Detect the local extrema of x.
Push the maxima into maxes and the minima into mins.
"""
function localmaxmin!(y, maxes::Vector{Int}, mins::Vector{Int})
    empty!(maxes)
    empty!(mins)
    for i in 2:length(y)-1
        if y[i+1] < y[i] > y[i-1]
            push!(maxes,i)
        elseif y[i+1] > y[i] < y[i-1]
            push!(mins, i)
        end
    end
end

"""
    zerocrossing!(y, crosses)

Compute the indices of zero crossings of vector `y`.
It searches for elements which are either zero or near a signflip and pushes the
indices into crosses.
"""
function zerocrossing!(y, crosses)
    empty!(crosses)
    for i ∈ 1:(length(y)-1)
        if y[i] == zero(y[i]) || sign(y[i]) * sign(y[i+1]) == -1
            push!(crosses, i)
        end
    end
end

"""
    get_edgepoint(y, xvec, extremas, pos, comp)

Compute the edgepoint which should be used as the extrema on the edge for the
spline computation.
"""
function get_edgepoint(y, xvec, extremas, pos, comp)
    if pos == first
        index = [1,2]
    elseif pos == last
        index = [length(extremas) - 1, length(extremas)]
    end
    # the x values must be embedded into a tuple
    knots = (xvec[extremas[index]],)
    itp = Interpolations.interpolate(knots, y[extremas[index]], Gridded(Linear()))
    expf = extrapolate(itp, Line())
    edgepoint = expf(pos(xvec))
    if comp(edgepoint, pos(y))
        edgepoint
    else
        edgepoint = pos(y)
    end
end

"""
    ismonotonic(x::AbstractVector)

Check whether `x` is monotonic.
"""
ismonotonic(x::AbstractVector) = issorted(x) || issorted(x, rev=true)

abstract type InterpMethod end
struct DierckXInterp <: InterpMethod end
#immutable InterpolationsInterp <: InterpMethod end

function interpolate(knotxvals::Vector, knotyvals::Vector, predictxvals::AbstractVector, m::DierckXInterp, k=3)
    spl = Dierckx.Spline1D(knotxvals, knotyvals, k=k)
    #@show spl, predictxvals
    Dierckx.evaluate(spl,predictxvals)
end

struct HaltingIterable{I, F}
    iter::I
    fun::F
end

function iterate(iter::HaltingIterable)
    next = iterate(iter.iter)
    return dispatch(iter, next)
end

function iterate(iter::HaltingIterable, (instruction, state))
    if instruction == :halt return nothing end
    next = iterate(iter.iter, state)
    return dispatch(iter, next)
end

function dispatch(iter::HaltingIterable, next)
    if next === nothing return nothing end
    return next[1], (iter.fun(next[1]) ? :halt : :continue, next[2])
end

halt(iter::I, fun::F) where {I,F} = HaltingIterable{I,F}(iter, fun)


text:
function NumericalPropagator(arg0::ODEIntegrator)
    return NumericalPropagator((ODEIntegrator,), arg0)
end

function NumericalPropagator(arg0::ODEIntegrator, arg1::AttitudeProvider)
    return NumericalPropagator((ODEIntegrator, AttitudeProvider), arg0, arg1)
end

function add_force_model(obj::NumericalPropagator, arg0::ForceModel)
    return jcall(obj, "addForceModel", void, (ForceModel,), arg0)
end

function get_all_force_models(obj::NumericalPropagator)
    return jcall(obj, "getAllForceModels", List, ())
end

function get_orbit_type(obj::NumericalPropagator)
    return jcall(obj, "getOrbitType", OrbitType, ())
end

function get_position_angle_type(obj::NumericalPropagator)
    return jcall(obj, "getPositionAngleType", PositionAngle, ())
end

function get_pv_coordinates(obj::NumericalPropagator, arg0::AbsoluteDate, arg1::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedPVCoordinates, (AbsoluteDate, Frame), arg0, arg1)
end

function propagate(obj::AbstractIntegratedPropagator, arg0::AbsoluteDate)
    return jcall(obj, "propagate", SpacecraftState, (AbsoluteDate,), arg0)
end

function remove_force_models(obj::NumericalPropagator)
    return jcall(obj, "removeForceModels", void, ())
end

function reset_initial_state(obj::NumericalPropagator, arg0::SpacecraftState)
    return jcall(obj, "resetInitialState", void, (SpacecraftState,), arg0)
end

function set_ignore_central_attraction(obj::NumericalPropagator, arg0::jboolean)
    return jcall(obj, "setIgnoreCentralAttraction", void, (jboolean,), arg0)
end

function set_initial_state(obj::NumericalPropagator, arg0::SpacecraftState)
    return jcall(obj, "setInitialState", void, (SpacecraftState,), arg0)
end

function set_mu(obj::NumericalPropagator, arg0::jdouble)
    return jcall(obj, "setMu", void, (jdouble,), arg0)
end

function set_orbit_type(obj::NumericalPropagator, arg0::OrbitType)
    return jcall(obj, "setOrbitType", void, (OrbitType,), arg0)
end

function set_position_angle_type(obj::NumericalPropagator, arg0::PositionAngle)
    return jcall(obj, "setPositionAngleType", void, (PositionAngle,), arg0)
end

function tolerances(::Type{NumericalPropagator}, arg0::jdouble, arg1::AbsolutePVCoordinates)
    return jcall(NumericalPropagator, "tolerances", Vector{Vector{jdouble}}, (jdouble, AbsolutePVCoordinates), arg0, arg1)
end

function tolerances(::Type{NumericalPropagator}, arg0::jdouble, arg1::Orbit, arg2::OrbitType)
    return jcall(NumericalPropagator, "tolerances", Vector{Vector{jdouble}}, (jdouble, Orbit, OrbitType), arg0, arg1, arg2)
end


text:
using Random
using Statistics
using StatsBase
using Plots
import Base.+
import Base.-

include("Point_Struct.jl")
include("fun.jl")
include("constants.jl")

# This script reproduce the Figure 3 in the Errard's paper

r = 0:20000

punti = [Point{Float64}(i, π / 2 , 0., undef, undef, undef)()   for i in r] #remember the functor!() else it is a cartesian poin
z = [i.z for i in punti]
w = [ω(o) for o in punti]
T = [Tₚ(o) for o in punti]
X = [χ₂(o,o) for o in punti]

plot( w[2:length(r)] ./ w[2]  , z[2:length(r)], yaxis=:log, xlims=(0,2), ylims = (1,2E4))
plot!( T[2:length(r)] ./ T[2], z[2:length(r)], yaxis=:log, xlims=(0,2), ylims = (1,2E4))
plot!( X[2:length(r)] ./ X[2], z[2:length(r)], yaxis=:log, xlims=(0,2), ylims = (1,2E4))

max_stars_repo_name: gruffoc/AFO

max_issues_repo_name: gruffoc/AFO

max_forks_repo_name: gruffoc/AFO


text:
"Nested dictionary to store indexed clauses."
const ClauseSubtable{T} = Dict{Symbol,Vector{T}} where  T <: AbstractClause
const ClauseTable{T} = Dict{Symbol,ClauseSubtable{T}}

function insert_clause!(table::ClauseTable{T}, c::Clause) where {T <: AbstractClause}
    subtable = get!(table, c.head.name, Dict{Symbol,Vector{T}}())
    if isa(c.head, Compound) && length(c.head.args) >= 1
        arg = c.head.args[1]
        if isa(arg, Var)
            push!(get!(subtable, :__var__, Vector{T}()), c)
        else
            push!(get!(subtable, Symbol(arg.name), Vector{T}()), c)
        end
        push!(get!(subtable, :__all__, Vector{T}()), c)
    else
        push!(get!(subtable, :__no_args__, Vector{T}()), c)
    end
end

"Insert clauses into indexed table for efficient look-up."
function insert_clauses!(table::ClauseTable{T}, clauses::Vector{T})  where {T <: AbstractClause}
    # Ensure no duplicates are added
    clauses = unique(clauses)
    if (length(table) > 0) setdiff!(clauses, deindex_clauses(table)) end
    # Iterate over clauses and insert into table
    for c in clauses
        insert_clause!(table, c)
    end
    return table
end

"Insert clauses into indexed table and return a new table."
function insert_clauses(table::ClauseTable{T}, clauses::Vector{T})  where {T <: AbstractClause}
    return insert_clauses!(deepcopy(table), clauses)
end

"Index clauses by functor name and first argument for efficient look-up."
function index_clauses(clauses::Vector{T})  where {T <: AbstractClause}
    return insert_clauses!(ClauseTable{T}(), clauses)
end

"Convert indexed clause table to flat list of clauses."
function deindex_clauses(table::ClauseTable{T})  where {T <: AbstractClause}
    clauses = Vector{T}()
    for (functor, subtable) in table
        if :__no_args__ in keys(subtable)
            append!(clauses, subtable[:__no_args__])
        else
            append!(clauses, subtable[:__all__])
        end
    end
    return clauses
end

"Retrieve matching clauses from indexed clause table."
function retrieve_clauses(table::ClauseTable{T}, term::Term, funcs::Dict=Dict())  where {T <: AbstractClause}
    clauses = Vector{T}()
    funcs = length(funcs) > 0 ? merge(default_funcs, funcs) : default_funcs
    subtable = get(table, term.name, nothing)
    if subtable === nothing
        return clauses
    end
    if isa(term, Compound) && length(term.args) >= 1
        arg = term.args[1]
        if isa(arg, Var) || arg.name in keys(funcs)
            clauses = get(subtable, :__all__, clauses)
        else
            clauses = [get(subtable, Symbol(arg.name), clauses);
                       get(subtable, :__var__, clauses)]
        end
    else
        clauses = get(subtable, :__no_args__, clauses)
    end
    return clauses
end

"Subtract one clause table from another (in-place)."
function subtract_clauses!(table1::ClauseTable{T}, table2::ClauseTable{T})  where {T <: AbstractClause}
    for (functor, subtable2) in table2
        if !(functor in keys(table1)) continue end
        subtable1 = table1[functor]
        for arg in keys(subtable2)
            if !(arg in keys(subtable1)) continue end
            setdiff!(subtable1[arg], subtable2[arg])
        end
    end
    return table1
end

"Subtract one clause table from another (returns new copy)."
function subtract_clauses(table1::ClauseTable{T}, table2::ClauseTable{T})  where {T <: AbstractClause}
    return subtract_clauses!(deepcopy(table1), table2)
end

"Subtract clauses from a indexed clause table (in-place)."
function subtract_clauses!(table::ClauseTable{T}, clauses::Vector{T})  where {T <: AbstractClause}
    return subtract_clauses!(table, index_clauses(clauses))
end

"Subtract clauses from a indexed clause table (returns new copy)."
function subtract_clauses(table::ClauseTable{T}, clauses::Vector{T})  where {T <: AbstractClause}
    return subtract_clauses(table, index_clauses(clauses))
end

"Return number of clauses in indexed clause table."
function num_clauses(table::ClauseTable{T})  where {T <: AbstractClause}
    n = 0
    for (functor, subtable) in table
        if :__no_args__ in keys(subtable)
            n += length(subtable[:__no_args__])
        else
            n += length(subtable[:__all__])
        end
    end
    return n
end


text:
"""
```
Network(id)
```

Construct an empty `Network` representing a system of interacting particles. A particle is represented by a `Species`, and an interaction is represented by a `Reaction`.

Add a `Species` or `Reaction` using `<=`:

```
m <= Species("X", 100)
m <= Reaction("birth", 2.0, "X --> X + X")
```
"""
struct Network
  identifier :: Symbol

  species_list  :: OrderedDict{Symbol,Species}
  reaction_list :: OrderedDict{Symbol,Reaction}
  dep_graph :: OrderedDict{Symbol,Vector{Symbol}}

  function Network(id)
    s  = OrderedDict{Symbol,Species}()
    r  = OrderedDict{Symbol,Reaction}()
    dep_graph = OrderedDict{Symbol,Vector{Symbol}}()

    return new(Symbol(id), s, r, dep_graph)
  end
end

"""
```
species_list(x::Network)
```

Retrieve `Species` from a `Network`.
"""
species_list(x::Network) = x.species_list

"""
```
reaction_list(x::Network)
```

Retrieve `Reaction`s from a `Network`.
"""
reaction_list(x::Network) = x.reaction_list

get_species(m::Network, key::AbstractString) = get_species(m, Symbol(key))

function get_species(m::Network, key::Symbol)
  s = get(m.species_list, key, nothing)

  if isnothing(s)
    error("""
    Species $(key) not defined.
    Add it with `m <= Species("$(key)")`.
    """)
  end

  return s
end

get_reaction(m::Network, key::AbstractString) = get_reaction(m, Symbol(key))

function get_reaction(m::Network, key::Symbol)
  r = get(m.reaction_list, key, nothing)

  if isnothing(r)
    error("""
    Reaction $(key) not defined.
    Add it with `m <= Reaction("$(key)", rate, formula)`.
    """)
  end

  return r
end

number_species(x::Network)   = length(species_list(x))
number_reactions(x::Network) = length(reaction_list(x))

function Base.show(io::IO, x::Network)
  print(io, "[ Model: $(x.identifier) ]\n")
  print(io, "  no. species:    $(number_species(x))\n")
  print(io, "  no. reactions:  $(number_reactions(x))")
end

function add_object!(model, object, fieldname)
  dict = getfield(model, fieldname)
  id   = object.identifier
  setindex!(dict, object, Symbol(id))
  return model
end

function (<=)(model::Network, object::Species)
  add_object!(model, object, :species_list)
end

function (<=)(model::Network, object::Reaction)
  reactants = object.reactants
  products  = object.products

  try
    validate(reactants, model)
    validate(products,  model)

    # add the reaction to the list of objects
    add_object!(model, object, :reaction_list)

    id = Symbol(object.identifier)
    dep_graph = model.dep_graph

    if haskey(dep_graph, id)
      # if the key already exists, empty the list of dependents
      empty!(dep_graph[id])

      # and remove the reaction as a dependency
      for (key, r) in reaction_list(model)
        list = dep_graph[key]
        ix = findfirst(list, id)
        ix > 0 && deleteat!(list, ix)
      end
    else
      # if the key does not exist, make a new vector for dependents
      dep_graph[id] = Symbol[]
    end

    # update the dependency graph
    update_dep_graph!(model, object)

  catch ex
    rethrow(ex)
  end
end

function validate(participants, model)
  species_dict = species_list(model)
  for species in keys(participants)
    if !haskey(species_dict, Symbol(species))
      error("$(species) is not defined.")
    end
  end
end

function rmv_object!(model, key, fieldname)
  dict = getfield(model, fieldname)
  if haskey(dict, key)
    delete!(dict, key)
  else
    error("$(key) not found in $(fieldname)")
  end
  return model
end

function update_dep_graph!(model, rj)
  id = Symbol(rj.identifier)

  rj_affects    = rj.affects
  rj_affectedby = rj.affectedby

  for (key, rk) in reaction_list(model)
    rk_affects    = rk.affects
    rk_affectedby = rk.affectedby

    # add edge R_j --> R_k
    if !isempty(intersect(rj_affectedby, rk_affects))
      push!(model.dep_graph[id], key)
    end

    # add edge R_k --> R_j
    if rj != rk && !isempty(intersect(rk_affectedby, rj_affects))
      push!(model.dep_graph[key], id)
    end
  end

  return model
end


text:
module ch14

using Printf
# -------------------
# Reading and Writing
# -------------------
fout = open("output.txt", "w")
line = "This here's the wattle,\n"
write(fout, line)
line = "the emblem of our land.\n"
write(fout, line)
close(fout)

# ----------
# Formatting
# ----------
fout = open("output.txt", "w")
write(fout, string(150))
camels = 42
#println(fout, "I have spotted $camels camels.")
@printf(fout, "I have spotted %i camels", camels)
close(fout)

# Filenames and Paths
cwd = pwd()
abspath("output.txt")
ispath("output.txt")
isdir("output.txt")
isfile("output.txt")
readdir(cwd)
# joinpath
# walkdir

# Command Objects
cmd = `echo hello`
run(cmd)
ret = read(cmd, String)

# repr()
# dump()
end


text:
# PNG, SVG, matrix import

# PNG

"""
    readpng(pathname)

Read a PNG file.

This returns a image object suitable for placing on the current drawing with `placeimage()`.
You can access its `width` and `height` fields:

    image = readpng("/tmp/test-image.png")
    w = image.width
    h = image.height
"""
function readpng(pathname)
    im = Cairo.read_from_png(pathname)
    if iszero(im.width) || iszero(im.height)
        throw(error("can't read this PNG image from $(pathname). Either it's not a valid PNG file, or the format is different from what Cairo is expecting."))
    end
    return im
end

function paint_with_alpha(ctx::Cairo.CairoContext, a = 0.5)
    Cairo.paint_with_alpha(get_current_cr(), a)
end

"""
    placeimage(pngimg, pos=O; centered=false)
    placeimage(pngimg, xpos, ypos; centered=false)

Place the PNG image on the drawing at `pos`, or (`xpos`/`ypos`). The image `img` has been previously
read using `readpng()`.

Use keyword `centered=true` to place the center of the image at the position.
"""
function placeimage(img::Cairo.CairoSurface, xpos, ypos; centered=false)
    if centered == true
        w, h = img.width, img.height
        xpos, ypos = xpos - (w/2), ypos - (h/2)
    end
    Cairo.set_source_surface(get_current_cr(), img, xpos, ypos)
    # no alpha
    Cairo.paint(get_current_cr())
end

placeimage(img::Cairo.CairoSurface, pt::Point=O; kwargs...) = placeimage(img, pt.x, pt.y; kwargs...)

"""
    placeimage(img, pt::Point=O, alpha; centered=false)
    placeimage(pngimg, xpos, ypos, alpha; centered=false)

Place a PNG image `pngimg` on the drawing at `pt` or `Point(xpos, ypos)`
with opacity/transparency `alpha`. The image has been
previously loaded using `readpng()`.

Use keyword `centered=true` to place the center of the image at the position.
"""
function placeimage(img::Cairo.CairoSurface, xpos, ypos, alpha; centered=false)
    if centered == true
        w, h = img.width, img.height
        xpos, ypos = xpos - (w/2), ypos - (h/2)
    end
    Cairo.set_source_surface(get_current_cr(), img, xpos, ypos)
    paint_with_alpha(get_current_cr(), alpha)
end

placeimage(img::Cairo.CairoSurface, pt::Point, alpha; kwargs...) =
  placeimage(img::Cairo.CairoSurface, pt.x, pt.y, alpha; kwargs...)

# SVG

struct SVGimage
    im::RsvgHandle
    xpos::Float64
    ypos::Float64
    width::Float64
    height::Float64
end

function _readsvgfile(fname)
    if Base.stat(fname).size == 0
         throw(error("readsvg: file $fname not found"))
    end
    r = Rsvg.handle_new_from_file(fname)

    if r.ptr == C_NULL
        throw(error("readsvg: couldn't read the file $fname"))
    end

    d = Rsvg.handle_get_dimensions(r)
    if iszero(d.width) || iszero(d.height)
        throw(error("readsvg(): can't get dimensions of this SVG image in $(pathname). Perhaps it's not a valid SVG file."))
    end

    return Luxor.SVGimage(r, d.em, d.ex, d.width, d.height)
end

function _readsvgstring(str)
    r = Rsvg.handle_new_from_data(str)

    if r.ptr == C_NULL
        throw(error("readsvg: can't read SVG from string starting $(str[1:20])."))
    end

    d = Rsvg.handle_get_dimensions(r)
    if iszero(d.width) || iszero(d.height)
        throw(error("readsvg_string(): can't get dimensions of this SVG image. Perhaps it's not valid SVG."))
    end
    return SVGimage(r, d.em, d.ex, d.width, d.height)
end

"""
    readsvg(str)

Read an SVG image. `str` is either pathname or pure SVG
code. This returns an SVG image object suitable for placing
on the current drawing with `placeimage()`.

Placing an SVG file:

```
@draw begin
    mycoollogo = readsvg("/tmp/mylogo.svg")
    placeimage(mycoollogo)
end
```

Placing SVG code:

```
# from https://github.com/edent/SuperTinyIcons
julialogocode = \"\"\"<svg xmlns="http://www.w3.org/2000/svg"
    aria-label="Julia" role="img"
    viewBox="0 0 512 512">
    <rect width="512" height="512" rx="15%" fill="#fff"/>
    <circle fill="#389826" cx="256" cy="137" r="83"/>
    <circle fill="#cb3c33" cx="145" cy="329" r="83"/>
    <circle fill="#9558b2" cx="367" cy="329" r="83"/>
</svg>\"\"\"

@draw begin
    julia_logo = readsvg(julialogocode)
    placeimage(julia_logo, centered=true)
end
```
"""
function readsvg(str)
    # str is either pathname or pure SVG code
    # unfortunately ispath fails on Mac if the string is longer than 255 characters
     if length(str) > 255 # don't check ispath, assume SVG string
         _readsvgstring(str)
     elseif ispath(str) # check is a file
         _readsvgfile(str)
     else
         _readsvgstring(str)
     end
 end

"""
    placeimage(svgimg, pos=O; centered=false)

Place an SVG image stored in `svgimg` on the drawing at `pos`. Use `readsvg()`
to read an SVG image from file, or from SVG code.

Use keyword `centered=true` to place the center of the image
at the position.
"""
function placeimage(im::SVGimage, pos=O; centered=false)
    if centered == true
        w, h = im.width, im.height
        pos = pos - ((w/2), (h/2))
    end
    @layer begin
        translate(pos)
        Rsvg.handle_render_cairo(Luxor.get_current_cr(), im.im)
    end
end

"""
    placeimage(matrix::AbstractMatrix{UInt32}, pos=O;
        alpha=1, centered=false)

Place an image matrix on the drawing at `pos` with opacity/transparency `alpha`.

Use keyword `centered=true` to place the center of the image
at the position.
"""
function placeimage(buffer::AbstractMatrix{UInt32}, pt=O;
        alpha=1, centered=false)
    if centered == true
        w, h = size(buffer)
        pt = Point(pt.x - (w/2), pt.y - (h/2))
    end
    Cairo.set_source_surface(Luxor.get_current_cr(), Cairo.CairoImageSurface(buffer, Cairo.FORMAT_ARGB32), pt.x, pt.y)
    paint_with_alpha(get_current_cr(), alpha)
end
placeimage(buffer::AbstractMatrix{ARGB32}, args...; kargs...) = placeimage(collect(reinterpret(UInt32, buffer)), args...; kargs...)
placeimage(buffer::AbstractMatrix{<:Colorant}, args...; kargs...) = placeimage(convert.(ARGB32, buffer), args...; kargs...)

function placeimage(buffer::Drawing, args...;
        kargs...)
    if buffer.surfacetype == :svg
        placeimage(_readsvgstring(String(copy(buffer.bufferdata))), args...; kargs...)
    else
        throw(error("surfacetype `$(buffer.surfacetype)` is not supported. Use `image` instead."))
    end
end


text:
############
# TimeData #
############
mutable struct TimeData
    ncalls::Int
    time::Int64
    allocs::Int64
    firstexec::Int64
end
TimeData(ncalls, time, allocs) = TimeData(ncalls, time, allocs, time)
Base.copy(td::TimeData) = TimeData(td.ncalls, td.time, td.allocs)
TimeData() = TimeData(0, 0, 0, time_ns())

function Base.:+(self::TimeData, other::TimeData)
    TimeData(self.ncalls + other.ncalls,
             self.time + other.time,
             self.allocs + other.allocs,
             min(self.firstexec, other.firstexec))
end

###############
# TimerOutput #
###############
mutable struct TimerOutput
    start_data::TimeData
    accumulated_data::TimeData
    inner_timers::Dict{String,TimerOutput}
    timer_stack::Vector{TimerOutput}
    name::String
    flattened::Bool
    enabled::Bool
    totmeasured::Tuple{Int64,Int64}
    prev_timer_label::String
    prev_timer::Union{TimerOutput,Nothing}

    function TimerOutput(label::String = "root")
        start_data = TimeData(0, time_ns(), gc_bytes())
        accumulated_data = TimeData()
        inner_timers = Dict{String,TimerOutput}()
        timer_stack = TimerOutput[]
        return new(start_data, accumulated_data, inner_timers, timer_stack, label, false, true, (0, 0), "", nothing)
    end

    # Jeez...
    TimerOutput(start_data, accumulated_data, inner_timers, timer_stack, name, flattened, enabled, totmeasured, prev_timer_label,
    prev_timer) = new(start_data, accumulated_data, inner_timers, timer_stack, name, flattened, enabled, totmeasured, prev_timer_label,
    prev_timer)

end

Base.copy(to::TimerOutput) = TimerOutput(copy(to.start_data), copy(to.accumulated_data), copy(to.inner_timers),
                                         copy(to.timer_stack), to.name, to.flattened, to.enabled, to.totmeasured, "", nothing)

const DEFAULT_TIMER = TimerOutput()

# push! and pop!
function Base.push!(to::TimerOutput, label::String)
    if length(to.timer_stack) == 0 # Root section
        current_timer = to
    else # Not a root section
        current_timer = to.timer_stack[end]
    end
    # Fast path
    if current_timer.prev_timer_label == label
        timer = current_timer.prev_timer
    else
        maybe_timer = get(current_timer.inner_timers, label, nothing)
        # this could be implemented more elegant using
        # get!(() -> TimerOutput(label), current_timer.inner_timers, label)
        # however that causes lots of allocations in
        # julia v1.3
        if maybe_timer === nothing
            timer = TimerOutput(label)
            current_timer.inner_timers[label] = timer
        else
            timer = maybe_timer
        end
    end
    current_timer.prev_timer_label = label
    current_timer.prev_timer = timer

    push!(to.timer_stack, timer)
    return timer.accumulated_data
end

Base.pop!(to::TimerOutput) = pop!(to.timer_stack)

# Only sum the highest parents
function totmeasured(to::TimerOutput)
    t, b = Int64(0), Int64(0)
    for section in values(to.inner_timers)
        timedata = section.accumulated_data
        t += timedata.time
        b += timedata.allocs
    end
    return t, b
end

function longest_name(to::TimerOutput, indent = 0)
    m = textwidth(to.name) + indent
    for inner_timer in values(to.inner_timers)
        m = max(m, longest_name(inner_timer, indent + 2))
    end
    return m
end

# merging timer outputs
const merge_lock = ReentrantLock() # needed for merges of objects on different threads

Base.merge(others::TimerOutput...) = merge!(TimerOutput(), others...)
function Base.merge!(self::TimerOutput, others::TimerOutput...; tree_point = String[])
    lock(merge_lock) do
        for other in others
            self.accumulated_data += other.accumulated_data
            its = self.inner_timers
            for point in tree_point
                its = its[point].inner_timers
            end
            _merge(its, other.inner_timers)
        end
        return self
    end
end
function _merge(self::Dict{String,TimerOutput}, other::Dict{String,TimerOutput})
    for key in keys(other)
        if haskey(self, key)
            self[key].accumulated_data += other[key].accumulated_data
            _merge(self[key].inner_timers, other[key].inner_timers)
        else
            self[key] = deepcopy(other[key])
        end
    end
end

#######
# API #
#######

# Accessors
ncalls(to::TimerOutput)    = to.accumulated_data.ncalls
allocated(to::TimerOutput) = to.accumulated_data.allocs
time(to::TimerOutput) = to.accumulated_data.time
totallocated(to::TimerOutput) = totmeasured(to)[2]
tottime(to::TimerOutput) = totmeasured(to)[1]

time() = time(DEFAULT_TIMER)
ncalls() = ncalls(DEFAULT_TIMER)
allocated() = allocated(DEFAULT_TIMER)
totallocated() = totmeasured(DEFAULT_TIMER)[2]
tottime() = totmeasured(DEFAULT_TIMER)[1]

get_defaulttimer() = DEFAULT_TIMER
Base.@deprecate get_defaultimer get_defaulttimer

# Macro
macro timeit(args...)
    return timer_expr(__module__, false, args...)
end

macro timeit_debug(args...)
    if !isdefined(__module__, :timeit_debug_enabled)
        Core.eval(__module__, :(timeit_debug_enabled() = false))
    end

    return timer_expr(__module__, true, args...)
end

function enable_debug_timings(m::Module)
    if !getfield(m, :timeit_debug_enabled)()
        Core.eval(m, :(timeit_debug_enabled() = true))
    end
end
function disable_debug_timings(m::Module)
    if getfield(m, :timeit_debug_enabled)()
        Core.eval(m, :(timeit_debug_enabled() = false))
    end
end

timer_expr(args...) = throw(ArgumentError("invalid macro usage for @timeit, use as @timeit [to] label codeblock"))

function is_func_def(f)
    if isa(f, Expr) && (f.head === :function || Base.is_short_function_def(f))
        return true
    else
        return false
    end
end

function timer_expr(m::Module, is_debug::Bool, ex::Expr)
    is_func_def(ex) && return timer_expr_func(m, is_debug, :($(TimerOutputs.DEFAULT_TIMER)), ex)
    return esc(_timer_expr(m, is_debug, :($(TimerOutputs).DEFAULT_TIMER), ex))
end

function timer_expr(m::Module, is_debug::Bool, label_or_to, ex::Expr)
    is_func_def(ex) && return timer_expr_func(m, is_debug, label_or_to, ex)
    return esc(_timer_expr(m, is_debug, :($(TimerOutputs).DEFAULT_TIMER), label_or_to, ex))
end

function timer_expr(m::Module, is_debug::Bool, label::String, ex::Expr)
    is_func_def(ex) && return timer_expr_func(m, is_debug, :($(TimerOutputs).DEFAULT_TIMER), ex, label)
    return esc(_timer_expr(m, is_debug, :($(TimerOutputs).DEFAULT_TIMER), label, ex))
end

function timer_expr(m::Module, is_debug::Bool, to, label, ex::Expr)
    is_func_def(ex) && return timer_expr_func(m, is_debug, to, ex, label)
    return esc(_timer_expr(m, is_debug, to, label, ex))
end

function _timer_expr(m::Module, is_debug::Bool, to::Union{Symbol, Expr, TimerOutput}, label, ex::Expr)
    @gensym local_to enabled accumulated_data b₀ t₀ val
    timeit_block = quote
        $local_to = $to
        $enabled = $local_to.enabled
        if $enabled
            $accumulated_data = $(push!)($local_to, $label)
        end
        $b₀ = $(gc_bytes)()
        $t₀ = $(time_ns)()
        $(Expr(:tryfinally,
            :($val = $ex),
            quote
                if $enabled
                    $(do_accumulate!)($accumulated_data, $t₀, $b₀)
                    $(pop!)($local_to)
                end
            end))
        $val
    end

    if is_debug
        return quote
            if $m.timeit_debug_enabled()
                $timeit_block
            else
                $ex
            end
        end
    else
        return timeit_block
    end
end

function timer_expr_func(m::Module, is_debug::Bool, to, expr::Expr, label=nothing)
    expr = macroexpand(m, expr)
    def = splitdef(expr)

    label === nothing && (label = string(def[:name]))

    def[:body] = if is_debug
        quote
            @inline function inner()
                $(def[:body])
            end
            $(_timer_expr(m, is_debug, to, label, :(inner())))
        end
    else
        _timer_expr(m, is_debug, to, label, def[:body])
    end

    return esc(combinedef(def))
end

function do_accumulate!(accumulated_data, t₀, b₀)
    accumulated_data.time += time_ns() - t₀
    accumulated_data.allocs += gc_bytes() - b₀
    accumulated_data.ncalls += 1
end

reset_timer!() = reset_timer!(DEFAULT_TIMER)
function reset_timer!(to::TimerOutput)
    to.inner_timers = Dict{String,TimerOutput}()
    to.start_data = TimeData(0, time_ns(), gc_bytes())
    to.accumulated_data = TimeData()
    to.prev_timer_label = ""
    to.prev_timer = nothing
    resize!(to.timer_stack, 0)
    return to
end

# We can remove this now that the @timeit macro is exception safe.
# Doesn't hurt to keep it for a while though
timeit(f::Function, label::String) = timeit(f, DEFAULT_TIMER, label)
function timeit(f::Function, to::TimerOutput, label::String)
    accumulated_data = push!(to, label)
    b₀ = gc_bytes()
    t₀ = time_ns()
    local val
    try
        val = f()
    finally
        accumulated_data.time += time_ns() - t₀
        accumulated_data.allocs += gc_bytes() - b₀
        accumulated_data.ncalls += 1
        pop!(to)
    end
    return val
end

Base.haskey(to::TimerOutput, name::String) = haskey(to.inner_timers, name)
Base.getindex(to::TimerOutput, name::String) = to.inner_timers[name]

function flatten(to::TimerOutput)
    t, b = totmeasured(to)
    inner_timers = Dict{String,TimerOutput}()
    for inner_timer in values(to.inner_timers)
        _flatten!(inner_timer, inner_timers)
    end
    toc = copy(to)
    return TimerOutput(toc.start_data, toc.accumulated_data, inner_timers, TimerOutput[], "Flattened", true, true, (t, b), "", to)
end

function _flatten!(to::TimerOutput, inner_timers::Dict{String,TimerOutput})
    for inner_timer in values(to.inner_timers)
        _flatten!(inner_timer, inner_timers)
    end

    if haskey(inner_timers, to.name)
        timer = inner_timers[to.name]
        timer.accumulated_data += to.accumulated_data
    else
        toc = copy(to)
        toc.inner_timers = Dict{String,TimerOutput}()
        inner_timers[toc.name] = toc
    end
end

enable_timer!(to::TimerOutput=DEFAULT_TIMER) = to.enabled = true
disable_timer!(to::TimerOutput=DEFAULT_TIMER) = to.enabled = false

# Macro to selectively disable timer for expression
macro notimeit(args...)
    notimeit_expr(args...)
end

# Default function throws an error for the benefit of the user
notimeit_expr(args...) = throw(ArgumentError("invalid macro usage for @notimeit, use as @notimeit [to] codeblock"))

complement!() = complement!(DEFAULT_TIMER)
function complement!(to::TimerOutput)
    if length(to.inner_timers) == 0
        return nothing
    end
    tot_time = to.accumulated_data.time
    tot_allocs = to.accumulated_data.allocs
    for timer in values(to.inner_timers)
        tot_time -= timer.accumulated_data.time
        tot_allocs -= timer.accumulated_data.allocs
        complement!(timer)
    end
    tot_time = max(tot_time, 0)
    tot_allocs = max(tot_allocs, 0)
    if !(to.name in ["root", "Flattened"])
        name = string("~", to.name, "~")
        timer = TimerOutput(to.start_data, TimeData(max(1,to.accumulated_data.ncalls), tot_time, tot_allocs), Dict{String,TimerOutput}(), TimerOutput[], name, false, true, (tot_time, tot_allocs), to.name, to)
        to.inner_timers[name] = timer
    end
    return to
end

# If @notimeit was called without a TimerOutput instance, use default timer
notimeit_expr(ex::Expr) = notimeit_expr(:($(TimerOutputs.DEFAULT_TIMER)), ex)

# Disable timer, evaluate expression, restore timer to previous value, and return expression result
function notimeit_expr(to, ex::Expr)
    return quote
        local to = $(esc(to))
        local enabled = to.enabled
        $(disable_timer!)(to)
        local val
        $(Expr(:tryfinally,
            :(val = $(esc(ex))),
            quote
                if enabled
                    $(enable_timer!)(to)
                end
            end))
        val
    end
end


text:
include("parameter.jl")
include("generate.jl")
include("standard.jl")


text:
InputFolder = './Images/EyeVesselsSpots/';
OutputFolder = './Results/Images/EyeVesselsSpotsFilam/';
Lbl = 1;
Fill = -1;

@iA = '*.tif';		% Image filter

@fxg_kIsoScanFilam [iA] > [M1];
params.GRad = 1;
params.NAngles = 8;
params.Len = 7;
params.Contrast1 = 25;
params.Contrast2 = 12;
params.Anis = 0.35;
params.MaxHoleArea = 3;
/endf

@fxkg_kSklClean [M1, iA] > [M2];
params.MinBrchLgth = 3;
params.SearchRad = 9;
params.MinMean = 30;
params.MinArea = 25;
/endf

/show iA > M2;
/keep M2 > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
struct DummyPotential <: AbstractPotential; end

#function addforce!(dt::Float64,system::AbstractSystem,model::AbstractPotential)
#  system.x .+= dt .* model.f
#end


text:
using LinearAlgebra
using SpecialFunctions: erfc
using StaticArrays
using Test
using Zygote
using FiniteDiff

const Mat3{T} = SMatrix{3, 3, T, 9} where T
const Vec3{T} = SVector{3, T} where T
const MVec3{T} = MVector{3, T} where T

"""
A modified ewald summation suitable for AD with Zygote.
Original docstring below.

Compute the electrostatic interaction energy per unit cell between point
charges in a uniform background of compensating charge to yield net
neutrality. the `lattice` and `recip_lattice` should contain the
lattice and reciprocal lattice vectors as columns. `charges` and
`positions` are the point charges and their positions (as an array of
arrays) in fractional coordinates.
"""
function energy_ewald_zygote_faster(lattice, charges, positions; η=nothing)
    T = eltype(lattice)

    for i=1:3
        if norm(lattice[:,i]) == 0
            ## TODO should something more clever be done here? For now
            ## we assume that we are not interested in the Ewald
            ## energy of non-3D systems
            return T(0)
        end
    end
    energy_ewald_zygote_faster(lattice, T(2π) * inv(lattice'), charges, positions; η=η)
end

function energy_ewald_zygote_faster(lattice, recip_lattice, charges, positions; η=nothing)
    T = eltype(lattice)
    @assert T == eltype(recip_lattice)
    @assert length(charges) == length(positions)
    if η === nothing
        # Balance between reciprocal summation and real-space summation
        # with a slight bias towards reciprocal summation
        η = sqrt(sqrt(T(1.69) * norm(recip_lattice ./ 2T(π)) / norm(lattice))) / 2
    end

    #
    # Numerical cutoffs
    #
    # The largest argument to the exp(-x) function to obtain a numerically
    # meaningful contribution. The +5 is for safety.
    max_exponent = -log(eps(T)) + 5

    # The largest argument to the erfc function for various precisions.
    # To get an idea:
    #   erfc(5) ≈ 1e-12,  erfc(8) ≈ 1e-29,  erfc(10) ≈ 2e-45,  erfc(14) ≈ 3e-87
    max_erfc_arg = get(
        Dict(Float32 => 5, Float64 => 8, BigFloat => 14),
        T,
        something(findfirst(arg -> 100 * erfc(arg) < eps(T), 1:100), 100) # fallback for not yet implemented cutoffs
    )

    #
    # Reciprocal space sum
    #
    # Initialize reciprocal sum with correction term for charge neutrality
    sum_recip = - (sum(charges)^2 / 4η^2)

    # Function to return the indices corresponding
    # to a particular shell
    # TODO switch to an O(N) implementation
    function shell_indices(ish)
        (Vec3(i,j,k) for i in -ish:ish for j in -ish:ish for k in -ish:ish
        if maximum(abs.((i,j,k))) == ish)
    end

    # Loop over reciprocal-space shells
    gsh = 1 # Exclude G == 0
    any_term_contributes = true
    while any_term_contributes
        any_term_contributes = false

        # Compute G vectors and moduli squared for this shell patch
        for G in shell_indices(gsh)
            Gsq = sum(abs2, recip_lattice * Zygote.dropgrad(G))

            # Check if the Gaussian exponent is small enough
            # for this term to contribute to the reciprocal sum
            exponent = Gsq / 4η^2
            if exponent > max_exponent
                continue
            end

            cos_strucfac = sum(Z * cos(2T(π) * dot(r, Zygote.dropgrad(G))) for (r, Z) in zip(positions, charges))
            sin_strucfac = sum(Z * sin(2T(π) * dot(r, Zygote.dropgrad(G))) for (r, Z) in zip(positions, charges))
            sum_strucfac = cos_strucfac^2 + sin_strucfac^2

            any_term_contributes = true
            sum_recip += sum_strucfac * exp(-exponent) / Gsq

        end
        gsh += 1
    end
    # Amend sum_recip by proper scaling factors:
    sum_recip *= 4T(π) / abs(det(lattice))

    #
    # Real-space sum
    #
    # Initialize real-space sum with correction term for uniform background
    sum_real = Zygote.dropgrad(-2η / sqrt(T(π)) * sum(abs2, charges))

    # Loop over real-space shells
    rsh = 0 # Include R = 0
    any_term_contributes = true
    while any_term_contributes || rsh <= 1
        any_term_contributes = false

        # Loop over R vectors for this shell patch
        for R in shell_indices(rsh)
            for i = 1:length(positions), j = 1:length(positions)
                ti = positions[i]
                Zi = charges[i]
                tj = positions[j]
                Zj = charges[j]

                # Avoid self-interaction
                if rsh == 0 && ti == tj
                    continue
                end

                dist = norm(lattice * (ti - tj - Zygote.dropgrad(R)))

                # erfc decays very quickly, so cut off at some point
                if η * dist > max_erfc_arg
                    continue
                end

                any_term_contributes = true
                sum_real += Zi * Zj * erfc(η * dist) / dist
            end # i,j
        end # R
        rsh += 1
    end
    energy = (sum_recip + sum_real) / 2  # Divide by 2 (because of double counting)
    energy
end

#=======#

function compute_forces_zygote(positions, lattice, charges)
    _positions = reduce(hcat, positions)
    forces = first(
        Zygote.gradient(_positions -> begin
                positions = collect(eachcol(_positions))
                -energy_ewald_zygote_faster(lattice, charges, positions)
            end,
            _positions
        )
    )
    return forces
end

function compute_forces_finitediff(positions, lattice, charges)
    _positions = reduce(hcat, positions)
    forces = FiniteDiff.finite_difference_gradient(_positions -> begin
            positions = collect(eachcol(_positions))
            -energy_ewald_zygote_faster(lattice, charges, positions)
        end,
        _positions
    )
    return forces
end

#===#

function benchmark_ewald_zygote(num_atoms=20)
    lattice = Mat3(
        [0.0  5.131570667152971 5.131570667152971;
        5.131570667152971 0.0 5.131570667152971;
        5.131570667152971 5.131570667152971  0.0]
    )
    # perturb positions away from equilibrium to get nonzero force
    positions = [ones(3)/8+rand(3)/20 for _ in 1:num_atoms]
    charges = [14 for _ in 1:num_atoms]

    forces_finitediff = compute_forces_finitediff(positions, lattice, charges)
    forces_zygote = compute_forces_zygote(positions, lattice, charges)

    # check consistency of forces
    # @test forces_zygote ≈ forces_finitediff atol=1e-6
    @show sum(abs.(forces_finitediff - forces_zygote)) / (3num_atoms)

    println("Timings:")
    println("Energy without forces")
    @time energy_ewald_zygote(lattice, charges, positions)
    println("FiniteDiff forces")
    @time compute_forces_finitediff(positions, lattice, charges)
    println("Zygote forces")
    @time compute_forces_zygote(positions, lattice, charges)

    return nothing
end

# benchmark_ewald_zygote(2)
    # Timings:
    # Energy without forces
    #   0.004692 seconds (113.77 k allocations: 4.967 MiB)
    # FiniteDiff forces
    #   0.075556 seconds (1.37 M allocations: 59.604 MiB)
    # Zygote forces
    #   0.318643 seconds (1.03 M allocations: 48.197 MiB, 11.16% gc time)

# benchmark_ewald_zygote(10)
    # Timings:
    # Energy without forces
    #   0.285002 seconds (2.46 M allocations: 105.066 MiB, 57.46% gc time)
    # FiniteDiff forces
    #   8.808988 seconds (147.72 M allocations: 6.156 GiB, 18.05% gc time)
    # Zygote forces
    #   4.530250 seconds (15.03 M allocations: 717.221 MiB, 13.84% gc time)

# Zygote-rewrite required us to
# 1. do not use try-catch
# 2. use Zygote.@ignore on shell_indices
# 3. use Zygote.dropgrad on shell int vectors

const lattice = Mat3(
    [0.0  5.131570667152971 5.131570667152971;
    5.131570667152971 0.0 5.131570667152971;
    5.131570667152971 5.131570667152971  0.0]
)
# perturb positions away from equilibrium to get nonzero force
const positions = [Vec3(ones(3)/8+rand(3)/20) for _ in 1:2]
const charges = [14 for _ in 1:2]

@time energy_ewald_zygote_faster(lattice, charges, positions) # 0.319368 seconds (470.33 k allocations: 24.824 MiB, 6.28% gc time, 99.67% compilation time)

using BenchmarkTools
# @btime energy_ewald_zygote_faster($lattice, $charges, $positions) # 555.540 μs (3378 allocations: 53.33 KiB)

# @time energy_ewald(lattice, charges, positions)
# @time energy_ewald(lattice, charges, positions)

#========#
# A larger physical system with pymatgen
using DFTK

function benchmark_ewald_zygote_pymatgen(num_atoms=20)
    # lattice = Mat3(
    #     [0.0  5.131570667152971 5.131570667152971;
    #     5.131570667152971 0.0 5.131570667152971;
    #     5.131570667152971 5.131570667152971  0.0]
    # )
    # # perturb positions away from equilibrium to get nonzero force
    # positions = [ones(3)/8+rand(3)/20 for _ in 1:num_atoms]
    # charges = [14 for _ in 1:num_atoms]

    a = 10.263141334305942  # Lattice constant in Bohr
    lattice = a / 2 .* [[0 1 1.]; [1 0 1.]; [1 1 0.]]
    Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
    atoms = [Si => [ones(3)/8, -ones(3)/8]];
    pystruct = pymatgen_structure(lattice, atoms)
    pystruct.make_supercell([1, 1, 1])
    lattice = load_lattice(pystruct)
    positions = [s.frac_coords for s in pystruct.sites];
    charges = [14 for _ in 1:length(positions)]

    # forces_finitediff = compute_forces_finitediff(positions, lattice, charges)
    forces_zygote = compute_forces_zygote(positions, lattice, charges)

    # check consistency of forces
    # @test forces_zygote ≈ forces_finitediff atol=1e-6
    # @show sum(abs.(forces_finitediff - forces_zygote)) / (3num_atoms)

    println("Timings:")
    println("Energy without forces")
    @time energy_ewald_zygote(lattice, charges, positions)
    println("FiniteDiff forces")
    # @time compute_forces_finitediff(positions, lattice, charges)
    println("Zygote forces")
    @time compute_forces_zygote(positions, lattice, charges)

    return nothing
end

benchmark_ewald_zygote_pymatgen()

max_stars_repo_name: niklasschmitz/DifferentiableDFTK

max_issues_repo_name: niklasschmitz/DifferentiableDFTK

max_forks_repo_name: niklasschmitz/DifferentiableDFTK


text:
module MinFEM

import WriteVTK
using SparseArrays, LinearAlgebra

export PDESystem,
       solve,
       refresh,
       evaluateMeshFunction,
       asmBoundarySource,
       getDirichletProjection,
       asmLaplacian,
       asmMassMatrix,
       asmBoundaryMassMatrix,
       asmCubicSecondDerivativeMatrix,
       asmCubicDerivativeMatrix,
       asmCubicTerm,
       asmElasticity,
       asmDirichletCondition,
       asmIncompressibility,
       asmVectorLaplacian,
       L2norm,
       asmGradient,
       computeGradient,
       Mesh,
       import_mesh,
       unit_square,
       getCellVolumes,
       open_vtk_file,
       save_vtk_file,
       write_point_data,
       write_cell_data

include("fe_types.jl")
include("mesh_methods.jl")
include("fe_methods.jl")

end # module


text:
# Compute first num_components basis vectors for PCA, after subtracting projection onto fixed_comp
function fit_gen_pca_rv{T}(X::Array{T,2}, fixed_comp::Array{T,1}; quad_term::Array{T,1} = zeros(T,0), num_components::Integer=4, tol::Float64=1e-12, max_it::Int64 = 10)
  num_lambda = size(X,1)
  num_spectra = size(X,2)
  M = rand(T, (num_lambda,num_components) )  # random initialization is part of algorithm (i.e., not zeros)
  s = zeros(T,num_lambda)                    # pre-allocated memory for compute_pca_component
  scores = zeros(num_spectra,num_components)
  mu = vec(mean(X,2))
  Xtmp = X.-mu                               # perform PCA after subtracting off mean
  totalvar = sum(abs2,Xtmp)
  fracvar = zeros(num_components)
  M[:,1] = fixed_comp                        # Force fixed (i.e., Doppler) component to replace first PCA component
  fixed_comp_norm = 1/sum(abs2,fixed_comp)
  for i in 1:num_spectra
    scores[i,1] = z = (dot(view(Xtmp,:,i),fixed_comp)*fixed_comp_norm)  # Normalize differently, so scores are z (i.e., doppler shift)
	  Xtmp[:,i] -= z*fixed_comp
    if length(quad_term) == length(fixed_comp)  # Optionally subtract off quadratic term from Doppler shift
      Xtmp[:,i] -= 0.5*z*z*quad_term   # subtract off projection onto fixed component and quadratic term separately
    end
  end
  fracvar[1] = sum(abs2,Xtmp)/totalvar
  #println("# sum(abs2,Xtmp) = ", sum(abs2,Xtmp) )
  #println("# j = ", 1, " sum(abs2,Xtmp) = ", sum(abs2,Xtmp), " frac_var_remain= ", fracvar[1] )
  for j in 2:num_components
    compute_pca_component!(Xtmp,view(M,:,j),s,tol=tol,max_it=max_it)
    for i in 1:num_spectra
       scores[i,j] = dot(view(Xtmp,:,i),view(M,:,j)) #/sum(abs2,view(M,:,j-1))
       Xtmp[:,i] .-= scores[i,j]*view(M,:,j)
	  end
    fracvar[j] = sum(abs2,Xtmp)/totalvar
    println("# j = ", j, " sum(abs2,Xtmp) = ", sum(abs2,Xtmp), " frac_var_remain= ", fracvar[j] )
  end
  return (mu,M,scores)
end

function est_rvs_from_pc{T}(X::Array{T,2}, mu::Vector{T}, doppler_comp::Vector{T})
  num_lambda = size(X,1)
  num_spectra = size(X,2)
  speed_of_light = 3.0e8
  rv = zeros(num_spectra)
  Xtmp = X.-mu
  fixed_comp_norm = 1/sum(abs2,doppler_comp)
  for i in 1:num_spectra
      rv[i] = dot(Xtmp[:,i],doppler_comp)*fixed_comp_norm*speed_of_light
  end
  rv
end

max_stars_repo_name: eford/RVSpectraKitLearn

max_issues_repo_name: eford/RVSpectraKitLearn

max_forks_repo_name: eford/RVSpectraKitLearn


text:
include("runtests.jl")
include("parallel_matrix.jl")
include("parallel_mult.jl")
include("parallel_latent_full.jl")
include("parallel_latent_tensor.jl")
include("heavy_copyto.jl")


text:
struct Figure3a
    Hs::Vector{Float64}
    dates::Vector{Float64}
    function Figure3a(waveseries,wavedates)
        @assert length(waveseries) == length(wavedates)
        new(4*std.(waveseries),wavedates)
    end
end

struct Figure3b
    speed::Vector{Float64}
    dates::Vector{Float64}
end

struct Figure3c
    power::Matrix{Float64}
    Ω::Vector{Float64}
    dates::Vector{Float64}
    function Figure3c(waveseries,wavedates)
        spectrogram, Ω, specdates = PaperAnalysis.compute_spectrogram(waveseries,wavedates)
        new(spectrogram, Ω, specdates)
    end
end

struct Figure3d
    direction::Matrix{Float64}
    Ω::Vector{Float64}
    dates::Vector{Float64}
    function Figure3d(waveseries,wavedates)
        wavedir = PaperAnalysis.mt_wavedirection(permutedims(reduce(vcat, wave for wave in waveseries)),2304,2304÷2)
        specdates = wavedates[1]:(wavedates[2]-wavedates[1])/2:wavedates[end]
        Ω = collect((2π*1.28).*wavedir.freq)
        dir = wavedir.power
        new(mod2pi.(dir),Ω,specdates)
    end
end

struct Figure3e
    direction::Vector{Float64}
    dates::Vector{Float64}
    function Figure3e(direction,dates)
        @assert length(direction) == length(dates)
        new(mod2pi_nan.(deg2rad.(direction)), dates)
    end
end

struct Figure3
    fig3a::Figure3a
    fig3b::Figure3b
    fig3c::Figure3c
    fig3d::Figure3d
    fig3e::Figure3e
    function Figure3()
        @load "dat/anonomous_data.jld"
        new(
            Figure3a(waveseries, wavetime./24 ./3600),
            Figure3b(windspeed, windtime./24 ./3600),
            Figure3c(waveseries, wavetime./24 ./3600),
            Figure3d(waveseries, wavetime./24 ./3600),
            Figure3e(winddirection, windtime./24 ./3600)
        )
    end
end

## plotting
@recipe function f(f::Figure3a)
    yguide := L"$H_s$ (m)"
    @series begin
        seriestype := :heatmap
        xlims := (f.dates[1],f.dates[end])
        xticks := (1:5,map(i->"day $i",1:5))
        f.dates[end] .+ [2/24,3/24], 2:3, (x,y) -> 1
    end
    @series begin
        label := false
        seriescolor := 1
        xlims := (f.dates[1],f.dates[end])
        xticks := (1:5,map(i->"day $i",1:5))
        ylims := (1.3,5.2)
        f.dates, f.Hs
    end
end
@recipe function f(f::Figure3b)
    yguide := "wind speed (m/s)"
    @series begin
        seriestype := :heatmap
        xlims := (f.dates[1],f.dates[end])
        xticks := (1:5,map(i->"day $i",1:5))
        f.dates[end] .+ [2/24,3/24], 2:3, (x,y) -> 1
    end
    @series begin
        label := false
        xlims := (f.dates[1],f.dates[end])
        xticks := (1:5,map(i->"day $i",1:5))
        seriescolor := 2
        f.dates, f.speed
    end
end
@recipe function f(f::Figure3c)
    yguide := "angular frequency"
    seriestype := :heatmap
    xticks := (1:5,map(i->"day $i",1:5))
    f.dates, f.Ω, 10.0 .* log10.(f.power)
end
@recipe function f(f::Figure3d)
    yguide := "angular frequency"
    seriestype := :heatmap
    seriescolor := :hsv
    xticks := (1:5,map(i->"day $i",1:5))
    f.dates, f.Ω, f.direction
end

@recipe function f(f::Figure3e)
    @series begin
        seriestype := :heatmap
        seriescolor := :hsv
        f.dates, range(0,2π,length=100),(x,y)->y
    end
    @series begin
        yguide := "wind direction"
        label := false
        seriescolor := :black
        xlims := (f.dates[1], f.dates[end])
        xticks := (1:5,map(i->"day $i",1:5))
        f.dates, f.direction
    end
end

@recipe function f(f::Figure3)
    layout := (fieldcount(typeof(f)),1)
    for (i,field) in enumerate(fieldnames(typeof(f)))
        @series begin
            subplot := i
            getfield(f,field)
        end
    end
end


text:
# Used in the version of LineSearches that was interfaced.
# Those types are depricated in the current state of LineSearches

# A cache for results from line search methods (to avoid recomputation)
mutable struct LineSearchResults{T}
    alpha::Vector{T}
    value::Vector{T}
    slope::Vector{T}
    nfailures::Int
end

LineSearchResults(::Type{T}) where {T} = LineSearchResults(T[], T[], T[], 0)

Base.length(lsr::LineSearchResults) = length(lsr.alpha)

function Base.push!(lsr::LineSearchResults{T}, a, v, d) where T
    push!(lsr.alpha, convert(T, a))
    push!(lsr.value, convert(T, v))
    push!(lsr.slope, convert(T, d))
    return
end

function clear!(lsr::LineSearchResults)
    empty!(lsr.alpha)
    empty!(lsr.value)
    empty!(lsr.slope)
    return
    # nfailures is deliberately not set to 0
end

mutable struct LineSearchException{T<:Real} <: Exception
    message::AbstractString
    alpha::T
    lsr::LineSearchResults
end


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

"""
    CompositeTask(tasks, [adjacency])

A learning task composed of multiple sub-`tasks`.
The tasks are assumed to be independent, unless
an `adjacency` matrix is provided with the list
of dependencies.

## Examples

A composite task composed of three independent sub-tasks
`t₁`, `t₂` and `t₃`:

```julia
julia> CompositeTask([t₁,t₂,t₃])
```

A composite sequential task `t₁ → t₂ → t₃`:

```julia
julia> CompositeTask([t₁,t₂,t₃], [0 1 0; 0 0 1; 0 0 0])
```

### Notes

Only directed acyclic graphs (DAGs) are permitted.
"""
struct CompositeTask <: AbstractLearningTask
  tasks::Vector{AbstractLearningTask}
  adjacency::Matrix{Int}

  function CompositeTask(tasks, adjacency)
    @assert all(length(tasks) .== size(adjacency)) "invalid adjacency for tasks"
    @assert isacyclic(adjacency) "dependencies must be acyclic"
    new(tasks, adjacency)
  end
end

CompositeTask(tasks::AbstractVector) =
  CompositeTask(tasks, zeros(Int, length(tasks), length(tasks)))

iscomposite(task::CompositeTask) = true