text:
using IntervalArithmetic, IntervalRootFinding
using ForwardDiff
using Base.Test

const D = IntervalRootFinding.derivative

include("wilkinson.jl")

function generate_wilkinson(n)#, T=BigFloat)   # SLOW

    p = poly(collect(1:n))
    p_prime = polyder(p)

    coeffs = p.a #[convert(T, coeff) for coeff in p.a]
    # coeffs_prime = [convert(T, coeff) for coeff in p_prime.a]

    function f(x)
        total = coeffs[1]
        for i in 2:length(coeffs)
            total += coeffs[i]*x^(i-1)
        end
        total
    end

    f
end

setprecision(Interval, Float64)
float_pi = @interval(pi)

setprecision(Interval, 10000)
big_pi = @interval(pi)
# Using precision "only" 256 leads to overestimation of the true roots for `cos`
# i.e the Newton method gives more accurate results!

half_pi = big_pi / 2
three_halves_pi = 3*big_pi/2

# Format:  (function, derivative, lower_bound, upper_bound, [true_roots])
function_list = [
                    (sin, cos,    -5,  5,    [ -big_pi, @interval(0), big_pi ] ) ,
                    (cos, x->D(cos, x), -7.5, 7.5, [ -three_halves_pi, -half_pi, half_pi, three_halves_pi ] ),
                    (W₃, x->D(W₃, x),   -10, 10,   [ @interval(1), @interval(2), @interval(3) ] ),
                    (W₇, x->D(W₇, x),   -10, 10,   [ @interval(i) for i in 1:7 ] ),
                    (x->exp(x)-2, y->D(x->exp(x),y),  -20, 20,  [log(@biginterval(2))] ),
                    (x->asin(sin(x)) - 0.1, y->D(x->asin(sin(x)),y), 0, 1, [@biginterval(0.1)])
                ]

@testset "Testing root finding" begin

    for rounding_type in (:wide, :narrow)
        @testset "Interval rounding: $rounding_type" begin
            # setrounding(Interval, rounding_type)

            for prec in ( (BigFloat,53), (BigFloat,256), (Float64,64) )
                @testset "Precision: $prec" begin
                    setprecision(Interval, prec)

                    for method in (IntervalRootFinding.newton, IntervalRootFinding.krawczyk)
                        @testset "Method $method" begin

                            for func in function_list

                                f, f_prime, a_lower, a_upper, true_roots = func
                                a = @interval(a_lower, a_upper)

                                @testset "Function $f; interval $a" begin

                                    for autodiff in (false, true)

                                        @testset "With autodiff=$autodiff" begin

                                            if autodiff
                                                rts = find_roots(f, a, method)
                                            else
                                                rts = find_roots(f, f_prime, a, method)
                                            end

                                            @test length(rts) == length(true_roots)

                                            for i in 1:length(rts)
                                                root = rts[i]

                                                @test isa(root, Root)
                                                @test is_unique(root)
                                                @test true_roots[i] ⊆ root.interval
                                            end
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

setprecision(Interval, Float64)

@testset "find_roots tests" begin
    f(x) = x^2 - 2

    rts = IntervalRootFinding.newton(f, @interval(10, 11))
    @test length(rts) == 0

    rts = IntervalRootFinding.find_roots_midpoint(f, -5, 5)
    @test length(rts) == 3
    @test length(rts[1]) == 2

    rts = find_roots(f, -5, 5)
    @test length(rts) == 2

    setprecision(Interval, 256)

    for method in (IntervalRootFinding.newton, IntervalRootFinding.krawczyk)
        new_roots = method(f, rts)
        @test length(new_roots) == length(rts)

        for (root1, root2) in zip(new_roots, rts)
            @test root1 ⊆ root2
        end
    end
end

@testset "Multiple roots" begin
    setprecision(Interval, Float64)
    let
        f(x) = (x-1) * (x^2 - 2)^3 * (x^3 - 2)^4

        rts = IntervalRootFinding.newton(f, -5..5.1, maxlevel=1000)

        @test length(rts) == 4
        @test rts[1].status == :unknown
        @test rts[1].interval == Interval(-1.4142135623730954, -1.414213562373095)
        @test rts[3].interval == Interval(1.259921049894873, 1.2599210498948734)

    end
end

@testset "Iterated logistic function fixed points" begin
    ∘(f, g) = x -> f(g(x))
    iterate(f, n) = n == 1 ? f : f ∘ iterate(f, n-1)

    f(x) = 4x*(1-x)  # logistic map

    for n in 1:10
        g = x -> iterate(f, n)(x) - x  # look for fixed points of f^n

        rts = IntervalRootFinding.newton(g, 0..1)
        @test length(rts) == 2^n

        rts = IntervalRootFinding.krawczyk(g, 0..1)
        @test length(rts) == 2^n

    end
end

# Example of a function with a double root at 0 from Burden & Faires, 9th ed, p.84
# exp(x) - x - 1


text:
# # Time series forecasting
# This tutorial is inspired by a [Tensorflow](https://www.tensorflow.org/tutorials/structured_data/time_series) Time Series Forecasting Tutorial.
# It implements the data analysis, plotting, and model training from scratch using Julia tools such as Plots.jl, Flux, and DataFrames.jl.

# ## Setup
# Importing the usual suspects. Two interesting notes:
# * StatsPlots re-exports Plots.jl functions, so we don't need that separately
# * MLDataPattern is useful for preprocessing data for ML in a lazy (read: non-memory intensive) way. We will be seeing more of it
using ZipFile
using CSV
using DataFrames
using StatsPlots
using Dates
using FFTW
using CUDA
using Flux
using MLDataPattern

using Random: seed!
using Statistics: mean, std
using Flux: unsqueeze

# Set a seed to make this reproducible.
seed!(4231)

ENV["LINES"] = 20;

# Disallow scalar indexing for CUDA - prevents slow accesses to GPU memory
CUDA.allowscalar(false)

# ## The weather dataset
# We are going to be using weather data for 7 years from the Max Planck Institute: https://www.bgc-jena.mpg.de/wetter/

function download_data(; fname="jena_climate_2009_2016.zip")
    DATA_PATH = "https://storage.googleapis.com/tensorflow/tf-keras-datasets/jena_climate_2009_2016.csv.zip"
    isfile(fname) || download(DATA_PATH, fname)

    zip = ZipFile.Reader(fname)
    csv = zip.files[1]
    df = CSV.File(csv) |> DataFrame
    close(zip)
    return df
end

# More discussion on reading zipped CSVs here: https://juliadata.github.io/CSV.jl/stable/index.html#Example:-reading-from-a-zip-file-1

df = download_data();

df = df[6:6:end, :]

df[!,"Date Time"] = Dates.DateTime.(df[:,"Date Time"], "dd.mm.yyyy HH:MM:SS");
# More info on working with dates and times: https://en.wikibooks.org/wiki/Introducing_Julia/Working_with_dates_and_times

col = ["Date Time", "T (degC)", "p (mbar)", "rho (g/m**3)"]

df[:,col]

# One easy way to plot data from a datafrom is with the StatsPlots `@df` macro
@df df plot(cols(1), cols(2:4); layout=(3, 1))

@df df[1:480,:] plot(cols(1), cols(2:4); layout=(3, 1))

# ## Inspect and cleanup
@show describe(df)

df."p (mbar)"; #cool that this works

# Need to replace bad database values
replace!(df[!,"wv (m/s)"], -9999.0 => 0);
replace!(df[!,"max. wv (m/s)"], -9999.0 => 0);

@show describe(df)

# ## Feature engineering
# Bucket wind direction and speed
histogram2d(df[!,"wd (deg)"], df[!,"wv (m/s)"], bins=(75,75), xlabel="wd (deg)", ylabel="wv (m/s)")

# It would be better to decompose this into X and Y velocities
wd_rad = df[!,"wd (deg)"] * π / 180
df.Wx = df[!,"wv (m/s)"] .* cos.(wd_rad)
df.Wy = df[!,"wv (m/s)"] .* sin.(wd_rad)
df."max Wx" = df[!,"max. wv (m/s)"] .* cos.(wd_rad)
df."max Wy" = df[!,"max. wv (m/s)"] .* sin.(wd_rad);

histogram2d(df.Wx, df.Wy, bins=(75,75), xlabel="Wind X [m/s]", ylabel="Wind Y [m/s]")

# Modify timestamp into a continuous "time of day" signal to deal with periodicity.
timestamp_s = Dates.datetime2unix.(df."Date Time");

day = 24*60*60 #seconds in a day
year = 365.2425 * day #seconds in a year

df[!,"Day sin"] = sin.(timestamp_s * (2 * π / day))
df[!,"Day cos"] = cos.(timestamp_s * (2 * π / day))
df[!,"Year sin"] = sin.(timestamp_s * (2 * π / year))
df[!,"Year cos"] = cos.(timestamp_s * (2 * π / year));

plot(df[1:25,"Day sin"], legend=false)
plot!(df[1:25,"Day cos"])
xlabel!("Time [h]")
title!("Time of Day Signal")

# If you don't know the frequency to choose for your periodicity, you can take the FFT.
# Here you can see that the yearly and daily periodicities are the most prominent.
fftrans = FFTW.rfft(df[!,"T (degC)"])
f_per_dataset = 1:size(fftrans)[1]

n_samples_h = size(df[!,"T (degC)"])[1]
hours_per_year = 24 * 365.2524
years_per_dataset = n_samples_h / hours_per_year
f_per_year = f_per_dataset / years_per_dataset;

plot(f_per_year, abs.(fftrans), xscale=:log10, ylim=(0, 400000), xlim=(0.3,Inf), leg=false)
xticks!([1, 365.2524], ["1/Year", "1/Day"])
xlabel!("Frequency (log scale)")

# ## Split the data
# Drop columns that won't be used further.
select!(df, Not([:("wv (m/s)"),:("max. wv (m/s)"), :("wd (deg)"), :("Date Time")]));

column_indices = pairs(names(df))
indices_columns = Dict(value => key for (key, value) in column_indices)
df = convert.(Float32, df) # Don't need high precision; reduces errors later on when using Params - gradients are Float32

# Split the data into a training and validation set
n = size(df)[1]
train_df = df[1:round(Int,n*0.7, RoundDown),:]
valid_df = df[round(Int,n*0.7, RoundUp):round(Int,n*0.9, RoundDown),:]
test_df = df[round(Int,n*0.9, RoundUp):end,:]; # matching TF tutorial exactly, can also use partition

num_features = size(df,2)

# ## Normalize the data
train_mean = mean.(eachcol(train_df))
train_std = std.(eachcol(train_df))

train_df = (train_df .- train_mean') ./ train_std'
valid_df = (valid_df .- train_mean') ./ train_std'
test_df = (test_df .- train_mean') ./ train_std'

df_std = (df .- train_mean') ./ train_std'
df_std = stack(df_std)

violin(df_std.variable, df_std.value, xrotation=30.0, legend=false, xticks=:all) # use plotattr() to learn about keywords
boxplot!(df_std.variable, df_std.value, fillalpha=0.75, outliers=false)

# ## Data Windowing
# We will define our own WindowGenerator, some constructors, and plotting functions. The data from the WindowGenerator will be used in training.
include("window_generator.jl")

h = 6; # historical window length
f = 1; # future window length

# WindowGenerator makes use of MLDataPattern's `slidingwindow` to generate the windows. It is good at flexibly generating sequences data, though requires a bit of mind-bending to fully understand how it generates sequences.
# [slidingwindow docs](https://mldatapatternjl.readthedocs.io/en/latest/documentation/dataview.html?highlight=slidingwindow#labeled-windows)
slidingwindow(i -> i+h:i+h+f-1, Array(train_df)', h, stride=1)

# WindowGenerator is implemented to make sequence-generating simpler. It uses `slidingwindow` behind the scenes to create sequences.
@doc WindowGenerator

# Create a WindowGenerator
wg = WindowGenerator(6, 1, train_df, valid_df, label_columns="T (degC)");

# Plotting methods are implemented for raw data and with predictions.
plot(wg)

# We will also make use of the utility function for batching lazily-evaluated timeseries from slidingwindow
include("batch_ts.jl")

# To understand better what all the above code does, let's look at an imitation training loop and look at everything's dimensions
practice_df = train_df[1:10,:]
a = slidingwindow(i -> i+h:i+h+f-1, Array(practice_df)', h, stride=1)

for i in eachbatch(shuffleobs(a), size=2)
    (x,y) = batch_ts(i)
    @show x
    @show y
    println()
end
#=
That took the data at 10 timestamps (`train_df[1:10,:]`) and threw it into 4 'observations' using all the data it could while requesting 6 historical and 1 target point. This works out to the following:
* (1,2,3,4,5,6), 7
* (2,3,4,5,6,7), 8
* (3,4,5,6,7,8), 9
* (4,5,6,7,8,9), 10

It then shuffles these obserations, and lazily batches them into 2 batches each of size 2.
=#

# # Single Step Models
loss(x,y) = Flux.Losses.mse(x, y)

# ### Baseline - 1h
struct Baseline
    label_index::Int
end

(m::Baseline)(x) = x[m.label_index,:,:]

target = "T (degC)"

# Since this model repeats the last point, we make slidingwindows with 1 historical and 1 target point.
single_step_1h = WindowGenerator(1, 1, train_df, valid_df, label_columns=target);

baseline_model = Baseline(wg.target_idx[1])

# Demonstrate the 'training error' for this model, although there is no actual training.
function run_single_step_baseline(wg, model)
    preds = Float32[]
    reals = Float32[]
    for (x,y) in wg.train
        val = model(x)[1]
        push!(preds, val)
        push!(reals, y[model.label_index])
    end

    l = loss(preds, reals)
    return l
end

run_single_step_baseline(single_step_1h, baseline_model)

# ### Baseline - 24h
# Let's try to predict the next hour's value for 24 hours
single_step_24h = WindowGenerator(24, 24, train_df, valid_df, label_columns=target; offset=1);

plot(single_step_24h, baseline_model)

# ### Linear Models
# ##### 1 hour
# Set up a fairly generic training loop for the rest of the models. It takes in a model, and a WindowGenerator object, and uses the ADAM optimizer to train the model.
opt = Flux.Optimise.ADAM(0.01)

function train_model!(model, wg::WindowGenerator, opt; epochs=20, bs=16, dev=Flux.gpu)
    model = model |> dev
    ps = params(model)
    t = shuffleobs(wg.train)
    v = batch_ts(getobs(wg.valid))

    local l
    vl_prev = Inf
    for e in 1:epochs
        for d in eachbatch(t, size=bs)
            x, y = batch_ts(d)
            y = y[wg.target_idx,:,:]
            x, y = x |> dev, y |> dev
            gs = gradient(ps) do
                l = loss(model(x),y)
            end
            Flux.update!(opt, ps, gs)
        end
        l = round(l;digits=4)
        vl = round(loss(model(v[1]),v[2][wg.target_idx,:,:]); digits=4)
        println("Epoch $e/$epochs - train loss: $l, valid loss: $vl")
        #=
        crude early-stopping
        vl_prev < (vl - 0.001) && break
        vl_prev = vl
        =#
    end
    model = model |> cpu
end

# Now, train the linear model
linear = Dense(size(single_step_1h.train[1][1],1), 1; initW=Flux.glorot_uniform, initb=Flux.zeros)

@time linear = train_model!(linear, single_step_1h, opt; bs=32, epochs=20)

# ##### 24 hours
# Show the model's predictions for each point in the window.
plot(single_step_24h, linear)

# Visualize the weights of the linear regression. Naturally, since we are trying to predict future temperatures, we would expect that the current temperature is on of the most useful features.
# This is shown by it having a high positive weight.
bar(names(train_df), linear.W[:], xrotation=30.0, legend=false, xticks=:all, tickfontsize=6)

# ### Dense
# A model with 2 Fully-connected layers and relu non-linearities.
dense = Chain(
    Dense(19, 64, relu),
    Dense(64, 64, relu),
    Dense(64, 1)
)

@time dense = train_model!(dense, single_step_1h, opt; bs=32, epochs=20)

# Dense model's predicions
plot(single_step_24h, dense)

# ### Multi-step Dense
# Now we are going to use 3 historical hours to predict 1 hour in the future.
single_step_3h = WindowGenerator(3, 1, train_df, valid_df, label_columns=target);

plot(single_step_3h)

multi_step_dense = Chain(
    i -> reshape(i, :, 1, size(i)[end]), # flatten first two dimensions, but preserve batch dimension
    Dense(19*3, 32, relu),
    Dense(32, 32, relu),
    Dense(32, 1)
)

@time multi_step_dense = train_model!(multi_step_dense, single_step_3h, opt; bs=32, epochs=20)

# Show the model's 1 hour predicion using 3 hours of historical data.
plot(single_step_3h, multi_step_dense)

# ### Convolutional Neural Network

# conv_model = Chain(
#     Conv((19,), 3=>32, relu), # need to explain why this conv pattern
#     x -> Flux.flatten(x),
#     Dense(32, 32, relu),
#     Dense(32, 1),
#     x -> unsqueeze(x, 1)
# )

conv_model = Chain(
    x -> permutedims(x, [2,1,3]), # put data in NTime x NCovariates X NBatch https://github.com/FluxML/Flux.jl/issues/1465
    Conv((3,), 19=>32, relu), # convolve over 3 inputs - 19 variables -> 32 filters
    x -> Flux.flatten(x),
    Dense(32, 32, relu),
    Dense(32, 1),
    x -> unsqueeze(x, 1)
)

single_step_3h = WindowGenerator(3, 1, train_df, valid_df, label_columns=target);

@time conv_model = train_model!(conv_model, single_step_3h, opt; bs=32, epochs=20)

# Show the model's 1 hour predicion using 3 hours of historical data.
plot(single_step_3h, conv_model)

# ### Recurrent Neural Network

rnn_single_24h = Chain(
    LSTM(19, 32), # input features * output dimensionality
    Dense(32, 1),
    x->reshape(x,:) # lstm output * single number
)

# Flux uses vectors of observations for RNNs, rather than the bulk tensor
a = [rand(Float32, 19) for i in 1:24]
b = rnn_single_24h.(a)
a[1]
rnn_single_24h(a[1])

# ### Performance

# ### Multi-output Models

# #### Baseline

# #### Dense

# #### RNN

# #### Advanced: Residual Connections

# #### Performance

# # Multi-Step Models

# ### Baselines

# ## Single-shot Models

# ### Linear

# ### Dense

# ### CNN

# ### RNN

# ## Advanced Autoregressive model
# ### RNN

# ## Performance

# Next Steps

ext: jl

max_stars_repo_name: adinhobl/model-zoo

max_issues_repo_name: adinhobl/model-zoo

max_forks_repo_name: adinhobl/model-zoo


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

using Markdown
using InteractiveUtils

# ╔═╡ 6be49f08-6894-11eb-307f-6d65bf5f9130
using Unitful, UnitfulAtomic

# ╔═╡ fa72a230-6612-11eb-2c8c-41f072db62c8
using DifferentialEquations

# ╔═╡ 19d62798-6613-11eb-3432-fd10668c8447
using StaticArrays

# ╔═╡ dc2f69a6-6613-11eb-39a0-456ca3cfcafe
using Plots; gr()

# ╔═╡ 5f050fe8-6614-11eb-38e1-3b878533a01c
using LaTeXStrings

# ╔═╡ ab527e9c-66a2-11eb-3531-8b666bb82451
using DataFrames

# ╔═╡ 3c0653ec-67ba-11eb-1bce-e5ec4c92989e
using UnitfulAstro

# ╔═╡ c7ffa252-660e-11eb-381b-2f2ceec373c1
md"## white dwarfs"

# ╔═╡ de5358ce-67c0-11eb-1c87-49efb42f32d7
md"### Course on Compact Objects: Assignment 4

Instructor: Prof. P. Ajith
Author: Md Arif Shaikh, Postdoc, ICTS-TIFR
Date: Feb 5, 2021"

# ╔═╡ 73e9d7ec-6793-11eb-3bca-bd8d40ef03d6
md"### Problem 1

Show that the pressure exerted by a gas of particle with isotropic momentum distribution $n(p)$ is given by

$P = \frac{1}{3}\int_0^\infty p v_p n(p) dp$

where $v_p$ is the velocity associated with momentum $p$.
"

# ╔═╡ b71be1a4-6793-11eb-0012-dba495d16df2
md"Let us consider the pressure due to a particle moving at an angle $\theta$ with the $z$ direction. The component of velocity along $z$ would be given by $v \cos\theta$

$v_z = v \cos\theta$

and therefore,

$p_z = p \cos\theta.$

The rate of change of momentum is given by

$\frac{p_z}{\Delta z/v_z} = \frac{p_z v_z}{\Delta z} = \frac{p v \cos^2\theta}{\Delta z}$

So pressure per particle moving at an angle $\theta$ is

$\frac{p v \cos^2\theta}{\Delta z \Delta x \Delta y} = \frac{p v \cos^2\theta}{V}.$

The total pressure for all particles moving at an angle $\theta$ would be given by

$P = \int_0^\infty \frac{p v \cos^2\theta}{V} N(p, \theta) dp$

where $N(p, \theta)$ is the number of particles with momentum between $p$ and $p+dp$ and with angle $\theta$ to $\theta + d\theta$. This number would be

$N(p,\theta) = \frac{d^3p}{h^3}g_s f(p) = g_s f(p)\frac{2\pi p^2\sin\theta d\theta dp}{h^3},$

where $g_s$ is the degeracy of the state and $f(p)$ is the probabilty that the state would be occupied.

Thus,

$P = \int_0^\infty \frac{p v \cos^2\theta}{V} g_s f(p)\frac{2\pi p^2 \sin\theta d\theta dp}{h^3} dp.$

To get the average total pressure we integrate over $\theta$

$P = \int_0^\pi d\theta \cos^2\theta \sin\theta \int_0^\infty \frac{p v }{V} \frac{g_s f(p) 2\pi p^2}{h^3} dp = \frac{1}{3} \int_0^\infty \frac{p v }{V} \frac{g_s f(p) 4\pi p^2}{h^3} dp = \frac{1}{3}\int_0^\infty p v n(p) dp$

where $n(p) = g_s f(p) 4\pi p^2/V.$
"

# ╔═╡ 6dab1184-6798-11eb-01e2-af1bc7dbdd1d
md"### Problem 2

Argue why we are justified in using a 'cold' degenerate equation of state to describe a white dwarf with a temperature
$T \sim 10^4$K (Hint: Show that the degeneracy parameter $\mu/kT \gg 0$, where $\mu$ is the chemical potential and $k$ the Boltzmann
constant. The density of the white dwarf is $\sim 10^6 g/cm^3$ and the chemical potential $\sim$ the Fermi energy. Assume that
$\mu_e = 2$)."

# ╔═╡ b8e416ba-6891-11eb-2c92-cdc819988514
md"We have to basically show that the thermal energy of the electron gas is negligible compared to Fermi energy"

# ╔═╡ 33138330-6892-11eb-36fa-fdb893d1a245
kB = u"k_au"

# ╔═╡ bd2e2f6e-6894-11eb-1246-db600b48a079
Temp = 1e4 * u"K"

# ╔═╡ 244ed64c-6892-11eb-3f02-4f264c247878
ThermalEnergy = uconvert(u"erg", (3/2) * kB * Temp)

# ╔═╡ ae72f48a-6896-11eb-000d-a75c99a66d0c
h = 2 * pi * u"ħ_au"

# ╔═╡ 9af314b6-6897-11eb-1016-9bcbf7e7a28d
μe = 2

# ╔═╡ c2d91a0c-6897-11eb-12d9-05df32fef1d9
amu = 1.66054e-24 * u"g"

# ╔═╡ b9c69186-6896-11eb-3118-3b8a5ed36b65
ne = 1e6 * u"g/cm^3" / (μe * amu)

# ╔═╡ e29a39f2-6897-11eb-0515-db38cc2614ed
me = 9.109 * 1e-31 * u"kg"

# ╔═╡ 019fc0f6-6898-11eb-3936-fbfca4d9f508
c = u"c"

# ╔═╡ f6d01dbc-6895-11eb-2ded-e536398b6f60
function FermiEnergy(ne, me, h)
	pF = h * (3 * ne / (8 * pi))^(1/3)
	sqrt(pF^2 * c^2 + me^2 * c^4)
end

# ╔═╡ d001849e-6897-11eb-3267-f7fa03a09e80
Fenergy = uconvert(u"erg", FermiEnergy(ne, me, h))

# ╔═╡ 4a6ca9b8-6899-11eb-390a-8d15b0fe2912
Fenergy/ThermalEnergy

# ╔═╡ 93d8e972-6799-11eb-0149-ed6bd5c87857
md"### Problem 3

Derive Lane-Emden equation"

# ╔═╡ b98c5b72-6611-11eb-3001-71f3dddb7f23
md"### Lane-Emden equation"

# ╔═╡ e1480d1c-660e-11eb-1c29-5decbd0775cd
md"$\frac{dm(r)}{dr} = 4\pi r^2\rho(r), \frac{dp(r)}{dr}=-\frac{Gm(r)\rho(r)}{r^2}$. These two equations can be combined to give

$\frac{1}{r^2}\frac{d}{dr}\left(\frac{r^2}{\rho(r)}\frac{dp(r)}{dr}\right) = -4\pi G \rho(r)$

This can be further written in terms of a dimensionless form using

$\rho = \rho_c \theta^n, r = a\xi, \Gamma = 1 + \frac{1}{n}$

Then we get the following equation

$\frac{1}{\xi^2}\frac{d}{d\xi}\left(\xi^2 \frac{d\theta}{d\xi}\right) = -\theta^n$

This the Lane-Emden equation.
"

# ╔═╡ add3d718-6799-11eb-0a12-398c88b47440
md"### Problem 4"

# ╔═╡ ce805c54-6611-11eb-1c6b-ff532c7973b6
md"### Solve Lane-Emden Equation"

# ╔═╡ dad519b6-6611-11eb-20f6-df1d57ee4fcb
md"To solve the lane-emden equation we write it as two first order odes by defining $\psi = \frac{d\theta}{d\xi}$

$\frac{d\theta}{d\xi} = \psi, \frac{d\psi}{d\xi} = - \frac{2\psi}{\xi} - \theta^n$

With the following initial values

$\theta(0) = 1, \psi(0) = 0$"

# ╔═╡ 9ca0c614-6612-11eb-2dce-1b8e502da696
function laneEmden(u, n, ξ)
	θ, ψ = u
	dθ = ψ
	dψ = - (2 * ψ / ξ) - θ^n
	@SVector [dθ, dψ]
end

# ╔═╡ 38699210-6613-11eb-2fae-1502e73646d3
u0 = @SVector [1., 0.]

# ╔═╡ 59092ad0-6613-11eb-3cca-0be9b3b034f6
ξspan1 = (1e-10, 10)

# ╔═╡ 71d03ac2-6613-11eb-31a9-376c57d17fbe
problem1 = ODEProblem(laneEmden, u0, ξspan1, 1.)

# ╔═╡ b568cdc6-6613-11eb-3bf6-0f0fce82c1d2
sol1 = solve(problem1)

# ╔═╡ 033df3a0-6614-11eb-17e8-d175880db41d
plot(sol1, vars = (0, 1), xlabel = L"\xi", ylabel = L"\theta", label = L"n = 1")

# ╔═╡ b386900e-6629-11eb-1a5e-7b5c82483b88
md"Solutions for different integer $n$"

# ╔═╡ f0587620-662a-11eb-1f28-5f6a0938d300

# ╔═╡ 3eb68340-662b-11eb-1465-1b4e43a9750d
condition(u, ξ, integrator) =  u[1];

# ╔═╡ c22a1caa-662b-11eb-2fe0-99791c25951d
affect!(integrator) = terminate!(integrator);

# ╔═╡ cdfd4ed8-662b-11eb-28e3-21f8f05045e8
cb = ContinuousCallback(condition, affect!);

# ╔═╡ eb305878-6629-11eb-224c-0775dbfbbe01
ξspan = (1e-10, 20)

# ╔═╡ c2483324-6629-11eb-096f-892e2660ddbe
begin
	p = plot();
	for n in 1:4
		problem = ODEProblem(laneEmden, u0, ξspan, n)
		sol = solve(problem, callback = cb)
		plot!(p, sol, vars = (0, 1), lw = 2, label = L"n = %$n")
	end
	plot!(xlabel = L"\xi", ylabel = L"\theta", ylim = (0, 1))
	p
end

# ╔═╡ 75295168-6640-11eb-2f4d-79dbafd7e74c
md"for fractional $n$ the solver would throw error as negative value of $\theta$ would give complex value for $\theta^n$. To avoid this problem we can just change $\theta$ in our set of odes to abs($\theta$)."

# ╔═╡ e245dec6-7002-11eb-339f-c590a8afc985
function laneEmden2(u, n, ξ)
	θ, ψ = u
	dθ = ψ
	dψ = - (2 * ψ / ξ) - abs(θ)^n
	@SVector [dθ, dψ]
end

# ╔═╡ c57492b8-6640-11eb-1f31-29cc2440c6e0
begin
	pfrac = plot();
	ns = 0.5:1.0:4.5
	ξspan2 = [1.0e-10, 30.0]
	ξ1s = []
	ψs = []
	for n in ns
		problem = ODEProblem(laneEmden2, u0, ξspan2, n)
		sol = solve(problem, callback = cb)
		push!(ξ1s, sol.t[end])
		push!(ψs, sol[2, :][end])
		plot!(pfrac, sol.t, sol[1, :], lw = 2, xlabel = L"\xi", ylabel = L"\theta", label = L"n = %$n")
	end
end

# ╔═╡ 1d4b843c-66a1-11eb-188c-871b357bfdd7
ξ1s

# ╔═╡ 72c6e67c-66a1-11eb-2171-852cad67d86f
pfrac

# ╔═╡ 913a247c-66b8-11eb-296a-b92827dff6c8
data = DataFrame(n = ns, ξ1 = ξ1s, ψ1 = ψs)

# ╔═╡ 774b990a-66b9-11eb-2f71-ef6e4b08ec6f
md"### Problem 5

Radius and Mass of white dwarf"

# ╔═╡ 27a64b22-679d-11eb-0763-6d5931f2d0bb
md"Mass of the star is given by

$M = \int_0^\infty 4\pi\rho r^2 dr$

now,

$r = \xi  \left(\frac{4 \pi  G {\rho_c}^{1-\frac{1}{n}}}{K (n+1)}\right)^{-\frac{1}{2}}$

$\rho=\rho_c \theta^n$

$\theta^n = - \frac{1}{\xi^2}\frac{d}{d\xi}\left(\xi^2 \frac{d\theta}{d\xi}\right)$

Thus,

$M = 4\pi \rho_c\left(\frac{4 \pi  G {\rho_c}^{1-\frac{1}{n}}}{K (n+1)}\right)^{-3/2}\int_0^{\xi_1} d\left(\xi^2 \frac{d\theta}{d\xi}\right) = 4\pi \rho_c\left(\frac{4 \pi  G {\rho_c}^{1-\frac{1}{n}}}{K (n+1)}\right)^{-3/2}\left(\xi^2 \frac{d\theta}{d\xi}\right)_{\xi_1}$

or

$M = 4\pi \rho_c^{(3-n)/2n}\left(\frac{4 \pi G}{K (n+1)}\right)^{-3/2}\left(\xi^2 \frac{d\theta}{d\xi}\right)_{\xi_1}$

Now, the radius would be given by value of $r$ at $\xi = \xi_1$ which is

$R_\star = \xi_1  \left(\frac{4 \pi  G {\rho_c}^{1-\frac{1}{n}}}{K (n+1)}\right)^{-\frac{1}{2}}$

From this we can express $\rho_c$ in terms of $R_\star$ and $\xi_1$ as

$\rho_c = \left(\frac{K (n+1) \xi_1 ^2}{4\pi G R_\star^2}\right)^{\frac{n}{n-1}}$

So,

$\boxed{M = 4\pi R_\star^{(3-n)/(1-n)} \left(\frac{K (n+1)}{4\pi G}\right)^{n/(n-1)} \xi_1^{-(3-n)/(1-n)} \xi_1^2\left.\frac{d\theta}{d\xi}\right|_{\xi_1}}$
"

# ╔═╡ d299d78a-67b7-11eb-2bf5-192e423764d7
md"### Problem 6

Compute mass and radius of white dwarfs
"

# ╔═╡ ea8073d8-67b7-11eb-3a59-99769f0574f9
md"From problem 4, we get that for $n=3/2$, $\xi_1 \approx 3.655$. Now given that

$\rho_c = 10^6 - 10^9 g/cm^3, K \approx 10^{13} \mu_e^{-5/3}, \mu_e = 2.$"

# ╔═╡ a1b993b2-66bf-11eb-3df9-ddd3a3363f1b
function MassWhiteDwarf(R, K, G, ξ1, ψ1, n)
	4 * pi * R^((3. -n)/(1. - n))*((K*(n+1))/(4*pi*G))^(n/(n-1))*ξ1^(-(3. - n)/(1. -n))*ξ1^2 * abs(ψ1)
end

# ╔═╡ bb649190-67b9-11eb-012d-2bc5d384ddba
function Rstar(G, ξ1, ρc, K, n)
	return ξ1 * (4 * pi * G * ρc^(1.0 - (1.0 / n))/(K * (n + 1)))^(-1.0 / 2.0)
end

# ╔═╡ 056ca994-67ba-11eb-2459-5bf59f7db5c6
ξ1 = data.ξ1[2]

# ╔═╡ 2f4ab814-67ba-11eb-34e6-db97da22b6ad
ψ1 = data.ψ1[2]

# ╔═╡ 47d0d1c0-67ba-11eb-2d35-f32177efc1ab
G = u"GMsun"/u"Msun"

# ╔═╡ 6b1ef04e-67ba-11eb-3c96-d15015059f76
μ = 2.0

# ╔═╡ dbb73808-67be-11eb-12d0-916378cf0c7f
n = 3/2

# ╔═╡ 58b6b5ea-67ba-11eb-3bca-1d89747d69bd
K = 1e13 * μ^(-5.0/3.0) * u"dyn/cm^2"/(u"g/cm^3")^(1+(1/n))

# ╔═╡ 7205031c-67ba-11eb-1e45-75ac2b20c76f
ρcs = [1e6, 1e7, 1e8, 1e9] * u"g/cm^3"

# ╔═╡ 7233ae8c-67bb-11eb-13dc-716bb11add02
ρcs[2]

# ╔═╡ be271b86-67ba-11eb-2828-517b586bc4ad
r1 = uconvert(u"km", Rstar(G, ξ1, ρcs[1], K, n))

# ╔═╡ 3e871518-67bd-11eb-100b-9dce3872fa40
mass1 = MassWhiteDwarf(r1, K, G, ξ1, ψ1, n)

# ╔═╡ 93838510-67bd-11eb-1ea0-7fcb8db20c6f
uconvert(u"Msun", mass1)

# ╔═╡ 171b5cca-67c0-11eb-3c51-c525076f9966
begin
	masses = []
	radii = []

	for ρc in ρcs
		radius = uconvert(u"km", Rstar(G, ξ1, ρc, K, n))
		mass = uconvert(u"Msun", MassWhiteDwarf(radius, K, G, ξ1, ψ1, n))
		push!(masses, mass)
		push!(radii, radius)
	end
end

# ╔═╡ 7299b18c-67c0-11eb-1feb-2533c617d02d
massRadius = DataFrame(central_density = ρcs, radius = radii, mass = masses)

# ╔═╡ Cell order:
# ╟─c7ffa252-660e-11eb-381b-2f2ceec373c1
# ╟─de5358ce-67c0-11eb-1c87-49efb42f32d7
# ╟─73e9d7ec-6793-11eb-3bca-bd8d40ef03d6
# ╟─b71be1a4-6793-11eb-0012-dba495d16df2
# ╟─6dab1184-6798-11eb-01e2-af1bc7dbdd1d
# ╟─b8e416ba-6891-11eb-2c92-cdc819988514
# ╠═6be49f08-6894-11eb-307f-6d65bf5f9130
# ╠═33138330-6892-11eb-36fa-fdb893d1a245
# ╠═bd2e2f6e-6894-11eb-1246-db600b48a079
# ╠═244ed64c-6892-11eb-3f02-4f264c247878
# ╠═f6d01dbc-6895-11eb-2ded-e536398b6f60
# ╠═ae72f48a-6896-11eb-000d-a75c99a66d0c
# ╠═9af314b6-6897-11eb-1016-9bcbf7e7a28d
# ╠═c2d91a0c-6897-11eb-12d9-05df32fef1d9
# ╠═b9c69186-6896-11eb-3118-3b8a5ed36b65
# ╠═e29a39f2-6897-11eb-0515-db38cc2614ed
# ╠═019fc0f6-6898-11eb-3936-fbfca4d9f508
# ╠═d001849e-6897-11eb-3267-f7fa03a09e80
# ╠═4a6ca9b8-6899-11eb-390a-8d15b0fe2912
# ╟─93d8e972-6799-11eb-0149-ed6bd5c87857
# ╟─b98c5b72-6611-11eb-3001-71f3dddb7f23
# ╟─e1480d1c-660e-11eb-1c29-5decbd0775cd
# ╟─add3d718-6799-11eb-0a12-398c88b47440
# ╟─ce805c54-6611-11eb-1c6b-ff532c7973b6
# ╟─dad519b6-6611-11eb-20f6-df1d57ee4fcb
# ╠═fa72a230-6612-11eb-2c8c-41f072db62c8
# ╠═19d62798-6613-11eb-3432-fd10668c8447
# ╠═9ca0c614-6612-11eb-2dce-1b8e502da696
# ╠═38699210-6613-11eb-2fae-1502e73646d3
# ╠═59092ad0-6613-11eb-3cca-0be9b3b034f6
# ╠═71d03ac2-6613-11eb-31a9-376c57d17fbe
# ╠═b568cdc6-6613-11eb-3bf6-0f0fce82c1d2
# ╠═dc2f69a6-6613-11eb-39a0-456ca3cfcafe
# ╠═5f050fe8-6614-11eb-38e1-3b878533a01c
# ╠═033df3a0-6614-11eb-17e8-d175880db41d
# ╟─b386900e-6629-11eb-1a5e-7b5c82483b88
# ╠═f0587620-662a-11eb-1f28-5f6a0938d300
# ╠═3eb68340-662b-11eb-1465-1b4e43a9750d
# ╠═c22a1caa-662b-11eb-2fe0-99791c25951d
# ╠═cdfd4ed8-662b-11eb-28e3-21f8f05045e8
# ╠═eb305878-6629-11eb-224c-0775dbfbbe01
# ╠═c2483324-6629-11eb-096f-892e2660ddbe
# ╟─75295168-6640-11eb-2f4d-79dbafd7e74c
# ╠═e245dec6-7002-11eb-339f-c590a8afc985
# ╠═c57492b8-6640-11eb-1f31-29cc2440c6e0
# ╠═1d4b843c-66a1-11eb-188c-871b357bfdd7
# ╠═72c6e67c-66a1-11eb-2171-852cad67d86f
# ╠═ab527e9c-66a2-11eb-3531-8b666bb82451
# ╠═913a247c-66b8-11eb-296a-b92827dff6c8
# ╟─774b990a-66b9-11eb-2f71-ef6e4b08ec6f
# ╟─27a64b22-679d-11eb-0763-6d5931f2d0bb
# ╟─d299d78a-67b7-11eb-2bf5-192e423764d7
# ╟─ea8073d8-67b7-11eb-3a59-99769f0574f9
# ╠═a1b993b2-66bf-11eb-3df9-ddd3a3363f1b
# ╠═bb649190-67b9-11eb-012d-2bc5d384ddba
# ╠═056ca994-67ba-11eb-2459-5bf59f7db5c6
# ╠═2f4ab814-67ba-11eb-34e6-db97da22b6ad
# ╠═3c0653ec-67ba-11eb-1bce-e5ec4c92989e
# ╠═47d0d1c0-67ba-11eb-2d35-f32177efc1ab
# ╠═6b1ef04e-67ba-11eb-3c96-d15015059f76
# ╠═dbb73808-67be-11eb-12d0-916378cf0c7f
# ╠═58b6b5ea-67ba-11eb-3bca-1d89747d69bd
# ╠═7205031c-67ba-11eb-1e45-75ac2b20c76f
# ╠═7233ae8c-67bb-11eb-13dc-716bb11add02
# ╠═be271b86-67ba-11eb-2828-517b586bc4ad
# ╠═3e871518-67bd-11eb-100b-9dce3872fa40
# ╠═93838510-67bd-11eb-1ea0-7fcb8db20c6f
# ╠═171b5cca-67c0-11eb-3c51-c525076f9966
# ╠═7299b18c-67c0-11eb-1feb-2533c617d02d


text:
cd(@__DIR__); include("setups/simpletree.jl")
gr(dpi = 200)

## build rNGWF
distROT = natural_eigdist(𝚽, 𝛌, Q; α = 1.0, input_format = :pmf1, distance = :ROT)
rNGWF, dic_l2x = rngwf_all_vectors(distROT, 𝚽; σ = 0.1 * maximum(distROT), thres = 0.15)
Γ = rngwf_lx(dic_l2x)

just_eigenvec = [2, 16]
important_idx = sortperm((rNGWF' * f).^2; rev = true)
plot(layout = Plots.grid(3, 7), size = (1400, 1500))
for i = 1:21
    l, x = Γ[important_idx[i]]
    w = rNGWF[:, important_idx[i]]
    w ./= norm(w, 2)
    scatter!(X[:, 1], X[:, 2], marker_z = w, ms = 5, c = :viridis, subplot = i, mswidth = 0)
    plot!(frame = :none, cbar = false, grid = false, legend = false, subplot = i)
    if l in just_eigenvec
        plot!(title = latexstring("\\bar{\\psi}_{", l, ",", x, "} \\equiv \\mathbf{\\phi}_{$(l)}"), titlefont = 18, subplot = i)
    else
        plot!(title = latexstring("\\bar{\\psi}_{", l, ",", x, "}"), titlefont = 18, subplot = i)
    end
end
plt = current()
savefig(plt, "../figs/simpletree_f_rngwf_top21.png")


text:
## Sequential Search / Linear Search (sorted)

function search(list, num)
for i in sort(list)
        if num == i
            println("Found it")
            break
        elseif i > num
            println("Nope not here!")
            break
        end
    end
end

search([1,2,3,1,2,2,3,1,56,43,43,4,3,4,3,6,36,43,], 50)


text:
# Astronomical forcing term from
# de Saedeleer, B., Crucifix, M., & Wieczorek, S. (2013).
# Is the astronomical forcing a reliable and unique pacemaker for climate?
# A conceptual model study. Climate Dynamics, 40(1–2), 273–294.
# https://doi.org/10.1007/s00382-012-1316-1

function astronomical_forcing(A, t)
    # notice that mean value (494.2447 W/m²) of forcing has been removed
    F = 0.0
    @inbounds for i in 1:size(A,1)
        ωt = A[i, 1]*t
        F += A[i, 2]*sin(ωt) + A[i, 3]*cos(ωt)
    end
    return F/11.77 # value of α_ε_1 in W/m^2
end

A = [ # Notice: all these harmonics can be reduced massively to just 6 (See Fig 2)
      # so if the code becomes too slow, I can do this.
    # Obliquity terms
    0.153249478547167 -11.2287376815124 3.51682075211241
    0.158148666238883 -3.82499371467540 -0.761851750263805
    0.117190147169570 2.28814805956066 1.80233702684623
    0.155061775112933 -1.29770081956440 -0.635152963728496
    0.217333905941751 0.380973541305497 -1.46301711999210
    0.150162587421217 1.54904176353302 -0.0883941912769817
    0.211709630908568 -0.810768209286259 -0.577980646565494
    0.156336369673117 -0.918358442095885 0.196083726889428
    0.148350290855451 0.256895610735773 -0.524697312305024
    0.206924898030688 -0.335783913402678 -0.0194792150128644
    0.212525165090383 0.267659228540196 0.128915417116900
    0.229992875969202 0.0696189733188958 0.0746231714061285
    0.306498957094334 0.0247349748169616 0.0140464395340974
    0.311398144786051 0.0138353727621181 0.0304736668840422
    0.004899187691716 -0.160479848721994 0.0594077968934257
    # Precession terms
    0.264933601588513 -15.5490493322904 -9.70406287110532
    0.280151350350945 15.4319556361701 4.75247271131525
    0.331110950251899 9.0992249352734 -10.6115244887390
    0.328024059125949 -7.87065384013669 6.61544246063503
    0.326211762560183 0.813786144754451 -4.52641408099246
    0.269742342439881 0.0690448504314857 -3.31639260969558
    0.332923246817665 1.44050770785967 1.06339286050120
    0.371638925683567 0.925324276580528 -1.02066758672154
    0.275366617473065 0.997628846513796 -0.362906496840039
    0.323124871434233 -0.378637986107629 0.527217891742183
    0.259396912994958 0.339477750517033 -0.560509461538342
    0.324937167999999 -0.576082669762308 1.18669572739338
    0.334197841377850 0.346906064369828 -0.648189701487285
    0.274551083291250 -0.441772417569753 0.289576210423804
    0.418183080135680 -0.0184884064645011 0.109632390175297
    0.111684123041346 -0.428006728186239 0.357006342316690
    0.433400828898112 -0.0049199219454561 -0.106148873639336
    0.126901871803777 0.257509918217341 -0.377639794223366
    0.336010137943616 -0.421809264016129 0.324327509437558
    0.177861471704732 -0.161827722328271 -0.362683869407858
]

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
struct DensityAnalysis
    options::Dict{Symbol, Any}
end
DensityAnalysis(; kwargs...) = DensityAnalysis(Dict{Symbol, Any}(kwargs))

# Work around lack of length 1 tuple method
# TODO: also add weights support here
_kde(data::NTuple{1, Any}; kwargs...) = kde(data...; kwargs...)
_kde(data::Tuple; kwargs...) = kde(data; kwargs...)

function _density(data...; extrema, npoints=200, kwargs...)
    k = _kde(data; kwargs...)
    rgs = map(e -> range(e...; length=npoints), extrema)
    res = pdf(k, rgs...)
    return (rgs..., res)
end

function (d::DensityAnalysis)(le::Entry)
    summaries = map(summary∘getvalue, le.mappings.positional)
    extrema = get(d.options, :extrema, Tuple(summaries))
    options = merge(d.options, pairs((; extrema)))
    return splitapply(le) do entry
        labels, mappings = map(getlabel, entry.mappings), map(getvalue, entry.mappings)
        result = _density(mappings.positional...; mappings.named..., options...)
        labeled_result = map(Labeled, vcat(labels.positional, "pdf"), collect(result))
        plottypes = [LinesFill, Heatmap, Volume]
        default_plottype = plottypes[length(mappings.positional)]
        return Entry(
            AbstractPlotting.plottype(entry.plottype, default_plottype),
            Arguments(labeled_result),
            entry.attributes
        )
    end
end

"""
    density(; extrema, npoints, kernel, bandwidth)

Fit a kernel density estimation of `data`.
"""
density(; kwargs...) = Layer((DensityAnalysis(; kwargs...),))


text:
@everywhere module TestDecisionVariable

using StochasticPrograms
using Test

function test_decision_no_bound(Structure)
    ξ₁ = @scenario a = 1 probability = 0.5
    ξ₂ = @scenario a = 2 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(model, x)
    end
    @second_stage sp = begin
        @recourse(model, y)
    end
    # First-stage
    x = sp[1,:x]
    @test StochasticPrograms.stage(x) == 1
    @test !JuMP.has_lower_bound(x)
    @test !JuMP.has_upper_bound(x)
    @test !JuMP.is_fixed(x)
    @test "x" == @inferred JuMP.name(x)
    @test x == decision_by_name(owner_model(x), 1, "x")
    # Scenario-dependent
    y = sp[2,:y]
    @test StochasticPrograms.stage(y) == 2
    @test_throws ErrorException JuMP.has_lower_bound(y)
    @test !JuMP.has_lower_bound(y, 1)
    @test !JuMP.has_lower_bound(y, 2)
    @test_throws ErrorException JuMP.has_upper_bound(y)
    @test !JuMP.has_upper_bound(y, 1)
    @test !JuMP.has_upper_bound(y, 2)
    @test_throws ErrorException JuMP.is_fixed(y)
    @test !JuMP.is_fixed(y, 1)
    @test !JuMP.is_fixed(y, 2)
    @test_throws ErrorException JuMP.is_integer(y)
    @test !JuMP.is_integer(y, 1)
    @test !JuMP.is_integer(y, 2)
    @test_throws ErrorException JuMP.is_binary(y)
    @test !JuMP.is_binary(y, 1)
    @test !JuMP.is_binary(y, 2)
    @test "y" == @inferred JuMP.name(y)
    @test y == decision_by_name(owner_model(y), 2, "y")
end

function test_decision_lower_bound(Structure)
    ξ₁ = @scenario a = 1 probability = 0.5
    ξ₂ = @scenario a = 2 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(model, x >= 0, Bin)
    end
    @second_stage sp = begin
        @uncertain a
        @recourse(model, y >= a, Bin)
    end
    # First-stage
    x = sp[1,:x]
    @test JuMP.has_lower_bound(x)
    @test 0.0 == @inferred JuMP.lower_bound(x)
    @test !JuMP.has_upper_bound(x)
    @test !JuMP.is_fixed(x)
    @test JuMP.is_binary(x)
    @test !JuMP.is_integer(x)
    JuMP.delete_lower_bound(x)
    @test !JuMP.has_lower_bound(x)
    # Scenario-dependent
    y = sp[2,:y]
    @test JuMP.has_lower_bound(y, 1)
    @test JuMP.has_lower_bound(y, 2)
    @test 1.0 == @inferred JuMP.lower_bound(y, 1)
    @test 2.0 == @inferred JuMP.lower_bound(y, 2)
    @test !JuMP.has_upper_bound(y, 1)
    @test !JuMP.has_upper_bound(y, 2)
    @test !JuMP.is_fixed(y, 1)
    @test !JuMP.is_fixed(y, 2)
    @test JuMP.is_binary(y, 1)
    @test JuMP.is_binary(y, 2)
    @test !JuMP.is_integer(y, 1)
    @test !JuMP.is_integer(y, 2)
    JuMP.delete_lower_bound(y, 1)
    @test !JuMP.has_lower_bound(y, 1)
    JuMP.delete_lower_bound(y, 2)
    @test !JuMP.has_lower_bound(y, 2)
end

function test_decision_upper_bound(Structure)
    ξ₁ = @scenario a = 1 probability = 0.5
    ξ₂ = @scenario a = 2 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(model, x <= 1.0, Int)
    end
    @second_stage sp = begin
        @uncertain a
        @recourse(model, y <= a, Int)
    end
    # First-stage
    x = sp[1,:x]
    @test !JuMP.has_lower_bound(x)
    @test JuMP.has_upper_bound(x)
    @test 1.0 == @inferred JuMP.upper_bound(x)
    @test !JuMP.is_fixed(x)
    @test !JuMP.is_binary(x)
    @test JuMP.is_integer(x)
    JuMP.delete_upper_bound(x)
    @test !JuMP.has_upper_bound(x)
    # Scenario-dependent
    y = sp[2,:y]
    @test !JuMP.has_lower_bound(y, 1)
    @test !JuMP.has_lower_bound(y, 2)
    @test JuMP.has_upper_bound(y, 1)
    @test JuMP.has_upper_bound(y, 2)
    @test 1.0 == @inferred JuMP.upper_bound(y, 1)
    @test 2.0 == @inferred JuMP.upper_bound(y, 2)
    @test !JuMP.is_fixed(y, 1)
    @test !JuMP.is_fixed(y, 2)
    @test !JuMP.is_binary(y, 1)
    @test !JuMP.is_binary(y, 2)
    @test JuMP.is_integer(y, 1)
    @test JuMP.is_integer(y, 2)
    JuMP.delete_upper_bound(y, 1)
    @test !JuMP.has_upper_bound(y, 1)
    JuMP.delete_upper_bound(y, 2)
    @test !JuMP.has_upper_bound(y, 2)
end

function test_decision_fix(Structure)
    ξ₁ = @scenario a = 1 probability = 0.5
    ξ₂ = @scenario a = 2 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(model, x == 1)
    end
    @second_stage sp = begin
        @uncertain a
        @recourse(model, y == a)
    end
    # First-stage
    x = sp[1,:x]
    @test !JuMP.has_lower_bound(x)
    @test !JuMP.has_upper_bound(x)
    @test JuMP.is_fixed(x)
    @test 1.0 == @inferred JuMP.value(x)
    JuMP.unfix(x)
    @test !JuMP.is_fixed(x)
    JuMP.fix(x, 2.0)
    @test !JuMP.has_lower_bound(x)
    @test !JuMP.has_upper_bound(x)
    @test JuMP.is_fixed(x)
    @test 2.0 == @inferred JuMP.value(x)
    # Scenario-dependent
    y = sp[2,:y]
    @test !JuMP.has_lower_bound(y, 1)
    @test !JuMP.has_lower_bound(y, 2)
    @test !JuMP.has_upper_bound(y, 1)
    @test !JuMP.has_upper_bound(y, 2)
    @test JuMP.is_fixed(y, 1)
    @test JuMP.is_fixed(y, 2)
    @test 1.0 == @inferred JuMP.value(y, 1)
    @test 2.0 == @inferred JuMP.value(y, 2)
    JuMP.unfix(y, 1)
    @test !JuMP.is_fixed(y, 1)
    @test JuMP.is_fixed(y, 2)
    JuMP.unfix(y, 2)
    @test !JuMP.is_fixed(y, 2)
    JuMP.fix(y, 1, 2.0)
    @test !JuMP.has_lower_bound(y, 1)
    @test !JuMP.has_upper_bound(y, 1)
    @test JuMP.is_fixed(y, 1)
    @test !JuMP.is_fixed(y, 2)
    @test 2.0 == @inferred JuMP.value(y, 1)
    JuMP.fix(y, 2, 1.0)
    @test !JuMP.has_lower_bound(y, 2)
    @test !JuMP.has_upper_bound(y, 2)
    @test JuMP.is_fixed(y, 2)
    @test 1.0 == @inferred JuMP.value(y, 2)
end

function test_decision_custom_index_sets(Structure)
    ξ₁ = @scenario a = 1 probability = 0.5
    ξ₂ = @scenario a = 2 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(model, x[0:1, 10:20, 1:1] >= 2)
    end
    @second_stage sp = begin
        @uncertain a
        @recourse(model, y[0:1, 10:20, 1:1] >= a)
    end
    # First-stage
    x = sp[1,:x]
    @test JuMP.has_lower_bound(x[0, 15, 1])
    @test 2 == @inferred JuMP.lower_bound(x[0, 15, 1])
    @test !JuMP.has_upper_bound(x[0, 15, 1])
    @first_stage sp = begin
        @decision(model, x[i in -10:10, s in [:a,:b]] <= 5.5, Int)
    end
    generate!(sp)
    x = sp[1,:x]
    @test 5.5 == @inferred JuMP.upper_bound(x[-4, :a])
    @test "x[-10,a]" == @inferred JuMP.name(x[-10,:a])
    @test x[-10,:a] == decision_by_name(sp, 1, "x[-10,a]")
    # Scenario-dependent
    y = sp[2,:y]
    @test JuMP.has_lower_bound(y[0, 15, 1], 1)
    @test JuMP.has_lower_bound(y[0, 15, 1], 2)
    @test 1 == @inferred JuMP.lower_bound(y[0, 15, 1], 1)
    @test 2 == @inferred JuMP.lower_bound(y[0, 15, 1], 2)
    @test !JuMP.has_upper_bound(y[0, 15, 1], 1)
    @test !JuMP.has_upper_bound(y[0, 15, 1], 2)
    @second_stage sp = begin
        @recourse(model, y[i in -10:10, s in [:a,:b]] <= 5.5, Int)
    end
    y = sp[2,:y]
    @test 5.5 == @inferred JuMP.upper_bound(y[-4, :a], 1)
    @test 5.5 == @inferred JuMP.upper_bound(y[-4, :a], 2)
    @test "y[-10,a]" == @inferred JuMP.name(y[-10,:a])
    @test nothing === decision_by_name(sp, 1, "y[-10,a]")
    @test y[-10,:a] == decision_by_name(sp, 2, "y[-10,a]")
end

function test_variable_is_valid_delete(Structure)
    ξ₁ = @scenario a = 1 probability = 0.5
    ξ₂ = @scenario a = 2 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(model, x)
    end
    @second_stage sp = begin
        @recourse(model, y)
    end
    # First-stage
    x = sp[1,:x]
    @test JuMP.is_valid(sp, x)
    JuMP.delete(sp, x)
    @test !JuMP.is_valid(sp, x)
    @test_throws Exception JuMP.delete(sp, x)
    # Scenario-dependent
    y = sp[2,:y]
    @test_throws ErrorException JuMP.is_valid(sp, y)
    @test JuMP.is_valid(sp, y, 1)
    @test JuMP.is_valid(sp, y, 2)
    @test_throws ErrorException JuMP.delete(sp, y)
    JuMP.delete(sp, y, 1)
    @test !JuMP.is_valid(sp, y, 1)
    @test JuMP.is_valid(sp, y, 2)
    JuMP.delete(sp, y, 2)
    @test !JuMP.is_valid(sp, y, 2)
    @test_throws Exception JuMP.delete(sp, y, 1)
    @first_stage sp = begin
        @decision(model, x[1:3] >= 1)
    end
    @second_stage sp = begin
        @recourse(model, y[1:3] >= 1)
    end
    # First-stage
    x = sp[1,:x]
    @test all(is_valid.(sp, x))
    delete(sp, x)
    @test all((!is_valid).(sp, x))
    second_sp = StochasticProgram([Scenario()], Structure...)
    @test_throws Exception JuMP.delete(second_sp, x[1])
    @test_throws Exception JuMP.delete(second_sp, x)
    # Scenario-dependent
    y = sp[2,:y]
    @test_throws ErrorException is_valid.(sp, y)
    @test all(is_valid.(sp, y, 1))
    @test all(is_valid.(sp, y, 2))
    @test_throws ErrorException delete.(sp, y)
    delete(sp, y, 1)
    @test all((!is_valid).(sp, y, 1))
    @test all(is_valid.(sp, y, 2))
    delete(sp, y, 2)
    @test all((!is_valid).(sp, y, 2))
end

function test_variable_bounds_set_get(Structure)
    ξ₁ = @scenario a = 1 probability = 0.5
    ξ₂ = @scenario a = 2 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(model, 0 <= x <= 2)
    end
    @second_stage sp = begin
        @uncertain a
        @recourse(model, 0 <= y <= a)
    end
    # First-stage
    x = sp[1,:x]
    @test 0 == @inferred JuMP.lower_bound(x)
    @test 2 == @inferred JuMP.upper_bound(x)
    set_lower_bound(x, 1.)
    @test 1. == @inferred JuMP.lower_bound(x)
    set_upper_bound(x, 3.)
    @test 3. == @inferred JuMP.upper_bound(x)
    # Scenario-dependent
    y = sp[2,:y]
    @test_throws ErrorException JuMP.lower_bound(y)
    @test 0 == @inferred JuMP.lower_bound(y, 1)
    @test 0 == @inferred JuMP.lower_bound(y, 2)
    @test_throws ErrorException JuMP.upper_bound(y)
    @test 1 == @inferred JuMP.upper_bound(y, 1)
    @test 2 == @inferred JuMP.upper_bound(y, 2)
    @test_throws ErrorException JuMP.set_lower_bound(y, 1.)
    set_lower_bound(y, 1, 1.)
    @test 1. == @inferred JuMP.lower_bound(y, 1)
    set_lower_bound(y, 2, 1.)
    @test 1. == @inferred JuMP.lower_bound(y, 2)
    @test_throws ErrorException JuMP.set_upper_bound(y, 3.)
    set_upper_bound(y, 1, 3.)
    @test 3. == @inferred JuMP.upper_bound(y, 1)
    set_upper_bound(y, 2, 3.)
    @test 3. == @inferred JuMP.upper_bound(y, 2)
end

function test_variable_starts_set_get(Structure)
    ξ₁ = @scenario a = 1 probability = 0.5
    ξ₂ = @scenario a = 2 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(model, x[1:3])
    end
    @second_stage sp = begin
        @recourse(model, y[1:3])
    end
    # First-stage
    x = sp[1,:x]
    x0 = collect(1:3)
    JuMP.set_start_value.(x, x0)
    @test JuMP.start_value.(x) == x0
    @test JuMP.start_value.([x[1],x[2],x[3]]) == x0
    # First-stage
    y = sp[2,:y]
    @test_throws ErrorException JuMP.set_start_value.(y, x0)
    JuMP.set_start_value.(y, 1, x0)
    JuMP.set_start_value.(y, 2, x0)
    @test_throws ErrorException JuMP.start_value.(y)
    @test JuMP.start_value.(y, 1) == x0
    @test JuMP.start_value.(y, 2) == x0
    @test_throws ErrorException JuMP.start_value.([y[1],y[2],y[3]])
    @test JuMP.start_value.([y[1],y[2],y[3]], 1) == x0
    @test JuMP.start_value.([y[1],y[2],y[3]], 2) == x0
end

function test_variable_integrality_set_get(Structure)
    ξ₁ = @scenario a = 1 probability = 0.5
    ξ₂ = @scenario a = 2 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(model, x[1:3])
    end
    @second_stage sp = begin
        @recourse(model, y[1:3])
    end
    # First-stage
    x = sp[1,:x]
    JuMP.set_integer(x[2])
    JuMP.set_integer(x[2])
    @test JuMP.is_integer(x[2])
    JuMP.unset_integer(x[2])
    @test !JuMP.is_integer(x[2])
    JuMP.set_binary(x[1])
    JuMP.set_binary(x[1])
    @test JuMP.is_binary(x[1])
    @test_throws Exception JuMP.set_integer(x[1])
    JuMP.unset_binary(x[1])
    @test !JuMP.is_binary(x[1])
    # Scenario-dependent
    y = sp[2,:y]
    @test_throws ErrorException JuMP.set_integer(y[2])
    JuMP.set_integer(y[2], 1)
    JuMP.set_integer(y[2], 1)
    JuMP.set_integer(y[3], 2)
    JuMP.set_integer(y[3], 2)
    @test_throws ErrorException JuMP.is_integer(y[2])
    @test JuMP.is_integer(y[2], 1)
    @test !JuMP.is_integer(y[2], 2)
    @test JuMP.is_integer(y[3], 2)
    @test !JuMP.is_integer(y[3], 1)
    @test_throws ErrorException JuMP.unset_integer(y[2])
    JuMP.unset_integer(y[2], 1)
    JuMP.unset_integer(y[3], 2)
    @test !JuMP.is_integer(y[2], 1)
    @test_throws ErrorException JuMP.set_binary(y[2])
    JuMP.set_binary(y[1], 1)
    JuMP.set_binary(y[1], 1)
    JuMP.set_binary(y[1], 2)
    JuMP.set_binary(y[1], 2)
    @test_throws ErrorException JuMP.is_binary(y[1])
    @test JuMP.is_binary(y[1], 1)
    @test JuMP.is_binary(y[1], 2)
    @test_throws Exception JuMP.set_integer(y[1], 1)
    @test_throws Exception JuMP.set_integer(y[1], 2)
    @test_throws ErrorException JuMP.unset_binary(y[1])
    JuMP.unset_binary(y[1], 1)
    @test !JuMP.is_binary(y[1], 1)
    JuMP.unset_binary(y[1], 2)
    @test !JuMP.is_binary(y[1], 2)
end

function test_variables_constrained_on_creation(Structure)
    ξ₁ = @scenario a = 1 probability = 0.5
    ξ₂ = @scenario a = 2 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(model, x[1:2] in SecondOrderCone())
    end
    @second_stage sp = begin
        @recourse(model, y[1:2] in SecondOrderCone())
    end
    x = sp[1,:x]
    y = sp[2,:y]
    @test num_constraints(sp, 1, typeof(x), MOI.SecondOrderCone) == 1
    @test num_constraints(sp, 2, typeof(y), MOI.SecondOrderCone) == 1
    @test name(x[1]) ==  "x[1]"
    @test name(x[2]) ==  "x[2]"
    @test name(y[1]) == "y[1]"
    @test name(y[2]) == "y[2]"
    @first_stage sp = begin
        @decision(model, x[1:2] in SecondOrderCone())
        @decision(model, [1:2] in SecondOrderCone())
    end
    @second_stage sp = begin
        @recourse(model, y[1:2] in SecondOrderCone())
        @recourse(model, [1:2] in SecondOrderCone())
    end
    x = sp[1,:x]
    y = sp[2,:y]
    @test num_constraints(sp, 1, typeof(x), MOI.SecondOrderCone) == 2
    @test num_constraints(sp, 2, typeof(y), MOI.SecondOrderCone) == 2
    @first_stage sp = begin
        @decision(model, x[1:2] in SecondOrderCone())
        @decision(model, [1:2] in SecondOrderCone())
        @decision(model, [1:3] in SecondOrderCone())
    end
    @second_stage sp = begin
        @recourse(model, y[1:2] in SecondOrderCone())
        @recourse(model, [1:2] in SecondOrderCone())
        @recourse(model, [1:3] in SecondOrderCone())
    end
    x = sp[1,:x]
    y = sp[2,:y]
    @test num_constraints(sp, 1, typeof(x), MOI.SecondOrderCone) == 3
    @test num_constraints(sp, 2, typeof(y), MOI.SecondOrderCone) == 3
    @first_stage sp = begin
        @decision(model, z in MOI.Semiinteger(1.0, 2.0))
    end
    generate!(sp)
    z = sp[1,:z]
    @test num_constraints(sp, 1, typeof(z), MOI.Semiinteger{Float64}) == 1
    @first_stage sp = begin
        @decision(model, z in MOI.Semiinteger(1.0, 2.0))
        @decision(model, set = MOI.Semiinteger(1.0, 2.0))
    end
    generate!(sp)
    z = sp[1,:z]
    @test num_constraints(sp, 1, typeof(z), MOI.Semiinteger{Float64}) == 2
end

function test_all_decision_variables(Structure)
    ξ₁ = @scenario a = 1 probability = 0.5
    ξ₂ = @scenario a = 2 probability = 0.5
    sp = StochasticProgram([ξ₁,ξ₂], Structure...)
    @first_stage sp = begin
        @decision(model, x₁)
        @decision(model, x₂)
    end
    @second_stage sp = begin
        @recourse(model, y₁)
        @recourse(model, y₂)
    end
    # First-stage
    x₁ = sp[1,:x₁]
    x₂ = sp[1,:x₂]
    @test [x₁, x₂] == @inferred all_decision_variables(sp, 1)
    # Scenario-dependent
    y₁ = sp[2,:y₁]
    y₂ = sp[2,:y₂]
    @test [y₁, y₂] == @inferred all_decision_variables(sp, 2)
end

function runtests()
    @testset "DecisionVariable" begin
        for structure in [(Deterministic(),),
                          (Vertical(),),
                          (Horizontal(),),
                          (Deterministic(), () -> MOIU.MockOptimizer(MOIU.Model{Float64}()))]
            name = length(structure) == 1 ? "$(structure[1])" : "$(structure[1]) with decision bridges"
            @testset "$name" begin
                for name in names(@__MODULE__; all = true)
                    if !startswith("$(name)", "test_")
                        continue
                    end
                    f = getfield(@__MODULE__, name)
                    @testset "$(name)" begin
                        f(structure)
                    end
                end
            end
        end
    end
end

function run_dtests()
    @testset "DecisionVariable" begin
        for structure in [(DistributedVertical(),), (DistributedHorizontal(),)]
            @testset "$(structure[1])" begin
                for name in names(@__MODULE__; all = true)
                    if !startswith("$(name)", "test_")
                        continue
                    end
                    f = getfield(@__MODULE__, name)
                    @testset "$(name)" begin
                        f(structure)
                    end
                end
            end
        end
    end
end

end


text:
using Printf

include("gen_rmesh_exp.jl")
include("lagrange_interp.jl")
include("radial_interp.jl")
include("rsch_integ_rk4.jl")
include("sch_rk4_step.jl")
include("get_min_idx.jl")
include("get_n_nodes.jl")

function main()
    # Build radial function
    r_min = 1e-9
    r_max = 50.0
    a = 1e9
    Nr = 1000
    rmesh = gen_rmesh_exp(r_min, r_max, a, Nr)

    Z = 1
    l = 0

    P = zeros(Float64,Nr)
    Q = zeros(Float64,Nr)
    V = -Z ./ rmesh

    E = -0.125
    imax = rsch_integ_rk4!( E, Z, l, rmesh, V, P, Q )
    println("imax = ", imax)

    println("Some P and Q")
    for i in 1:5
        @printf("%8d %20.10e %20.10e\n", i, P[i], Q[i])
    end
    @printf("...\n")
    for i in imax-5:imax
        @printf("%8d %20.10e %20.10e\n", i, P[i], Q[i])
    end

    minidx = get_min_idx(imax, P)
    println("minidx = ", minidx)
    println("imax   = ", imax)

    Nnodes = get_n_nodes( minidx-1, P )
    println("Nnodes = ", Nnodes)

    println("Pass here")
end

@time main()
@time main()


text:
###########################
# FunctionCell Interface  #
###########################
"""
	sample([f,] opts)

Returns a cell that samples data piped to it according to the sampling options specified
in `opts`. If the dataset is labeled, an additional function `f` can be specified to
obtain the labels calling `LearnBase.targets(f,data)`. The argument `opts` can be a
number specifying the number of samples or fraction desired or, a `Dict(class=>samples)`
if different sampling is to be applied for different classes. In this case,

* `class` has to have the same type as the labels
* `samples` is a number. If it is and `Int`, it specifies the exact number of samples
from its corresponding class to sample. If it is a `Float`, a fraction of samples is
generated.

When data is passed into the sampling cell, a datasubset of the original data and
corresponding indices in the dataset are returned.

# Examples
```
julia> a = datacell([1 2 3; 4 5 6],["a","a","b"])
DataCell, 3 obs, 2 vars, 1 target(s)/obs, 2 distinct values: "b"(1),"a"(2)

julia> +a
2×3 Array{Int64,2}:
 1  2  3
 4  5  6

julia> a |> sample(7)
([*]DataCell, 7 obs, 2 vars, 1 target(s)/obs, 2 distinct values: "b"(1),"a"(6), [1, 2, 2, 3, 2, 1, 2])

julia> +(ans[1])
2×7 SubArray{Int64,2,Array{Int64,2},Tuple{Base.Slice{Base.OneTo{Int64}},Array{Int64,1}},false}:
 1  2  2  3  2  1  2
 4  5  5  6  5  4  5

julia> a = datacell([1 2 3 4 5; 6 7 8 9 0],["a","a","b","b","b"]);

julia> +a
2×5 Array{Int64,2}:
 1  2  3  4  5
 6  7  8  9  0

julia> b = sample(a,Dict("a"=>4, "b"=>2.3))
([*]DataCell, 11 obs, 2 vars, 1 target(s)/obs, 2 distinct values: "b"(7),"a"(4), [4, 4, 4, 3, 4, 5, 4, 1, 2, 2, 2])

julia> +b[1]
2×11 SubArray{Int64,2,Array{Int64,2},Tuple{Base.Slice{Base.OneTo{Int64}},Array{Int64,1}},false}:
 4  4  4  3  4  5  4  1  2  2  2
 9  9  9  8  9  0  9  6  7  7  7
```
"""
sample(f::Function, opts::S where S<:Dict{T} where T<:DataElement) = FunctionCell(sample, (f,opts), "Labeled data sampler")
sample(opts::S where S<:Dict{T} where T<:DataElement) = sample(identity, opts)

sample(f::Function, opts::S where S<:Real) = FunctionCell(sample, (f,opts), "Data sampler")
sample(opts::S where S<:Real) = sample(identity, opts)

############################
# DataCell/Array Interface #
############################
"""
	sample(data, [f,] opts)

Samples `data` according to the sampling options specified in `opts`. If the dataset is labeled,
an additional function `f` can be specified to obtain the labels calling `LearnBase.targets(data,f)`.
Returns a `datasubset` and the indices corresponding to the sampled observations.
"""
# Sampling based on number/percentage of samples desired (f is ignored)
sample(x::T where T<:CellData, opts) = begin
	d, idx = sample(strip(x), identity, opts)
	return (datacell(d), idx)
end

sample(x::T where T<:CellData, f::Function, opts::S where S<:Real) = begin
	d, idx = sample(strip(x), f, opts)
	return (datacell(d), idx)
end

sample(x::T where T<:Union{AbstractArray, Tuple{<:AbstractArray, <:AbstractArray}}, f::Function, opts::S where S<:Real) = begin
	idx = _sample_(opts)(collect(1:nobs(x)))
	return (datasubset(x,idx), idx)
end

# Class-based sampling
sample(x::T where T<:Union{CellDataL, CellDataLL}, opts::Dict) = sample(x, identity, opts)

sample(x::T where T<:Union{CellDataL, CellDataLL}, f::Function, opts::Dict) = begin
	d, idx = sample(strip(x), f, opts)
	return (datacell(d), idx)
end

sample(x::T where T<:Tuple{<:AbstractArray, <:AbstractArray}, f::Function, opts::Dict) = begin
	lm = labelmap(targets(f,x[2]))
	idx = reduce(vcat, (k in keys(opts) ? _sample_(opts[k])(lm[k]) : lm[k] for k in keys(lm)))
	return (datasubset(x, idx), idx)
end

# Main sampling function: return the elements of a vector in proportion to a integer (exact number of returned samples)
# or float (fraction of original number of samples)
_sample_(m::T where T<:Int) = (x)->(m <= length(x)) ? x[randperm(length(x))[1:m]] : rand(x,m)
_sample_(m::T where T<:AbstractFloat) =
	(x)->begin
		mi = round(eltype(1), m*length(x))
		n = length(x)
		(mi <= n) ? x[randperm(n)[1:mi]] : rand(x,mi)
	end


text:
# All function related to the different Delay pre-selection statistics. Currently
# there is only the continuity statistic [`pecora`](@ref) and the range-function
# implemented (range function means that there is actually no preselection performed)

# Constructors:

## Constructors for DelayPreSelection Functions
"""
    Continuity_function <: AbstractDelayPreselection

    Constructor for the continuity function `⟨ε★⟩` by Pecora et al.[^Pecora2007],
    see [`pecora`](@ref).

    ## Fieldnames
    * `K::Int`: the amount of nearest neighbors in the δ-ball. Must be at
      least 8 (in order to gurantee a valid statistic). `⟨ε★⟩` is computed taking
      the minimum result over all `k ∈ K` (read algorithm description in [`pecora`](@ref)).
    * `samplesize::Real`: determine the fraction of all phase space points
      to be considered (fiducial points v) to average ε★ to produce `⟨ε★⟩`.
    * `α::Real = 0.05`: The significance level for obtaining the continuity statistic
    * `p::Real = 0.5`: The p-parameter for the binomial distribution used for the
      computation of the continuity statistic ⟨ε★⟩.

    ## Defaults
    * When calling `Continuity_function()` a Continuity_function-object is constructed
      with `K=13`, `samplesize=1.`, `α=0.05` and `p=0.5`.

"""
struct Continuity_function <: AbstractDelayPreselection
    K::Int
    samplesize::Real
    α::Real
    p::Real
    # Constraints and Defaults
    Continuity_function(k=13,x=1.,y=0.05,z=0.5) = begin
        @assert k > 7 "At least 8 nearest neighbors must be in the δ-ball, in order to gurantee a valid statistic."
        @assert 0. < x ≤ 1. "Please select a valid `samplesize`, which denotes a fraction of considered fiducial points, i.e. `samplesize` ∈ (0 1]"
        @assert 1. > y > 0.
        @assert 1. > z > 0.
        new(k,x,y,z)
    end
end

"""
    Range_function <: AbstractDelayPreselection

    Constructor for a range of possible delay values. In this case there is
    actually no "pre-selection" of delay, but rather all possible delays, given
    in the input `τs` (see, [`mcdts_embedding`](@ref)) are considered. This can significantly affect the
    computation time. There are no fieldnames, simply construct by typing
    `RangeFunction()`.
"""
struct Range_function <: AbstractDelayPreselection end

# Functions:

"""
    get_max_idx(Λ::AbstractDelayPreselection, dps::Vector, τ_vals, ts_vals) → max_idx

    Compute the candidate delay values from the given delay pre-selection statistic
    `dps` with respect to `Λ`, which determined how `dps` was obtained and how
    to select the candidates (e.g. pick the maxima of `dps` in case of the
    `Λ` being the Continuity function). See [`Continuity_function`](@ref) and
    [`Range_function`](@ref).
"""
function get_max_idx(Λ::Range_function, dps::Vector{T}, τ_vals, ts_vals, ts) where {T}
    max_idx = Vector(dps[2:end].+1)
    ts_idx = findall(e->e==ts, ts_vals) # do not consider already taken delays
    filter!(e->e∉(τ_vals[ts_idx] .+ 2), max_idx) # do not consider already taken delays
    return max_idx
end
function get_max_idx(Λ::Continuity_function, dps::Vector{T}, τ_vals, ts_vals, ts) where {T}
    _, max_idx = get_maxima(dps) # determine local maxima in delay_pre_selection_statistic
    ts_idx = findall(e->e==ts, ts_vals) # do not consider already taken delays
    filter!(e->e∉(τ_vals[ts_idx] .+ 2), max_idx) # do not consider already taken delays
    return max_idx
end

## Methods for delay preselection stats

"""
    get_delay_statistic(optimalg.Λ<: AbstractDelayPreselection, Ys, τs, w, τ_vals, ts_vals; kwargs... )

    Compute the delay statistic according to the chosen method in `optimalg.Λ` (see [`MCDTSOptimGoal`](@ref))
"""
function get_delay_statistic(Λ::Continuity_function, Ys, τs, w, τ_vals, ts_vals; metric = Euclidean(), kwargs... )

    ε★, _ = DelayEmbeddings.pecora(Ys, Tuple(τ_vals), Tuple(ts_vals); delays = τs, w = w,
            samplesize = Λ.samplesize, K = Λ.K, metric = metric, α = Λ.α,
            p = Λ.p)
    return ε★
end
function get_delay_statistic(Λ::Range_function, Ys, τs, w, τ_vals, ts_vals; kwargs... )
    return repeat(Vector(1:length(τs)), outer = [1,size(Ys,2)])
end


text:
# Load simulation engine
# import FLOWFVS
reload("FLOWUnsteady")
uns = FLOWUnsteady
vlm = uns.vlm
gt = vlm.vtk

# ------------ GENERIC MODULES -------------------------------------------------
using Base.Test
using PyPlot
using JLD

# ------------ HEADERS ---------------------------------------------------------
# Load modules
for module_name in ["bertinswing"]
    include("test_"*module_name*".jl")
end

# ------------ TESTS -----------------------------------------------------------

# Test VLM solver: Isolated wing
@test bertin_VLM(; wake_coupled=false, nsteps=0, verbose=true, disp_plot=true)

# Test VLM regularization: Isolated wing
@test bertin_VLM(; wake_coupled=false, vlm_fsgm=0.00125, nsteps=0, verbose=true, disp_plot=true)

# Test VPM+VLM coupling: Isolated wing
@test bertin_VLM(; wake_coupled=true, nsteps=150, verbose=true, disp_plot=true)

# Test VPM+VLM solver on kinematic velocity: Isolated wing
@test bertin_kinematic(; nsteps=150, verbose=true, disp_plot=true)

# Test simulation stability: Isolated wing with a very fine time stepping and shedding
@test bertin_kinematic(; nsteps=300, p_per_step=4, vlm_rlx=0.75, verbose=true, disp_plot=true)

# Test quasisteady solver: No VPM wake shedding
@test bertin_kinematic(; nsteps=10, VehicleType=uns.QVLMVehicle, verbose=true, disp_plot=true)

# @test bertin_kinematic(; nsteps=150, p_per_step=4, vlm_rlx=0.75, save_path="temps/bertinswing11/", verbose=true, disp_plot=true)

max_stars_repo_name: dendisuhubdy/FLOWUnsteady

max_issues_repo_name: dendisuhubdy/FLOWUnsteady

max_forks_repo_name: dendisuhubdy/FLOWUnsteady


text:
# # Linear regression

# Import the necessary packages:

using Distributions
using MLJBase
using Soss
using SossMLJ
using Statistics

# In this example, we fit a Bayesian linear regression model with the
# canonical link function.

# Suppose that we are given a matrix of features `X` and a column vector of
# labels `y`. `X` has `n` rows and `p` columns. `y` has `n` elements. We assume
# that our observation vector `y` is a realization of a random variable `Y`.
# We define `μ` (mu) as the expected value of `Y`, i.e. `μ := E[Y]`. Our model
# comprises three components:
#
# 1. The probability distribution of `Y`: for linear regression, we assume that each `Yᵢ` follows a normal distribution with mean `μᵢ` and variance `σ²`.
# 2. The systematic component, which consists of linear predictor `η` (eta), which we define as `η := Xβ`, where `β` is the column vector of `p` coefficients.
# 3. The link function `g`, which provides the following relationship: `g(E[Y]) = g(μ) = η = Xβ`. It follows that `μ = g⁻¹(η)`, where `g⁻¹` denotes the inverse of `g`. For linear regression, the canonical link function is the identity function. Therefore, when using the canonical link function, `μ = g⁻¹(η) = η`.
#
# In this model, the parameters that we want to estimate are `β` and `σ`.
# We need to select prior distributions for these parameters. For each `βᵢ`
# we choose a normal distribution with zero mean and variance `s²`. Here, `βᵢ`
# denotes the `i`th component of `β`. For `σ`, we will choose a half-normal
# distribution with variance `t²`. `s` and `t` are hyperparameters that we will
# need to choose.

# We define this model using the Soss probabilistic programming library:

m = @model X, s, t begin
    p = size(X, 2) # number of features
    β ~ Normal(0, s) |> iid(p) # coefficients
    σ ~ HalfNormal(t) # dispersion
    η = X * β # linear predictor
    μ = η # `μ = g⁻¹(η) = η`
    y ~ For(eachindex(μ)) do j
        Normal(μ[j], σ) # `Yᵢ ~ Normal(mean=μᵢ, variance=σ²)`
    end
end;

# Generate some synthetic features. Let us generate two continuous
# features and two binary categorical features.

num_rows = 1_000
x1 = randn(num_rows)
x2 = randn(num_rows)
x3 = Int.(rand(num_rows) .> 0.5)
x4 = Int.(rand(num_rows) .> 0.5)
X = (x1 = x1, x2 = x2, x3 = x3, x4 = x4)

# Define the hyperparameters of our prior distributions:

hyperparams = (s=0.1, t=0.1)

# Convert the Soss model into a `SossMLJModel`:

model = SossMLJModel(;
    model       = m,
    hyperparams = hyperparams,
    infer       = dynamicHMC,
    response    = :y,
);

# Generate some synthetic labels:

args = merge(model.transform(X), hyperparams)
truth = rand(m(args))

# Create an MLJ machine for fitting our model:

mach = MLJBase.machine(model, X, truth.y)

# Fit the machine. This may take several minutes.

fit!(mach)

# Construct the posterior distribution and the joint posterior predictive distribution:

predictor_joint = predict_joint(mach, X)
typeof(predictor_joint)

# Draw a single sample from the joint posterior predictive distribution:

single_sample = rand(predictor_joint; response = :y)

# Evaluate the logpdf of the joint posterior predictive distribution at this sample:

logpdf(predictor_joint, single_sample)

# True `β`:

truth.β

# Posterior distribution of `β`

predict_particles(mach, X; response = :β)

# Difference between the posterior distribution of `β` to the true values:

truth.β - predict_particles(mach, X; response = :β)

# Compare the joint posterior predictive distribution of `μ` to the true values:

truth.μ - predict_particles(mach, X; response = :μ)

# Compare the joint posterior predictive distribution of `y` to the true values:

truth.y - predict_particles(mach, X; response = :y)

# Construct each of the marginal posterior predictive distributions:

predictor_marginal = MLJBase.predict(mach, X)
typeof(predictor_marginal)

# `predictor_marginal` has one element for each row in `X`

size(predictor_marginal)

# Draw a single sample from each of the marginal posterior predictive distributions:

only.(rand.(predictor_marginal))

# Use cross-validation to evaluate the model with respect to the expected value of the root mean square error (RMSE)

evaluate!(mach, resampling=CV(; nfolds = 6, shuffle = true), measure=rms_expected, operation=predict_particles)

# Use cross-validation to evaluate the model with respect to the median of the root mean square error (RMSE)

evaluate!(mach, resampling=CV(; nfolds = 6, shuffle = true), measure=rms_median, operation=predict_particles)


text:
export Orthogonal
import LinearAlgebra: svd

doc"""
    Orthogonal([gain=1.0])

# References
* Saxe et al., [Exact solutions to the nonlinear dynamics of learning in deep linear neural networks](http://arxiv.org/abs/1312.6120)
"""
struct Orthogonal
    gain::Float64
end
Orthogonal() = Orthogonal(1.0)

function (o::Orthogonal)(::Type{T}, dim1::Int, dim2::Int) where T
    a = randn(T, dim1, dim2)
    F = svd(a)
    # u, _, v = svd(a)
    q = size(F.U) == (dim1,dim2) ? F.U : F.Vt
    q .* T(o.gain)
end


text:
adjacency_mat = [
    false true true true true true true true true false true true true true false false false true false true false true false false false false false false false false false true false false
    true false true true false false false true false false false false false true false false false true false true false true false false false false false false false false true false false false
    true true false true false false false true true true false false false true false false false false false false false false false false false false false true true false false false true false
    true true true false false false false true false false false false true true false false false false false false false false false false false false false false false false false false false false
    true false false false false false true false false false true false false false false false false false false false false false false false false false false false false false false false false false
    true false false false false false true false false false true false false false false false true false false false false false false false false false false false false false false false false false
    true false false false true true false false false false false false false false false false true false false false false false false false false false false false false false false false false false
    true true true true false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false
    true false true false false false false false false false false false false false false false false false false false false false false false false false false false false false true false true true
    false false true false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false true
    true false false false true true false false false false false false false false false false false false false false false false false false false false false false false false false false false false
    true false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false
    true false false true false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false
    true true true true false false false false false false false false false false false false false false false false false false false false false false false false false false false false false true
    false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false true true
    false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false true true
    false false false false false true true false false false false false false false false false false false false false false false false false false false false false false false false false false false
    true true false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false
    false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false true true
    true true false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false true
    false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false true true
    true true false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false
    false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false true true
    false false false false false false false false false false false false false false false false false false false false false false false false false true false true false true false false true true
    false false false false false false false false false false false false false false false false false false false false false false false false false true false true false false false true false false
    false false false false false false false false false false false false false false false false false false false false false false false true true false false false false false false true false false
    false false false false false false false false false false false false false false false false false false false false false false false false false false false false false true false false false true
    false false true false false false false false false false false false false false false false false false false false false false false true true false false false false false false false false true
    false false true false false false false false false false false false false false false false false false false false false false false false false false false false false false false true false true
    false false false false false false false false false false false false false false false false false false false false false false false true false false true false false false false false true true
    false true false false false false false false true false false false false false false false false false false false false false false false false false false false false false false false true true
    true false false false false false false false false false false false false false false false false false false false false false false false true true false false true false false false true true
    false false true false false false false false true false false false false false true true false false true false true false true true false false false false false true true true false true
    false false false false false false false false true true false false false true true true false false true true true false true true false false true true true true true true true false
   ]

max_stars_repo_name: alok/GCN

max_issues_repo_name: alok/GCN

max_forks_repo_name: alok/GCN


text:
#=
HDF5 TYPES SUPPORTED
signed and unsigned integers of 8, 16, 32, and 64 bits, Float32, Float64;
Complex versions of these numeric types
Arrays of these numeric types (including complex versions)

 ASCIIString and UTF8String; and Arrays of these two string types.

Note: can't really do String support for interoperability reasons.
=#
const HDF5Type = Union{ Float32, Float64,
                        Complex{Float32}, Complex{Float64},
                        UInt8, UInt16, UInt32, UInt64,
                        Complex{UInt8}, Complex{UInt16},
                        Complex{UInt32}, Complex{UInt64},
                        Int8, Int16, Int32, Int64,
                        Complex{Int8}, Complex{Int16},
                        Complex{Int32}, Complex{Int64} }

const HDF5Array = Union{Array{Float64, N},
                        Array{Float32, N},
                        Array{Complex{Float64}, N},
                        Array{Complex{Float32}, N},
                        Array{UInt8, N},
                        Array{UInt16, N},
                        Array{UInt32, N},
                        Array{UInt64, N},
                        Array{Complex{UInt8}, N},
                        Array{Complex{UInt16}, N},
                        Array{Complex{UInt32}, N},
                        Array{Complex{UInt64}, N},
                        Array{Int8, N},
                        Array{Int16, N},
                        Array{Int32, N},
                        Array{Int64, N},
                        Array{Complex{Int8}, N},
                        Array{Complex{Int16}, N},
                        Array{Complex{Int32}, N},
                        Array{Complex{Int64}, N}} where N

# These are adapted from unix2datetime.(1.0e-9.*[typemin(Int64), typemax(Int64)])
const unset_s = "1677-09-21T00:12:44"
const unset_t = "2262-04-11T23:47:16"


text:
"""
    MaternKernel(; ν::Real=1.5, metric=Euclidean())

Matérn kernel of order `ν` with respect to the `metric`.

# Definition

For inputs ``x, x'`` and metric ``d(\\cdot, \\cdot)``, the Matérn kernel of order
``\\nu > 0`` is defined as
```math
k(x,x';\\nu) = \\frac{2^{1-\\nu}}{\\Gamma(\\nu)}\\big(\\sqrt{2\\nu} d(x, x')\\big) K_\\nu\\big(\\sqrt{2\\nu} d(x, x')\\big),
```
where ``\\Gamma`` is the Gamma function and ``K_{\\nu}`` is the modified Bessel function of
the second kind of order ``\\nu``.
By default, ``d`` is the Euclidean metric ``d(x, x') = \\|x - x'\\|_2``.

A Gaussian process with a Matérn kernel is ``\\lceil \\nu \\rceil - 1``-times
differentiable in the mean-square sense.

!!! note

    Differentiation with respect to the order ν is not currently supported.

See also: [`Matern12Kernel`](@ref), [`Matern32Kernel`](@ref), [`Matern52Kernel`](@ref)
"""
struct MaternKernel{Tν<:Real,M} <: SimpleKernel
    ν::Vector{Tν}
    metric::M

    function MaternKernel(ν::Real, metric)
        @check_args(MaternKernel, ν, ν > zero(ν), "ν > 0")
        return new{typeof(ν),typeof(metric)}([ν], metric)
    end
end

MaternKernel(; nu::Real=1.5, ν::Real=nu, metric=Euclidean()) = MaternKernel(ν, metric)

@functor MaternKernel

@inline _get_ν(k::MaternKernel) = only(k.ν)
ChainRulesCore.@non_differentiable _get_ν(k)  # work-around; should be "NotImplemented" rather than NoTangent

@inline function kappa(k::MaternKernel, d::Real)
    result = _matern(_get_ν(k), d)
    return ifelse(iszero(d), one(result), result)
end

function _matern(ν::Real, d::Real)
    y = sqrt(2ν) * d
    return exp((one(d) - ν) * logtwo - loggamma(ν) + ν * log(y) + log(besselk(ν, y)))
end

metric(k::MaternKernel) = k.metric

function Base.show(io::IO, κ::MaternKernel)
    return print(io, "Matern Kernel (ν = ", only(κ.ν), ", metric = ", κ.metric, ")")
end

## Matern12Kernel = ExponentialKernel aliased in exponential.jl

"""
    Matern32Kernel(; metric=Euclidean())

Matérn kernel of order ``3/2`` with respect to the `metric`.

# Definition

For inputs ``x, x'`` and metric ``d(\\cdot, \\cdot)``, the Matérn kernel of order ``3/2`` is
 given by
```math
k(x, x') = \\big(1 + \\sqrt{3} d(x, x') \\big) \\exp\\big(- \\sqrt{3} d(x, x') \\big).
```
By default, ``d`` is the Euclidean metric ``d(x, x') = \\|x - x'\\|_2``.

See also: [`MaternKernel`](@ref)
"""
struct Matern32Kernel{M} <: SimpleKernel
    metric::M
end

Matern32Kernel(; metric=Euclidean()) = Matern32Kernel(metric)

kappa(::Matern32Kernel, d::Real) = (1 + sqrt(3) * d) * exp(-sqrt(3) * d)

metric(k::Matern32Kernel) = k.metric

function Base.show(io::IO, k::Matern32Kernel)
    return print(io, "Matern 3/2 Kernel (metric = ", k.metric, ")")
end

"""
    Matern52Kernel(; metric=Euclidean())

Matérn kernel of order ``5/2`` with respect to the `metric`.

# Definition

For inputs ``x, x'`` and metric ``d(\\cdot, \\cdot)``, the Matérn kernel of order ``5/2`` is
given by
```math
k(x, x') = \\bigg(1 + \\sqrt{5} d(x, x') + \\frac{5}{3} d(x, x')^2\\bigg)
           \\exp\\big(- \\sqrt{5} d(x, x') \\big).
```
By default, ``d`` is the Euclidean metric ``d(x, x') = \\|x - x'\\|_2``.

See also: [`MaternKernel`](@ref)
"""
struct Matern52Kernel{M} <: SimpleKernel
    metric::M
end

Matern52Kernel(; metric=Euclidean()) = Matern52Kernel(metric)

kappa(::Matern52Kernel, d::Real) = (1 + sqrt(5) * d + 5 * d^2 / 3) * exp(-sqrt(5) * d)

metric(k::Matern52Kernel) = k.metric

function Base.show(io::IO, k::Matern52Kernel)
    return print(io, "Matern 5/2 Kernel (metric = ", k.metric, ")")
end


text:
using FundamentalsNumericalComputation,LinearAlgebra,DifferentialEquations
using Test

@testset "Chapter 1" begin
	@test FNC.horner([-1,3,-3,1],1.6) ≈ 0.6^3
end

@testset "Chapter 2" begin
	A = [ 1 2 3 0; -1 1 2 -1; 3 1 2 4; 1 1 1 1 ]
	L,U = FNC.lufact(A)
	@test norm(L*U - A) < 100eps()
	@test norm(U - triu(U)) < 100eps()
	@test norm(L - tril(L)) < 100eps()
	b = [1,10,0,-1] / 5;
	@test norm(L\b - FNC.forwardsub(L,b)) < 100eps()
	@test norm(U\b - FNC.backsub(U,b)) < 100eps()
end

@testset "Chapter 3" begin
	A = [3 4 5;-1 0 1;4 2 0; 1 1 2; 3 -4 1]
	b = 5:-1:1
	@test FNC.lsnormal(A,b) ≈ A\b
	@test FNC.lsqrfact(A,b) ≈ A\b
	Q,R = qr(A)
	QQ,RR = FNC.qrfact(A)
	@test Q ≈ QQ
	@test R ≈ RR[1:3,:]
end

@testset "Chapter 4" begin

	for c = [2,4,7.5,11]
		f = x -> exp(x) - x - c;
		dfdx = x -> exp(x) - 1;
		x = FNC.newton(f,dfdx,1.0);  r = x[end];
		@test abs(f(r)) < 100eps()
	end

	for c = [2,4,7.5,11]
		f = x -> exp(x) - x - c;
		dfdx = x -> exp(x) - 1;
		x = FNC.secant(f,3,0.5);  r = x[end];
		@test abs(f(r)) < 100eps()
	end

	function nlfun(x)
		f = zeros(3)
		f[1] = exp(x[2]-x[1]) - 2;
		f[2] = x[1]*x[2] + x[3];
		f[3] = x[2]*x[3] + x[1]^2 - x[2];
		return f
	end
	function nljac(x)
		J = zeros(3,3)
		J[1,:] = [-exp(x[2]-x[1]),exp(x[2]-x[1]), 0]
		J[2,:] = [x[2], x[1], 1]
		J[3,:] = [2*x[1], x[3]-1, x[2]]
		return J
	end

	x = FNC.newtonsys(nlfun,nljac,[0,0,0]);
	@test norm(nlfun(x[end])) < 100eps()
	x = FNC.newtonsys(nlfun,nljac,[1,2,3]);
	@test norm(nlfun(x[end])) < 100eps()

	x = FNC.levenberg(nlfun,[10,-4,-3])
	@test norm(nlfun(x[end])) < 1e-12

end

@testset "Chapter 5" begin
	f = t->cos(5t)
	Q,t = FNC.intadapt(f,-1,3,1e-8)
	@test Q ≈ (sin(15)+sin(5))/5 rtol = 1e-5
	T,_ = FNC.trapezoid(f,-1,3,820)
	@test T ≈ (sin(15)+sin(5))/5 rtol = 1e-4

	t = [-2,-0.5,0,1,1.5,3.5,4]/10
	S = FNC.spinterp(t,exp.(t))
	@test S(0.33) ≈ exp(0.33) rtol = 1e-5
	w = FNC.fdweights(t.-0.12,2)
	f = x->cos(3x)
	@test dot(w,f.(t)) ≈ -9cos(0.36) rtol = 1e-3
	y = FNC.hatfun(0.22,t,5)
	@test y ≈ (0.22-t[5])/(t[6]-t[5])
	@test FNC.hatfun(0.6,t,5)==0
	p = FNC.plinterp(t,f.(t))
	@test p(0.22) ≈ f(t[5]) + (f(t[6])-f(t[5]))*(0.22-t[5])/(t[6]-t[5])
end

@testset "Chapter 6" begin
	f = (u,p,t) -> u + p*t^2
	û = exp(1.5) - 2*(-2 + 2*exp(1.5) - 2*1.5 - 1.5^2)
	ivp = ODEProblem(f,1,(0,1.5),-2)
	t,u = FNC.euler(ivp,4000)
	@test û ≈ u[end] rtol = 0.005
	t,u = FNC.am2(ivp,4000)
	@test û ≈ u[end] rtol = 0.005

	g = (u,p,t) -> [t+p-sin(u[2]),u[1]]
	ivp = ODEProblem(g,[-1.,4],(1.,2.),-6)
	sol = solve(ivp)
	t,u = FNC.euler(ivp,4000)
	@test u[end] ≈ sol.u[end] rtol=0.004
	t,u = FNC.ie2(ivp,4000)
	@test u[end] ≈ sol.u[end] rtol=0.0005
	t,u = FNC.rk4(ivp,800)
	@test u[end] ≈ sol.u[end] rtol=0.0005
	t,u = FNC.ab4(ivp,800)
	@test u[end] ≈ sol.u[end] rtol=0.0005
	t,u = FNC.rk23(ivp,1e-4)
	@test u[end] ≈ sol.u[end] rtol=0.0005
	t,u = FNC.am2(ivp,2000)
	@test u[end] ≈ sol.u[end] rtol=0.0005
end


text:
module GalacticMetaheuristics

using Reexport, GalacticOptim, GalacticOptim.SciMLBase
@reexport using Metaheuristics

function initial_population!(opt, prob, bounds, f)
    opt_init = deepcopy(opt)
    opt_init.options.iterations = 2
    Metaheuristics.optimize(f, bounds, opt_init)

    pop_size = opt_init.parameters.N
    population_rand = [bounds[1, :] + rand(length(prob.u0)) .* (bounds[2, :] - bounds[1, :]) for i in 1:(pop_size-1)]
    push!(population_rand, prob.u0)
    population_init = [Metaheuristics.create_child(x, f(x)) for x in population_rand]
    prev_status = Metaheuristics.State(Metaheuristics.get_best(population_init), population_init)
    opt.parameters.N = pop_size
    opt.status = prev_status
    return nothing
end

function __map_optimizer_args!(prob::OptimizationProblem, opt::Metaheuristics.AbstractAlgorithm;
    callback=nothing,
    maxiters::Union{Number,Nothing}=nothing,
    maxtime::Union{Number,Nothing}=nothing,
    abstol::Union{Number,Nothing}=nothing,
    reltol::Union{Number,Nothing}=nothing,
    kwargs...)

    for j in kwargs
        if j.first .∈ Ref(propertynames(Metaheuristics.Information()))
            error("Set $(j.first) by directly passing it to Information Structure which is passed to $(typeof(opt)) algorithms when calling solve().")
        elseif j.first .∈ Ref(propertynames(Metaheuristics.Options()))
            setproperty!(opt.options, j.first, j.second)
        else
            error("$(j.first) keyword is not a valid option for $(typeof(opt).super) algorithm.")
        end
    end

    if !isnothing(callback)
        @warn "Callback argument is currently not used by $(typeof(opt).super)"
    end

    if !isnothing(maxiters)
        opt.options.iterations = maxiters
    end

    if !isnothing(maxtime)
        opt.options.time_limit = maxtime
    end

    if !isnothing(abstol)
        opt.options.f_tol = abstol
    end

    if !isnothing(reltol)
        @warn "common reltol is currently not used by $(typeof(opt).super)"
    end
    return nothing
end

function SciMLBase.__solve(prob::OptimizationProblem, opt::Metaheuristics.AbstractAlgorithm;
    callback=nothing,
    maxiters::Union{Number,Nothing}=nothing,
    maxtime::Union{Number,Nothing}=nothing,
    abstol::Union{Number,Nothing}=nothing,
    reltol::Union{Number,Nothing}=nothing,
    progress=false,
    use_initial::Bool=false,
    kwargs...)

    local x

    maxiters = GalacticOptim._check_and_convert_maxiters(maxiters)
    maxtime = GalacticOptim._check_and_convert_maxtime(maxtime)

    _loss = function (θ)
        x = prob.f(θ, prob.p)
        return first(x)
    end

    if !isnothing(prob.lb) & !isnothing(prob.ub)
        opt_bounds = [prob.lb prob.ub]'
    else
        error("$(opt) requires lower and upper bounds to be defined.")
    end

    if !isnothing(prob.f.cons)
        @warn "Equality constraints are current not passed on by GalacticOptim"
    end

    if !isnothing(prob.lcons)
        @warn "Inequality constraints are current not passed on by GalacticOptim"
    end

    if !isnothing(prob.ucons)
        @warn "Inequality constraints are current not passed on by GalacticOptim"
    end

    __map_optimizer_args!(prob, opt, callback=callback, maxiters=maxiters, maxtime=maxtime, abstol=abstol, reltol=reltol; kwargs...)

    if use_initial
        initial_population!(opt, prob, opt_bounds, _loss)
    end

    t0 = time()
    opt_res = Metaheuristics.optimize(_loss, opt_bounds, opt)
    t1 = time()

    SciMLBase.build_solution(prob, opt, Metaheuristics.minimizer(opt_res), Metaheuristics.minimum(opt_res); original=opt_res)
end

end


text:
const stats_path = joinpath(dirname(@__DIR__),"stats")
const deltas_cols = [:arena_place,:arena_victories,:arena_defences,:total_power,:maxed_heroes]
const z_cols = [:clan_tokens,:raid_points,:portal_stones]

const curr_members = ["AngelaPapi", "Annanana",
    "Berserkresse","Breyn","D_Mac", "Displace",
    "DocStrange", "Dom","Drizzt4Ever", "Fisichella",
    "Foxgloves", "Glory", "Imper0s", "IrishWrath",
    "KARMA", "Kropman", "Lux", "Mbad", "Mikaeus",
    "Minimoose27", "MrMilanna", "Psykosis", "Rhellik",
    "TRUNG", "The Wiza", "Tristanslav", "WhoTFCars",
    "bnv", "matrixttd", "sillien", "skinwalker"]

const past_members = ["Rhellik", "Tristanslav",
    "Fryseboksen"]


text:
graph = ARGS[1]
idx = parse(Int, ARGS[2])
nx = parse(Int, ARGS[3])
ny = parse(Int, ARGS[4])
rep = parse(Int, ARGS[5])

scratch_dir = ARGS[6]

if isfile("$(scratch_dir)/res_$(graph)_$(rep).jld")
    println("the file already exists!")
    exit()
end

using KLIEPInference
using ProximalBase, CoordinateDescent
using LinearAlgebra, SparseArrays, Statistics, Random
using Distributions, StatsBase, JLD

include("SymKLIEP.jl")

println("importing parameters from $(graph).jld")
file = jldopen("graphs/$(graph).jld", "r")
γx = pack(read(file, "Θx"))
γy = pack(read(file, "Θy"))
close(file)

p = length(γx)

println("generating samples")
Random.seed!(123 + rep)
spl = IsingSampler(γx; thin=2000)
X = rand(spl, nx)
spl = IsingSampler(γy; thin=2000)
Y = rand(spl, ny)

Ψx = Ψising(X)
Ψy = Ψising(Y)

res = zeros(3, 2, 3, 5)

println("original procedure with ny > nx ...")
@show λ1 = exp.(range(log(sqrt(16. * log(p) / nx)), stop=log(sqrt(2. * log(p) / nx)), length=5))
λ2 = sqrt(2. * log(p) / ny)
for t = 1:5
    # step 1
    θ = spKLIEP(Ψx, Ψy, λ1[t], CD_KLIEP(); loadings=false)

    # step 2
    H = KLIEP_Hessian(θ, Ψy)
    ω = Hinv_row(H, idx, λ2)
    supp = KLIEPInference._find_supp(idx, ω)
    ω[supp] = view(H, supp, supp)\(supp .=== idx)

    # SparKLIE+1
    θ1 = KLIEPInference._debias1(Ψx, Ψy, θ, ω, idx)
    σ = stderr_SparKLIE(Ψx, Ψy, θ, ω)

    res[1, 1, 1, t] = θ1 - γx[idx] + γy[idx]
    res[1, 1, 2, t] = σ
    res[1, 1, 3, t] = (θ1 - γx[idx] + γy[idx]) / σ

    # SparKLIE+2
    θ2 = KLIEPInference._debias2(Ψx, Ψy, θ, ω, idx)

    supp = KLIEPInference._find_supp(idx, θ, ω)
    ω = KLIEP_Hessian(θ[supp], Ψy[supp, :])\(supp .=== idx)
    σ = stderr_SparKLIE(Ψx[supp, :], Ψy[supp, :], θ[supp], ω)

    res[1, 2, 1, t] = θ2 - γx[idx] + γy[idx]
    res[1, 2, 2, t] = σ
    res[1, 2, 3, t] = (θ2 - γx[idx] + γy[idx]) / σ
end
println("... done!")

println("reverse KL ...")
@show λ1 = exp.(range(log(sqrt(16. * log(p) / nx)), stop=log(sqrt(2. * log(p) / nx)), length=5))
λ2 = sqrt(2. * log(p) / nx)
for t = 1:5
    # step 1
    θ = spKLIEP(Ψy, Ψx, λ1[t], CD_KLIEP(); loadings=false)

    # step 2
    H = KLIEP_Hessian(θ, Ψx)
    ω = Hinv_row(H, idx, λ2)
    supp = KLIEPInference._find_supp(idx, ω)
    ω[supp] = view(H, supp, supp)\(supp .=== idx)

    # SparKLIE+1
    θ1 = KLIEPInference._debias1(Ψy, Ψx, θ, ω, idx)
    σ = stderr_SparKLIE(Ψy, Ψx, θ, ω)

    res[2, 1, 1, t] = -θ1 + γx[idx] - γy[idx]
    res[2, 1, 2, t] = σ
    res[2, 1, 3, t] = (-θ1 + γx[idx] - γy[idx]) / σ

    # SparKLIE+2
    θ2 = KLIEPInference._debias2(Ψy, Ψx, θ, ω, idx)

    supp = KLIEPInference._find_supp(idx, θ, ω)
    ω = KLIEP_Hessian(θ[supp], Ψx[supp, :])\(supp .=== idx)
    σ = stderr_SparKLIE(Ψy[supp, :], Ψx[supp, :], θ[supp], ω)

    res[2, 2, 1, t] = -θ2 + γx[idx] - γy[idx]
    res[2, 2, 2, t] = σ
    res[2, 2, 3, t] = (-θ2 + γx[idx] - γy[idx]) / σ
end
println("... done!")

println("symmetric KL ...")
@show λ1 = exp.(range(log(sqrt(16. * log(p) / nx)), stop=log(sqrt(2. * log(p) / nx)), length=5))
λ2 = .5 * (sqrt(2. * log(p) / nx) + sqrt(2. * log(p) / ny))
for t = 1:5
    # step 1
    θ = spSymKLIEP(Ψx, Ψy, λ1[t], CD_SymKLIEP())

    # step 2
    H = SymKLIEP_Hessian(θ, Ψx, Ψy)
    ω = Hinv_row(H, idx, λ2)
    supp = KLIEPInference._find_supp(idx, ω)
    ω[supp] = view(H, supp, supp)\(supp .=== idx)

    # SparKLIE+1
    θ1 = SymKLIEP_debias1(Ψx, Ψy, θ, ω, idx)
    σ = SymKLIEP_stderr(Ψx, Ψy, θ, ω)

    res[3, 1, 1, t] = θ1 - γx[idx] + γy[idx]
    res[3, 1, 2, t] = σ
    res[3, 1, 3, t] = (θ1 - γx[idx] + γy[idx]) / σ

    # SparKLIE+2
    θ2 = SymKLIEP_debias2(Ψx, Ψy, θ, ω, idx)

    supp = KLIEPInference._find_supp(idx, θ, ω)
    ω = SymKLIEP_Hessian(θ[supp], Ψx[supp, :], Ψy[supp, :])\(supp .=== idx)
    σ = SymKLIEP_stderr(Ψx[supp, :], Ψy[supp, :], θ[supp], ω)

    res[3, 2, 1, t] = θ2 - γx[idx] + γy[idx]
    res[3, 2, 2, t] = σ
    res[3, 2, 3, t] = (θ2 - γx[idx] + γy[idx]) / σ
end
println("... done!")

println("saving results to $(scratch_dir)/res_$(graph)_$(rep).jld ...")
@save "$(scratch_dir)/res_$(graph)_$(rep).jld" res
println("... done!")


text:
# Preamble ########################################################################################
using Pkg
Pkg.activate(Base.current_project())

cd("..") # Go to the parent directory
###################################################################################################
# Packages needed #################################################################################
using PolyGibbs
using DelimitedFiles
using ImageCore
using ImageIO
using ImageMagick
using FileIO
using PyPlot
###################################################################################################
# Angle of rotation ###############################################################################

const θ = π/4

###################################################################################################
# Retrieve of data from image #####################################################################

@info "Retrieving the image channels and extracting the two-dimensional fields"
input = Float64.(channelview(RGB.(load("data/origin/wmap/wmap50.png"))))

rchn = input[1, :, :]
gchn = input[2, :, :]
bchn = input[3, :, :]

###################################################################################################
# Identification of the RGB fields dimensions ################################################

j = (size(input)[2] - 1)/2;

N = Integer(2*j + 1);

@info "The size of the field is $N"

###################################################################################################
# Rotation of the image channels ###########################################################################

@info "Performing rotation of R channel at θ = $θ"
routput = Emodes(rchn, θ)

@info "Performing rotation of G channel at θ = $θ"
goutput = Emodes(gchn, θ)

@info "Performing rotation of B channel at θ = $θ"
boutput = Emodes(bchn, θ)

###################################################################################################
# Saving results ##################################################################################

@info "Writting the results as separate data files for each channel"
writedlm("data/target/wmap/wmap50R_rpi4", real(routput))
writedlm("data/target/wmap/wmap50G_rpi4", real(goutput))
writedlm("data/target/wmap/wmap50B_rpi4", real(boutput))

###################################################################################################

begin
    imgrep = permutedims(StackedView(transpose(real(routput)), transpose(real(goutput)), transpose(real(boutput))), [3, 2, 1])
	fig = figure()
	ax = gca()
	ax.imshow(imgrep, cmap = "gray")
	tight_layout()
	savefig("data/target/wmap/wmap50_rpi4", dpi = 300)
end

max_stars_repo_name: rurz/PolyGibbs

max_issues_repo_name: rurz/PolyGibbs

max_forks_repo_name: rurz/PolyGibbs


text:
# # Foundations of Bilevel Programming: Example 6
# This example is from the book Princeton Handbook of Test Problems in Local and Global Optimization
# Dempe, Chapter 9.3.2 -parg 221 [url](https://www.springer.com/gp/book/9780792358015)

# Here, only the second level is described

# Model of the problem
# First level
# ```math
# \min (x-5)^2+(2y+1)^2,\\
# \notag s.t.\\
# x \geq 0,\\
# y \geq 0,\\
# ```
# Second level
# ```math
# \min (y-1)^2-1.5xy,\\
# \notag s.t.\\
# -3x+y \leq -3,\\
# x-0.5y \leq 4,\\
# x+y \leq 7,\\
# ```

using BilevelJuMP
using Ipopt
using JuMP
using Test

model = BilevelModel(Ipopt.Optimizer, mode = BilevelJuMP.ProductMode(1e-9))

# Global variables
atol = 1e-3

# First we need to create all of the variables in the upper and lower problems:

# Upper level variables
@variable(Upper(model), x)

#Lower level variables
@variable(Lower(model), y)

# Then we can add the objective and constraints of the upper problem:

# Upper level objecive function
@objective(Upper(model), Min, (x-5)^2 + (2y+1)^2)

# Upper level constraints
@constraint(Upper(model), x >= 0)
@constraint(Upper(model), y >= 0) # only in lowrrin GAMS

# Followed by the objective and constraints of the lower problem:

# Lower objective function
@objective(Lower(model), Min, (y-1)^2 -1.5*x*y)

# Lower constraints
@constraint(Lower(model), -3x +    y <= -3)
@constraint(Lower(model),   x - 0.5y <= 4)
@constraint(Lower(model),   x +    y <= 7)

# Initial Starting conditions  #src

# Now we can solve the problem and verify the solution again that reported by
# Dempe.

optimize!(model)
primal_status(model)
termination_status(model)

@test value(x) ≈ 1 atol=atol
@test value(y) ≈ 0 atol=atol


text:
fixprob(p::Float, eps_val::Float = eps(Float)) = min(max(p, eps_val), Float(1)-eps_val)

fixnan(v::Float) = isnan(v) ? Float(0) : v

fixdiv(x::Float) = x == 0 ? Float(1) : x

logprob(p::Float, eps_val::Float = eps(Float)) = log(fixprob(p, eps_val))

tau(min_samps, min_weight::Float) = (Float(1) - exp(log(Float(1)-min_weight)/min_samps))

function logsumexp{N}(arr::Array{Float, N})
  add_val = maximum(arr, 1)
  out = log.(sum(exp.(arr .- add_val), 1)) .+ add_val
  return out
end

function probmat{N}(arr::Array{Float, N})
  p = exp.(arr)
  p ./= sum(p,1)
  return p
end

function set!(rr::RemoteChannel, val)
  if isready(rr)
    take!(rr)
  end
  put!(rr, val)
end

function fetchcall(rr_func::RemoteChannel, args...)
  return fetch(rr_func)(args...)
end


text:
using LOLTools.SummonerV4

api_key = get(ENV, "RIOT_TOKEN", "")
platform = "kr"
summonerdto = SummonerV4.summoner_by_name(api_key, platform, "오뚜기3분미트볼")
@info SummonerV4.summoner_by_id(api_key, platform, summonerdto.id) == summonerdto

using LOLTools.ChampionMasteryV4
masteries = ChampionMasteryV4.by_summoner_id(api_key, platform, summonerdto.id)
@info length(masteries)
for mastery in masteries[1:3]
    @info mastery
end