text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

const Mbprec =
    SMatrix{6,6}(+0.9999256795,      -0.0111814828,      -0.0048590040,
                 -0.000551,          -0.238560,          +0.435730,
                 +0.0111814828,      +0.9999374849,      -0.0000271557,
                 +0.238509,          -0.002667,          -0.008541,
                 +0.0048590039,      -0.0000271771,      +0.9999881946,
                 -0.435614,          +0.012254,          +0.002117,
                 -0.00000242389840,  +0.00000002710544,  +0.00000001177742,
                 +0.99990432,        -0.01118145,        -0.00485852,
                 -0.00000002710544,  -0.00000242392702,  +0.00000000006585,
                 +0.01118145,        +0.99991613,        -0.00002716,
                 -0.00000001177742,  +0.00000000006585,  -0.00000242404995,
                 +0.00485852,        -0.00002717,        +0.99996684)

# Note: IDL version of `bprecess' changes in-place "muradec", "parallax" and
# "radvel".  We don't do anything like this, but calculations are below,
# commented, in case someone is interested.
function _bprecess(ra::T, dec::T, parallax::T, radvel::T,
                   epoch::T, muradec::Vector{T}) where {T<:AbstractFloat}
    @assert length(muradec) == 2
    sinra,  cosra  = sincos(deg2rad(ra))
    sindec, cosdec = sincos(deg2rad(dec))
    r0 = SVector(cosra*cosdec,  sinra*cosdec,  sindec)
    r0_dot = SVector(-muradec[1]*sinra*cosdec - muradec[2]*cosra*sindec,
                     muradec[1]*cosra*cosdec - muradec[2]*sinra*sindec,
                     muradec[2]*cosdec) .+ 21.095 .* radvel * parallax * r0
    R_1 = Mbprec * vcat(r0, r0_dot)
    r1 = R_1[1:3]
    if isfinite(epoch)
        r1_dot = @view R_1[4:6]
        r1 .+= deg2rad.(r1_dot .* (epoch .- 1950) ./ 360000)
        A = A_precess .+ deg2rad.(A_dot_precess .* (epoch .- 1950) ./ 360000)
    else
        A = A_precess
    end
    rmag = norm(r1)
    s1 = r1 ./ rmag
    r = Array{T}(undef, 3)
    s = copy(s1)
    for j = 0:2
        r .= s1 .+ A .- dot(s, A) .* s
        s .= r ./ rmag
    end
    rmag = norm(r)
    # s1_dot = r1_dot ./ rmag
    # r_dot = s1_dot + A_dot_precess - sum(s.*A_dot_precess)*s
    x = r[1]
    y = r[2]
    z = r[3]
    # x_dot = r_dot[1]
    # y_dot = r_dot[2]
    # z_dot = r_dot[3]
    # muradec[1] = (x*y_dot - y*x_dot)/(x*x + y*y)
    # muradec[2] = (z_dot*(x*x + y*y) - z*(x*x_dot + y*y_dot))/(rmag*rmag*hypot(x, y))
    dec1950 = asin(z / rmag)
    ra1950  = atan(y, x)
    # if parallax > 0
    #     radvel = (x*x_dot + y*y_dot + z*z_dot)/(21.095*parallax*rmag)
    #     parallax = parallax / rmag
    # end
    if ra1950 < 0
        ra1950 += 2 * T(pi)
    end
    ra1950  = rad2deg(ra1950)
    dec1950 = rad2deg(dec1950)
    return ra1950, dec1950
end

# Main interface.
bprecess(ra::Real, dec::Real, muradec::Vector{<:Real}; parallax::Real=0.0, radvel::Real=0.0) =
    _bprecess(promote(float(ra), float(dec), float(parallax), float(radvel), NaN)...,
              float(muradec))

bprecess(ra::Real, dec::Real, epoch::Real=2000.0) =
    _bprecess(promote(float(ra), float(dec), 0.0, 0.0, float(epoch))...,
              zeros(typeof(float(ra)), 2))

# Tuple arguments.
bprecess(radec::Tuple{Real,Real}, muradec::Vector{<:Real};
         parallax::Real=0.0, radvel::Real=0.0) =
             bprecess(radec..., muradec, parallax=parallax, radvel=radvel)

bprecess(radec::Tuple{Real,Real}, epoch::Real=2000.0) =
    bprecess(radec..., epoch)

# Vectorial arguments.
function bprecess(ra::AbstractArray{R}, dec::AbstractArray{<:Real},
                  muradec::AbstractArray{<:Real};
                  parallax::AbstractArray{<:Real}=zeros(R, length(ra)),
                  radvel::AbstractArray{<:Real}=zeros(R, length(ra))) where {R<:Real}
    @assert length(ra) == length(dec) == size(muradec)[2] == length(parallax) == length(radvel)
    typer = float(R)
    ra1950  = similar(ra, typer)
    dec1950 = similar(dec, typer)
    for i in eachindex(ra)
        ra1950[i], dec1950[i] = bprecess(ra[i], dec[i], muradec[:,i],
                                         parallax=parallax[i], radvel=radvel[i])
    end
    return ra1950, dec1950
end

function bprecess(ra::AbstractArray{R}, dec::AbstractArray{D},
                  epoch::Real=2000.0) where {R<:Real,D<:Real}
    @assert length(ra) == length(dec)
    typer = float(R)
    ra1950  = similar(ra, typer)
    dec1950 = similar(dec, typer)
    for i in eachindex(ra)
        ra1950[i], dec1950[i] = bprecess(ra[i], dec[i], epoch[i])
    end
    return ra1950, dec1950
end

"""
    bprecess(ra, dec[, epoch]) -> ra1950, dec1950
    bprecess(ra, dec, muradec[, parallax=parallax, radvel=radvel]) -> ra1950, dec1950

### Purpose ###

Precess positions from J2000.0 (FK5) to B1950.0 (FK4).

### Explanation ###

Calculates the mean place of a star at B1950.0 on the FK4 system from the mean
place at J2000.0 on the FK5 system.

`bprecess` function has two methods, one for each of the following cases:

* the proper motion is known and non-zero
* the proper motion is unknown or known to be exactly zero (i.e. extragalactic
  radio sources).  Better precision can be achieved in this case by inputting
  the epoch of the original observations.

### Arguments ###

The function has 2 methods.  The common mandatory arguments are:

* `ra`: input J2000 right ascension, in degrees.
* `dec`: input J2000 declination, in degrees.

The two methods have a different third argument (see "Explanation" section for
more details).  It can be one of the following:

* `muradec`: 2-element vector containing the proper motion in seconds of arc
  per tropical *century* in right ascension and declination.
* `epoch`: scalar giving epoch of original observations.

If none of these two arguments is provided (so `bprecess` is fed only with right
ascension and declination), it is assumed that proper motion is exactly zero and
`epoch = 2000`.

If it is used the method involving `muradec` argument, the following keywords
are available:

* `parallax` (optional numerical keyword): stellar parallax, in seconds of arc.
* `radvel` (optional numerical keyword): radial velocity in km/s.

Right ascension and declination can be passed as the 2-tuple `(ra, dec)`.  You
can also pass `ra`, `dec`, `parallax`, and `radvel` as arrays, all of the same
length N.  In that case, `muradec` should be a matrix 2×N.

### Output ###

The 2-tuple of right ascension and declination in 1950, in degrees, of input
coordinates is returned.  If `ra` and `dec` (and other possible optional
arguments) are arrays, the 2-tuple of arrays `(ra1950, dec1950)` of the same
length as the input coordinates is returned.

### Method ###

The algorithm is taken from the Explanatory Supplement to the Astronomical
Almanac 1992, page 186.  See also Aoki et al (1983), A&A, 128, 263.  URL:
http://adsabs.harvard.edu/abs/1983A%26A...128..263A.

### Example ###

The SAO2000 catalogue gives the J2000 position and proper motion for the star HD
119288.  Find the B1950 position.

* RA(2000) = 13h 42m 12.740s
* Dec(2000) = 8d 23' 17.69''
* Mu(RA) = -.0257 s/yr
* Mu(Dec) = -.090 ''/yr

```jldoctest
julia> using AstroLib

julia> muradec = 100*[-15*0.0257, -0.090]; # convert to century proper motion

julia> ra = ten(13, 42, 12.74)*15;

julia> decl = ten(8, 23, 17.69);

julia> adstring(bprecess(ra, decl, muradec), precision=2)
" 13 39 44.526  +08 38 28.63"
```

### Notes ###

"When transferring individual observations, as opposed to catalog mean place,
the safest method is to transform the observations back to the epoch of the
observation, on the FK4 system (or in the system that was used to to produce the
observed mean place), convert to the FK5 system, and transform to the the epoch
and equinox of J2000.0" -- from the Explanatory Supplement (1992), p. 180

`jprecess` performs the precession to J2000 coordinates.

Code of this function is based on IDL Astronomy User's Library.
"""
bprecess


text:
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

using NLPModels
using OptimizationProblems
using Stopping
using LineSearch
using LinearAlgebra
using NLPModelsJuMP
using JuMP

# algorithms which reach a solution on that problem
# algorithms_good = [:armijo_ls, :Biss_Cub_ls, :Biss_SecA_ls, :Biss_Sec_ls, :Biss_Nwt_ls, :TR_Nwt_ls, :TR_Sec_ls, :TR_SecA_ls, :ARC_Nwt_ls]
algorithms_good = [:hagerzhang]

for algo in algorithms_good
   println(" ")
   println("$algo")
   nlp = MathProgNLPModel(arwhead());
   h = LineModel(nlp, nlp.meta.x0, -grad(nlp,nlp.meta.x0));
   lsatx = LSAtT(0.0, h₀ = obj(h, 0.0), g₀ = grad(h, 0.0),
                      ht = obj(h, 0.0), gt = grad(h, 0.0))
   stop_meta = StoppingMeta(max_iter = 50)
   ls_meta = LS_Function_Meta()
   stop_ls = LS_Stopping(h, (x,y)-> armijo(x, y), lsatx, meta = stop_meta);

   state, optimality = eval(algo)(h, stop_ls, f_meta = ls_meta, verboseLS = true)
   @test optimality == true
   reset!(nlp)
end

# algorithms which cannot reach a solution on that problem
# algorithms_bad = [:ARC_Sec_ls, :ARC_SecA_ls]
#
# for algo in algorithms_bad
#    println(" ")
#    println("$algo")
#    nlp = MathProgNLPModel(arwhead());
#    h = LineModel(nlp, nlp.meta.x0, -grad(nlp,nlp.meta.x0));
#    lsatx = LSAtT(0.0, h₀ = obj(h, 0.0), g₀ = grad(h, 0.0),
#                       ht = obj(h, 0.0), gt = grad(h, 0.0))
#    stop_ls = LS_Stopping(h, (x,y)-> armijo(x, y), lsatx);
#    stop_meta = StoppingMeta(max_iter = 50)
#
#    state, optimality = eval(algo)(h, stop_ls, verboseLS = true)
#    @test optimality == false
#    reset!(nlp)
# end

ext: jl


text:
using Chain
using Utils

inputfile = replace(@__FILE__, r"\.jl$" => ".txt")
load(text) = @chain text begin
    per_split_parse(_, ',')
end
fload() = inputfile |> slurp |> load
input = fload()

const COUNT = 4

function part1(memory, replace=nothing)
    try
        memory = zerobased(copy(memory))
        if !isnothing(replace)
            memory[1] = replace[1]
            memory[2] = replace[2]
        end
        ip = 0
        while true
            opc = memory[ip]
            if opc == 99
                break
            else
                op = opc == 1 ? Base.:+ : Base.:*
                memory[memory[ip+3]] = op(memory[memory[ip+1]], memory[memory[ip+2]])
            end
            ip += 4
        end
        memory[0]
    catch BoundsError
        typemin(eltype(memory))
    end
end

function part2()
    tuples = zerobased(collect(Iterators.product(0:99, 0:99)))
    @chain tuples begin
        findfirst(t -> part1(input, t) == 19690720, _)
        get(tuples, _, 0)
        100 * _[1] + _[2]
    end
end

println("Last state before the fire: ", part1(input, (12, 2)))
println("100 * noun + verb: ", part2())

max_stars_repo_name: CmdQ/AoC


text:
using Test
using Revise
using Deeplearning

@testset "basic.jl" begin
    Deeplearning.df1()
    Deeplearning.df_arg_matrix()
    Deeplearning.df_m()
    Deeplearning.simple_models()
    Deeplearning.build_layers()
    Deeplearning.build_layers_o1()
    Deeplearning.build_layers_dense()
end


text:
###############################################################################
#
# Photosynthesis mode
#
###############################################################################
"""
    abstract type AbstractCalculationMode

Hierarchy of AbstractCalculationMode
- [`GCO₂Mode`](@ref)
- [`PCO₂Mode`](@ref)
"""
abstract type AbstractCalculationMode end

"""
    struct GCO₂Mode <: AbstractCalculationMode

Calculate leaf photosynthesis using leaf internal CO₂ partial pressure
"""
struct GCO₂Mode <: AbstractCalculationMode end

"""
    struct PCO₂Mode <: AbstractCalculationMode

Calculate leaf photosynthesis using total leaf diffusive conductance to CO₂
"""
struct PCO₂Mode <: AbstractCalculationMode end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
#__precompile__()
module GtkIDE

global const HOMEDIR = @__DIR__()
global const REDIRECT_STDOUT = true
pkgdir(pkg::Module) = abspath(joinpath(dirname(Base.pathof(pkg)), ".."))

using Gtk, GtkSourceWidget, GtkUtilities, GtkMarkdownTextView
using Immerse
using GtkREPL
using JSON, ConfParser, Refactoring
using Distributed, REPL, InteractiveUtils, Sockets, Markdown, Pkg
using GtkREPL.GtkTextUtils

import GtkREPL: ConsoleManager, Console, current_console, print_to_console, new_prompt,
worker, add_remote_console_cb, ConsoleCommand, on_command_done, index, style_css, get_tab,
PROPAGATE, INTERRUPT, nonmutable, offset

include("Options.jl")

import Gtk: GtkTextIter, get_default_mod_mask, GdkKeySyms, selected, hasselection, mutable
import Gadfly.Colors
import Immerse: Cairo, Compose

import REPL: REPLCompletions.completions
import Base: push!, search
import Cairo.text

export image, plot, figure, rprint
export GtkREPL, Pkg #This gets called by remote consoles

function method_filename(m)
    tv, decls, file, line = Base.arg_decl_parts(m.ms[1])
    return file, line
end

#export add_console, figure

#Order matters
include("gtk_utils.jl")
include("PlotWindow.jl")
include("StyleAndLanguageManager.jl")
include("MainWindow.jl")
include("Project.jl")
include("Console.jl")
include("ConsoleCommands.jl")
include("Editor.jl")
include("NtbookUtils.jl")
include("PathDisplay.jl")
include("MainMenu.jl")
include("SidePanels.jl")
include(joinpath("sidepanels", "FilesPanel.jl"))
include(joinpath("sidepanels", "WorkspacePanel.jl"))
include(joinpath("sidepanels", "ProjectsPanel.jl"))
include("Logo.jl")

include("init.jl")

#__init__()
end#module
#


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/PlaneElasticity.jl/blob/master/LICENSE

""" Plane elasticity implementation for JuliaFEM. """
module PlaneElasticity

using FEMBase

import FEMBase: get_unknown_field_name,
                get_formulation_type,
                assemble_elements!

type PlaneStress <: FieldProblem
end

type PlaneStrain <: FieldProblem
end

function get_unknown_field_name(::Problem{PlaneStress})
    return "displacement"
end

function get_unknown_field_name(::Problem{PlaneStrain})
    return "displacement"
end

function assemble_elements!{B}(::Problem{PlaneStress}, ::Assembly,
                               elements::Vector{Element{B}}, ::Float64)

    for element in elements
        info("Not doing anything useful right now.")
    end

end

function assemble_elements!{B}(::Problem{PlaneStrain}, ::Assembly,
                               elements::Vector{Element{B}}, ::Float64)

    for element in elements
        info("Not doing anything useful right now.")
    end

end

export PlaneStress, PlaneStrain

end


text:
using Swirlflow
using Test

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


text:
module ParameterHandling

using Compat: only
using ChainRulesCore
using InverseFunctions: inverse
using LogExpFunctions: logit, logistic
using LinearAlgebra
using SparseArrays

export flatten,
    value_flatten, positive, bounded, fixed, deferred, orthogonal, positive_definite

include("flatten.jl")
include("parameters.jl")

include("test_utils.jl")

end # module


text:
"""
    CoreForwardModel(forwarddates, forwardprices, yieldcurve)

Specifies the value of a `SingleStock` in terms of its forward prices.

* `forwarddates`: ordered dates for which forward prices are available
* `forwardprices`: forward prices of the asset
* `yieldcurve`: risk-free discount curve
"""
struct CoreForwardModel{DT <: AbstractVector{Date}, T, S<:YieldTermStructure} <: AbstractCoreModel
    forwarddates::DT
    forwardprices::Vector{T}
    yieldcurve::S
end

CoreForwardModel(startdate::Date, forwarddates::AbstractVector{Date}, forwardprices::Vector{T}, yieldrate::Float64) where {T} =
    CoreForwardModel(forwarddates, forwardprices,
                     ConstantContinuousYieldCurve(Actual365(), yieldrate, startdate))

CoreForwardModel(startdate::Date, forwarddate::Date, forwardprice, yieldrate::Float64) =
    CoreForwardModel(startdate, [forwarddate], [forwardprice], yieldrate)

numeraire(m::CoreForwardModel) = unit(m.forwardprices[1])
startdate(m::CoreForwardModel) = startdate(m.yieldcurve)

yearfractionto(m::CoreForwardModel, dt::Date) =
    yearfraction(daycount(m.yieldcurve),  startdate(m),  dt)

function value(m::CoreForwardModel, c::WhenAt{Receive{T}}) where T
    value(m, c.c) * discount(m.yieldcurve, maturitydate(c))
end

function date2index(m::CoreForwardModel, dt::Date)
    ii = searchsorted(m.forwarddates, dt)
    isempty(ii) && throw(DomainError())
    return ii[1]
end

forwardprice(m::CoreForwardModel, ::SingleStock, dt::Date) =
    m.forwardprices[date2index(m, dt)]

function value(m::CoreForwardModel, c::SingleStock)
    forwardprice(m, c, startdate(m))
end
function value(m::CoreForwardModel, c::WhenAt{SingleStock})
    dt = maturitydate(c)
    forwardprice(m, c.c, dt) * discount(m.yieldcurve, dt)
end


text:
using Test, ForwardDiff, TensorNetworkEvolve
using Cracker, OMEinsum

@testset "OMEinsum" begin
    T = Float64
    A, B, C = track(rand(T, 2, 2)), track(rand(T, 2, 2)), track(rand(T, 2))
    Z = ein"ij, jk, k->ik"(A, B, C)
    @test Z.record.f === einsum

    ret = abs2(sum(Z))
    Cracker.backpropagate!(ret, 1.0)
    a, b, c = untrack.((A, B, C))
    gA = ForwardDiff.gradient(a->abs2(sum(ein"ij, jk, k->ik"(a, b, c))), a)
    gB = ForwardDiff.gradient(b->abs2(sum(ein"ij, jk, k->ik"(a, b, c))), b)
    gC = ForwardDiff.gradient(c->abs2(sum(ein"ij, jk, k->ik"(a, b, c))), c)
    @test A.record.grad ≈ gA
    @test B.record.grad ≈ gB
    @test C.record.grad ≈ gC
end


text:
"""
    A = JopDct(T, n...[;dims=()])

where `A` is a type II discrete cosine transfrom (DCT), T is the element type of
the domain and range if A, and n is the size of A.  If `dims` is not set, then `A`
is the N-dimensional DCT and where `N=length(n)`.  Otherwise, `A` is a DCT
over the dimensions specified by `dims`.

# Examples

## 1D
```julia
A = JopDct(Float64, 128)
m = rand(domain(A))
d = A*m
```

## 2D
```julia
A = JopDct(Float64, 128, 64)
m = rand(domain(A))
d = A*m
```

## 2D transform over 3D array
```julia
A = JopDct(Float64,128,64,32;dims=(1,2))
m = rand(domain(A))
d = A*m
```

# Notes
* the adjoint pair is achieved by application `sqrt(2/N)` and `sqrt(2)` factors.  See:
    * https://en.wikipedia.org/wiki/Discrete_cosine_transform
    * http://www.fftw.org/fftw3_doc/1d-Real_002deven-DFTs-_0028DCTs_0029.html#g_t1d-Real_002deven-DFTs-_0028DCTs_0029
"""
function JopDct(T, n...; dims=())
    dims = dims == () ? ntuple(i->i, length(n)) : dims
    JopLn(df! = JopDct_df!, df′! = JopDct_df′!, dom = JetSpace(T, n), rng = JetSpace(T, n), s = (dims=dims,))
end
export JopDct

function JopDct_df!(d::AbstractArray{T}, m::AbstractArray{T}; dims, kwargs...) where {T}
    d .= m
    dct!(d, dims)
    sc1, sc = dctscale(dims, size(m), T)
    d[1] *= sc1
    d .*= sc
end

function JopDct_df′!(m::AbstractArray{T}, d::AbstractArray{T}; dims, kwargs...) where {T}
    m .= d
    sc1, sc = dctscale(dims, size(m), T)
    m[1] *= sc1
    idct!(m, dims)
    m .*= sc
end

dctscale(dims, n, ::Type{T}) where {T} = length(dims)*sqrt(T(2)), T(sqrt(1.0/mapreduce(dim->n[dim], *, dims)))


text:
using Test
using Altro
using BenchmarkTools
using TrajectoryOptimization
using RobotDynamics
using RobotZoo
using StaticArrays, LinearAlgebra
using JLD2
using ForwardDiff
using FiniteDiff
using Random
using SparseArrays
const RD = RobotDynamics
const TO = TrajectoryOptimization

TEST_TIME = false
Random.seed!(1)

##
# @testset "Initialization" begin
#     include("large_init.jl")
# end

if TEST_TIME
    @testset "Benchmark Suite" begin
        include("benchmark_problems.jl")
    end
end

@testset "Benchmark Problems" begin
    include("nl_cartpole.jl")
    include("cartpole.jl")
    include("quadrotor.jl")
    include("escape_solve.jl")
    include("rocket_test.jl")
    include(joinpath(@__DIR__,"..","examples","quickstart.jl"))
end

@testset "Solvers" begin
    include("constructors.jl")
    include("augmented_lagrangian_tests.jl")
    include("alilqr_test.jl")
    include("ilqr_test.jl")
    include("alcon_test.jl")
    include("alconset_test.jl")
    include("solve_tests.jl")
    include("socp_test.jl")
    include("projected_newton_test.jl")
    include("expansion_test.jl")
    include("infeasible_problem.jl")
end

@testset "Solver Options" begin
    include("solver_opts_test.jl")
    include("solver_stats_test.jl")
end

@testset "Quickstart" begin
    include(joinpath(@__DIR__, "..", "examples", "quickstart.jl"))
    @test max_violation(altro) ≈ 5.8962e-7 rtol=1e-4
    @test cost(altro) ≈ 1.539 rtol = 1e-4
    @test iterations(altro) == 44
end

@testset "Utils" begin
    @testset "QDLDL" begin
        include("qdldl_test.jl")
    end
    @testset "Sparseblocks" begin
        include("sparseblocks_test.jl")
    end
end


text:
push!(LOAD_PATH, pwd())
include("startup.jl")

using HDF5
#using PyCall
#pygui(:tk) # :tk, :gtk3, :gtk, :qt5, :qt4, :qt, or :wx
using PyPlot
import .InterpolationFunctions
import  .MasslistFunctions
import .ResultFileFunctions
import .PlotFunctions

#fp = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/H3ORF400Vpp_final/DryAtmo/results/"
#fn = "_result_withDMS.hdf5"
fp = "/media/wiebke/Extreme SSD/DMS_TME_HPMTF_IBK2020/data/O3ramp/results/"
#fp = "/media/wiebke/Extreme SSD/DMS_TME_HPMTF_IBK2020/data/DMSramp/results/"
#fp = "/media/wiebke/Extreme SSD/DMS_TME_HPMTF_IBK2020/data/RO2ramp/results/"
#fp = "/media/wiebke/Extreme SSD/DMS_TME_HPMTF_IBK2020/data/OCSexperiment/results/"
fn = "_result.hdf5"
file = string(fp,fn)

plotstages = false
stagesfile = "/media/wiebke/Extreme SSD/DMS_TME_HPMTF_IBK2020/DMS_TME_BAG_List.csv"		# add filename, when plotstages = true

plotHighTimeRes = true # Plot every datapoint or only file averages
plotFittedInsteadOfSummed = true # Use multi peak fitted data instead of raw
smoothing = 300 # Average n samples, 1 for raw
plotsymbol = "-"
isobarToPlot = 0
timedelay = Dates.Hour(0) # CLOUD12, ...
#timedelay = Dates.Hour(1) # CLOUDX, CLOUD11
#timedelay = Dates.Hour(4) # SALTENA

backgroundSubstractionMode = 2 # 0=no correction, 1=minimum found in all data,  2=mean of range given below
backgroundStart = DateTime(2020,5,28,8,02)
backgroundEnd = DateTime(2020,5,28,8,20)

massrange2plot = false
plotWronglyAssigned = false

include("manualMassLibrary.jl")

massesToPlot = [
# Examples for selecting what to plot:

#H3O[1]
#H3OH2O[1]
#H3OH2OH2O[1]
#MasslistFunctions.createCompound(H=8, O=4)[1]
#NH3H2O[1]
#NH3[1]
#NH3NH3[1]
#MasslistFunctions.createCompound(Hplus=0, O=2)[1]
#NO[1]

#MasslistFunctions.massFromComposition(C=10,H=16,O=2)
#APINENE[1]
#MasslistFunctions.massFromComposition(C=9,H=8,O=1)	# Cinnamyl-Alkohol?
#MasslistFunctions.massFromComposition(C=8,H=14,O=1)	#
#MasslistFunctions.massFromComposition(C=7,H=8,O=1)	#
#MasslistFunctions.massFromComposition(C=7,H=8,O=2)	#
#MasslistFunctions.massFromComposition(C=7,H=8,O=3)	#
#MasslistFunctions.massFromComposition(C=7,H=6,O=1)	#
#MasslistFunctions.massFromComposition(C=7,H=6,O=2)	#
#MasslistFunctions.massFromComposition(C=7,H=6,O=3)	#
#MasslistFunctions.massFromComposition(C=7,H=10,O=3)	#
#MasslistFunctions.massFromComposition(C=7,H=10,O=4)	#
#MasslistFunctions.massFromComposition(C=16,H=18,O=10)	# D5
#MasslistFunctions.massFromComposition(C=15,H=24,O=15)	# D6
#MasslistFunctions.massFromComposition(C=21,H=26,O=15)	# D7
#MasslistFunctions.massFromComposition(C=20,H=32,O=20)	# D8
#MasslistFunctions.massFromComposition(C=20,H=38,O=7)	# continuosly increases during afternoon - temperature effect?
#MasslistFunctions.massFromComposition(C=17,H=30,O=14)	# ?! calib?
#MasslistFunctions.massFromComposition(C=19,H=32,O=17)	# ?! calib?
#MasslistFunctions.massFromComposition(C=21,H=28,O=16)	#
#MasslistFunctions.massFromComposition(C=22,H=22,O=17)	# daily mini movement
#MasslistFunctions.massFromComposition(C=21,H=22,O=19)	# daily mini movement

#MasslistFunctions.massFromComposition(C=3,H=8,O=3,N=1)

#MasslistFunctions.massFromComposition(C=2,H=8,S=1,N=1,O=4)
#MasslistFunctions.createMassList(; C=1:2, O=0:10, N=0, S=1, nHplus=1, allowRadicals=true)[1]
#DMS[1] #both days    left + right
#MasslistFunctions.massFromComposition(C=2,H=9,N=1,S=1)
#DMSO[1] #both days  left
#DMSO2[1] #11-05    left + right
#MSIA[1]
#MSA[1] #both days    left?  + right
#MasslistFunctions.massFromComposition(C=2,H=4,S=1,O=3)
#MasslistFunctions.massFromComposition(C=2,H=4,S=1,O=1)
#MasslistFunctions.massFromComposition(C=2,H=6,S=1,O=3)
#MasslistFunctions.massFromComposition(C=1, S=1,O=1, Hplus = 0)

#MasslistFunctions.massFromComposition(C=1,H=2,O=1)
#MasslistFunctions.massFromComposition(C=1,H=7,S=1,O=1)
#MasslistFunctions.massFromComposition(C=2,H=7,S=1,O=1) #11-05   left
#MasslistFunctions.massFromComposition(C=2,H=7,S=1,O=2) #both days    left
#MasslistFunctions.massFromComposition(C=1,H=4,S=1,O=2) #    left
#MasslistFunctions.massFromComposition(C=1,H=4,S=1)
#MasslistFunctions.massFromComposition(C=2,H=5,S=1) #both days     right
#MasslistFunctions.massFromComposition(C=2,H=5,S=1,O=2) #11-05 az    right
#MasslistFunctions.massFromComposition(C=2,H=5,S=1,O=4,N=1)
#MasslistFunctions.massFromComposition(C=2,H=5,S=1,O=1)
#MasslistFunctions.massFromComposition(C=1,H=3,S=1)
#MasslistFunctions.massFromComposition(C=1,H=3,S=1,O=1)
#MasslistFunctions.massFromComposition(C=1,H=3,S=1,O=2) #20-05    right
#MasslistFunctions.massFromComposition(C=1,H=3,S=1,O=3)    #  right
#MasslistFunctions.massFromComposition(S=1,O=2)

# DMS + TME radicals
#MasslistFunctions.massFromComposition(C=1,H=3,S=1,O=2)
#MasslistFunctions.massFromComposition(C=1,H=3,O=2)
MasslistFunctions.massFromComposition(C=3,H=5,O=3)
  #MasslistFunctions.massFromComposition(C=6,H=13,O=3)
MasslistFunctions.massFromComposition(C=2,H=5,S=1,O=2)
#MasslistFunctions.massFromComposition(C=1,H=3,S=1,O=3)
#MasslistFunctions.massFromComposition(C=1,H=3,S=1,O=4)
#MasslistFunctions.massFromComposition(C=2,H=7,S=1,O=3)
#MasslistFunctions.massFromComposition(C=2,H=5,S=1,O=4)

# DMS + TME dimers

MasslistFunctions.massFromComposition(C=2,H=6,O=2)
MasslistFunctions.massFromComposition(C=4,H=8,O=3)
MasslistFunctions.massFromComposition(C=6,H=10,O=4)
  #MasslistFunctions.massFromComposition(C=7,H=16,O=3)

MasslistFunctions.massFromComposition(C=9,H=18,O=4)
MasslistFunctions.massFromComposition(C=12,H=26,O=4)
MasslistFunctions.massFromComposition(C=2,H=6,S=1,O=2)

MasslistFunctions.massFromComposition(C=2,H=6,S=1,O=3)
MasslistFunctions.massFromComposition(C=2,H=6,S=1,O=4)
MasslistFunctions.massFromComposition(C=3,H=8,S=1,O=2)
MasslistFunctions.massFromComposition(C=3,H=8,S=1,O=4)
  #MasslistFunctions.massFromComposition(C=3,H=10,S=1,O=3)
MasslistFunctions.massFromComposition(C=4,H=8,S=1,O=3)
MasslistFunctions.massFromComposition(C=4,H=8,S=1,O=4)
MasslistFunctions.massFromComposition(C=4,H=8,S=1,O=5)
MasslistFunctions.massFromComposition(C=5,H=10,S=1,O=3)
MasslistFunctions.massFromComposition(C=5,H=10,S=1,O=5)
MasslistFunctions.massFromComposition(C=5,H=12,S=1,O=4)
MasslistFunctions.massFromComposition(C=7,H=16,S=1,O=3)
MasslistFunctions.massFromComposition(C=7,H=16,S=1,O=4)
MasslistFunctions.massFromComposition(C=7,H=16,S=1,O=5)
MasslistFunctions.massFromComposition(C=8,H=18,S=1,O=3)
MasslistFunctions.massFromComposition(C=8,H=18,S=1,O=5)
MasslistFunctions.massFromComposition(C=8,H=20,S=1,O=4)
MasslistFunctions.massFromComposition(C=3,H=8,S=2,O=2)
MasslistFunctions.massFromComposition(C=3,H=8,S=2,O=3)
MasslistFunctions.massFromComposition(C=3,H=8,S=2,O=4)
MasslistFunctions.massFromComposition(C=4,H=10,S=2,O=2)
MasslistFunctions.massFromComposition(C=4,H=10,S=2,O=4)
MasslistFunctions.massFromComposition(C=4,H=12,S=2,O=3)

]

if massrange2plot
	result4masses = ResultFileFunctions.loadResults(file; useAveragesOnly= true, masslistOnly = true)
	massesToPlot = result4masses.MasslistMasses[59.8 .< result4masses.MasslistMasses .< 60.0]
end

measResult = ResultFileFunctions.loadResults(file, massesToLoad=massesToPlot, useAveragesOnly=!plotHighTimeRes, raw=!plotFittedInsteadOfSummed, massMatchTolerance=0.001)
if isobarToPlot != 0
  isobarResult = ResultFileFunctions.loadResults(file, massesToLoad=[isobarToPlot+0.3], massMatchTolerance=0.5, useAveragesOnly=!plotHighTimeRes, raw=!plotFittedInsteadOfSummed)
  measResult=ResultFileFunctions.joinResultsMasses(measResult, isobarResult)
end

measResult.Times = measResult.Times .- timedelay

fig=figure()
ax = subplot(111)

if (backgroundSubstractionMode == 0)
  background=0
elseif (backgroundSubstractionMode == 1)
  background = minimum(InterpolationFunctions.averageSamples(measResult.Traces,smoothing),dims=1)
elseif backgroundSubstractionMode == 2
  import Statistics
  background = Statistics.mean(measResult.Traces[(measResult.Times.>backgroundStart) .& (measResult.Times.<backgroundEnd),:],dims=1)
end

bgCorrectedTraces = measResult.Traces .- background

ax.semilogy(Dates.unix2datetime.(InterpolationFunctions.averageSamples(Dates.datetime2unix.(measResult.Times), smoothing)),InterpolationFunctions.averageSamples(bgCorrectedTraces[:,1],smoothing), ".-") #linewidth=1)
ax.semilogy(Dates.unix2datetime.(InterpolationFunctions.averageSamples(Dates.datetime2unix.(measResult.Times), smoothing)),InterpolationFunctions.averageSamples(bgCorrectedTraces[:,2:end],smoothing), plotsymbol) #linewidth=1)

#semilogy(Dates.unix2datetime(InterpolationFunctions.averageSamples(Dates.datetime2unix(measResult.Times), smoothing)),InterpolationFunctions.averageSamples(bgCorrectedTraces,smoothing), plotsymbol) #linewidth=1)

if plotstages == true
	PlotFunctions.plotStages(stagesfile, ax, measResult.Times[1], measResult.Times[end]; headerrow = 9)
end

startTimeString = Dates.format(measResult.Times[1],"yyyy/mm/dd")
endTimeString = Dates.format(measResult.Times[end],"yyyy/mm/dd")
title("$startTimeString - $endTimeString")
xlabel("Time [UTC]")
ylabel("Signal [CPS]")

legStrings = []
for i = 1:length(measResult.MasslistMasses)
  push!(legStrings,"m/z $(round(measResult.MasslistMasses[i],digits=3)) - $(MasslistFunctions.sumFormulaStringFromCompositionArray(measResult.MasslistCompositions[:,i]))")
end

if plotWronglyAssigned
	legStrings = []
	for i = 1:length(measResult.MasslistMasses)
	  push!(legStrings,"m/z $(round(measResult.MasslistMasses[i],digits=3))")
	end
end
box = ax.get_position()
cols = 1

majorformatter = matplotlib.dates.DateFormatter("%m/%d %H:%M")
ax.xaxis.set_major_formatter(majorformatter)
legend(legStrings)
#grid()

#tight_layout()


text:
using Test
using SparseArrays
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

@testset "interface.jl file Tests" begin
	@testset "characteristicMatrix Tests" begin
		V,(VV,EV,FV,CV) = Lar.cuboid([1.,1.,1.], true);
		@test Matrix(Lar.characteristicMatrix(FV)) == [
		 1  1  1  1  0  0  0  0;
		 0  0  0  0  1  1  1  1;
		 1  1  0  0  1  1  0  0;
		 0  0  1  1  0  0  1  1;
		 1  0  1  0  1  0  1  0;
		 0  1  0  1  0  1  0  1]
		@test size(Lar.characteristicMatrix(CV))==(1,8)
		@test typeof(Lar.characteristicMatrix(CV))==SparseMatrixCSC{Int8,Int64}
		@test Matrix(Lar.characteristicMatrix(EV)) == [
		 1  1  0  0  0  0  0  0;
		 0  0  1  1  0  0  0  0;
		 0  0  0  0  1  1  0  0;
		 0  0  0  0  0  0  1  1;
		 1  0  1  0  0  0  0  0;
		 0  1  0  1  0  0  0  0;
		 0  0  0  0  1  0  1  0;
		 0  0  0  0  0  1  0  1;
		 1  0  0  0  1  0  0  0;
		 0  1  0  0  0  1  0  0;
		 0  0  1  0  0  0  1  0;
		 0  0  0  1  0  0  0  1]
		@test size(Lar.characteristicMatrix(EV))==(12,8)
		@test typeof(Lar.characteristicMatrix(EV))==SparseMatrixCSC{Int8,Int64}
	end

	@testset "signed_boundary_1 Tests" begin
		V,(VV,EV,FV,CV) = Lar.cuboid([1.,1.,1.], true);
		signed_boundary_1 = Lar.boundary_1( EV::Lar.Cells )

		@test EV==[[1, 2], [3, 4], [5, 6], [7, 8], [1, 3], [2, 4], [5, 7],
			[6, 8], [1, 5], [2, 6], [3, 7], [4, 8]]
		@test length(EV)==12
		@test typeof(EV)==Array{Array{Int64,1},1}
		@test size(signed_boundary_1)==(8,12)
		@test typeof(signed_boundary_1)==SparseMatrixCSC{Int8,Int64}
		@test nnz(signed_boundary_1)==24

		@test Matrix(Lar.boundary_1(EV::Lar.Cells))==[
		 -1   0   0   0  -1   0   0   0  -1   0   0   0;
		  1   0   0   0   0  -1   0   0   0  -1   0   0;
		  0  -1   0   0   1   0   0   0   0   0  -1   0;
		  0   1   0   0   0   1   0   0   0   0   0  -1;
		  0   0  -1   0   0   0  -1   0   1   0   0   0;
		  0   0   1   0   0   0   0  -1   0   1   0   0;
		  0   0   0  -1   0   0   1   0   0   0   1   0;
		  0   0   0   1   0   0   0   1   0   0   0   1]
	end

	@testset "unsigned_coboundary_1 Tests" begin
		V,(VV,EV,FV,CV) = Lar.cuboid([1.,1.,1.], true);
		unsigned_coboundary_1 = Lar.u_coboundary_1(FV,EV)

		@test size(unsigned_coboundary_1)==(6,12)
		@test nnz(unsigned_coboundary_1)==24
		@test typeof(unsigned_coboundary_1)==SparseMatrixCSC{Int8,Int64}

		@test Matrix(unsigned_coboundary_1)==[
		 1  1  0  0  1  1  0  0  0  0  0  0;
		 0  0  1  1  0  0  1  1  0  0  0  0;
		 1  0  1  0  0  0  0  0  1  1  0  0;
		 0  1  0  1  0  0  0  0  0  0  1  1;
		 0  0  0  0  1  0  1  0  1  0  1  0;
		 0  0  0  0  0  1  0  1  0  1  0  1]
	 end

	@testset "signed_coboundary_1 Tests" begin
		FV = [[1,2,3,4,5,17,16,12],
		[1,2,3,4,6,7,8,9,10,11,12,13,14,15],
		[4,5,9,11,12,13,14,15,16,17],
		[2,3,6,7], [8,9,10,11]]

		EV = [[1,2],[2,3],[3,4],[4,5],[1,12],[2,6],[3,7],[4,9],[5,17],[6,7],[8,9],
		[8,10],[9,11],[10,11],[11,15],[12,13],[12,16],[13,14],[14,15],[16,17]]

		V = [0   2   5   7  10   2   5   3   7  3  7  0  3  3  7  0  10;
			16  16  16  16  16  13  13  11  11  8  8  5  5  2  2  0   0]
		signed_coboundary_1 = Lar.coboundary_1( V,FV,EV,false,true );

		@test size(Matrix(signed_coboundary_1))==(5,20)
		@test typeof(signed_coboundary_1)==SparseArrays.SparseMatrixCSC{Int8,Int64}
		@test nnz(signed_coboundary_1)==40
		@test signed_coboundary_1[1,1]==-1
		@test signed_coboundary_1[3,20]==-1
		@test Matrix(signed_coboundary_1)==
		[-1  -1  -1  -1   1   0   0   0  -1   0   0   0   0   0   0   0   1   0   0   1;
		  1   0   1   0  -1   1  -1   1   0   1  -1   1   0   1   1  -1   0  -1  -1   0;
		  0   0   0   1   0   0   0  -1   1   0   0   0  -1   0  -1   1  -1   1   1  -1;
		  0   1   0   0   0  -1   1   0   0  -1   0   0   0   0   0   0   0   0   0   0;
		  0   0   0   0   0   0   0   0   0   0   1  -1   1  -1   0   0   0   0   0   0]
	end

	@testset "coboundary_1 with non-convex cells Tests" begin
		V = [5  0  6  9  5  3  0  5  8  10  10  10  0  0  3  9  8  6  5;
			 9  9  5  5  5  6  6  0  0   9   5   0  0  3  3  8  3  8  3]

		EV = [[1,2],[3,4],[3,5],[1,10],[16,18],[7,14],[6,15],[8,13],[13,14],[10,11],
		[11,12],[8,19],[9,17],[1,5],[17,19],[6,7],[14,15],[4,11],[3,18],[4,16],[8,9],[9,12],[2,7]]

		FV = [[1,2,7,6,15,14,13,8,19,17,9,12,11,4,3,5],
		[3,4,16,18],[6,7,14,15],[1,3,4,5,10,11,16,18],
		[8,9,17,19],[1,2,7,14,13,8,9,12,11,10]]

		copFE = Lar.coboundary_1( V::Lar.Points, FV::Lar.Cells, EV::Lar.Cells, false,true );

		@test Matrix(copFE) ==
		[-1  0  0  1  0 -1  0  1  1  1  1  0  0  0  0  0  0  0  0  0 -1 -1 -1;
		  1 -1  1  0  0  0  1 -1 -1  0 -1  1 -1 -1 -1 -1 -1 -1  0  0  0  1  1;
		  0  1  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0 -1  1  0  0  0;
		  0  0  0  0  0  1 -1  0  0  0  0  0  0  0  0  1  1  0  0  0  0  0  0;
		  0  0 -1 -1 -1  0  0  0  0 -1  0  0  0  1  0  0  0  1  1 -1  0  0  0;
		  0  0  0  0  0  0  0  0  0  0  0 -1  1  0  1  0  0  0  0  0  1  0  0]
	end

	@testset "chaincomplex 2D Tests" begin
		W =
		 [0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0  2.0  2.0  2.0  2.0  3.0  3.0  3.0  3.0;
		  0.0  1.0  2.0  3.0  0.0  1.0  2.0  3.0  0.0  1.0  2.0  3.0  0.0  1.0  2.0  3.0]
		EW =
		[[1, 2],[2, 3],[3, 4],[5, 6],[6, 7],[7, 8],[9, 10],[10, 11],[11, 12],[13, 14],
		 [14, 15],[15, 16],[1, 5],[2, 6],[3, 7],[4, 8],[5, 9],[6, 10],[7, 11],[8, 12],
		 [9, 13],[10, 14],[11, 15],[12, 16]]
		V,bases,coboundaries = Lar.chaincomplex(W,EW)

		@test length(bases[1])==24	# edges
		@test typeof(bases[1])==Array{Array{Int64,1},1}	# edges

		@test length(bases[2])==9 # faces -- previously unknown !!
		@test typeof(bases[2])==Array{Array{Int64,1},1} # faces

		@test size(coboundaries[1])==(24, 16) # coboundary_1
		@test typeof(coboundaries[1])==SparseArrays.SparseMatrixCSC{Int8,Int64}
		# coboundary_1
		@test SparseArrays.nnz(coboundaries[1])==48 # coboundary_1

		@test size(coboundaries[2])==(9,24) # coboundary_2: oriented 2-cycles of faces
		@test typeof(coboundaries[2])==SparseArrays.SparseMatrixCSC{Int8,Int64}
		@test SparseArrays.nnz(coboundaries[2])==36
		# coboundary_2: oriented 2-cycles of faces
		@test Matrix(coboundaries[2]) ==
		[-1  0  0  1  0  0  0  0  0  0  0  0  1 -1  0  0  0  0  0  0  0  0  0  0;
		  0 -1  0  0  1  0  0  0  0  0  0  0  0  1 -1  0  0  0  0  0  0  0  0  0;
		  0  0 -1  0  0  1  0  0  0  0  0  0  0  0  1 -1  0  0  0  0  0  0  0  0;
		  0  0  0 -1  0  0  1  0  0  0  0  0  0  0  0  0  1 -1  0  0  0  0  0  0;
		  0  0  0  0 -1  0  0  1  0  0  0  0  0  0  0  0  0  1 -1  0  0  0  0  0;
		  0  0  0  0  0 -1  0  0  1  0  0  0  0  0  0  0  0  0  1 -1  0  0  0  0;
		  0  0  0  0  0  0  0 -1  0  0  1  0  0  0  0  0  0  0  0  0  0  1 -1  0;
		  0  0  0  0  0  0 -1  0  0  1  0  0  0  0  0  0  0  0  0  0  1 -1  0  0;
		  0  0  0  0  0  0  0  0 -1  0  0  1  0  0  0  0  0  0  0  0  0  0  1 -1]
	end

#	@testset "chaincomplex 3D Tests" begin
#		cube_1 = ([0 0 0 0 1 1 1 1; 0 0 1 1 0 0 1 1; 0 1 0 1 0 1 0 1],
#		[[1,2,3,4],[5,6,7,8],[1,2,5,6],[3,4,7,8],[1,3,5,7],[2,4,6,8]],
#		[[1,2],[3,4],[5,6],[7,8],[1,3],[2,4],[5,7],[6,8],[1,5],[2,6],[3,7],[4,8]] )
#
#		cube_2 = Lar.Struct(
#			[Lar.t(
#			0,0,0.5), Lar.r(0,0,pi/3), cube_1])
#		V,FV,EV = Lar.struct2lar(
#			Lar.Struct([ cube_1, cube_2 ]))
#
#		V,bases,coboundaries = Lar.chaincomplex(V,FV,EV)
#		(EV, FV, CV), (cscEV, cscFE, cscCF) = bases,coboundaries
#
#		@test typeof(FV)==Array{Array{Int64,1},1} # bases[2]
#		@test length(FV)==18 # bases[2]
#		@test FV==[[1,3,4,6],[2,3,5,6],[7,8,9,10],[1,2,3,7,8],[4,6,9,10,11,12],[5,6,11,12],
#		[1,4,7,9],[2,5,11,13],[2,8,10,11,13],[2,3,14,15,16],[11,12,13,17],
#		[11,12,13,18,19,20],[2,3,13,17],[2,13,14,18],[15,16,19,20],[3,6,12,15,19],
#		[3,6,12,17],[14,16,18,20]]
#		@test typeof(CV)==Array{Array{Int64,1},1} # bases[3]
#		@test length(CV)==3 # bases[3]
#		@test CV==[[2, 3, 5, 6, 11, 12, 13, 14, 15, 16, 18, 19, 20],
#		[2, 3, 5, 6, 11, 12, 13, 17],
#		[1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 13, 17]]
#
#		@test typeof(cscEV)==SparseArrays.SparseMatrixCSC{Int8,Int64} # coboundaries[1]
#		@test size(cscEV)==(34,20) # coboundaries[1]
#		@test nnz(cscEV)==68 # coboundaries[1]
#
#		@test typeof(cscFE)==SparseArrays.SparseMatrixCSC{Int8,Int64} # coboundaries[2]
#		@test size(cscFE)==(18,34) # coboundaries[2]
#		@test nnz(cscFE)==80 # coboundaries[2]
#
#		@test typeof(cscCF)==SparseArrays.SparseMatrixCSC{Int8,Int64} # coboundaries[3]
#		@test size(cscCF)==(4,18) # coboundaries[3]
#		@test nnz(cscCF)==36 # coboundaries[3]
#	end
end


text:
#=
module NonsmoothFwdAD
=====================
A quick implementation of:

- the vector forward mode of automatic differentiation (AD) for generalized
  derivative evaluation, developed in the article:
  KA Khan and PI Barton (2015), https://doi.org/10.1080/10556788.2015.1025400

- the "compass difference" rule for generalized derivatives of scalar-valued bivariate
  functions, developed in the article:
  KA Khan and Y Yuan (2020), https://doi.org/10.46298/jnsao-2020-6061

This implementation applies generalized differentiation rules via operator overloading,
and is modeled on the standard "smooth" vector forward AD mode implementation described
in Chapter 6 of "Evaluating Derivatives (2nd ed.)" by Griewank and Walther (2008).
The "AFloat" struct here is analogous to Griewank and Walther's "adouble" class.

The following operators have been overloaded:
x+y, -x, x-y, x*y, inv(x), x/y, x^p, exp, log, sin, cos, abs, max, min, hypot
(where "p" is a fixed integer)

Additional operators may be overloaded in the same way; a new unary operation would be
overloaded just like "log" or "abs", and a new binary operation would be overloaded
just like "x*y" or "hypot".

Written by Kamil Khan on February 5, 2022
Edited by Maha Chaudhry in April 2022
=#
module NonsmoothFwdAD

using LinearAlgebra

export AFloat

export eval_ld_derivative,
    eval_dir_derivative,
    eval_gen_derivative,
    eval_gen_gradient,
    eval_compass_difference

# default tolerances; feel free to change.
const DEFAULT_ZTOL = 1e-8   # used to decide if we're at the kink of "abs" or "hypot"

struct AFloat
    val::Float64
    dot::Vector{Float64}
    ztol::Float64  # used in "abs" and "hypot" to decide if quantities are 0
end

# outer constructors for when this.dot and/or this.ztol aren't provided explicitly
AFloat(val::Float64, dot::Vector{Float64}) = AFloat(val, dot, DEFAULT_ZTOL)
AFloat(val::Float64, n::Int, ztol::Float64 = DEFAULT_ZTOL) = AFloat(val, zeros(n), ztol)

# display e.g. with "println"
Base.show(io::IO, x::AFloat) = print(io, "(", x.val, "; ", x.dot, ")")

# define promotion from Float64 to AFloat.
# Base.convert can't be used, because it doesn't know the intended dimension
# of the new "dot" vector.
Base.promote(uA::AFloat, uB::Float64) = (uA, AFloat(uB, length(uA.dot), uA.ztol))
Base.promote(uA::Float64, uB::AFloat) = reverse(promote(uB, uA))
Base.promote(uA::AFloat, uB::AFloat) = (uA, uB)

# for u::AFloat, define "u[1]" to mean "u". Helps handle vector/scalar outputs.
function Base.getindex(u::AFloat, i::Int)
    return (i == 1) ? u : throw(DomainError("i: must be 1"))
end

## define high-level generalized differentiation operations, given a mathematical
## function f composed from supported elemental operations, and written as though
## its input is a Vector

# compute:
#  y = the function value f(x), as a Vector, and
#  yDot = the LD-derivative f'(x; xDot)
#  yDot is the same type as xDot, which is either a Matrix or a Vector{Vector}
function eval_ld_derivative(
    f::Function,
    x::Vector{Float64},
    xDot::Matrix{Float64};
    ztol::Float64 = DEFAULT_ZTOL
)
    # use operator overloading to compute f(x) and f'(x; xDot)
    yLD = eval_AFloat_vec_output(f, x, xDot, ztol = ztol)

    # recover outputs from AFloats
    y = [vLD.val for vLD in yLD]
    yDot = reduce(vcat, [vLD.dot for vLD in yLD]')

    return y, yDot
end

function eval_ld_derivative(
    f::Function,
    x::Vector{Float64},
    xDot::Vector{Vector{Float64}};
    ztol::Float64 = DEFAULT_ZTOL
)
    # use operator overloading to compute f(x) and f'(x; xDot)
    yLD = eval_AFloat_vec_output(f, x, xDot, ztol = ztol)

    # recover outputs from AFloats
    y = [vLD.val for vLD in yLD]
    yDot = [vLD.dot for vLD in yLD]

    return y, yDot
end

function eval_AFloat_vec_output(
    f::Function,
    x::Vector{Float64},
    xDot::Matrix{Float64};
    ztol::Float64 = DEFAULT_ZTOL
)
    # express inputs as AFloats
    xLD = [AFloat(v, Vector(vDot), ztol) for (v, vDot) in zip(x, eachrow(xDot))]

    # use operator overloading to compute f(x) and f'(x; xDot)
    yLD = f(xLD)
    if !(yLD isa Vector)
        yLD = [yLD]
    end

    return yLD
end

function eval_AFloat_vec_output(
    f::Function,
    x::Vector{Float64},
    xDot::Vector{Vector{Float64}};
    ztol::Float64 = DEFAULT_ZTOL
)
    # express inputs as AFloats
    xLD = map(AFloat, x, xDot, ztol)

    # use operator overloading to compute f(x) and f'(x; xDot)
    yLD = f(xLD)
    if !(yLD isa Vector)
        yLD = [yLD]
    end

    return yLD
end

# compute:
#   y = the function value f(x), as a vector, and
#   yDot = the directional derivative f'(x; xDot::Vector), as a vector
function eval_dir_derivative(
    f::Function,
    x::Vector{Float64},
    xDot::Vector{Float64};
    ztol::Float64 = DEFAULT_ZTOL
)
    xDotMatrix = reshape(xDot, :, 1)
    (y, yDotMatrix) = eval_ld_derivative(f, x, xDotMatrix, ztol = ztol)
    yDot = vec(yDotMatrix)
    return y, yDot
end

# compute:
#   y = the function value f(x), and
#   yDeriv = the lexicographic derivative D_L f(x; xDot)
#   (xDot defaults to I if not provided, skipping the final linear solve.)
function eval_gen_derivative(
    f::Function,
    x::Vector{Float64},
    xDot::Matrix{Float64};
    ztol::Float64 = DEFAULT_ZTOL
)
    (y, yDot) = eval_ld_derivative(f, x, xDot, ztol = ztol)
    yDeriv = yDot / xDot
    return y, yDeriv
end

function eval_gen_derivative(
    f::Function,
    x::Vector{Float64};
    ztol::Float64 = DEFAULT_ZTOL
)
    (y, yDot) = eval_ld_derivative(f, x, Matrix{Float64}(I(length(x))), ztol = ztol)
    return y, yDot
end

# for scalar-valued f (ordinarily returning a Float64), compute:
#   y = the function value f(x), and
#   yGrad = the lexicographic gradient grad_L f(x; xDot),
#           which is a vector (and the transpose of the lex. derivative)
#   (xDot defaults to I if not provided, skipping the final linear solve.)
function eval_gen_gradient(
    f::Function,
    x::Vector{Float64},
    xDot::Matrix{Float64};
    ztol::Float64 = DEFAULT_ZTOL
)
    # use operator overloading to compute f(x) and f'(x; xDot)
    yLD = eval_AFloat_vec_output(f, x, xDot, ztol = ztol)

    # compute generalized gradient element
    return yLD[1].val, (yLD[1].dot'/ xDot)'
end

function eval_gen_gradient(
    f::Function,
    x::Vector{Float64};
    ztol::Float64 = DEFAULT_ZTOL
)
    # use operator overloading to compute f(x) and f'(x; xDot)
    yLD = eval_AFloat_vec_output(f, x, Matrix{Float64}(I(length(x))), ztol = ztol)

    return yLD[1].val, yLD[1].dot
end

# for a scalar-valued function f, compute:
#   y = the function value f(x), and
#   yCompass = the compass difference of f at x.
# x must be either a scalar or vector, and yCompass is the same type as x.
function eval_compass_difference(
    f::Function,
    x::Vector{Float64};
    ztol::Float64 = DEFAULT_ZTOL
)
    y = f(x)

    (length(y) == 1) ||
        throw(DomainError("f; this function is not scalar-valued"))

    fVec(u) = [f(u)[1]] # account for f returning either a Float64 or Vector{Float64}

    yCompass = copy(x)
    coordVec = zeros(length(x))
    for i in eachindex(x)
        coordVec[i] = 1.0
        _, yDotPlus = eval_dir_derivative(fVec, x, coordVec, ztol = ztol)
        _, yDotMinus = eval_dir_derivative(fVec, x, -coordVec, ztol = ztol)
        yCompass[i] = 0.5*(yDotPlus[1] - yDotMinus[1])
        coordVec[i] = 0.0
    end
    return y, yCompass
end

function eval_compass_difference(
    f::Function,
    x::Float64;
    ztol::Float64 = DEFAULT_ZTOL
)
    (y, yCompassVec) = eval_compass_difference(f, [x], ztol = ztol)
    return y, yCompassVec[1]
end

## define generalized differentiation rules for the simplest elemental operations

# macro to define (::AFloat, ::Float64) and (::Float64, ::AFloat) variants
# of a defined bivariate operation(::AFloat, ::AFloat).
macro define_mixed_input_variants(op)
    return quote
        Base.$op(uA::AFloat, uB::Float64) = $op(promote(uA, uB)...)
        Base.$op(uA::Float64, uB::AFloat) = $op(promote(uA, uB)...)
    end
end

# addition
function Base.:+(uA::AFloat, uB::AFloat)
    vVal = uA.val + uB.val
    vDot = uA.dot + uB.dot
    return AFloat(vVal, vDot, uA.ztol)
end
@define_mixed_input_variants +

# negative and subtraction
Base.:-(u::AFloat) = AFloat(-u.val, -u.dot, u.ztol)

Base.:-(uA::AFloat, uB::AFloat) = uA + (-uB)
@define_mixed_input_variants -

# multiplication
function Base.:*(uA::AFloat, uB::AFloat)
    vVal = uA.val * uB.val
    vDot = (uB.val * uA.dot) + (uA.val * uB.dot)
    return AFloat(vVal, vDot, uA.ztol)
end
@define_mixed_input_variants *

# reciprocal and division
function Base.inv(u::AFloat)
    vVal = inv(u.val)
    vDot = -u.dot / ((u.val)^2)
    return AFloat(vVal, vDot, u.ztol)
end

Base.:/(uA::AFloat, uB::AFloat) = uA * inv(uB)
@define_mixed_input_variants /

# integer powers
function Base.:^(u::AFloat, p::Int)
    if p == 0
        return 0.0*u + 1.0
    else
        vVal = (u.val)^p
        vDot = p * (u.val)^(p-1) * u.dot
        return AFloat(vVal, vDot, u.ztol)
    end
end

# exponential, logarithm, sine, and cosine
Base.exp(u::AFloat) = AFloat(exp(u.val), exp(u.val) * u.dot, u.ztol)

Base.log(u::AFloat) = AFloat(log(u.val), u.dot / u.val, u.ztol)

Base.sin(u::AFloat) = AFloat(sin(u.val), cos(u.val) * u.dot, u.ztol)

Base.cos(u::AFloat) = AFloat(cos(u.val), -sin(u.val) * u.dot, u.ztol)

# absolute value.
# When evaluating zDot, if any quantity "q" has abs(q) < u.ztol,
# then "q" is considered to be 0.
function Base.abs(u::AFloat)
    vVal = abs(u.val)
    if vVal > u.ztol
        s = sign(u.val)
    else
        s = 0.0
        for d in u.dot
            if abs(d) > u.ztol
                s = sign(d)
                break
            end
        end
    end
    vDot = s * u.dot
    return AFloat(vVal, vDot, u.ztol)
end

# bivariate max and min.
# Uses the identity max(uA,uB) = 0.5*(uA+uB+abs(uA-uB)), and similarly for min
Base.max(uA::AFloat, uB::AFloat) = 0.5*(uA + uB + abs(uA - uB))
@define_mixed_input_variants max

Base.min(uA::AFloat, uB::AFloat) = 0.5*(uA + uB - abs(uA - uB))
@define_mixed_input_variants min

# sqrt(uA^2 + uB^2); this is LinearAlgebra.hypot in Julia.
# uA.ztol is used as in abs(::AFloat)
function Base.hypot(uA::AFloat, uB::AFloat)
    vVal = hypot(uA.val, uB.val)
    if vVal > uA.ztol
        sA = uA.val / vVal
        sB = uB.val / vVal
    else
        sA = 0.0
        sB = 0.0
        for (dA, dB) in zip(uA.dot, uB.dot)
            sNorm = hypot(dA, dB)
            if sNorm > uA.ztol
                sA = dA / sNorm
                sB = dB / sNorm
                break
            end
        end
    end
    vDot = sA * uA.dot + sB * uB.dot
    return AFloat(vVal, vDot, uA.ztol)
end
@define_mixed_input_variants hypot

## overload comparisons to permit simple if-statements, but these conditions' values
## should not change under small perturbations in inputs.
Base.:<(uA::AFloat, uB::AFloat) = (uA.val < uB.val)
@define_mixed_input_variants <

Base.:>(uA::AFloat, uB::AFloat) = (uA.val > uB.val)
@define_mixed_input_variants >

Base.:>=(uA::AFloat, uB::AFloat) = (uA.val >= uB.val)
@define_mixed_input_variants >=

Base.:<=(uA::AFloat, uB::AFloat) = (uA.val <= uB.val)
@define_mixed_input_variants <=

end # module

max_stars_repo_name: kamilkhanlab/nonsmooth-forward-ad

max_issues_repo_name: kamilkhanlab/nonsmooth-forward-ad

max_forks_repo_name: kamilkhanlab/nonsmooth-forward-ad


text:
using ReservoirComputing
using MLJLinearModels
#model parameters
const approx_res_size = 30
const radius = 1.2
const activation = tanh
const degree = 6
const sigma = 0.1
const beta = 0.0
const alpha = 1.0
const nla_type = NLADefault()
const in_size = 3
const out_size = 3
const extended_states = false
const delta = 0.5

const train_len = 50
const predict_len = 12
data = ones(Float64, in_size, 100)
train = data[:, 1:1+train_len-1]
test = data[:, train_len:train_len+predict_len-1]

#constructor 1
esn = ESN(approx_res_size,
    train,
    degree,
    radius,
    activation = activation,
    sigma = sigma,
    alpha = alpha,
    nla_type = nla_type,
    extended_states = extended_states)

#test train Ridge
t = Ridge(beta, Analytical())
@test isequal(t.lambda, beta)
@test isequal(t.solver, Analytical())
W_out = ESNtrain(t, esn)
@test size(W_out) == (out_size, esn.res_size)

#test train Lasso
t = Lasso(beta, ProxGrad())
@test isequal(t.lambda, beta)
@test isequal(t.solver, ProxGrad())
W_out = ESNtrain(t, esn)
@test size(W_out) == (out_size, esn.res_size)

#test train ElastNet
t = ElastNet(beta, beta, ProxGrad())
@test isequal(t.lambda, beta)
@test isequal(t.gamma, beta)
@test isequal(t.solver, ProxGrad())
W_out = ESNtrain(t, esn)
@test size(W_out) == (out_size, esn.res_size)

#test train RobustHuber
t = RobustHuber(delta, beta, 0.0, Newton())
@test isequal(t.delta, delta)
@test isequal(t.lambda, beta)
@test isequal(t.gamma, 0.0)
@test isequal(t.solver, Newton())
W_out = ESNtrain(t, esn)
@test size(W_out) == (out_size, esn.res_size)


text:
module Tower

include("nodes.jl")

end # module


text:
import Base: ==, !=, +, -, *, /, ^, <, >, |, <=, >=, ~, :, !, <<, >>, >>>, &,
             cld, fld, add12, mod, rem, div

function defbasics(::Type{ST}, ::Type{BT}) where {ST,BT}
    @eval begin
        Base.@pure Base.values(::$ST{V}) where V = V::$BT
        Base.@pure Base.values(::Type{$ST{V}}) where V = V::$BT

        (::Type{<:$ST})(val::Val{V}) where V = $ST{$BT(V)}()

        Base.eltype(::$ST) = $BT
        Base.eltype(::Type{<:$ST}) = $BT

        Base.typemax(::$ST) = $ST{Base.typemax($BT)}()
        Base.typemax(::Type{$ST}) = $ST{Base.typemax($BT)}()

        Base.typemin(::$ST) = $ST{typemin($BT)}()
        Base.typemin(::Type{$ST}) = $ST{typemin($BT)}()

        promote_rule(::Type{<:$ST}, ::Type{$BT}) = $BT

        (::Type{$BT})(::$ST{X}) where X = X::$BT
        (::Type{<:$ST{<:Any}})(x::$ST{X}) where X = x::$ST{X}
        (::Type{<:$ST{<:Any}})(x::$BT) = $ST{x}()
        (::Type{<:$ST{<:Any}})(x::BaseNumber) = $ST($BT(x))

        ofeltype(::Type{$BT}, val::$ST) = val
        ofeltype(::Type{<:$ST{<:Any}}, val::$BT) = val
        ofeltype(::Type{<:$ST{<:Any}}, val::$ST) = val
        ofeltype(::Type{$BT}, val::$BT) = val
        ofeltype(::$ST, val::$BT) = val
        ofeltype(::$BT, val::$BT) = val
        ofeltype(::$ST, val::$ST) = val
        ofeltype(::$BT, val::$ST) = val

        isstatic(::$ST) = true
        isstatic(::Type{<:$ST}) = true
        isstatic(::$BT) = false
        isstatic(::Type{$BT}) = false

        tostatic(x::$ST) = x
        tostatic(x::Type{<:$ST}) = x
        tostatic(x::$BT) = $ST{x}()
        tostatic(x::Type{$BT}) = $ST

        seek_static_val(::Type{$BT}, val::Val{V}) where V = $ST{V}()
    end
end

function defmath(::Type{ST}, ::Type{BT}) where {ST,BT}
    STOne = ST(BT(1))
    STOneType = typeof(STOne)
    STZero = ST(BT(0))
    STZeroType = typeof(STZero)

    for f in (:-, :+, :*, :cld, :fld, :mod, :rem, :div)
        @eval begin
            $f(x::$ST{V1}, y::$ST{V2}) where {V1,V2} = $ST{$f(V1::$BT, V2::$BT)}()
        end
    end

    @eval begin
        (-)(::$ST{V}) where V = $ST{-V::$BT}()
        function add12(x::$ST, y::$ST)
            x, y = ifelse(abs(y) > abs(x), (y, x), (x, y))
            Base.canonicalize2(x, y)
        end
        function Base.muladd(x::$ST{X}, y::$ST{Y}, z::$ST{Z}) where {X,Y,Z}
            $ST{muladd(X::$BT,Y::$BT,Z::$BT)}()
        end

        Base.max(x::$ST{X}, y::$ST{Y}) where {X,Y} = X::$BT > Y::$BT ? x : y
        Base.min(x::$ST{X}, y::$ST{Y}) where {X,Y} = X::$BT > Y::$BT ? y : x
        Base.minmax(x::$ST{X}, y::$ST{Y}) where {X,Y} = X::$BT > Y::$BT ? (y, x) : (x, y)

        @pure Base.one(::$ST) = $STOne
        @pure Base.one(::Type{<:$ST}) = $STOne
        @pure Base.oneunit(::$ST) = $STOne
        @pure Base.oneunit(::Type{<:$ST}) = $STOne

        @pure Base.isone(::$STOneType) = true
        @pure Base.isone(::$ST{T}) where T = false

        @pure Base.zero(::$ST) = $STZero
        @pure Base.zero(::Type{<:$ST}) = $STZero

        @pure Base.iszero(::$STZeroType) = true
        @pure Base.iszero(::$ST{T}) where T = false
    end

    return STOne, STZero
end

function defbool(::Type{ST}, ::Type{BT}) where {ST,BT}
    @eval begin
        @pure Base.:(==)(::$ST{V}, ::$ST{V}) where V = true
        @pure Base.:(==)(::$ST{V1}, ::$ST{V2}) where {V1,V2} = false

        @pure Base.:(!=)(::$ST{V}, ::$ST{V}) where V = false
        @pure Base.:(!=)(::$ST{V1}, ::$ST{V2}) where {V1,V2} = true

        Base.:(>)(::$ST{V1}, ::$ST{V2}) where {V1,V2} = V1::$BT > V2::$BT
        Base.:(>=)(::$ST{V1}, ::$ST{V2}) where {V1,V2} = V1::$BT >= V2::$BT

        Base.:(<)(::$ST{V1}, ::$ST{V2}) where {V1,V2} = V1::$BT < V2::$BT
        Base.:(<=)(::$ST{V1}, ::$ST{V2}) where {V1,V2} = V1::$BT <= V2::$BT

        Base.round(::$ST{X}, r::RoundingMode) where X = $ST{round(X::$BT, r)}()
        Base.isless(::$ST{V1}, ::$ST{V2}) where {V1,V2} = isles(V1::$BT, V2::$BT)
    end
end

#=
function defsigned(::Type{ST}, ::Type{BT}) where {ST,BT}
    flipsign(::$ST{X}, ::$ST{Y}) where {X,Y} = flipsign_int(x, y)
end
flipsign(x::T, y::T)

# these require SInt to be set
=#


text:
using AxisArrayConversion: to, AxisArrayConversion
const AAC = AxisArrayConversion
using Test
import AxisKeys; const AK = AxisKeys
import AxisArrays; const AA = AxisArrays
import DimensionalData; const DD = DimensionalData

using Documenter: doctest
doctest(AxisArrayConversion)

@testset "name_axes" begin
    ax = (a=[1,2], b=3:4)
    @test @inferred(AAC.name_axes(ax)) === ax
    ax = (1:2, [3,4], [5,6,10])
    expected = NamedTuple{(:x1,:x2,:x3)}(ax)
    @test @inferred(AAC.name_axes(ax)) === expected
end

@testset "NamedTuple" begin

    nt = (axes=(a=1:3,), values=[10,20,30])
    @test AAC.to(NamedTuple, nt) === nt

    nt = (axes=(1:3,), values=[10,20,30])
    @test AAC.to(NamedTuple, nt) === (axes=(x1=nt.axes[1],), values=nt.values)
end

@testset "Array" begin

    @testset "Vector" begin
        vec = [1,2,3]
        nt = @inferred to(NamedTuple, vec)
        @test propertynames(nt) === (:axes, :values)
        @test nt.values === vec
        @test nt.axes   === (x1 = axes(vec,1),)
        vec2 = @inferred to(Array, nt)
        @test vec2 === vec
        vec3 = @inferred to(Vector{Int}, nt)
        @test vec3 === vec
        vec_f64 = @inferred to(Array{Float64}, nt)
        @test vec_f64 == vec
        @test typeof(vec_f64) === Vector{Float64}
    end

    @testset "Matrix" begin
        mat = rand(Int, (2,3))

        nt = @inferred to(NamedTuple, mat)
        @test nt === (
            axes = NamedTuple{(:x1, :x2)}(axes(mat)),
            values = mat
        )
        @test @inferred(to(Matrix, nt)) === mat
        @test @inferred(to(Array,  nt)) === mat
    end

    @testset "lax conversion" begin
        # we might want to disable lax conversion by default?
        nt = (axes=(a=[2,3,4],), values=[1,2,3])
        @test to(Vector, nt) === nt.values

        ka::AK.KeyedArray = to(AK.KeyedArray, nt)
        @test to(Vector, ka) === nt.values

    end
end

@testset "check_consistency" begin
    AAC.check_consistency((axes=(a=1:2,), values=[1,20]))
    AAC.check_consistency((axes=(1:2,), values=[1,20]))
    AAC.check_consistency((values=[1,20], axes=(a=1:2,)))
    @test_throws ArgumentError AAC.check_consistency((values=[1,20], axes=(a=1:3,)))
    @test_throws ArgumentError AAC.check_consistency((values=[1,20], axes=(a=1:2,b=3:4)))
    @test_throws ArgumentError AAC.check_consistency((values=[1,20], axs=(a=1:2,)))
end

@testset "to" begin
    axs = (a=1:2, b=[10,20, 30])
    nt = (axes=axs, values=(randn(2,3)))
    @test nt === to(NamedTuple, nt)

    aa = to(AA.AxisArray, nt)
    @test aa isa AA.AxisArray

    da = to(DD.DimArray, nt)
    @test da isa DD.DimArray

    ka = to(AK.KeyedArray, nt)
    @test ka isa AK.KeyedArray

    @test aa == to(AA.AxisArray, aa)
    @test aa == to(AA.AxisArray, da)
    @test aa == to(AA.AxisArray, ka)
    @test aa == to(AA.AxisArray, nt)

    @test da == to(DD.DimArray, aa)
    @test da == to(DD.DimArray, da)
    @test da == to(DD.DimArray, ka)
    @test da == to(DD.DimArray, nt)

    @test ka == to(AK.KeyedArray, aa)
    @test ka == to(AK.KeyedArray, da)
    @test ka == to(AK.KeyedArray, ka)
    @test ka == to(AK.KeyedArray, nt)

    @test nt == to(NamedTuple, aa)
    @test nt == to(NamedTuple, da)
    @test nt == to(NamedTuple, ka)
    @test nt == to(NamedTuple, nt)

    @inferred to(NamedTuple, aa)
    @inferred to(NamedTuple, da)
    @inferred to(NamedTuple, ka)
    @inferred to(NamedTuple, nt)

    saa = to(AAC.SimpleAxisArray, aa)

    arrs = [aa, da, ka, nt, saa]
    for a1 in arrs
        for a2 in arrs
            T = typeof(a1)
            b1 = to(T, a2)
            @test b1 isa T
            @test b1 == a1

            TUpper = @inferred AAC.roottype(T)
            @test T <: TUpper
            @test T !== TUpper
            @test to(TUpper, a2) == to(T, a2)
            @test typeof(to(TUpper, a2)) === typeof(to(T, a2))

            AAC.check_consistency(a1)

            @inferred to(T, a2)
            # try
            #     @inferred to(T, a2)
            # catch err
            #     @warn """Inference of `to(T, o)` broken for
            #     T = $T
            #     typeof(o) = $(typeof(o))
            #     """
            # end
        end
    end
end


text:
using LowRankApprox, FillArrays, Test

@testset "LowRankMatrix" begin
    @testset "Constructors" begin
        @test Matrix(LowRankMatrix(Zeros(10,5))) == zeros(10,5)

        @test LowRankMatrix{Float64}(Zeros(10,5)) == LowRankMatrix(Zeros(10,5)) ==
                    LowRankMatrix{Float64}(Zeros(10,5),1) == LowRankMatrix{Float64}(Zeros{Int}(10,5),1) ==
                    LowRankMatrix{Float64}(zeros(10,5)) == LowRankMatrix(zeros(10,5))  ==
                    LowRankMatrix{Float64}(zeros(Int,10,5))

        @test isempty(LowRankMatrix{Float64}(zeros(10,5)).U)
        @test isempty(LowRankMatrix{Float64}(zeros(10,5)).V)
        @test isempty(LowRankMatrix(Zeros(10,5)).U)
        @test isempty(LowRankMatrix(Zeros(10,5)).V)

        @test rank(LowRankMatrix(Zeros(10,5))) == 0

        @test Matrix(LowRankMatrix(Ones(10,5))) == fill(1.0,10,5)
        @test LowRankMatrix{Float64}(Ones(10,5)) == LowRankMatrix(Ones(10,5)) ==
                    LowRankMatrix{Float64}(Ones{Int}(10,5))
        @test LowRankMatrix{Float64}(fill(1.0,10,5)) == LowRankMatrix(fill(1.0,10,5))  ==
                    LowRankMatrix{Float64}(fill(1,10,5))
        @test rank(LowRankMatrix(Ones(10,5))) == 1

        @test LowRankMatrix(Ones(10,5)) ≈ LowRankMatrix(fill(1.0,10,5))

        @test Matrix(LowRankMatrix(Ones(10,5))) == fill(1.0,10,5)
        @test LowRankMatrix{Float64}(Ones(10,5)) == LowRankMatrix(Ones(10,5))
        @test rank(LowRankMatrix(fill(1.0,10,5))) == 1

        x = 2
        @test Matrix(LowRankMatrix(Fill(x,10,5))) ≈ fill(x,10,5)
        @test LowRankMatrix{Float64}(Fill(x,10,5)) == LowRankMatrix(Fill(x,10,5)) ==
                    LowRankMatrix{Float64}(Fill{Float64}(x,10,5))
        @test LowRankMatrix{Float64}(fill(x,10,5)) == LowRankMatrix(fill(x,10,5))  ==
                    LowRankMatrix{Float64}(fill(x,10,5))
        @test rank(LowRankMatrix(Fill(x,10,5))) == 1

    end

    @testset "LowRankMatrix algebra" begin
        A = LowRankApprox._LowRankMatrix(randn(20,4), randn(12,4))
        @test Matrix(2*A) ≈ Matrix(A*2) ≈ 2*Matrix(A)

        B = LowRankApprox._LowRankMatrix(randn(20,2), randn(12,2))
        @test Matrix(A+B) ≈ Matrix(A) + Matrix(B) ≈ Matrix(Matrix(A) + B) ≈
                    Matrix(A + Matrix(B))
        @test rank(A+B) ≤ rank(A) + rank(B)

        @test Matrix(A-B) ≈ Matrix(A) - Matrix(B) ≈ Matrix(Matrix(A) - B) ≈
                    Matrix(A - Matrix(B))
        @test rank(A-B) ≤ rank(A) + rank(B)

        B = LowRankApprox._LowRankMatrix(randn(12,2), randn(14,2))

        @test A*B isa LowRankMatrix
        @test rank(A*B) == size((A*B).U,2) == 4

        @test Matrix(A)*Matrix(B) ≈ Matrix(A*Matrix(B)) ≈ Matrix(Matrix(A)*B) ≈ Matrix(A*B)

        B = LowRankApprox._LowRankMatrix(randn(10,2), randn(14,2))
        @test_throws DimensionMismatch A*B
        @test_throws DimensionMismatch B*A

        B = randn(12,14)
        @test A*B isa LowRankMatrix
        @test rank(A*B) == size((A*B).U,2) == 4
        @test Matrix(A)*B ≈ Matrix(A*B)

        B = randn(20,20)
        @test B*A isa LowRankMatrix
        @test rank(B*A) == size((B*A).U,2) == 4
        @test B*Matrix(A) ≈ Matrix(B*A)

        v = randn(12)
        @test all(LowRankApprox.mul!(randn(size(A,1)), A, v) .=== A*v )
        @test A*v ≈ Matrix(A)*v
    end
end


text:
include("random_sky.jl")
include("gen_noise_map.jl")
using ProgressMeter

e_des = Array{Float64, 1}(undef,100)
d_des = Array{Float64, 1}(undef,100)

for i in 1:10
    global dir=false
    mappa = zeros(NN, NN)
    hit   = zeros(NN, NN)
    N_baselines = 10
    dataset_baselines = Array{Main.strip2.TOD_fake_pointing, 1}(undef, N_baselines)

    for i = 1:N_baselines
        dataset = strip2.observe_sky(NN, conv_map; noise=true, direction_l_r=dir)
        m, h = strip2.binned_map(NN, dataset.time_order_data, dataset.pointing)
        mappa .+= m
        hit   .+= h
        dataset_baselines[i] = dataset
        global dir =  true ⊻ dir
    end

    map_des = strip2.destriper(NN, dataset_baselines)
    map_d = reshape(map_des, NN, NN)

    ell_max = 5000.0
    delta_ell = 50.0

    win_map_dest = strip2.windowing!(NN, (map_d .- avg_dest))

    a, b = strip2.get_power_spectrum(win_map_dest, win_map_dest,
                                            ell_max, delta_ell, pix_size, NN)
    e_des .+= a
    d_des .+= b
end

e_des ./= 10.0
d_des ./= 10.0

Multiply = DlTT[Int.(round.(e_des)) .- 1] ./ d_des

R = 0.0001:0.0001:0.1

mag_k = ((2.0 * π) ./ (R .+ 0.01 )).^(5.0 / 3.0)
pix_to_rad = pix_size / 60. * π  / 180.  # Scale factor between pix to rad
ell_scale_factor = 2. * π / pix_to_rad
ell_atm = R .* ell_scale_factor

(1 .+ ( max.(Freq, Freq[2]) ./ fknee).^(-alpha)) .* (sigma^2)

p(f)=(1 + f/f_knee )^(-α) * σ

pp=plot(ell, DlTT./maximum(DlTT), label="CMB Temp. normalized", yaxis=:log, xlims =(0.00001, 1000), ylims=(0.1,1), xlabel="l", ylabel="Dl/max(Dl)", size=(500,500))
pp=plot!(ell, DlBB./maximum(DlBB), label="CMB BB normalized", size=(500,500))
pp=plot!(ell_atm, mag_k ./ maximum(mag_k), label="Atmospheric normalized", size=(500,500))


text:
function simulate_paths(space, c_store, mparams, opts, smax, smin;
                        misspec_damages = :base, alt_damages = false, d1 = [], d2 = [])
"""
Simulate model, outer function
"""

    # Simulating
    opts[:sim] = true

    # Boolean option for doing alternative (but correctly specified) damage parameters
    opts[:alt_damages] = alt_damages
    if opts[:alt_damages]
        mparams[:d1] = d1
        mparams[:d2] = d2
    end

    # Correctly specified or misspecified damages
    mparams[:damages] = misspec_damages

    # Initialize state array
    if opts[:learn_e]
        states = zeros(mparams[:sim_length] + 1, 5)
    else
        states = zeros(mparams[:sim_length] + 1, 3)
    end

    ## Initialize storage arrays for simulations
    # State array
    s_store = Vector{Array}(undef, mparams[:num_runs])
    # d_2 mean quadrature
    loc_store = Vector{Vector{Array}}(undef, mparams[:num_runs])
    # Carbon tax channels
    channels = Vector{Dict{Symbol,Array}}(undef, mparams[:num_runs])
    # Fractional net output quadrature
    dam_store = Vector{Vector{Array}}(undef, mparams[:num_runs])
    # Controls
    controls_store = Vector{Vector}(undef, mparams[:num_runs])
    # Investment rate
    inv_rate = Vector{Vector}(undef, mparams[:num_runs])
    # Jumps up array
    jumps_up = Vector{Array}(undef, mparams[:num_runs])
    # Jumps down array
    jumps_down = Vector{Array}(undef, mparams[:num_runs])

    # Initial state
    if opts[:learn_e]
        states[1, :] = [mparams[:kstart], mparams[:Mstart], mparams[:expstart], mparams[:varstart], mparams[:lossstart]]
    else
        states[1, :] = [mparams[:kstart], mparams[:Mstart], mparams[:lossstart]]
    end

    # Exogenous variables
    exog_states = simulate_exog_vars(mparams)

    # Load shocks, pre-made so common shocks across different frameworks in the paper
    if mparams[:damages] == :base && mparams[:num_runs] > 1000
        temp = BSON.load("results/loaded_sim_values_large.bson")
    elseif mparams[:damages] == :base && mparams[:num_runs] <= 1000
        temp = BSON.load("results/loaded_sim_values_large_adapt.bson")
    else
        println("Simulating with misspecified damage function.")
        temp = BSON.load("results/loaded_sim_values_alt.bson")
    end
    sim_values = temp["sim_values"]

    # How to generate lognormal random variables:
    # exp.(rand(Normal(location_parameter, sqrt(scale_parameter), number_of_draws))
    # shocks = [exp.(rand(Normal(-0.5885394796543533,sqrt(1.1770789593087065)), 101)) for i = 1:50000]
    # bson("results/loaded_sim_values_large.bson", Dict{String, Array{Array{Float64}}}("shocks" => shocks))
    if opts[:alt_damages]
        shocks_dict = BSON.load("results/loaded_sim_values_shocks.bson")
        extract(shocks_dict)
        sim_values[2] = shocks
    end

    # Initialize array of variables to send to simulator
    sim_parameters = Array{Array}(undef, mparams[:num_runs], 1)
    sim_parameters = [simParameters(states, exog_states, mparams, space, opts, c_store,
        sim_values[1][i], sim_values[2][i], smax, smin) for i = 1:mparams[:num_runs]]

    # Simulate
    sim_results = pmap(simulate_inner, sim_parameters)

    # Pull out and store results
    for run = 1:mparams[:num_runs]
        # Controls
        controls_store[run] = sim_results[run][1]
        # States
        s_store[run] = sim_results[run][2]
        # d_2 mean quadrature
        loc_store[run] = sim_results[run][3]
        # Fractional net output
        dam_store[run] = sim_results[run][4]
        # Carbon tax channels
        channels[run] = sim_results[run][5]
        # Investment rate
        inv_rate[run] = sim_results[run][6]
    end

    # Functions to pull out and separate states and controls
    states_func(M) = [[M[i][:, j] for i = 1:length(M)] for j = 1:size(s_store[1], 2)]
    controls_func(M) = [[[M[i][j][k] for j = 1:mparams[:sim_length]] for i = 1:length(M)] for k = 1:2]

    # Pull out states
    if opts[:learn_e]
        capital, co2, prior_loc, prior_scale, damage = states_func(s_store)
    else
        capital, co2, damage = states_func(s_store)
    end

    # Pull out controls
    consumption, abatement = controls_func(controls_store)

    # Initialize dictionary for mean tax channels
    channels_out = Dict{Symbol,Array}()
    channels_out_samples = Dict{Symbol,Array}()

    # Get 90% confidence bounds on carbon tax
    channels_out[:tax_lower], channels_out[:tax_upper] = get_bounds([channels[i][:scc_direct] for i = 1:mparams[:num_runs]], .05)

    # Pull out individual channels to average out
    for key in keys(channels[1])
        channels_out[key] = [channels[i][key] for i = 1:mparams[:num_runs]]
        # Keep individual run tax channels
        channels_out_samples[key] = [channels[i][key] for i = 1:mparams[:num_runs]]
        # Average channels
        channels_out[key] = mean(hcat(channels_out[key]...), dims = 2)
    end

    if ~opts[:learn_e]
        prior_loc = 0.
        prior_scale = 0.
    end

    return consumption, abatement, inv_rate, capital, co2,
        prior_loc, prior_scale, damage, sim_values, loc_store, dam_store,
        channels_out, channels_out_samples, exog_states

end

function simulate_paths_adapt(space, c_store, mparams, opts, smax, smin;
    alt_damages = false)
"""
Simulate model, outer function for adapted grid
"""

    include("initialization.jl")

    # Boolean option for doing alternative (but correctly specified) damage parameters
    opts[:alt_damages] = alt_damages

    # Simulating
    opts[:sim] = false

    # Initialize state array
    if opts[:learn_e]
        states = zeros(mparams[:sim_length] + 1, 5)
    else
        states = zeros(mparams[:sim_length] + 1, 3)
    end

    ## Initialize storage arrays for simulations
    # State array
    s_store = Vector{Array}(undef, mparams[:num_runs])
    # d_2 mean quadrature
    loc_store = Vector{Vector{Array}}(undef, mparams[:num_runs])
    # Carbon tax channels
    channels = Vector{Dict{Symbol,Array}}(undef, mparams[:num_runs])
    # Fractional net output quadrature
    dam_store = Vector{Vector{Array}}(undef, mparams[:num_runs])
    # Controls
    controls_store = Vector{Vector}(undef, mparams[:num_runs])
    # Investment rate
    inv_rate = Vector{Vector}(undef, mparams[:num_runs])

    # Initial state
    if opts[:learn_e]
        states[1, :] = [mparams[:kstart], mparams[:Mstart], mparams[:expstart], mparams[:varstart], mparams[:lossstart]]
    else
        states[1, :] = [mparams[:kstart], mparams[:Mstart], mparams[:lossstart]]
    end

    # Exogenous variables
    exog_states = simulate_exog_vars(mparams)

    # Load shocks, pre-made so common shocks across different frameworks in the paper
    temp = BSON.load("results/loaded_sim_values_large_adapt.bson")

    sim_values = temp["sim_values"]

    # Initialize array of variables to send to simulator
    sim_parameters = Array{simParameters}(undef, mparams[:num_runs], 1)
    sim_parameters =
        [simParameters(states, exog_states, mparams, space, opts, c_store,
            sim_values[1][i], sim_values[2][i], smax, smin) for i = 1:mparams[:num_runs]]

    # Simulate
    sim_results = pmap(simulate_inner, sim_parameters)

    # Pull out and store results
    for run = 1:mparams[:num_runs]
        # Controls
        controls_store[run] = sim_results[run][1]
        # States
        s_store[run] = sim_results[run][2]
        # d_2 mean quadrature
        loc_store[run] = sim_results[run][3]
        # Fractional net output
        dam_store[run] = sim_results[run][4]
    end

    # Functions to pull out and separate states and controls
    states(M) = [[M[i][:, j] for i = 1:length(M)] for j = 1:size(s_store[1], 2)]
    controls(M) = [[[M[i][j][k] for j = 1:mparams[:sim_length]] for i = 1:length(M)] for k = 1:2]

    # Pull out states
    if opts[:learn_e]
        capital, co2, prior_loc, prior_scale, damage = states(s_store)
    else
        capital, co2, damage = states(s_store)
    end

    # Pull out controls
    consumption, abatement = controls(controls_store)

    if ~opts[:learn_e]
        prior_loc = 0.
        prior_scale = 0.
    end

    return consumption, abatement, capital, co2,
        prior_loc, prior_scale, damage, loc_store, dam_store

end

function simulate_paths_adapt(space, c_store, mparams, opts;
     alt_damages = false)
"""
Simulate model for adaptive grid on first iteration
"""

        include("initialization.jl")

        # Switch for doing alternative (but correctly specified) damage parameters
        opts[:alt_damages] = alt_damages

        # Simulating
        opts[:sim] = false

        # Initialize state array
        if opts[:learn_e]
            states = zeros(mparams[:sim_length] + 1, 5)
        else
            states = zeros(mparams[:sim_length] + 1, 3)
        end

        ## Initialize storage arrays for simulations
        loc_next = Vector{Array}(undef, mparams[:sim_length] + 1)
        dam_next = Vector{Array}(undef, mparams[:sim_length] + 1)
        loc_next[1] = mparams[:expstart]*ones(mparams[:nqnodes], mparams[:nqnodes])
        dam_next[1] = mparams[:lossstart]*ones(mparams[:nqnodes], mparams[:nqnodes])
        s_store = Vector{Array}(mparams[:num_runs_adapt])
        loc_store = Vector{Vector{Array}}(undef, mparams[:num_runs_adapt])
        dam_store = Vector{Vector{Array}}(undef, mparams[:num_runs_adapt])
        controls = Vector{Vector}(undef, mparams[:sim_length])
        controls_store = Vector{Vector}(undef, mparams[:num_runs_adapt])
        loc_next[1] = mparams[:expstart]*ones(mparams[:nqnodes], mparams[:nqnodes])
        dam_next[1] = mparams[:lossstart]*ones(mparams[:nqnodes], mparams[:nqnodes])

        # Initial state
        if opts[:learn_e]
            states[1, :] = [mparams[:kstart], mparams[:Mstart], mparams[:expstart], mparams[:varstart], mparams[:lossstart]]
        else
            states[1, :] = [mparams[:kstart], mparams[:Mstart], mparams[:lossstart]]
        end

        # Exogenous variables
        exog_states = simulate_exog_vars(mparams)

        # Load shocks, pre-made so common shocks across different frameworks in the paper
        temp = BSON.load("results/loaded_sim_values_large_adapt.bson")

        sim_values = temp["sim_values"]

        # Initialize array of variables to send to simulator
        sim_parameters = Array{simParameters}(undef, mparams[:num_runs], 1)
        sim_parameters =
            [simParameters(states, exog_states, mparams, space, opts,
            c_store, sim_values[1][i], sim_values[2][i], smax, smin) for i = 1:mparams[:num_runs]]

        # Simulate
        sim_results = pmap(simulate_inner, sim_parameters)

        # Pull out and store results
        for run = 1:mparams[:num_runs]
            # Controls
            controls_store[run] = sim_results[run][1]
            # States
            s_store[run] = sim_results[run][2]
            # d_2 mean quadrature
            loc_store[run] = sim_results[run][3]
            # Fractional net output
            dam_store[run] = sim_results[run][4]
        end

        # Functions to pull out and separate states and controls
        states(M) = [[M[i][:, j] for i = 1:length(M)] for j = 1:size(s_store[1], 2)]
        controls(M) = [[[M[i][j][k] for j = 1:mparams[:sim_length]] for i = 1:length(M)] for k = 1:2]

        # Pull out states
        if opts[:learn_e]
            capital, co2, prior_loc, prior_scale, damage = states(s_store)
        else
            capital, co2, damage = states(s_store)
        end

        # Pull out controls
        consumption, abatement = controls(controls_store)

        if ~opts[:learn_e]
            prior_loc = 0.
            prior_scale = 0.
        end

        return consumption, abatement, capital, co2,
            prior_loc, prior_scale, damage, loc_store, dam_store

    end

function simulate_inner(sim_parameters_in)
"""
Simulate model, inner function
"""

    # Current state array
    states = sim_parameters_in.states
    # Exogenous variables
    exog_states = sim_parameters_in.exog_states
    # Model parameters
    mparams = sim_parameters_in.mparams
    # Approximation space
    space = sim_parameters_in.space
    # Boolean options
    opts = sim_parameters_in.opts
    # Coefficient vector
    c_store = sim_parameters_in.coefficients

    # Randomly sampled damage parameters
    sim_values = Vector{Array}(undef, 2)
    sim_values[1] = [sim_parameters_in.dam_exponent]
    sim_values[2] = sim_parameters_in.dam_coefficient

    if opts[:adapted]
        smax = sim_parameters_in.smax
        smin = sim_parameters_in.smin
    end

    # Initialize d_2 mean and fractional net output quadrature
    loc_next = Vector{Array}(undef, mparams[:sim_length])
    dam_next = Vector{Array}(undef, mparams[:sim_length])

    # Initialize control vector
    controls = Vector{Vector}(undef, mparams[:sim_length])

    # Initial control guess
    qiguess = [.2, exog_states[:Psi][1]/2]
    cost = 0.
    cons = 0.

    # Initialize investment rate
    inv_rate = Vector{Float64}(undef, mparams[:sim_length])

    # Initialize abatement rate
    mu = Vector{Float64}(undef, mparams[:sim_length])

    # Simulate forward
    for t in 1:mparams[:sim_length]

        if opts[:adapted]
            space = fundefn(:cheb, mparams[:approximation_level], smin[t + 1, :], smax[t + 1, :])
        end

        # Load coefficient vector for continuation value function
        space[:c_cont] = c_store[t + 1]

        # Load year
        space[:year] = t - 1
        mparams[:year] = t - 1

        # Time t exogenous variables
        value_vars = reduce_dict(exog_states, t)

        # Maximize Bellman
        results = maximize_bellman(stateVariables(vec(states[t, :]), qiguess, value_vars, mparams, space, opts, space, 0))

        # Pull out consumption
        cons = results[1][1]

        # Abatement cost
        cost = results[1][2]

        # Abatement rate
        mu[t] = (cost/value_vars[:Psi]).^(1/mparams[:a2])

        # Next period's guess
        qiguess = [max(1, cons), cost]
        qi = [cons, cost]
        v_out = results[2]

        # Store controls in vector
        controls[t] = [qi[1], mu[t]]

        # Store investment rate
        if opts[:learn_e]
            inv_rate[t] = -(qi[1] .+ (qi[2] - 1)*states[t, 1].^mparams[:kappa].*states[t, 5])./(states[t, 1].^mparams[:kappa].*states[t, 5])
        else
            inv_rate[t] = -(qi[1] .+ (qi[2] - 1)*states[t, 1].^mparams[:kappa].*states[t, 3])./(states[t, 1].^mparams[:kappa].*states[t, 3])
        end

        # Compute actual next state and distributions for the next state
        states[t + 1, :], loc_next[t], dam_next[t] = transitions(states[t, :], qi, value_vars, mparams, opts, sim_values[1][1], sim_values[2][t])

    end

    # Compute tax and perform decomposition
    if opts[:sim]
        if opts[:adapted]
            tax_parameters = taxParameters(states, controls, exog_states, mparams, space, opts, c_store, smax, smin, loc_next, dam_next)
            channels = compute_tax(tax_parameters)
        else
            tax_parameters = taxParameters(states, controls, exog_states, mparams, space, opts, c_store, [0. 0.; 0. 0.], [0. 0.; 0. 0.], loc_next, dam_next)
            channels = compute_tax(tax_parameters)
        end

        return controls, states, loc_next, dam_next, channels, inv_rate
    else
        return controls, states, loc_next, dam_next
    end

end

function simulate_exog_vars(mparams)

    # Time periods vector
    t = 0:mparams[:sim_length]

    # Non-industrial emissions
    exog_states = Dict{Symbol,Array}(:Bs => mparams[:B0]*exp.(mparams[:gB]*t))

    # Labor
    exog_states[:Lt] = mparams[:L0] .+ (mparams[:Linf]-mparams[:L0])*(1.0 .- exp.(-mparams[:dL]*t))

    # Labor growth
    exog_states[:gL] = mparams[:dL]./((mparams[:Linf]/(mparams[:Linf]-mparams[:L0]))*exp.(mparams[:dL]*t) .- 1)

    # Technological growth rate
    exog_states[:gA] = mparams[:gA0]*exp.(-mparams[:dA]*t)

    # Technology
    exog_states[:At] = mparams[:A0]*exp.( (mparams[:gA0]/mparams[:dA])*(1 .- exp.(-mparams[:dA]*t)))

    # Growth rate for exogenous growth in labor and productivity
    exog_states[:effkgrowth] = exp.(-(exog_states[:gA]+exog_states[:gL]))

    # Effective discount factor
    exog_states[:betaeff] = exp.((-mparams[:delta] .+ exog_states[:gA]*mparams[:rho] + exog_states[:gL]))

    # Exogenous decarbonization of production
    exog_states[:sigma] = mparams[:sigma0]*exp.((mparams[:gsigma0]/mparams[:dsigma])*(1 .- exp.(-mparams[:dsigma]*t)))
    exog_states[:emint] = exog_states[:sigma].*exog_states[:At].*exog_states[:Lt]

    # Abatement cost function coefficient
    exog_states[:Psi] = (mparams[:a0]*exog_states[:sigma]./mparams[:a2]).*(1 .- (1 .- exp.(mparams[:gPsi]*t))./mparams[:a1] )

    # External forcing
    exog_states[:EF] = mparams[:EF0] .+ 0.01*(mparams[:EF1] .- mparams[:EF0])*min.(t, 100)

    # External forcing in next year
    exog_states[:EF_next] = mparams[:EF0] .+ 0.01*(mparams[:EF1] .- mparams[:EF0])*min.(t .+ 1, 100)

    return exog_states
end

max_stars_repo_name: irudik/optimal-climate-policy-aej

max_issues_repo_name: irudik/optimal-climate-policy-aej

max_forks_repo_name: irudik/optimal-climate-policy-aej


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

@recipe function f(point::Point{Dim,T}) where {Dim,T}
  seriestype --> :scatter
  seriescolor --> :black
  legend --> false

  [Tuple(coordinates(point))]
end

@recipe function f(points::AbstractVector{Point{Dim,T}}) where {Dim,T}
  seriestype --> :scatter
  seriescolor --> :black
  legend --> false

  Tuple.(coordinates.(points))
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "PDAL"
version = v"2.3.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/PDAL/PDAL/releases/download/$version/PDAL-$version-src.tar.gz", "5b0b92258874ef722b5027054d64c8b318b524e7a9b2b250d0330d76e19b8618"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""

cd $WORKSPACE/srcdir/PDAL-*/

if [[ "${target}" == x86_64-linux-musl ]]; then
    # Delete libexpat to prevent it from being picked up by mistake
    rm /usr/lib/libexpat.so*
    # Convince to link `pdal` to libcurl.  Honestly I don't understand why this
    # doesn't work out-of-the-box, like it does for other platforms, this is a
    # quick patch to make it work.
    atomic_patch -p1 ../patches/x86_64-linux-musl-pdal-curl.patch
fi

atomic_patch -p1 ${WORKSPACE}/srcdir/patches/relative_path_dimbuilder.patch
# We'll build `dimbuilder` separately.
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/cmake-disable-dimbuilder.patch
if [[ "${target}" == *-mingw* ]]; then
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/mingw-arbiter.patch
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/mingw-charbuf.patch
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/mingw-cmake-macros.patch
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/mingw-cmake-patches.patch
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/mingw-delete-japanese-unicode-test.patch
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/mingw-dynamiclibrary.patch
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/mingw-fileutils.patch
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/mingw-greedyprojection.patch
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/mingw-optech.patch
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/mingw-pointlayout.patch
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/mingw-portableendian.patch
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/mingw-sbetwriter.patch
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/mingw-sparsematrix.patch
    if [[ "${nbits}" -eq 32 ]]; then
        atomic_patch -p1 ../patches/mingw-trapoline-stdcall.patch
    fi
fi

mkdir -p build/dimbuilder && cd build/dimbuilder

# Build dimbuilder with the host compiler.
(
    # For some reason, CMake seems to ignore the toolchain file, let's force the
    # compiler with the CXX environment variable
    export CXX=${HOST_CXX}
    cmake ../../dimbuilder -G Ninja \
        -DCMAKE_TOOLCHAIN_FILE=${CMAKE_HOST_TOOLCHAIN} \
        -DNLOHMANN_INCLUDE_DIR="$(realpath ../../vendor/nlohmann)" \
        -DCMAKE_BUILD_TYPE=Release
    ninja -j${nproc}
    mkdir -p ../bin
    mv dimbuilder ../bin/.
)

cd ..
cmake .. -G Ninja \
    -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_LIBRARY_PATH:FILEPATH="${libdir}" \
    -DCMAKE_INCLUDE_PATH:FILEPATH="${includedir}" \
    -DCMAKE_BUILD_TYPE=Release \
    -DBUILD_PLUGIN_I3S=OFF \
    -DBUILD_PLUGIN_NITF=OFF \
    -DBUILD_PLUGIN_TILEDB=OFF \
    -DBUILD_PLUGIN_ICEBRIDGE=OFF \
    -DBUILD_PLUGIN_HDF=OFF \
    -DBUILD_PLUGIN_PGPOINTCLOUD=OFF \
    -DBUILD_PLUGIN_E57=OFF \
    -DBUILD_PGPOINTCLOUD_TESTS=OFF \
    -DWITH_LAZPERF=OFF \
    -DBUILD_PGPOINTCLOUD_TESTS=OFF \
    -DWITH_LASZIP=ON \
    -DWITH_ZSTD=ON \
    -DWITH_ZLIB=ON \
    -DWITH_TESTS=OFF

ninja -j${nproc}
ninja install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
#platforms = [
#    Platform("x86_64", "linux"; libc = "glibc")
#]

platforms = expand_cxxstring_abis(supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libpdal_util", :libpdal_util),
    LibraryProduct(["libpdal_base", "libpdalcpp"], :libpdal_base),
    LibraryProduct("libpdal_plugin_kernel_fauxplugin", :libpdal_plugin_kernel_fauxplugin),
    ExecutableProduct("pdal", :pdal),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="GDAL_jll", uuid="a7073274-a066-55f0-b90d-d619367d196c"))
    Dependency(PackageSpec(name="libgeotiff_jll", uuid="06c338fa-64ff-565b-ac2f-249532af990e"))
    Dependency(PackageSpec(name="LASzip_jll", uuid="8372b9c3-1e34-5cc3-bfab-1a98e101de11"))
]

# Build the tarballs, and possibly a `build.jl` as well.
# PDAL GitHub CI scripts currently run on GCC 7.5, so we'll match them in major version at least
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"7.1.0", julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# """
#     readimageheader(io::IO)

# Reads four 32 bit integers at the current position of `io` and
# interprets them as a MNIST-image-file header, which is described
# in detail in the table below

#             ║     First    │  Second  │  Third  │   Fourth
#     ════════╬══════════════╪══════════╪═════════╪════════════
#     offset  ║         0000 │     0004 │    0008 │       0012
#     descr   ║ magic number │ # images │  # rows │  # columns

# These four numbers are returned as a Tuple in the same storage order
# """
function readimageheader(io::IO)
    magic_number = bswap(read(io, UInt32))
    total_items  = bswap(read(io, UInt32))
    nrows = bswap(read(io, UInt32))
    ncols = bswap(read(io, UInt32))
    UInt32(magic_number), Int(total_items), Int(nrows), Int(ncols)
end

# """
#     readimageheader(file::AbstractString)

# Opens and reads the first four 32 bits values of `file` and
# returns them interpreted as an MNIST-image-file header
# """
function readimageheader(file::AbstractString)
    gzopen(readimageheader, file, "r")::Tuple{UInt32,Int,Int,Int}
end

# """
#     readlabelheader(io::IO)

# Reads two 32 bit integers at the current position of `io` and
# interprets them as a MNIST-label-file header, which consists of a
# *magic number* and the *total number of labels* stored in the
# file. These two numbers are returned as a Tuple in the same
# storage order.
# """
function readlabelheader(io::IO)
    magic_number = bswap(read(io, UInt32))
    total_items  = bswap(read(io, UInt32))
    UInt32(magic_number), Int(total_items)
end

# """
#     readlabelheader(file::AbstractString)

# Opens and reads the first two 32 bits values of `file` and
# returns them interpreted as an MNIST-label-file header
# """
function readlabelheader(file::AbstractString)
    gzopen(readlabelheader, file, "r")::Tuple{UInt32,Int}
end


text:
#!/usr/bin/env julia

using Luxor

function point_intersect(fname)
  Drawing(1000, 1000, fname)
  origin()
  setopacity(0.4)

  for i in 1:100
      randomhue()
      p1 = Point(rand(-300:300), rand(-300:300))
      p2 = Point(rand(-300:300), rand(-300:300))
      p3 = Point(rand(-300:300), rand(-300:300))
      p4 = Point(rand(-300:300), rand(-300:300))
      flag, intersection_point = intersection(p1, p2, p3, p4)
      if flag
          circle(intersection_point, 5, :fill)
          move(p1)
          line(p2)
          stroke()
          move(p3)
          line(p4)
          stroke()
      end
  end
  finish()
  println("finished test: output in $(fname)")
end

point_intersect("/tmp/line_intersection.pdf")


text:
function create(c::Client, ::Type{Member}, g::AbstractGuild, u::User; kwargs...)
    return add_guild_member(c, g.id, u.id; kwargs...)
end

function retrieve(c::Client, ::Type{Member}, g::AbstractGuild, u::User)
    return get_guild_member(c, g.id, u.id)
end
function retrieve(c::Client, ::Type{Member}, g::AbstractGuild)
    return list_guild_members(c, g.id)
end

function update(c::Client, m::Member, g::AbstractGuild; kwargs...)
    return modify_guild_member(c, g.id, m.user.id; kwargs...)
end
function update(c::Client, m::Member, r::Role, g::AbstractGuild)
    return add_guild_member_role(c, g.id, m.user.id, r.id)
end

function delete(c::Client, m::Member, g::AbstractGuild)
    return remove_guild_member(c, g.id, m.user.id)
end
function delete(c::Client, m::Member, r::Role, g::AbstractGuild)
    return remove_guild_member_role(c, g.id, m.user.id, r.id)
end


text:
using Base.Test
using DESPOT
using POMDPs

if is_windows()
    error("This test is only valid on Linux and OS X platforms at this time")
end

# Test on a simple RockSample problem
include("../src/problems/RockSample/runRockSample.jl")

include("pomdps_compatibility_tests.jl")
include("test_with_other_particle_filter.jl")
include("test_mersenne.jl")

# Common problem parameters
n_particles             = 500 # for both solver and belief updater
main_seed               = 42
discount                = 0.95
search_depth            = 90
time_per_move           = -1. #default: 1
pruning_constant        = 0.
eta                     = 0.95
sim_len                 = -1
max_trials              = 100 #default: -1
approximate_ubound      = false
debug                   = 0

# Standard batch (3 runs) RockSample(4,4) test
grid_size               = 4
num_rocks               = 4

# The return values below are batch averages
sim_steps, undiscounted_return, discounted_return, run_time =
                main(
                    grid_size = grid_size,
                    num_rocks = num_rocks,
                    n_particles = n_particles,
                    main_seed = main_seed,
                    discount = discount,
                    search_depth = search_depth,
                    time_per_move = time_per_move,
                    pruning_constant = pruning_constant,
                    eta = eta,
                    sim_len = sim_len,
                    max_trials = max_trials,
                    approximate_ubound = approximate_ubound,
                    debug = debug,
                    n_reps = 3
                    )

#println(sim_steps, undiscounted_return, discounted_return, run_time)
if is_linux()
    @test sim_steps == 11
    @test discounted_return ≈ 12.62 atol=1e-2
    @test undiscounted_return == 20.00
end
# osx tests
if is_apple()
    @test sim_steps == 11
    @test discounted_return ≈ 12.97 atol=1e-2
    @test undiscounted_return == 20.00
end

println("DESPOT/RockSample(4,4) batch test status: PASSED")

# Non-standard RockSample(5,6) test
# (non-standard means that the initial state space is generated programmatically)

grid_size               = 5
num_rocks               = 6

sim_steps, undiscounted_return, discounted_return, run_time =
            execute(grid_size = grid_size,
                    num_rocks = num_rocks,
                    n_particles = n_particles,
                    main_seed = main_seed,
                    discount = discount,
                    search_depth = search_depth,
                    time_per_move = time_per_move,
                    pruning_constant = pruning_constant,
                    eta = eta,
                    sim_len = sim_len,
                    max_trials = max_trials,
                    approximate_ubound = approximate_ubound,
                    debug = debug
                    )

println(sim_steps, undiscounted_return, discounted_return, run_time)
if is_linux()
    @test sim_steps == 21
    @test discounted_return ≈ 28.97 atol=1e-2
    @test undiscounted_return == 50.00
end
# osx tests
if is_apple()
    @test sim_steps == 20
    @test discounted_return ≈ 20.80 atol=1e-2
    @test undiscounted_return == 40.00
end

println("DESPOT/RockSample(5,6) test status: PASSED")

# Standard RockSample(7,8) test
grid_size               = 7
num_rocks               = 8

sim_steps, undiscounted_return, discounted_return, run_time =
            execute(
                    grid_size = grid_size,
                    num_rocks = num_rocks,
                    n_particles = n_particles,
                    main_seed = main_seed,
                    discount = discount,
                    search_depth = search_depth,
                    time_per_move = time_per_move,
                    pruning_constant = pruning_constant,
                    eta = eta,
                    sim_len = sim_len,
                    max_trials = max_trials,
                    approximate_ubound = approximate_ubound,
                    debug = debug
                    )

#println(sim_steps, undiscounted_return, discounted_return, run_time)

if is_linux()
    @test sim_steps == 32
    @test discounted_return ≈ 24.82 atol=1e-2
    @test undiscounted_return == 50.00
end
# osx tests
if is_apple()
    @test sim_steps == 30
    @test discounted_return ≈ 16.98 atol=1e-2
    @test undiscounted_return == 40.00
end

println("DESPOT/RockSample(7,8) test status: PASSED")

# Standard RockSample(11,11) test
grid_size               = 11
num_rocks               = 11

sim_steps, undiscounted_return, discounted_return, run_time =
            execute(
                    grid_size = grid_size,
                    num_rocks = num_rocks,
                    n_particles = n_particles,
                    main_seed = main_seed,
                    discount = discount,
                    search_depth = search_depth,
                    time_per_move = time_per_move,
                    pruning_constant = pruning_constant,
                    eta = eta,
                    sim_len = sim_len,
                    max_trials = max_trials,
                    approximate_ubound = approximate_ubound,
                    debug = debug
                    )

#println(sim_steps, undiscounted_return, discounted_return, run_time)
if is_linux()
    @test sim_steps == 50
    @test discounted_return ≈ 21.24 atol=1e-2
    @test undiscounted_return == 60.00
end
# osx tests
if is_apple()
    @test sim_steps == 45
    @test discounted_return ≈ 13.5 atol=1e-2
    @test undiscounted_return == 40.00
end

println("DESPOT/RockSample(11,11) test status: PASSED")


text:
#######################################################################################
# GPU-crutch solve method
# Makes the simplest possible adaptive method for GPU-compatibility
# Out of place only
#######################################################################################
struct GPUSimpleATsit5 end
export GPUSimpleATsit5

function DiffEqBase.solve(prob::ODEProblem,
                          alg::GPUSimpleATsit5;
                          dt = 0.1,
                          abstol = 1e-6, reltol = 1e-3)
  @assert !isinplace(prob)
  u0 = prob.u0
  tspan = prob.tspan
  f = prob.f
  p = prob.p
  ts = Vector{eltype(dt)}(undef,1)
  ts[1] = prob.tspan[1]
  t = tspan[1]
  tf = prob.tspan[2]
  us = Vector{typeof(u0)}(undef,0)
  push!(us,recursivecopy(u0))
  u = u0
  qold = qoldinit
  k7 = f(u, p, t)

  cs, as, btildes, rs = _build_atsit5_caches(eltype(u0))
  c1, c2, c3, c4, c5, c6 = cs
  a21, a31, a32, a41, a42, a43, a51, a52, a53, a54,
  a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = as
  btilde1, btilde2, btilde3, btilde4, btilde5, btilde6, btilde7 = btildes

  # FSAL
  while t < tspan[2]
      uprev = u
      k1 = k7
      EEst = Inf

      while EEst > 1
        dt < 1e-14 && error("dt<dtmin")

        tmp = uprev+dt*a21*k1
        k2 = f(tmp, p, t+c1*dt)
        tmp = uprev+dt*(a31*k1+a32*k2)
        k3 = f(tmp, p, t+c2*dt)
        tmp = uprev+dt*(a41*k1+a42*k2+a43*k3)
        k4 = f(tmp, p, t+c3*dt)
        tmp = uprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4)
        k5 = f(tmp, p, t+c4*dt)
        tmp = uprev+dt*(a61*k1+a62*k2+a63*k3+a64*k4+a65*k5)
        k6 = f(tmp, p, t+dt)
        u = uprev+dt*(a71*k1+a72*k2+a73*k3+a74*k4+a75*k5+a76*k6)
        k7 = f(u, p, t+dt)

        tmp = dt*(btilde1*k1+btilde2*k2+btilde3*k3+btilde4*k4+
                     btilde5*k5+btilde6*k6+btilde7*k7)
        tmp = tmp./(abstol+max.(abs.(uprev),abs.(u))*reltol)
        EEst = defaultnorm(tmp)

        if iszero(EEst)
          q = inv(qmax)
        else
          @fastmath q11 = EEst^beta1
          @fastmath q = q11/(qold^beta2)
        end

        if EEst > 1
          dt = dt/min(inv(qmin),q11/gamma)
        else # EEst <= 1
          @fastmath q = max(inv(qmax),min(inv(qmin),q/gamma))
          qold = max(EEst,qoldinit)
          dtold = dt
          dt = dt/q #dtnew
          dt = min(abs(dt),abs(tf-t-dtold))

          if (tf - t - dtold) < 1e-14
            t = tf
          else
            t += dtold
          end
        end
      end

      push!(us,recursivecopy(u))
      push!(ts,t)
  end
  sol = DiffEqBase.build_solution(prob,alg,ts,us,
                                  calculate_error = false)
  DiffEqBase.has_analytic(prob.f) && DiffEqBase.calculate_solution_errors!(sol;timeseries_errors=true,dense_errors=false)
  sol
end