text:
using GLPlot, GLVisualize, GLAbstraction, Colors, GeometryTypes, Plots, FileIO
using Reactive, GLWindow
import Quaternions: qrotation

# create a mandelbulb volume
N = 400
dims = (N, N, N)

function mandelbulb{T}(x0::T,y0::T,z0::T, n, iter)
    x,y,z = x0,y0,z0
    for i=1:iter
        r = sqrt(x*x + y*y + z*z)
        theta = atan2(sqrt(x*x + y*y) , z)
        phi = atan2(y,x)
        rn = r^n
        x1 = rn * sin(theta*n) * cos(phi*n) + x0
        y1 = rn * sin(theta*n) * sin(phi*n) + y0
        z1 = rn * cos(theta*n) + z0
        (x1 * x1 + y1 * y1 + z1 * z1) > n && return T(i)
        x, y, z = x1,y1,z1
    end
    T(iter)
end

x, y, z = ntuple(3) do i
    # linearly spaced array (not dense) from -1 to 1
    reshape(linspace(-1f0, 1f0, dims[i]), ntuple(j-> j == i ? dims[i] : 1, 3))
end
# change value for different results!
volume = mandelbulb.(x, y, z, 3f0, 22)
# create two sliders to interact with the 2 parameters of the mandelbulb function
maxi = maximum(volume)
const cmap = RGBA{Float32}.(colormap("RdBu", 100))
vol = map(volume) do level
    idx = floor(Int, ((level / maxi) * (length(cmap) - 1))) + 1
    cmap[idx]
end

# prepare the slices

slices = ntuple(3) do i
    range_s = play_widget(1:size(vol, i))
    slice = map(range_s) do j
        vol[ifelse(i == 1, j, :), ifelse(i == 2, j, :) , ifelse(i == 3, j, :)]
    end
    slice, range_s
end;

rotations = (
    Mat4f0(qrotation(unit(Vec3f0, 1), Float32(pi/2))),
    rotationmatrix_x(Float32(pi/2)) * Mat4f0(qrotation(unit(Vec3f0, 2), Float32(pi/2))),
    rotationmatrix_z(Float32(pi/2))
)
const ranges = (x, y, z)

plots = map(1:3) do i
    rot = rotations[i]
    slice, slider = slices[i]
    axis = unit(Vec3f0, 1)
    model = map(slider) do r
        trans = unit(Vec3f0, mod1(3-i, 3)) * ranges[i][r]
        translationmatrix(trans) * rot
    end
    rs = ranges[mod1(i + 1, 3)], ranges[mod1(i + 2, 3)]
    mini, maxi = extrema(ranges[i])
    glplot(
        slice,
        model = model,
        primitive = SimpleRectangle(-1, -1, 2, 2), # TODO, this is an inconsistency in the API and should also just accept ranges kw_arg
        preferred_camera = :perspective,
    )
end


text:
module Foo
    function parse(x)
        println(x)
    end
    include("parse.jl")
end

ext: jl


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct FlowLogPropertiesFormat <: SwaggerModel
    targetResourceId::Any # spec type: Union{ Nothing, String } # spec name: targetResourceId
    targetResourceGuid::Any # spec type: Union{ Nothing, String } # spec name: targetResourceGuid
    storageId::Any # spec type: Union{ Nothing, String } # spec name: storageId
    enabled::Any # spec type: Union{ Nothing, Bool } # spec name: enabled
    retentionPolicy::Any # spec type: Union{ Nothing, RetentionPolicyParameters } # spec name: retentionPolicy
    format::Any # spec type: Union{ Nothing, FlowLogFormatParameters } # spec name: format
    flowAnalyticsConfiguration::Any # spec type: Union{ Nothing, TrafficAnalyticsProperties } # spec name: flowAnalyticsConfiguration
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function FlowLogPropertiesFormat(;targetResourceId=nothing, targetResourceGuid=nothing, storageId=nothing, enabled=nothing, retentionPolicy=nothing, format=nothing, flowAnalyticsConfiguration=nothing, provisioningState=nothing)
        o = new()
        validate_property(FlowLogPropertiesFormat, Symbol("targetResourceId"), targetResourceId)
        setfield!(o, Symbol("targetResourceId"), targetResourceId)
        validate_property(FlowLogPropertiesFormat, Symbol("targetResourceGuid"), targetResourceGuid)
        setfield!(o, Symbol("targetResourceGuid"), targetResourceGuid)
        validate_property(FlowLogPropertiesFormat, Symbol("storageId"), storageId)
        setfield!(o, Symbol("storageId"), storageId)
        validate_property(FlowLogPropertiesFormat, Symbol("enabled"), enabled)
        setfield!(o, Symbol("enabled"), enabled)
        validate_property(FlowLogPropertiesFormat, Symbol("retentionPolicy"), retentionPolicy)
        setfield!(o, Symbol("retentionPolicy"), retentionPolicy)
        validate_property(FlowLogPropertiesFormat, Symbol("format"), format)
        setfield!(o, Symbol("format"), format)
        validate_property(FlowLogPropertiesFormat, Symbol("flowAnalyticsConfiguration"), flowAnalyticsConfiguration)
        setfield!(o, Symbol("flowAnalyticsConfiguration"), flowAnalyticsConfiguration)
        validate_property(FlowLogPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type FlowLogPropertiesFormat

const _property_map_FlowLogPropertiesFormat = Dict{Symbol,Symbol}(Symbol("targetResourceId")=>Symbol("targetResourceId"), Symbol("targetResourceGuid")=>Symbol("targetResourceGuid"), Symbol("storageId")=>Symbol("storageId"), Symbol("enabled")=>Symbol("enabled"), Symbol("retentionPolicy")=>Symbol("retentionPolicy"), Symbol("format")=>Symbol("format"), Symbol("flowAnalyticsConfiguration")=>Symbol("flowAnalyticsConfiguration"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_FlowLogPropertiesFormat = Dict{Symbol,String}(Symbol("targetResourceId")=>"String", Symbol("targetResourceGuid")=>"String", Symbol("storageId")=>"String", Symbol("enabled")=>"Bool", Symbol("retentionPolicy")=>"RetentionPolicyParameters", Symbol("format")=>"FlowLogFormatParameters", Symbol("flowAnalyticsConfiguration")=>"TrafficAnalyticsProperties", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ FlowLogPropertiesFormat }) = collect(keys(_property_map_FlowLogPropertiesFormat))
Swagger.property_type(::Type{ FlowLogPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_FlowLogPropertiesFormat[name]))}
Swagger.field_name(::Type{ FlowLogPropertiesFormat }, property_name::Symbol) =  _property_map_FlowLogPropertiesFormat[property_name]

function check_required(o::FlowLogPropertiesFormat)
    (getproperty(o, Symbol("targetResourceId")) === nothing) && (return false)
    (getproperty(o, Symbol("storageId")) === nothing) && (return false)
    true
end

function validate_property(::Type{ FlowLogPropertiesFormat }, name::Symbol, val)
end


text:
module Gumbo

using Compat

if isfile(joinpath(dirname(dirname(@__FILE__)),"deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("Gumbo not properly installed. Please run Pkg.build(\"Gumbo\")")
end

include("CGumbo.jl")

export HTMLElement,
       HTMLDocument,
       HTMLText,
       NullNode,
       HTMLNode,
       attrs,
       text,
       tag,
       children,
       getattr,
       setattr!,
       parsehtml,
       postorder,
       preorder,
       breadthfirst,
       prettyprint

include("htmltypes.jl")
include("manipulation.jl")
include("comparison.jl")
include("io.jl")
include("conversion.jl")

end


text:
# set up cores
using Distributed
addprocs(1)

@everywhere begin

    # READ IN CSV - i'll aproximate it with a vector
    # what I want to have here is vec_to_pass be read in as data from the local
    # machine that then is available to all workers
    vec_to_pass = [1,2,3,4,5]

    path = pwd()
    path = path*"/src/reprexs/def_functions.jl"
    include(path)

    # variable I want to be the same across the different cores
    a = 4
    b = 3
end

@distributed (+) for i in 1:10

    add_2!(a, b, vec_to_pass[i], i)
end

max_stars_repo_name: colebrookson/trait-based-rewiring

max_issues_repo_name: colebrookson/trait-based-rewiring

max_forks_repo_name: colebrookson/trait-based-rewiring


text:
using Plots
using JLD2

struct SimParam
    samplingrate::Float64
end

function sin_generator(freq, t)
    sin(2*pi*t*freq)
end

function main()
    simparam = SimParam(100*1000*1000)
    mod_ratio = 0.8

    f_t = open("time.dat","w")
    f_w_s = open("w_s.dat", "w");
    f_w_c = open("w_c.dat", "w");
    f_w_m = open("w_m.dat", "w");

    for step in 1:100*1000*1000
        t = step/simparam.samplingrate
        w_s = sin_generator(100, t)
        #w_s = 0
        w_c = sin_generator(954*1000, t)
        w_m = (1+mod_ratio*w_s)*w_c
        write(f_t, t)
        write(f_w_c, w_c)
        write(f_w_s, w_s)
        write(f_w_m, w_m)
    end
    close(f_w_m)
    close(f_w_c)
    close(f_w_s)
    close(f_t)

    times = []
    w_m_vals = []
    open("time.dat","r") do f_t
        open("w_m.dat","r") do f_w_m
            for step in 1:10*1000*1000
                t = read(f_t, Float64)
                w_m = read(f_w_m, Float64)
                push!(times, t)
                push!(w_m_vals, w_m)
            end
        end
    end

    plot_graph=plot(times, w_m_vals, fmt =:png)
    #plot!(plot_graph,times,w_c_vals)
    #plot!(plot_graph,times,w_m_vals)

    #gui(plot_graph)

    savefig(plot_graph,"signal.png")

    jldopen(ARGS[1],"w") do file
        file["samplingrate"] = simparam.samplingrate

        vals = []
        open("time.dat","r") do f
            while !eof(f_t)
                push!(vals, read(f, Float64))
            end
        end
        file["times"] = vals

        vals = []
        open("w_s.dat","r") do f
            while !eof(f)
                push!(vals, read(f, Float64))
            end
        end
        file["w_s"] = vals

        vals = []
        open("w_c.dat","r") do f
            while !eof(f)
                push!(vals, read(f, Float64))
            end
        end
        file["w_c"] = vals

        vals = []
        open("w_m.dat","r") do f
            while !eof(f)
                push!(vals, read(f, Float64))
            end
        end
        file["w_m"] = vals

    end

    rm("time.dat")
    rm("w_s.dat");
    rm("w_c.dat");
    rm("w_m.dat");
end

main()

max_stars_repo_name: seisyuu-hantatsushi/julia

max_issues_repo_name: seisyuu-hantatsushi/julia

max_forks_repo_name: seisyuu-hantatsushi/julia


text:
unset = Unset()

struct FuncArg
    name::Any
    type::Any
    default::Any
end

struct FuncHeader
    name::Any
    args::Any
    kwargs::Any
    ret::Any
    fresh::Any
end

FuncHeader() = FuncHeader(unset, unset, unset, unset, unset)

is_func_header(a::FuncHeader) = a.args != unset

function func_arg(@nospecialize(ex))::FuncArg
    @switch ex begin
        @case :(::$ty)
        @with func_arg(gensym("_")).type = ty
        @case :($var::$ty)
        @with func_arg(var).type = ty
        @case Expr(:kw, var, default)
        @with func_arg(var).default = default
        @case Expr(:(=), var, default)
        @with func_arg(var).default = default
        @case var::Symbol
        FuncArg(var, unset, unset)
        @case Expr(:..., _)
        error(
            "GG does not support variadic argument($ex) so far.\n" *
            "Try\n" *
            "  f(x...) = _f(x)\n" *
            "  @gg _f(x) = ...\n" *
            "See more at: https://github.com/JuliaStaging/GeneralizedGenerated.jl/issues/38",
        )
        @case _
        error("GG does not understand the argument $ex.")
    end
end

function func_header(@nospecialize(ex))::FuncHeader
    @switch ex begin
        @case :($hd::$ret)
        @with func_header(hd).ret = ret

        @case :($f($(args...); $(kwargs...)))
        inter = @with func_header(f).args = map(func_arg, args)
        @with inter.kwargs = map(func_arg, kwargs)

        @case :($f($(args...)))
        @with func_header(f).args = map(func_arg, args)

        @case :($f where {$(args...)})
        @with func_header(f).fresh = args

        @case Expr(:tuple, Expr(:parameters, kwargs...), args...)
        inter = @with FuncHeader().args = map(func_arg, args)
        @with inter.kwargs = map(func_arg, kwargs)

        @case Expr(:tuple, args...)
        @with FuncHeader().args = map(func_arg, args)

        @case f
        @with FuncHeader().name = f
    end
end

function of_args(::Unset)
    Argument[]
end

function of_args(args::AbstractArray{FuncArg})
    ret = Argument[]
    for (i, each) in enumerate(args)
        name = each.name === unset ? gensym("_$i") : each.name
        type = each.type === unset ? nothing : each.type
        arg = Argument(name, type, each.default)
        push!(ret, arg)
    end
    ret
end

extract_tvar(var::Union{Symbol,Expr})::Symbol = @match var begin
    :($a <: $_) => a
    :($a >: $_) => a
    :($_ >: $a >: $_) => a
    :($_ <: $a <: $_) => a
    a::Symbol => a
end


text:
include("temperature_render.jl")
using Dates
using Plots

size_multipliers = 1:20
times = []

for mult in size_multipliers
    start=now()
    temperature_render(;size_multiplier=mult)
    push!(times, now()-start)
end

plt = plot(size_multipliers, times, xlabel="Size Multiplier", ylabel="Code Execution Time")
display(plt)


text:
#BLAS.set_num_threads(1) # this one should have worked but did not
ENV["OPENBLAS_NUM_THREADS"] = 1 # prevents thrashing by pmap + blas

# travis-ci limits to 4 or so
const nw = 4# Sys.CPU_CORES-2
if nworkers() < nw
    if nworkers() > 1
        @info("Removing existing parallel workers for tests...")
        rmprocs(workers())
    end
    @info("Starting $nw parallel workers for tests...")
    addprocs(nw)
    @info("$(nworkers()) parallel workers started")
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 33bd93ec-77d9-11eb-2da8-87271f8b30c7
begin
  using Pkg
  Pkg.activate("MLJ_env", shared=true)

  # using MLJ
  using CSV
  using DataFrames
  using PlutoUI
  using Test
  using Printf
  # using Plots
  using LinearAlgebra
end

# ╔═╡ 6d69b148-77d9-11eb-04f5-034603897c31
begin
  using JSON

  important_words = open("../../ML_UW_Spec/C03/data/important_words.json", "r") do f
      JSON.parse(f)
  end;

  length(important_words)
end

# ╔═╡ 0f6bb770-77df-11eb-0f1f-4dddbd7e2498
include("./text_proc.jl");

# ╔═╡ 1c52f3b6-78d1-11eb-3077-bf6c471251c8
include("./utils.jl")

# ╔═╡ baa8ca1c-77d8-11eb-2e26-3b2445acbfa1
md"""
## Logistic Regression with L2 regularization

The goal of this second notebook is to implement your own logistic regression classifier with L2 regularization. We will do the following:

 - Extract features from Amazon product reviews.
 - Turn a DataFrame into a Julia array.
 - Write a function to compute the derivative of log likelihood function with an L2 penalty with respect to a single coefficient.
 - Implement gradient ascent with an L2 penalty.
 - Empirically explore how the L2 penalty can ameliorate overfitting.

"""

# ╔═╡ 3ad4770e-77d9-11eb-3a93-addf3e427025
md"""
### Load review dataset

For this assignment, we will use a subset of the Amazon product review dataset. The subset was chosen to contain similar numbers of positive and negative reviews, as the original dataset consisted primarily of positive reviews.
"""

# ╔═╡ 6e264178-77d9-11eb-33c1-19462dde66e2
begin
  products = train = CSV.File("../../ML_UW_Spec/C03/data/amazon_baby_subset.csv";
    header=true) |> DataFrame;

  first(products, 3)
end

# ╔═╡ 75f2e820-77de-11eb-19eb-a5a5977b5257
size(products)

# ╔═╡ 36f2c616-77de-11eb-03b4-3b92a532949a
## how many missing values for the column :review?
sum(ismissing.(products.review))

# ╔═╡ 6d4ef592-77d9-11eb-3844-1bea7dc71094
## in-place replace(ment) of missing values for review  with empty string
replace!(products.review, missing =>"");

# ╔═╡ d7b9c812-77de-11eb-0ebb-6bb75de9725b
## no more missing values for column :review
@test sum(ismissing.(products.review)) == 0

# ╔═╡ 6e09ca34-77d9-11eb-1a99-6578f9141004
md"""
One column of this dataset is 'sentiment', corresponding to the class label with +1 indicating a review with positive sentiment and -1 indicating one with negative sentiment.
"""

# ╔═╡ 6ded8978-77d9-11eb-3d94-e7507a85197e
products.sentiment[1:20]

# ╔═╡ 6dd12c56-77d9-11eb-0350-1dfc08d17dd6
md"""
Let us quickly explore more of this dataset. The 'name' column indicates the name of the product. Here we list the first 10 products in the dataset. We then count the number of positive and negative reviews.
"""

# ╔═╡ 6db87918-77d9-11eb-2082-edecdf8b7576
first(products, 10).name

# ╔═╡ 6d316f86-77d9-11eb-107b-1bf5708655d5
md"""
Now, we will perform 2 simple data transformations:

  1. Remove punctuation using Python's built-in string functionality.
  2. Compute word counts (only for important_words)

"""

# ╔═╡ 734ed4d4-77df-11eb-3dab-e5c10b4cda62
md"""
###### Step 1.
"""

# ╔═╡ 0f528d0e-77df-11eb-1415-2fb0e9d6792d
begin
	products[!, :review_clean] = remove_punctuation.(products[:, :review]);
	select!(products, Not(:review));
end

# ╔═╡ 0f362fec-77df-11eb-02a6-210a58c2839b
names(products)

# ╔═╡ 0f1c6b16-77df-11eb-2197-93c4a82cc135
md"""
###### Step 2.

For each word in `important_words`, we compute a count for the number of times the word occurs in the review. We will store this count in a separate column (one for each word). The result of this feature processing is a single column for each word in `important_words` which keeps a count of the number of times the respective word occurs in the review text.

"""

# ╔═╡ 0057a50e-77e0-11eb-29d6-85c2dac67743
wc(text::String, word::String) = length(findall(word, text))

# ╔═╡ 51a94e8e-77e3-11eb-04c8-d1389b1f5856
function addcols!(df, words)
  for word ∈ words
    df[!, word] = wc.(df.review_clean, word)
  end
end

# ╔═╡ 0f03e474-77df-11eb-2791-f165fd011737
@time addcols!(products, important_words)

# ╔═╡ 0ee68550-77df-11eb-0a6e-c5bf5e3fd9ac
length(names(products))

# ╔═╡ 0ecc0054-77df-11eb-21a4-3f677548ecf9
products.perfect, length(products.perfect)

# ╔═╡ 0eb50d9a-77df-11eb-0b7d-557ab3ea2c81
md"""
Now, write some code to compute the number of product reviews that contain the word perfect.
"""

# ╔═╡ d39040fe-78d0-11eb-236d-15b506adaffd
md"""
### Train-Validation split

We split the data into a train-validation split with 80% of the data in the training set and 20% of the data in the validation set.

**Note:** In previous assignments, we have called this a **train-test split**. However, the portion of data that we don't train on will be used to help **select model parameters**. Thus, this portion of data should be called a **validation set**. Recall that examining performance of various potential models (i.e. models with different parameters) should be on a validation set, while evaluation of selected model should always be on a test set.
"""

# ╔═╡ 005dc370-78d1-11eb-0fc0-039f1017af8e
train_data, validation_data = train_test_split(products; split=0.8, seed=42);

# ╔═╡ 004d22cc-78d1-11eb-10ce-09eb125a96b3
size(train_data), size(validation_data)

# ╔═╡ 682a78ba-78d1-11eb-1706-15ca470f7cf2
md"""
### Using Julia Arrays
"""

# ╔═╡ 0e383a04-77df-11eb-2ec9-c5aeffc4b131
function get_data(df::DF, features, output) where {DF <: Union{DataFrame, SubDataFrame}}
  s_features = [Symbol(f) for f ∈ features]
  df[:, :intercept] .= 1.0
  s_features = [:intercept, s_features...]
  X_matrix = convert(Matrix, select(df, s_features))
  y = df[!, output]

  (X_matrix, Vector{eltype(y)}(y))
end

# ╔═╡ 7ad4c752-78d1-11eb-267b-236c60163a33
begin
	# need a copy to add intercept, train_data is actually a view
	feature_matrix_train, sentiment_train = get_data(copy(train_data),
		important_words, :sentiment);
	# same remark as above concerning validation_data
	feature_matrix_valid, sentiment_valid = get_data(copy(validation_data),
		important_words, :sentiment);
end

# ╔═╡ 9f29e2b8-78d1-11eb-0389-5978263881b6
size(feature_matrix_train), size(sentiment_train)

# ╔═╡ b093b9ac-78d1-11eb-3a2e-ffb6e1102b1f
size(feature_matrix_valid), size(sentiment_valid)

# ╔═╡ fd23d180-77e6-11eb-2091-af1e84d8fbcf
md"""
#### Building on logistic regression with no L2 penalty assignment

Recall from lecture that the link function is given by:

$$P(y_i = +1 | \mathbf{x}_i,\mathbf{w}) = \frac{1}{1 + \exp(-\mathbf{w}^T h(\mathbf{x}_i))},$$

where the feature vector $h(\mathbf{x}_i)$ represents the word counts of `important_words` in the review  $\mathbf{x}_i$.

Implement the `predict_probability` function that implements the link function:
"""

# ╔═╡ fd120ad4-77e6-11eb-15ea-bb7d7ea35587
function predict_probability(feature_matrix::Matrix, coeffs)
  1. ./ (1. .+ exp.(-(feature_matrix * coeffs)))
end

# ╔═╡ fcf1dd2e-77e6-11eb-0619-e7fd2a622877
md"""
### Adding  L2 penalty

Let us now work on extending logistic regression with L2 regularization. As discussed in the lectures, the L2 regularization is particularly useful in preventing overfitting. In this assignment, we will explore L2 regularization in detail.

Recall from lecture and the previous assignment that for logistic regression without an L2 penalty, the derivative of the log likelihood function is:

$$\frac{\partial\ell}{\partial w_j} = \sum_{i=1}^N h_j(\mathbf{x}_i)\left(\mathbf{1}[y_i = +1] - P(y_i = +1 | \mathbf{x}_i, \mathbf{w})\right)$$

**Adding L2 penalty to the derivative**

It takes only a small modification to add a L2 penalty. All terms indicated in **red** refer to terms that were added due to an **L2 penalty**.

  - Recall from the lecture that the link function is still the sigmoid:
$$P(y_i = +1 | \mathbf{x}_i,\mathbf{w}) = \frac{1}{1 + \exp(-\mathbf{w}^T h(\mathbf{x}_i))},$$

  - We add the L2 penalty term to the per-coefficient derivative of log likelihood:
$$\frac{\partial\ell}{\partial w_j} = \sum_{i=1}^N h_j(\mathbf{x}_i)\left(\mathbf{1}[y_i = +1] - P(y_i = +1 | \mathbf{x}_i, \mathbf{w})\right) \color{red}{-2\lambda w_j }$$

The **per-coefficient derivative for logistic regression with an L2 penalty** is as follows:

$$\frac{\partial\ell}{\partial w_j} = \sum_{i=1}^N h_j(\mathbf{x}_i)\left(\mathbf{1}[y_i = +1] - P(y_i = +1 | \mathbf{x}_i, \mathbf{w})\right) \color{red}{-2\lambda w_j }$$

and for the intercept term, we have
$$\frac{\partial\ell}{\partial w_0} = \sum_{i=1}^N h_0(\mathbf{x}_i)\left(\mathbf{1}[y_i = +1] - P(y_i = +1 | \mathbf{x}_i, \mathbf{w})\right)$$

"""

# ╔═╡ d6f4542a-78d2-11eb-0209-6533fe570748
md"""
Note: As we did in the Regression course, we do not apply the L2 penalty on the intercept. A large intercept does not necessarily indicate overfitting because the intercept is not associated with any particular feature.

Write a function that computes the derivative of log likelihood with respect to a single coefficient $w_j$. Unlike its counterpart in the last assignment, the function accepts five arguments:

  - `errors` vector containing $(\mathbf{1}[y_i = +1] - P(y_i = +1 | \mathbf{x}_i, \mathbf{w}))$ for all $i$
  - `feature` vector containing $h_j(\mathbf{x}_i)$  for all $i$
  -  `coefficient` containing the current value of coefficient $w_j$.
  -  `l2_penalty` representing the L2 penalty constant $\lambda$
  -  `feature_is_constant` telling whether the $j$-th feature is constant or not.
"""

# ╔═╡ d6976230-78d2-11eb-1fb7-dd1259c9b9c8
function feature_derivative_with_l2(errors, feature, coeff, l2_penalty;
		is_intercept=false)
    deriv = dot(errors, feature)
    ## add L2 penalty term for any feature that isn't the intercept.
    !is_intercept && (deriv -= 2 * l2_penalty * coeff)
	deriv
end

# ╔═╡ 4870ab12-78d3-11eb-12f9-1f5e36d0706d
md"""
**Quiz Question:** In the code above, was the intercept term regularized?
 - No (because of flag `feature_is_constant`)
"""

# ╔═╡ 665516c2-78d3-11eb-13d1-8dafac6a1d6a
md"""
To verify the correctness of the gradient ascent algorithm, we provide a function for computing log likelihood (which we recall from the last assignment was a topic detailed in an advanced optional video, and used here for its numerical stability).

$$\ell\ell(\mathbf{w}) = \sum_{i=1}^N \Big( (\mathbf{1}[y_i = +1] - 1)\mathbf{w}^T h(\mathbf{x}_i) - \ln\left(1 + \exp(-\mathbf{w}^T h(\mathbf{x}_i))\right) \Big) \color{red}{-\lambda\|\mathbf{w}\|_2^2}$$
"""

# ╔═╡ 70a490b2-78d3-11eb-3acd-9bc21f53815e
function log_likelihood_with_l2(f_matrix, sentiment, coeffs, l2_penalty)
	indic = sentiment .== 1
 	scores = f_matrix * coeffs
  	log_exp = log.(1. .+ exp.(-scores))
 	# mask = isinf.(log_exp)
  	# log_exp[mask] = -scores[mask]
	sum((indic .- 1) .* scores .- log_exp) - l2_penalty * sum(coeffs[2:end] .^ 2)
end

# ╔═╡ 0cea789c-78d4-11eb-36c4-c794c5e1a816
md"""
**Quiz Question:** Does the term with L2 regularization increase or decrease $\ell\ell(\mathbf{w})$?
  - should decrease
"""

# ╔═╡ 1a7b28f8-78d4-11eb-3906-5ba3566b9ae9
md"""
The `logistic_regression` function looks almost like the one in the last assignment, with a minor modification to account for the L2 penalty. Fill in the code below to complete this modification.
"""

# ╔═╡ 970fba24-77ed-11eb-1a3e-e769d6049078
function logistic_regression_with_l2(f_matrix::Matrix, sentiment::Vector,
		init_coeffs, l2_penalty;
		η=1e-7, max_iter=200)
  coeffs = copy(init_coeffs)

  for ix ∈ 1:max_iter
    preds = predict_probability(f_matrix, coeffs)
    indic = sentiment .== 1
    errors = indic .- preds

    for jx ∈ 1:length(coeffs)
	  is_intercept = (jx == 1)
      ∂c = feature_derivative_with_l2(errors, f_matrix[:, jx], coeffs[jx],
				l2_penalty; is_intercept)
      coeffs[jx] += η * ∂c
    end

    if ix ≤ 15 || (ix ≤ 100 && ix % 10 == 0) || (ix ≤ 1000 && ix % 100 == 0) ||
      (ix ≤ 10000 && ix % 1000 == 0) || ix % 10000 == 0
      lp = log_likelihood_with_l2(f_matrix, sentiment, coeffs, l2_penalty)
      @printf("iter: %d / %d: log likelihood of observed labels: %.8f\n",
        ceil(Int, log10(max_iter)), ix, lp)
    end
  end

  coeffs
end

# ╔═╡ 972e0c72-77ed-11eb-0ae7-c127694faf89
md"""
### Explore effects of L2 regularization

Now that we have written up all the pieces needed for regularized logistic regression, let's explore the benefits of using **L2 regularization** in analyzing sentiment for product reviews. **As iterations pass, the log likelihood should increase**.

Below, we train models with increasing amounts of regularization, starting with no L2 penalty, which is equivalent to our previous logistic regression implementation.
"""

# ╔═╡ afa2bff4-78d4-11eb-10cc-0deed0adb371
function run_model(;matrix_tr=feature_matrix_train,
		s_train=sentiment_train,
		init_coeffs=zeros(Float64, (194, 1)),
		l2_penalty=0, η=5e-6, max_iter=501)

	logistic_regression_with_l2(matrix_tr, s_train, init_coeffs, l2_penalty; η,
		max_iter)
end

# ╔═╡ af858fea-78d4-11eb-1499-3558c5ea680c
begin
	hcoeffs = Dict{String, Matrix}()

	const KV_LIST =  [(0, "coefficients [L2=0]"),
		(4, "coefficients [L2=4]"),
		(10, "coefficients [L2=10]"),
		(1e2, "coefficients [L2=1e2]"),
		(1e3, "coefficients [L2=1e3]"),
		(1e5, "coefficients [L2=1e5]")]

	for (l2, label) ∈ KV_LIST
    	print("Running model $(label)")
    	hcoeffs[label] = run_model(;l2_penalty=l2)
	end
end

# ╔═╡ 6d428fe4-78da-11eb-2873-9fd54b4485ea
md"""
```julia

...

Running model coefficients [L2=4]iter: 3 / 1: log likelihood of observed labels: -29098.64917580
iter: 3 / 2: log likelihood of observed labels: -28817.48280975
iter: 3 / 3: log likelihood of observed labels: -28554.20221492
iter: 3 / 4: log likelihood of observed labels: -28306.77144605
iter: 3 / 5: log likelihood of observed labels: -28073.62724612
iter: 3 / 6: log likelihood of observed labels: -27853.41039332
iter: 3 / 7: log likelihood of observed labels: -27644.93882078
iter: 3 / 8: log likelihood of observed labels: -27447.18172819
iter: 3 / 9: log likelihood of observed labels: -27259.23614139
iter: 3 / 10: log likelihood of observed labels: -27080.30692363
iter: 3 / 11: log likelihood of observed labels: -26909.69016432
iter: 3 / 12: log likelihood of observed labels: -26746.75951399
iter: 3 / 13: log likelihood of observed labels: -26590.95496394
iter: 3 / 14: log likelihood of observed labels: -26441.77360865
iter: 3 / 15: log likelihood of observed labels: -26298.76200098
iter: 3 / 20: log likelihood of observed labels: -25663.11897138
iter: 3 / 30: log likelihood of observed labels: -24683.35211346
iter: 3 / 40: log likelihood of observed labels: -23956.08137461
iter: 3 / 50: log likelihood of observed labels: -23390.05514100
iter: 3 / 60: log likelihood of observed labels: -22934.42235545
iter: 3 / 70: log likelihood of observed labels: -22558.31038461
iter: 3 / 80: log likelihood of observed labels: -22241.72295664
iter: 3 / 90: log likelihood of observed labels: -21971.04270712
iter: 3 / 100: log likelihood of observed labels: -21736.62200251
iter: 3 / 200: log likelihood of observed labels: -20412.68742436
iter: 3 / 300: log likelihood of observed labels: -19829.12657619
iter: 3 / 400: log likelihood of observed labels: -19498.14242117
iter: 3 / 500: log likelihood of observed labels: -19285.97423306

Running model coefficients [L2=10]iter: 3 / 1: log likelihood of observed labels: -29098.66002261
iter: 3 / 2: log likelihood of observed labels: -28817.53741305
iter: 3 / 3: log likelihood of observed labels: -28554.32996469
iter: 3 / 4: log likelihood of observed labels: -28306.99840133
iter: 3 / 5: log likelihood of observed labels: -28073.97638005
iter: 3 / 6: log likelihood of observed labels: -27853.90205286
iter: 3 / 7: log likelihood of observed labels: -27645.59110470
iter: 3 / 8: log likelihood of observed labels: -27448.01079546
iter: 3 / 9: log likelihood of observed labels: -27260.25646522
iter: 3 / 10: log likelihood of observed labels: -27081.53150287
iter: 3 / 11: log likelihood of observed labels: -26911.13070122
iter: 3 / 12: log likelihood of observed labels: -26748.42656501
iter: 3 / 13: log likelihood of observed labels: -26592.85806880
iter: 3 / 14: log likelihood of observed labels: -26443.92140137
iter: 3 / 15: log likelihood of observed labels: -26301.16230617
iter: 3 / 20: log likelihood of observed labels: -25666.87619856
iter: 3 / 30: log likelihood of observed labels: -24690.11746503
iter: 3 / 40: log likelihood of observed labels: -23966.03330325
iter: 3 / 50: log likelihood of observed labels: -23403.24288550
iter: 3 / 60: log likelihood of observed labels: -22950.83216858
iter: 3 / 70: log likelihood of observed labels: -22577.89582847
iter: 3 / 80: log likelihood of observed labels: -22264.42031143
iter: 3 / 90: log likelihood of observed labels: -21996.77941010
iter: 3 / 100: log likelihood of observed labels: -21765.32150230
iter: 3 / 200: log likelihood of observed labels: -20467.04068096
iter: 3 / 300: log likelihood of observed labels: -19903.56356097
iter: 3 / 400: log likelihood of observed labels: -19588.90851897
iter: 3 / 500: log likelihood of observed labels: -19390.35674625

...
```
"""

# ╔═╡ af6691e6-78d4-11eb-1d8a-1f920b745781
md"""
#### Compare coefficients

We now compare the coefficients for each of the models that were trained above. We will create a table of features and learned coefficients associated with each of the different L2 penalty values.

Below is a simple helper function that will help us create this table.
"""

# ╔═╡ af2feeb6-78d4-11eb-2873-9fd54b4485ea
df = DataFrame(Dict{Symbol, Vector{String}}(:word => ["(intercept)", important_words...]));

# ╔═╡ d52b25f6-78d5-11eb-36ac-8d6c7b6cdaf9
function add_coeffs_to_df!(df, coeffs, colname)
    insertcols!(df, colname => coeffs)
end

# ╔═╡ 001765f6-78d6-11eb-27e4-55ec5250ae23
for (_l2, label) ∈ KV_LIST
    add_coeffs_to_df!(df, reshape(hcoeffs[label], length(hcoeffs[label])), label)
end

# ╔═╡ 2276fa08-78d6-11eb-3e76-49a4a2f4094f
first(df, 5)

# ╔═╡ 280cdf00-78d6-11eb-37ac-a320a56a628e
md"""
Using the coefficients trained with L2 penalty 0, find the 5 most positive words (with largest positive coefficients). Save them to positive_words. Similarly, find the 5 most negative words (with largest negative coefficients) and save them to negative_words.
"""

# ╔═╡ 323dbeb8-78d6-11eb-08c1-13392c483a14
function top_n(;words=important_words,
		coeffs=hcoeffs["coefficients [L2=0]"][2:end], key=:positive)
    sort(
      [(w, c) for (w, c) ∈ zip(words, coeffs)],
      by=t -> t[2],
      rev= key == :positive
  )
end

# ╔═╡ 7ecbea02-78d6-11eb-1f51-2729cf921bc5
md"""
**Quiz Question**. Which of the following is **not** listed in either **positive_words** or **negative_words**?"
"""

# ╔═╡ 32217b90-78d6-11eb-05fa-d953b60010ff
positive_words = map(t -> t[1], top_n()[1:5])

# ╔═╡ 3207dc58-78d6-11eb-2b8f-a9b22b5c46fe
negative_words = map(t -> t[1], top_n(;key=:negative)[1:5])

# ╔═╡ 31ef8180-78d6-11eb-218d-ff1d2947d397
md"""
Let us observe the effect of increasing L2 penalty on the 10 words just selected. We provide you with a utility function to plot the coefficient path.
"""

# ╔═╡ 95841904-78d6-11eb-3423-936b80c46605
# TODO plot

# ╔═╡ ab93ec4a-78d6-11eb-1706-15ca470f7cf2
md"""
**Quiz Question**: (True/False) All coefficients consistently get smaller in size as the L2 penalty is increased.
  - True

**Quiz Question**: (True/False) The relative order of coefficients is preserved as the L2 penalty is increased. (For example, if the coefficient for 'cat' was more positive than that for 'dog', this remains true as the L2 penalty increases.)
  - True
"""

# ╔═╡ 6bd80650-780f-11eb-20f9-8f99728a72f3
md"""
### Measuring accuracy

We will now measure the classification accuracy of the model. Recall from the lecture that the classification accuracy can be computed as follows:

```math
\mbox{accuracy} = \frac{{\mbox{\# correctly classified data points}}}{\mbox{\# total data points}}
```

Recall from lecture that that the class prediction is calculated using
```math
\hat{y}_i =
\left\{
\begin{array}{ll}
      +1 & h(\mathbf{x}_i)^T\mathbf{w} > 0 \\
      -1 & h(\mathbf{x}_i)^T\mathbf{w} \leq 0 \\
\end{array}
\right.
```

**Note**: It is important to know that the model prediction code doesn't change even with the addition of an L2 penalty. The only thing that changes is the estimated coefficients used in this prediction.

"""

# ╔═╡ dcef9dea-78d6-11eb-267b-236c60163a33
function get_classif_accuracy(feature_matrix, sentiment, coeffs)
	scores = feature_matrix * coeffs
  	preds = ifelse.(scores .> 0., 1, -1)
	num_correct = sum(preds .== sentiment)
    return num_correct / size(feature_matrix)[1]  ## == accuracy
end

# ╔═╡ 30861830-78d7-11eb-0389-5978263881b6
md"""
Below, we compare the accuracy on the training data and validation data for all the models that were trained in this assignment. We first calculate the accuracy values and then build a simple report summarizing the performance for the various models.
"""

# ╔═╡ 5b0e8fd0-78d7-11eb-13d1-8dafac6a1d6a
begin
	train_accuracy, validation_accuracy = Dict(), Dict()

	for (l2, label) ∈ KV_LIST
    	train_accuracy[l2] = get_classif_accuracy(feature_matrix_train,
			sentiment_train, hcoeffs[label])
    	validation_accuracy[l2] = get_classif_accuracy(feature_matrix_valid,
			sentiment_valid, hcoeffs[label])
	end
end

# ╔═╡ 5af4099c-78d7-11eb-12f9-1f5e36d0706d
begin
  with_terminal() do
	# Build a simple report
	for key ∈ keys(validation_accuracy) |> collect |> sort
    	@printf("L2 penalty = %g", key)
    	@printf(" train accuracy = %1.3f, validation_accuracy = %1.3f\n",
				train_accuracy[key], validation_accuracy[key])
        println("----------------------------------------------------")
	end
  end
end

# ╔═╡ 8cee74b0-78d8-11eb-3acd-9bc21f53815e
# TODO Plot

# ╔═╡ 5ad9fbba-78d7-11eb-0209-6533fe570748
md"""
* **Quiz Question**: Which model (L2 = 0, 4, 10, 100, 1e3, 1e5) has the **highest** accuracy on the **training** data?
  - model with L2 = 0 (first model)

* **Quiz Question**: Which model (L2 = 0, 4, 10, 100, 1e3, 1e5) has the **highest** accuracy on the **validation** data?
  - model with L2 = 10 (third model)

* **Quiz Question**: Does the **highest** accuracy on the **training** data imply that the model is the best one?
  - No (it can indicate that the model is overfitting the training data, and therefore won't generalize well, unless...)
"""

# ╔═╡ Cell order:
# ╟─baa8ca1c-77d8-11eb-2e26-3b2445acbfa1
# ╠═33bd93ec-77d9-11eb-2da8-87271f8b30c7
# ╟─3ad4770e-77d9-11eb-3a93-addf3e427025
# ╠═6e264178-77d9-11eb-33c1-19462dde66e2
# ╠═75f2e820-77de-11eb-19eb-a5a5977b5257
# ╠═36f2c616-77de-11eb-03b4-3b92a532949a
# ╠═6d4ef592-77d9-11eb-3844-1bea7dc71094
# ╠═d7b9c812-77de-11eb-0ebb-6bb75de9725b
# ╟─6e09ca34-77d9-11eb-1a99-6578f9141004
# ╠═6ded8978-77d9-11eb-3d94-e7507a85197e
# ╟─6dd12c56-77d9-11eb-0350-1dfc08d17dd6
# ╠═6db87918-77d9-11eb-2082-edecdf8b7576
# ╠═6d69b148-77d9-11eb-04f5-034603897c31
# ╟─6d316f86-77d9-11eb-107b-1bf5708655d5
# ╠═0f6bb770-77df-11eb-0f1f-4dddbd7e2498
# ╟─734ed4d4-77df-11eb-3dab-e5c10b4cda62
# ╠═0f528d0e-77df-11eb-1415-2fb0e9d6792d
# ╠═0f362fec-77df-11eb-02a6-210a58c2839b
# ╟─0f1c6b16-77df-11eb-2197-93c4a82cc135
# ╠═0057a50e-77e0-11eb-29d6-85c2dac67743
# ╠═51a94e8e-77e3-11eb-04c8-d1389b1f5856
# ╠═0f03e474-77df-11eb-2791-f165fd011737
# ╠═0ee68550-77df-11eb-0a6e-c5bf5e3fd9ac
# ╠═0ecc0054-77df-11eb-21a4-3f677548ecf9
# ╟─0eb50d9a-77df-11eb-0b7d-557ab3ea2c81
# ╟─d39040fe-78d0-11eb-236d-15b506adaffd
# ╠═1c52f3b6-78d1-11eb-3077-bf6c471251c8
# ╠═005dc370-78d1-11eb-0fc0-039f1017af8e
# ╠═004d22cc-78d1-11eb-10ce-09eb125a96b3
# ╟─682a78ba-78d1-11eb-1706-15ca470f7cf2
# ╠═0e383a04-77df-11eb-2ec9-c5aeffc4b131
# ╠═7ad4c752-78d1-11eb-267b-236c60163a33
# ╠═9f29e2b8-78d1-11eb-0389-5978263881b6
# ╠═b093b9ac-78d1-11eb-3a2e-ffb6e1102b1f
# ╟─fd23d180-77e6-11eb-2091-af1e84d8fbcf
# ╠═fd120ad4-77e6-11eb-15ea-bb7d7ea35587
# ╟─fcf1dd2e-77e6-11eb-0619-e7fd2a622877
# ╟─d6f4542a-78d2-11eb-0209-6533fe570748
# ╠═d6976230-78d2-11eb-1fb7-dd1259c9b9c8
# ╟─4870ab12-78d3-11eb-12f9-1f5e36d0706d
# ╟─665516c2-78d3-11eb-13d1-8dafac6a1d6a
# ╠═70a490b2-78d3-11eb-3acd-9bc21f53815e
# ╟─0cea789c-78d4-11eb-36c4-c794c5e1a816
# ╟─1a7b28f8-78d4-11eb-3906-5ba3566b9ae9
# ╠═970fba24-77ed-11eb-1a3e-e769d6049078
# ╟─972e0c72-77ed-11eb-0ae7-c127694faf89
# ╠═afa2bff4-78d4-11eb-10cc-0deed0adb371
# ╠═af858fea-78d4-11eb-1499-3558c5ea680c
# ╟─6d428fe4-78da-11eb-2873-9fd54b4485ea
# ╟─af6691e6-78d4-11eb-1d8a-1f920b745781
# ╠═af2feeb6-78d4-11eb-2873-9fd54b4485ea
# ╠═d52b25f6-78d5-11eb-36ac-8d6c7b6cdaf9
# ╠═001765f6-78d6-11eb-27e4-55ec5250ae23
# ╠═2276fa08-78d6-11eb-3e76-49a4a2f4094f
# ╟─280cdf00-78d6-11eb-37ac-a320a56a628e
# ╠═323dbeb8-78d6-11eb-08c1-13392c483a14
# ╟─7ecbea02-78d6-11eb-1f51-2729cf921bc5
# ╠═32217b90-78d6-11eb-05fa-d953b60010ff
# ╠═3207dc58-78d6-11eb-2b8f-a9b22b5c46fe
# ╟─31ef8180-78d6-11eb-218d-ff1d2947d397
# ╠═95841904-78d6-11eb-3423-936b80c46605
# ╟─ab93ec4a-78d6-11eb-1706-15ca470f7cf2
# ╟─6bd80650-780f-11eb-20f9-8f99728a72f3
# ╠═dcef9dea-78d6-11eb-267b-236c60163a33
# ╟─30861830-78d7-11eb-0389-5978263881b6
# ╠═5b0e8fd0-78d7-11eb-13d1-8dafac6a1d6a
# ╠═5af4099c-78d7-11eb-12f9-1f5e36d0706d
# ╠═8cee74b0-78d8-11eb-3acd-9bc21f53815e
# ╟─5ad9fbba-78d7-11eb-0209-6533fe570748


text:
############################################################################
### hessian based standard errors
############################################################################

function se_hessian(sem_fit::SemFit; hessian = :finitediff)

    c = H_scaling(sem_fit.model)

    if hessian == :analytic
        H = hessian!(sem_fit.model, sem_fit.solution)
    elseif hessian == :analytic_last
        H = hessian(sem_fit.model)
    elseif hessian == :finitediff
        H = FiniteDiff.finite_difference_hessian(
                x -> objective!(sem_fit.model, x)[1],
                sem_fit.solution
                )
    elseif hessian == :optimizer
        throw(ArgumentError("standard errors from the optimizer hessian are not implemented yet"))
    elseif hessian == :expected
        throw(ArgumentError("standard errors based on the expected hessian are not implemented yet"))
    else
        throw(ArgumentError("I dont know how to compute `$how` standard-errors"))
    end

    invH = c*inv(H)
    se = sqrt.(diag(invH))

    return se
end

# Addition functions -------------------------------------------------------------
H_scaling(model::AbstractSemSingle) =
    H_scaling(
        model,
        model.observed,
        model.imply,
        model.diff,
        model.loss.functions...)

H_scaling(model, obs, imp, diff, lossfun::Union{SemML, SemWLS}) =
    2/(n_obs(model)-1)

H_scaling(model, obs, imp, diff, lossfun::SemFIML) =
    2/n_obs(model)

H_scaling(model::SemEnsemble) =
    2/n_obs(model)


text:
using TrajectoryOptimization: state, control

@testset "Objectives" begin
    n, m = rand(10:20), rand(5:10)
    N = 11

    Q = Diagonal(@SVector fill(0.1, n))
    R = Diagonal(@SVector fill(0.01, m))
    H = rand(m, n)
    q = @SVector rand(n)
    r = @SVector rand(m)
    c = rand()
    Qf = Diagonal(@SVector fill(10.0, n))
    xf = @SVector ones(n)

    @testset "Constructors" begin
        # Simple copy objective constructor
        qcost = QuadraticCost(Q, R, H, q, r, c)
        obj = Objective(qcost, N)
        @test length(obj) == N
        @test state_dim(obj) == n
        @test control_dim(obj) == m
        @test TO.get_J(obj) === obj.J
        @test obj[end].terminal == false

        # Maintains associativity
        @test obj[1].Q === obj[2].Q
        @test obj[1].q === obj[2].q
        qcost.r[1] = 1
        @test obj[3].r[1] ≈ 1

        # Terminal cost constructor
        qterm = QuadraticCost(Qf, zero(R), q = q, c = c, terminal = true)
        dterm = DiagonalCost(Qf, zero(R), q = q, c = c, terminal = true)
        obj = Objective(qcost, qterm, N)
        @test eltype(obj) <: QuadraticCost
        @test obj[1].Q === obj[2].Q
        @test obj[1].q === obj[2].q
        @test obj[1].Q ≈ Q
        @test obj[end].Q ≈ Qf
        @test obj[end].R ≈ zero(R)
        @test obj[end].q ≈ q
        @test obj[end].terminal == true

        obj = Objective(qcost, dterm, N)
        @test eltype(obj) == QuadraticCost{
            n,
            m,
            Float64,
            Diagonal{Float64,SVector{n,Float64}},
            Diagonal{Float64,SVector{m,Float64}},
        }
        @test obj[1].Q === obj[2].Q
        @test obj[1].q === obj[2].q
        @test obj[1].Q ≈ Q
        @test obj[end].Q ≈ Qf
        @test obj[end].R ≈ zero(R)
        @test obj[end].q ≈ q
        @test obj[end].terminal == true

        dcost = DiagonalCost(Q, R, q = q, r = r)
        obj = Objective(dcost, dterm, N)
        @test eltype(obj) == DiagonalCost{n,m,Float64}
        @test obj[1].Q == Q
        @test obj[end].Q == Qf

        # Cost vector constructor
        costs = [copy(k < N ? qcost : qterm) for k = 1:N]
        obj = Objective(costs)
        @test length(obj) == N
        @test eltype(obj) == QuadraticCost{
            n,
            m,
            Float64,
            Diagonal{Float64,MVector{n,Float64}},
            Diagonal{Float64,MVector{m,Float64}},
        }
        @test obj[1] !== obj[2]
        @test obj[1].Q !== obj[2].Q
        @test obj[1].q !== obj[2].q

        # LQR constructors
        xf = @SVector rand(n)
        obj = LQRObjective(Q, R, Qf, xf, N)
        @test length(obj) == N
        @test eltype(obj) == DiagonalCost{n,m,Float64}
        @test obj[1].Q ≈ Q
        @test obj[1].q ≈ -Q * xf
        @test obj[2].r ≈ zero(r)
        @test obj[2].c ≈ 0.5 * xf'Q * xf
        @test obj[end].c ≈ 0.5 * xf'Qf * xf

        obj = LQRObjective(Matrix(Q), R, Qf, xf, N)
        @test eltype(obj) <: QuadraticCost{n,m,Float64,<:SizedMatrix,<:Diagonal}
        @test obj[1].Q ≈ Q
        @test obj[1].q ≈ -Q * xf
        @test obj[2].r ≈ zero(r)
        @test obj[2].c ≈ 0.5 * xf'Q * xf
        @test obj[end].c ≈ 0.5 * xf'Qf * xf

        obj = LQRObjective(Q.diag, R.diag, Qf, xf, N)
        @test eltype(obj) <: DiagonalCost{n,m}
        @test obj[1].Q ≈ Q
        @test obj[1].q ≈ -Q * xf
        @test obj[2].r ≈ zero(r)
        @test obj[2].c ≈ 0.5 * xf'Q * xf
        @test obj[end].c ≈ 0.5 * xf'Qf * xf

        obj = LQRObjective(Diagonal(Vector(Q.diag)), R.diag, Vector(Qf.diag), xf, N)
        @test eltype(obj) <: DiagonalCost{n,m}
    end

    @testset "Evaluation and expansion" begin
        # Evaluation and expansion functions
        dt = 0.01
        N = 101
        Z = Traj([KnotPoint(rand(SVector{n}), rand(SVector{m}), dt * (k < N)) for k = 1:N])
        uref = @SVector rand(m)
        obj = LQRObjective(Q, R, Qf, xf, N, uf = uref)
        J = sum([
                0.5 * (x - xf)'Q * (x - xf) + 0.5 * (u - uref)'R * (u - uref)
                for (x, u) in zip(states(Z)[1:N-1], controls(Z))
            ]) * dt + 0.5 * (state(Z[end]) - xf)'Qf * (state(Z[end]) - xf)
        @test cost(obj, Z) ≈ J

        TO.cost!(obj, Z)
        @test sum(TO.get_J(obj)) ≈ J

        E = TO.QuadraticObjective(n, m, N)
        TO.cost_gradient!(E, obj, Z)
        @test all([E[k].q ≈ obj[k].Q * (state(Z[k]) - xf) * (k < N ? dt : 1.0) for k = 1:N])
        @test all([E[k].r ≈ R * (control(Z[k]) - uref) * dt for k = 1:N-1])

        TO.cost_hessian!(E, obj, Z)
        @test all([E[k].Q ≈ obj[k].Q * (k < N ? dt : 1.0) for k = 1:N])
        @test all([E[k].R ≈ R * dt for k = 1:N-1])
    end
end


text:
module RootFinding

export bisect_root, secant_root

function bisect_root(f, a, b; rtol=1e-8, atol=1e-8)
    fa = f(a)
    fb = f(b)

    @assert fa*fb < zero(fa)

    while abs(a-b) > atol + rtol*(abs(a) + abs(b))/2
        m = (a + b)/2
        fm = f(m)

        if fa*fm < zero(fa)
            fb = fm
            b = m
        else
            fa = fm
            a = m
        end
    end

    return (a+b)/2
end

function secant_root(f, a, b; rtol=1e-8, atol=1e-8)
    fa = f(a)
    fb = f(b)

    @assert fa*fb < zero(fa)

    # Secant guesses and function values:
    sg0 = a
    sg1 = b
    sf0 = fa
    sf1 = fb

    while abs(sg1-sg0) > atol + rtol*(abs(sg0) + abs(sg1))/2
        sg = sg0 - sf0*(sg1-sg0)/(sf1-sf0)

        if sg < a || sg > b
            m = (a+b)/2
            fm = f(m)

            if fm*fa > zero(fa)
                a = m
                fa = fm
            else
                b = m
                fb = fm
            end

            sg0 = a
            sg1 = b
            sf0 = fa
            sf1 = fb
        else
            sg0 = sg1
            sf0 = sf1
            sg1 = sg
            sf1 = f(sg)
        end
    end

    return sg1
end

end # module


text:
include("debug.jl")
include("fpump.jl")

function create_root_model!(optimizer::MOI.AbstractOptimizer, jp::JuniperProblem)
    ps = jp.options.log_levels

    jp.model = Model()
    lb = jp.l_var
    ub = jp.u_var
    # all continuous we solve relaxation first
    @variable(jp.model, lb[i] <= x[i=1:jp.num_var] <= ub[i])

    for i=1:jp.num_var
        JuMP.set_start_value(x[i], jp.primal_start[i])
    end

    if jp.options.registered_functions !== nothing
        for reg_f in jp.options.registered_functions
            if reg_f.gradf === nothing
                JuMP.register(jp.model, reg_f.s, reg_f.dimension, reg_f.f; autodiff=reg_f.autodiff)
            elseif reg_f.grad2f === nothing
                JuMP.register(jp.model, reg_f.s, reg_f.dimension, reg_f.f, reg_f.gradf)
            else
                JuMP.register(jp.model, reg_f.s, reg_f.dimension, reg_f.f, reg_f.gradf, reg_f.grad2f)
            end
        end
    end

    # TODO check whether it is supported
    if optimizer.nlp_data.has_objective
        obj_expr = MOI.objective_expr(optimizer.nlp_data.evaluator)
        expr_dereferencing!(obj_expr, jp.model)
        JuMP.set_NL_objective(jp.model, optimizer.sense, obj_expr)
    elseif optimizer.objective !== nothing
        MOI.set(jp.model, MOI.ObjectiveFunction{typeof(optimizer.objective)}(), optimizer.objective)
        MOI.set(jp.model, MOI.ObjectiveSense(), optimizer.sense)
    end

    backend = JuMP.backend(jp.model);
    llc = optimizer.linear_le_constraints
    lgc = optimizer.linear_ge_constraints
    lec = optimizer.linear_eq_constraints
    qlc = optimizer.quadratic_le_constraints
    qgc = optimizer.quadratic_ge_constraints
    qec = optimizer.quadratic_eq_constraints
    for constr_type in [llc, lgc, lec, qlc, qgc, qec]
        for constr in constr_type
            MOI.add_constraint(backend, constr[1], constr[2])
        end
    end
    for i in 1:jp.num_nl_constr
        constr_expr = MOI.constraint_expr(optimizer.nlp_data.evaluator, i)
        expr_dereferencing!(constr_expr, jp.model)
        JuMP.add_NL_constraint(jp.model, constr_expr)
    end

    (:All in ps || :Info in ps) && print_info(jp)

    jp.x = x
end

function solve_root_model!(jp::JuniperProblem)
    status, backend = optimize_get_status_backend(jp.model; solver=jp.nl_solver)
    jp.relaxation_status = status
    restarts = 0
    max_restarts = jp.options.num_resolve_root_relaxation
    jp.options.debug && debug_init(jp.debugDict)
    while !state_is_optimal(jp.relaxation_status; allow_almost=jp.options.allow_almost_solved) &&
        restarts < max_restarts && time()-jp.start_time < jp.options.time_limit

        # TODO freemodel for Knitro
        restart_values = generate_random_restart(jp)
        jp.options.debug && debug_restart_values(jp.debugDict,restart_values)
        for i=1:jp.num_var
            JuMP.set_start_value(jp.x[i], restart_values[i])
        end
        JuMP.optimize!(jp.model)
        jp.relaxation_status = MOI.get(backend, MOI.TerminationStatus())
        restarts += 1
    end
    if only_almost_solved(jp.relaxation_status)
        @warn "The relaxation is only almost solved."
    end
    return restarts
end


text:
using Plots
include(joinpath(@__DIR__, "rocket_fuel_slosh.jl"))
include(joinpath(pwd(), "src/models/visualize.jl"))

X̄, Ū = unpack(Z̄, prob)

t_nominal = [0.0, [sum([Ū[i][end] for i = 1:t]) for t = 1:T-1]...]
plot(t_nominal[1:end-1], hcat(Ū...)[1:end-1, :]', linetype=:steppost)
@show sum([Ū[t][end] for t = 1:T-1]) # works when 2.72

vis = Visualizer()
open(vis)
visualize!(vis, model_nominal, X̄, Δt = Ū[1][end])

# COM traj
xx_nom = [X̄[t][1] for t = 1:T]
zz_nom = [X̄[t][2] for t = 1:T]

X̄_dpo, Ū_dpo = unpack(Z[prob_dpo.prob.idx.nom], prob_dpo.prob.prob.nom)
xx = [X̄_dpo[t][1] for t = 1:T]
zz = [X̄_dpo[t][2] for t = 1:T]

using Plots
plot(xx_nom, zz_nom, color = :purple)
plot!(xx, zz, color = :orange)

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
const openers = "([{<"
const closers = ")]}>"
const corrupted_scores = [3, 57, 1197, 25137]
const incomplete_scores = [1, 2, 3, 4]

part1(input) = sum(check_line.(eachline(input), part = 1))

function part2(input)
    scores = check_line.(eachline(input), part = 2)
    return sort(filter(>(0), scores))[(1 + end) ÷ 2]
end

function check_line(line; part)
    stack = Char[]
    for c in line
        if c in openers
            push!(stack, c)
        elseif c in closers
            if findfirst(c, closers) != findfirst(pop!(stack), openers)
                part == 1 && return corrupted_scores[findfirst(c, closers)]
                part == 2 && return 0
            end
        end
    end
    part == 1 && return 0
    part == 2 && return foldl((x, y) -> 5 * x + y,
                              (incomplete_scores[findfirst(c, openers)]
                               for c in reverse(stack)))
end


text:
abstract type TN93 <: NASM end

const _πA(mod::TN93) = mod.πA
const _πC(mod::TN93) = mod.πC
const _πG(mod::TN93) = mod.πG
const _πT(mod::TN93) = mod.πT

function P(mod::TN93, t::Array{Float64})
  return [P(mod, i) for i in t]
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct DataInput <: File
        name::String
        descriptor::String
        col::Array{Symbol,1}
    end

Read .csv file with specific column names

# Arguments
- `name::String`: input file name
- `descriptor::String`: file descriptor
- `col::Array{Symbol,1}`: data column names
"""
mutable struct DataInput <: File
    "input file name"
    name::String
    "file descriptor"
    descriptor::String
    "data column names"
    col::Array{Symbol,1}
end

function DataInput(; name, descriptor, col, )
    DataInput(name, descriptor, col, )
end

"""Get [`DataInput`](@ref) `name`."""
get_name(value::DataInput) = value.name
"""Get [`DataInput`](@ref) `descriptor`."""
get_descriptor(value::DataInput) = value.descriptor
"""Get [`DataInput`](@ref) `col`."""
get_col(value::DataInput) = value.col

"""Set [`DataInput`](@ref) `name`."""
set_name!(value::DataInput, val) = value.name = val
"""Set [`DataInput`](@ref) `descriptor`."""
set_descriptor!(value::DataInput, val) = value.descriptor = val
"""Set [`DataInput`](@ref) `col`."""
set_col!(value::DataInput, val) = value.col = val

max_stars_repo_name: NREL/SLiDE

max_issues_repo_name: NREL/SLiDE

max_forks_repo_name: NREL/SLiDE


text:
"""
    kernel = Skernel(k,σ=1.0)

Construct the scattering kernel functions in momentum space.
The kernel function is used to construct both the energy-damping term, and noise in the energy-damped
stochastic projected Gross-Pitaevskii equation.

`k`: momentum-space grid in where scattering is to be evaluated.

`σ`: width parameter for `k`-space degrees of freedom in the case of N=1,2 spatial dimensions.

External links

[Low-Dimensional Stochastic Projected Gross-Pitaevskii Equation, Bradley, Rooney, MacDonald, Physical Review A 92, 033631 (2015)](https://arxiv.org/abs/1507.02023)

"""
function Skernel1(k,σ=1.0)
    erfcx(abs(k)*σ/sqrt(2))/sqrt(8*π*σ^2)
end
function Skernel2(k,σ=1.0)
    k==0.0 ? 0.0 : exp(abs(k*σ/2)^2)*besselk(0,abs(k*σ/2)^2)/(2*π)
end
function Skernel3(k)
    k==0.0 ? 0.0 : 1/abs(k)
end

Skernel(k::Array{Float64,1},σ=1.0) = Skernel1.(k,σ)
Skernel(k::Array{Float64,2},σ=1.0) = Skernel2.(k,σ)
Skernel(k::Array{Float64,3},σ=1.0) = Skernel3.(k)


text:
using FinEtools
using FinEtoolsFlexBeams.FESetCorotBeamModule: FESetL2CorotBeam
using FinEtoolsFlexBeams.CrossSectionModule: CrossSectionCircle, CrossSectionRectangle
using FinEtoolsFlexBeams.MeshFrameMemberModule: frame_member, merge_members
using FinEtoolsFlexBeams.VisUtilModule: plot_nodes, plot_midline, render, plot_space_box, plot_solid, space_aspectratio, default_layout_3d, plot_from_json
using PlotlyJS
using JSON

L = 42
xyz = [0 0 0;
0 L/4 L*1/4;
L/4 L/4 L*2/4;
L/4 0 L*3/4;
0 0 L]
nL = 20

cs = CrossSectionCircle(s -> 5.9910, s -> [0.0, 0.0, 1.0])
fens, fes = frame_member(xyz, nL, cs)
plots = cat(plot_nodes(fens),
    plot_midline(fens, fes; color = "rgb(155, 155, 255)", lwidth = 4),
    dims = 1)
layout = default_layout_3d(width = 500, height = 500)
layout[:title] = ""

pl = plot(plots, layout;
    config=PlotConfig(showLink = true, toImageButtonOptions = Dict(:format=>"webp")))
display(pl)

f(frame) = if frame < 10; "0$frame"; else "$frame"; end

sleep(0.5)
let
    frame = 1
    for eyex in 1.0:-0.1:-1.0
        pl.plot.layout["scene"][:camera][:eye][:x] = eyex
        pl.plot.layout[:title] = "$(eyex)"
        react!(pl, pl.plot.data, pl.plot.layout)
        sleep(0.05)
        savefig(pl, "a$(f(frame)).png")
        frame += 1
    end
end

# Convert to a gif
# magick a*.png a.gif

# The below works when run interactively from the command line:
# savejson(pl, "plot.json")

# pl = plot_from_json("plot.json")

# sleep(0.5)
# for eyex in 1.0:-0.1:-1.0
#     pl.plot.layout["scene"][:camera][:eye][:x] = eyex
#     pl.plot.layout[:title] = "$(eyex)"
#     react!(pl, pl.plot.data, pl.plot.layout)
#     sleep(0.5)
#     savefig(pl, "a.png")
# end


text:
"""
	CLFQP <: Policy

Struct representing a Control Lyapunov Function based quadratic program (CLFQP).

# Fields
- `control`: function u = control(x) describing the control law.
"""
struct CLFQP <: Policy
    control
end

"""
	(k::CLFQP)(x)

Solve CLF-QP for control input at state x.
"""
function (κ::CLFQP)(x)
    return κ.control(x)
end

"""
	CLFQP(Σ::ControlAffineSystem, CLF::ControlLyapunovFunction)
    CLFQP(Σ::ControlAffineSystem, CLF::ControlLyapunovFunction, A, b)
    CLFQP(Σ::ControlAffineSystem, CLF::ControlLyapunovFunction, A, b, p::Float64)

Construct a CLF-QP from a ControlAffineSystem and ControlLyapunovFunction.
"""
function CLFQP(Σ::ControlAffineSystem, CLF::ControlLyapunovFunction)
    function control(x)
        model = Model(OSQP.Optimizer)
        set_silent(model)
        Σ.m == 1 ? @variable(model, u) : @variable(model, u[1:(Σ.m)])
        @objective(model, Min, (1 / 2)u'u)
        @constraint(model, clf, CLF.∇V(x) * (Σ.f(x) + Σ.g(x) * u) <= -CLF.α(CLF(x)))
        optimize!(model)

        return Σ.m == 1 ? value(u) : value.(u)
    end

    return CLFQP(control)
end

function CLFQP(Σ::ControlAffineSystem, CLF::ControlLyapunovFunction, A, b)
    function control(x)
        model = Model(OSQP.Optimizer)
        set_silent(model)
        Σ.m == 1 ? @variable(model, u) : @variable(model, u[1:(Σ.m)])
        @objective(model, Min, (1 / 2)u'u)
        @constraint(model, clf, CLF.∇V(x) * (Σ.f(x) + Σ.g(x) * u) <= -CLF.α(CLF(x)))
        @constraint(model, A * u .<= b)
        optimize!(model)

        return Σ.m == 1 ? value(u) : value.(u)
    end

    return CLFQP(control)
end

function CLFQP(Σ::ControlAffineSystem, CLF::ControlLyapunovFunction, A, b, p::Float64)
    function control(x)
        model = Model(OSQP.Optimizer)
        set_silent(model)
        Σ.m == 1 ? @variable(model, u) : @variable(model, u[1:(Σ.m)])
        @variable(model, δ)
        @objective(model, Min, (1 / 2)u'u + p * δ^2)
        @constraint(model, clf, CLF.∇V(x) * (Σ.f(x) + Σ.g(x) * u) <= -CLF.α(CLF(x)) + δ)
        @constraint(model, A * u .<= b)
        optimize!(model)

        return Σ.m == 1 ? value(u) : value.(u)
    end

    return CLFQP(control)
end


text:
# TODO: Fill me in
type VegaTransform
	x::Dict{Any, Any}
end

tojs(x::VegaTransform) = x.x

Base.copy(x::VegaTransform) = VegaTransform()

max_stars_repo_name: jpthor/vega


text:
# C. Cooperative Groups

export this_grid, sync_grid

"""
    this_grid()

Returns a `grid_handle` of the grid group this thread belongs to. Only available if a
cooperative kernel is launched.
"""
this_grid

if VERSION >= v"1.2.0-DEV.512"
    @inline this_grid() =
        ccall("extern cudaCGGetIntrinsicHandle", llvmcall, Culonglong, (Cuint,), UInt32(1))
else
    @eval @inline this_grid() = Base.llvmcall(
        $("declare i64 @cudaCGGetIntrinsicHandle(i32)",
          "%rv = call i64 @cudaCGGetIntrinsicHandle(i32 1)
           ret i64 %rv"), Culonglong,
        Tuple{})
end

"""
    sync_grid(grid_handle::Culonglong)

Waits until all threads in all blocks in the grid `grid_handle` have reached this point and
all global memory accesses made by these threads prior to `sync_grid()` are visible to all
threads in the grid. A 32-bit integer `cudaError_t` is returned.
"""
sync_grid

if VERSION >= v"1.2.0-DEV.512"
    @inline sync_grid(grid_handle::Culonglong) =
        ccall("extern cudaCGSynchronize", llvmcall, cudaError_t,
              (Culonglong, Cuint), grid_handle, UInt32(0))
else
    @eval @inline sync_grid(grid_handle::Culonglong) = Base.llvmcall(
        $("declare i32 @cudaCGSynchronize(i64, i32)",
         "%rv = call i32 @cudaCGSynchronize(i64 %0, i32 0)
           ret i32 %rv"), cudaError_t,
        Tuple{Culonglong}, grid_handle)
end


text:
"""
    KrylovKit

A Julia package collecting a number of Krylov-based algorithms for linear problems,
singular value and eigenvalue problems and the application of functions of linear maps or
operators to vectors.

KrylovKit accepts general functions or callable objects as linear maps, and general Julia
objects with vector like behavior as vectors.

The high level interface of KrylovKit is provided by the following functions:

  - [`linsolve`](@ref): solve linear systems
  - [`eigsolve`](@ref): find a few eigenvalues and corresponding eigenvectors
  - [`geneigsolve`](@ref): find a few generalized eigenvalues and corresponding vectors
  - [`svdsolve`](@ref): find a few singular values and corresponding left and right
    singular vectors
  - [`exponentiate`](@ref): apply the exponential of a linear map to a vector
  - [`expintegrator`](@ref): exponential integrator for a linear non-homogeneous ODE,
    computes a linear combination of the `ϕⱼ` functions which generalize `ϕ₀(z) = exp(z)`.
"""
module KrylovKit

using LinearAlgebra
using Printf
using ChainRulesCore
const IndexRange = AbstractRange{Int}

export linsolve, eigsolve, geneigsolve, svdsolve, schursolve, exponentiate, expintegrator
export orthogonalize, orthogonalize!, orthonormalize, orthonormalize!
export basis, rayleighquotient, residual, normres, rayleighextension
export initialize, initialize!, expand!, shrink!
export ClassicalGramSchmidt, ClassicalGramSchmidt2, ClassicalGramSchmidtIR
export ModifiedGramSchmidt, ModifiedGramSchmidt2, ModifiedGramSchmidtIR
export LanczosIterator, ArnoldiIterator, GKLIterator
export CG, GMRES, BiCGStab, Lanczos, Arnoldi, GKL, GolubYe
export KrylovDefaults, EigSorter
export RecursiveVec, InnerProductVec

# Multithreading
const _NTHREADS = Ref(1)
get_num_threads() = _NTHREADS[]

function set_num_threads(n::Int)
    N = Base.Threads.nthreads()
    if n > N
        n = N
        _set_num_threads_warn(n)
    end
    _NTHREADS[] = n
end
@noinline function _set_num_threads_warn(n)
    @warn "Maximal number of threads limited by number of Julia threads,
            setting number of threads equal to Threads.nthreads() = $n"
end

enable_threads() = set_num_threads(Base.Threads.nthreads())
disable_threads() = set_num_threads(1)

function __init__()
    set_num_threads(Base.Threads.nthreads())
end

struct SplitRange
    start::Int
    step::Int
    stop::Int
    innerlength::Int
    outerlength1::Int
    outerlength::Int
end
function splitrange(r::OrdinalRange, n::Integer)
    start = first(r)
    stp = step(r)
    stop = last(r)
    l = length(r)
    innerlength = div(l, n)
    outerlength1 = l - n * innerlength
    outerlength = n
    return SplitRange(start, stp, stop, innerlength, outerlength1, outerlength)
end
function Base.iterate(r::SplitRange, i = 1)
    step = r.step
    if i <= r.outerlength1
        offset = (i-1)*(r.innerlength+1)*step
        start = r.start + offset
        stop = start + step*r.innerlength
    elseif i <= r.outerlength
        offset = (r.outerlength1 + (i-1)*r.innerlength)*step
        start = r.start + offset
        stop = start + step*(r.innerlength-1)
    else
        return nothing
    end
    return StepRange(start, step, stop), i+1
end
Base.length(r::SplitRange) = r.outerlength

# Algorithm types
include("algorithms.jl")

# Structures to store a list of basis vectors
"""
    abstract type Basis{T} end

An abstract type to collect specific types for representing a basis of vectors of type `T`.

Implementations of `Basis{T}` behave in many ways like `Vector{T}` and should have a
`length`, can be indexed (`getindex` and `setindex!`), iterated over (`iterate`), and
support resizing (`resize!`, `pop!`, `push!`, `empty!`, `sizehint!`).

The type `T` denotes the type of the elements stored in an `Basis{T}` and can be any custom
type that has vector like behavior (as defined in the docs of KrylovKit).

See [`OrthonormalBasis`](@ref) for a specific implementation.
"""
abstract type Basis{T} end

include("orthonormal.jl")

# Dense linear algebra structures and functions used in the algorithms below
include("dense/givens.jl")
include("dense/linalg.jl")
include("dense/packedhessenberg.jl")
include("dense/reflector.jl")

# Simple coordinate basis vector, i.e. a vector of all zeros and a single one on position `k`:
"""
    SimpleBasisVector(m, k)

Construct a simple struct `SimpleBasisVector <: AbstractVector{Bool}` representing a
coordinate basis vector of length `m` in the direction of `k`, i.e. for
`e_k = SimpleBasisVector(m, k)` we have `length(e_k) = m` and `e_k[i] = (i == k)`.
"""
struct SimpleBasisVector <: AbstractVector{Bool}
    m::Int
    k::Int
end
Base.axes(e::SimpleBasisVector) = (Base.OneTo(e.m),)
Base.size(e::SimpleBasisVector) = (e.m,)
@inline function Base.getindex(e::SimpleBasisVector, i)
    @boundscheck Base.checkbounds(e, i)
    return e.k == i
end

# some often used tools
function checkposdef(z)
    r = checkhermitian(z)
    r > 0 || error("operator does not appear to be positive definite: diagonal element $z")
    return r
end
function checkhermitian(z, n = abs(z))
    imag(z) <= sqrt(max(eps(n), eps(one(n)))) ||
        error("operator does not appear to be hermitian: diagonal element $z")
    return real(z)
end

# apply operators
include("apply.jl")

# Krylov and related factorizations and their iterators
include("factorizations/krylov.jl")
include("factorizations/lanczos.jl")
include("factorizations/arnoldi.jl")
include("factorizations/gkl.jl")

# A general structure to pass on convergence information
"""
    struct ConvergenceInfo{S,T}
        converged::Int
        residual::T
        normres::S
        numiter::Int
        numops::Int
    end

Used to return information about the solution found by the iterative method.

  - `converged`: the number of solutions that have converged according to an appropriate
    error measure and requested tolerance for the problem. Its value can be zero or one for
    [`linsolve`](@ref), [`exponentiate`](@ref) and  [`expintegrator`](@ref), or any integer
    `>= 0` for [`eigsolve`](@ref), [`schursolve`](@ref) or [`svdsolve`](@ref).
  - `residual:` the (list of) residual(s) for the problem, or `nothing` for problems without
    the concept of a residual (i.e. `exponentiate`, `expintegrator`). This is a single
    vector (of the same type as the type of vectors used in the problem) for `linsolve`, or
    a `Vector` of such vectors for `eigsolve`, `schursolve` or `svdsolve`.
  - `normres`: the norm of the residual(s) (in the previous field) or the value of any other
    error measure that is appropriate for the problem. This is a `Real` for `linsolve` and
    `exponentiate`, and a `Vector{<:Real}` for `eigsolve`, `schursolve` and `svdsolve`. The
    number of values in `normres` that are smaller than a predefined tolerance corresponds
    to the number `converged` of solutions that have converged.
  - `numiter`: the number of iterations (sometimes called restarts) used by the algorithm.
  - `numops`: the number of times the linear map or operator was applied
"""
struct ConvergenceInfo{S,T}
    converged::Int # how many vectors have converged, 0 or 1 for linear systems, exponentiate, any integer for eigenvalue problems
    residual::T
    normres::S
    numiter::Int
    numops::Int
end

function Base.show(io::IO, info::ConvergenceInfo)
    print(io, "ConvergenceInfo: ")
    info.converged == 0 && print(io, "no converged values ")
    info.converged == 1 && print(io, "one converged value ")
    info.converged > 1 && print(io, "$(info.converged) converged values ")
    println(
        io,
        "after ",
        info.numiter,
        " iterations and ",
        info.numops,
        " applications of the linear map;"
    )
    return println(io, "norms of residuals are given by $((info.normres...,)).")
end

# eigsolve en schursolve
include("eigsolve/eigsolve.jl")
include("eigsolve/lanczos.jl")
include("eigsolve/arnoldi.jl")
include("eigsolve/geneigsolve.jl")
include("eigsolve/golubye.jl")
include("eigsolve/svdsolve.jl")

# linsolve
include("linsolve/linsolve.jl")
include("linsolve/cg.jl")
include("linsolve/gmres.jl")
include("linsolve/bicgstab.jl")

# exponentiate
include("matrixfun/exponentiate.jl")
include("matrixfun/expintegrator.jl")

# rules for automatic differentation
include("adrules/linsolve.jl")

# custom vector types
include("recursivevec.jl")
include("innerproductvec.jl")

end


text:
export  HermiteMap,
        setcoeff!,
        clearcoeff!,
        evaluate!,
        evaluate,
        log_pdf,
        grad_x_log_pdf!,
        grad_x_log_pdf,
        hess_x_log_pdf!,
        hess_x_log_pdf,
        reduced_hess_x_log_pdf!,
        reduced_hess_x_log_pdf,
        optimize,
        inverse!,
        hybridinverse!

import Base: @propagate_inbounds

"""
$(TYPEDEF)

## Fields

$(TYPEDFIELDS)
"""
struct HermiteMap
        m::Int64
        Nx::Int64
        L::LinearTransform
        C::Array{HermiteMapComponent,1}
end

function HermiteMap(L::LinearTransform, C::Array{HermiteMapComponent,1})
        m = C[end].m
        Nx = C[end].Nx

        return HermiteMap(m, Nx, L, C)
end

"""
$(TYPEDSIGNATURES)

Access the i-th `HermiteMapComponent` of `M`.
"""
@propagate_inbounds Base.getindex(M::HermiteMap, i::Int) = getindex(M.C,i)

"""
$(TYPEDSIGNATURES)

Set the i-th `HermiteMapComponent` of `M`.
"""
@propagate_inbounds Base.setindex!(M::HermiteMap, C::HermiteMapComponent, i::Int) = setindex!(M.C,C,i)

function Base.show(io::IO, M::HermiteMap)
        println(io,"Hermite map of dimension "*string(M.Nx)*":")
        for i=1:M.Nx
                show(io, M.C[i])
        end
end

function HermiteMap(m::Int64, X::Array{Float64,2}; diag::Bool=true, factor::Float64=1.0, α::Float64 = αreg, b::String="CstProHermiteBasis")
        L = deepcopy(LinearTransform(X; diag = diag, factor = factor))

        if b ∈ ["ProHermiteBasis"; "PhyHermiteBasis";
                "CstProHermiteBasis"; "CstPhyHermiteBasis";
                "CstLinProHermiteBasis"; "CstLinPhyHermiteBasis"]
            B = eval(Symbol(b))(m)
        else
            error("The basis "*b*" is not defined.")
        end

        Nx = size(X,1)
        Nψ = 1
        coeff = zeros(Nψ)
        C = HermiteMapComponent[]
        idx = zeros(Int, Nψ, Nx)
        @inbounds for i=1:Nx
                MultiB = MultiBasis(B, i)
                vidx = idx[:,1:i]
                push!(C, HermiteMapComponent(IntegratedFunction(ExpandedFunction(MultiB, vidx, coeff)); α = α))
        end

        return HermiteMap(m, Nx, L, C)
end

"""
$(TYPEDSIGNATURES)

Set all the coefficients of the HermiteMap to zero
"""
function clearcoeff!(M::HermiteMap; start::Int64 = 1)
        for i=start:M.Nx
                clearcoeff!(M.C[i])
        end
end

"""
$(TYPEDSIGNATURES)

Evaluates in-place the HermiteMap `M` for the ensemble matrix `X`.
"""
function evaluate!(out, M::HermiteMap, X; apply_rescaling::Bool=true, start::Int64=1, P::Parallel = serial)

        Nx = M.Nx
        NxX, Ne = size(X)
        @assert NxX == Nx
        @assert size(out, 2) == Ne

        # We can apply the rescaling to all the components once
        if apply_rescaling == true
                transform!(M.L, X)
        end

        if typeof(P) <: Serial
        # We can skip the evaluation of the map on the observations componentd
        @inbounds for k=start:Nx
                Xk = view(X,1:k,:)
                col = view(out,k,:)
                evaluate!(col, M.C[k], Xk)
        end
        elseif typeof(P) <: Thread
        @inbounds Threads.@threads for k=start:Nx
                Xk = view(X,1:k,:)
                col = view(out,k,:)
                evaluate!(col, M.C[k], Xk)
        end
        end

        if apply_rescaling == true
                itransform!(M.L, X)
        end

        return out
end

"""
$(TYPEDSIGNATURES)

Evaluates the HermiteMap `M` for the ensemble matrix `X`.
"""
evaluate(M::HermiteMap, X; apply_rescaling::Bool=true, start::Int64=1, P::Parallel = serial) =
         evaluate!(zero(X), M, X; apply_rescaling = apply_rescaling, start = start, P = P)

"""
$(TYPEDSIGNATURES)

Evaluates the logarithm of the pullback of the reference density (standard normal) by the map `M` along the `component` directions  for the ensemble matrix `X`.
"""
function log_pdf(M::HermiteMap, X, component::Union{Int64, Array{Int64,1}, UnitRange{Int64}};
                 apply_rescaling::Bool = true)

        @assert size(X,1) == M.Nx "Wrong dimension of the input vector"

        # We can apply the rescaling to all the components once
        if apply_rescaling == true
                transform!(M.L, X)
        end

        logπ = zeros(size(X,2))

        if apply_rescaling == true
                # We add a - in front of the logdet
                # since we are interested inthe logdet of the inverse of M.L
                @inbounds for i in component
                        if typeof(M.L.L) <: Diagonal
                                logπ += log_pdf(M.C[i], view(X,1:i,:)) .+ log(1/(M.L.L.diag[i]))
                        elseif typeof(M.L.L) <: LowerTriangular
                                # For a triangular matrix, the inverse matrix has the inverse
                                # coefficients on its diagonal
                                logπ += log_pdf(M.C[i], view(X,1:i,:)) .+ log(1/(M.L.L[i,i]))
                        end
                end

        else
                @inbounds for i in component
                        logπ += log_pdf(M.C[i], view(X,1:i,:)) #.+ log(1/(M.L.L.diag[i]))
                end
        end

        # We can apply the rescaling to all the components once
        if apply_rescaling == true
                itransform!(M.L, X)
        end
        return logπ
end

"""
$(TYPEDSIGNATURES)

Evaluates the logarithm of the pullback of the reference density (standard normal) by the map `M` for the ensemble matrix `X`.
"""
log_pdf(M::HermiteMap, X; apply_rescaling::Bool = true) = log_pdf(M, X, 1:M.Nx; apply_rescaling = apply_rescaling)

## Compute grad_x_log_pdf

"""
$(TYPEDSIGNATURES)

Evaluates in-place the gradient of the logarithm of the pullback of the reference density (standard normal) by the map `M` for the ensemble matrix `X`.
"""
function grad_x_log_pdf!(result, cache_grad, cache, M::HermiteMap, X; apply_rescaling::Bool = true)
        Nx = M.Nx
        NxX, Ne = size(X)
        @assert size(X,1) == Nx "Wrong dimension of the input vector"
        @assert size(result) == (Ne, Nx) "Wrong dimension of the result"
        @assert size(cache_grad) == (Ne, Nx) "Wrong dimension of cache_grad"
        @assert size(cache) == (Ne, ) "Wrong dimension of cache"

        # We can apply the rescaling to all the components once
        if apply_rescaling == true
                transform!(M.L, X)
        end

        # The rescaling doesn't appears in the gradient (nor hessian), log(xy) = log(x) + log(y)
        @inbounds for i=1:Nx
                Xi = view(X,1:i,:)
                resulti = view(result,:,1:i)
                cache_gradi = view(cache_grad,:,1:i)
                grad_x_log_pdf!(cache_gradi, cache, M.C[i], Xi)
                @avx @. resulti += cache_gradi
        end

        if apply_rescaling == true
                itransform!(M.L, X)
        end

        return result
end

"""
$(TYPEDSIGNATURES)

Evaluates the gradient of the logarithm of the pullback of the reference density (standard normal) by the map `M` for the ensemble matrix `X`.
"""
grad_x_log_pdf(M::HermiteMap, X; apply_rescaling::Bool = true) = grad_x_log_pdf!(zeros(size(X,2), size(X,1)), zeros(size(X,2), size(X,1)), zeros(size(X,2)),
                                                                 M, X; apply_rescaling = apply_rescaling)

## Compute hess_x_log_pdf

"""
$(TYPEDSIGNATURES)

Evaluates in-place the hessian of the logarithm of the pullback of the reference density (standard normal) by the map `M` for the ensemble matrix `X`.
"""
function hess_x_log_pdf!(result, M::HermiteMap, X; apply_rescaling::Bool = true)

        Nx = M.Nx
        NxX, Ne = size(X)
        @assert size(X,1) == Nx "Wrong dimension of the input vector"
        @assert size(result) == (Ne, Nx, Nx) "Wrong dimension of the result"
        # @assert size(cache_hess) == (Ne, Nx, Nx) "Wrong dimension of cache_hess"
        # @assert size(cache_grad) == (Ne, Nx) "Wrong dimension of cache_grad"
        # @assert size(cache) == (Ne, ) "Wrong dimension of cache"

        # We can apply the rescaling to all the components once
        if apply_rescaling == true
                transform!(M.L, X)
        end

        # The rescaling doesn't appears in the hessian, log(xy) = log(x) + log(y)
        @inbounds for i=1:Nx
                Xi = view(X,1:i,:)
                resulti = view(result,:,1:i,1:i)
                # cache_gradi = view(cache_grad,:,1:i)
                # cache_hessi = view(cache_hess,:,1:i,1:i)
                # hess_x_log_pdf!(cache_hessi, cache_gradi, cache, M.C[i], Xi)
                # @avx @. resulti += cache_hessi
                resulti .+= hess_x_log_pdf(M.C[i], Xi)

        end

        if apply_rescaling == true
                itransform!(M.L, X)
        end

        return result
end

"""
$(TYPEDSIGNATURES)

Evaluates the hessian of the logarithm of the pullback of the reference density (standard normal) by the map `M` for the ensemble matrix `X`.
"""
hess_x_log_pdf(M::HermiteMap, X; apply_rescaling::Bool = true) = hess_x_log_pdf!(zeros(size(X,2), size(X,1), size(X,1)), M, X; apply_rescaling = apply_rescaling)

# This forms output an array of dimension (Ne, Nx, Nx) but the intermediate computation have been done
# by only looking at the active dimensions

"""
$(TYPEDSIGNATURES)

Evaluates in-place the hessian of the logarithm of the pullback of the reference density (standard normal) by the map `M` for the ensemble matrix `X`.
This routine exploits the sparsity pattern of the map `M`.
"""
function reduced_hess_x_log_pdf!(M::HermiteMap, X; apply_rescaling::Bool = true)

        Nx = M.Nx
        NxX, Ne = size(X)

        # Find the maximal length of active dimensions among the different map components
        Nmax = maximum(i -> length(active_dim(M.C[i])), 1:Nx)

        @assert size(X,1)       ==  Nx              "Wrong dimension of the input vector"
        # @assert length(result)  ==  Ne              "Wrong dimension of the result"
        # @assert size(d2cache)   == (1, Nmax, Nmax)  "Wrong dimension of the d2cache"
        # @assert size(dcache)    == (1, Nmax)        "Wrong dimension of the dcache"
        # @assert size(cache)     == (1, )            "Wrong dimension of the cache"

        Nmax = maximum(i -> length(active_dim(M.C[i])), 1:M.Nx)
        NxX, Ne = size(X)
        Nx = M.Nx

        sparsity_pattern = spzeros(Nx, Nx)
        @inbounds for i=1:Nx
                dimi = active_dim(M.C[i])
                sparsity_pattern[dimi, dimi] .= 0.0
        end

        result = ntuple(x->copy(sparsity_pattern), Ne)

        # We can apply the rescaling to all the components once
        if apply_rescaling == true
                transform!(M.L, X)
        end

        # The rescaling doesn't appears in the hessian, log(xy) = log(x) + log(y)
        @inbounds for j=1:Ne
                for i=1:Nx
                dimi = active_dim(M.C[i])
                Xi = X[1:i,j:j]

                result[j][dimi, dimi] .+= view(reduced_hess_x_log_pdf(M.C[i], Xi),1,:,:)
                end
        end

        if apply_rescaling == true
                itransform!(M.L, X)
        end

        return result
end

"""
$(TYPEDSIGNATURES)

Evaluates in-place the hessian of the logarithm of the pullback of the reference density (standard normal) by the map `M` for the ensemble matrix `X`.
This routine exploits the sparsity pattern of the map `M`.
"""
function reduced_hess_x_log_pdf!(result, d2cache, dcache, cache,M::HermiteMap, X; apply_rescaling::Bool = true)

        Nx = M.Nx
        NxX, Ne = size(X)

        # Find the maximal length of active dimensions among the different map components
        Nmax = maximum(i -> length(active_dim(M.C[i])), 1:Nx)

        @assert size(X,1)       ==  Nx              "Wrong dimension of the input vector"
        @assert length(result)  ==  Ne              "Wrong dimension of the result"
        @assert size(d2cache)   == (1, Nmax, Nmax)  "Wrong dimension of the d2cache"
        @assert size(dcache)    == (1, Nmax)        "Wrong dimension of the dcache"
        @assert size(cache)     == (1, )            "Wrong dimension of the cache"

        # We can apply the rescaling to all the components once
        if apply_rescaling == true
                transform!(M.L, X)
        end

        # The rescaling doesn't appears in the hessian, log(xy) = log(x) + log(y)
        @inbounds for j=1:Ne
                for i=1:Nx
                dimi = active_dim(M.C[i])
                Xi = X[1:i,j:j]
                d2cachei = view(d2cache,:,1:length(dimi), 1:length(dimi))
                dcachei  = view(dcache,:,1:length(dimi))
                fill!(d2cachei, 0.0)
                fill!(dcachei, 0.0)
                fill!(cache, 0.0)

                reduced_hess_x_log_pdf!(d2cachei, dcachei, cache, M.C[i], Xi)
                result[j][dimi, dimi] .+= view(d2cachei,1,:,:)
                end
        end

        if apply_rescaling == true
                itransform!(M.L, X)
        end

        return result
end

# In this version the samples are treated one at a time for scalability
"""
$(TYPEDSIGNATURES)

Evaluates the hessian of the logarithm of the pullback of the reference density (standard normal) by the map `M` for the ensemble matrix `X`.
This routine exploits the sparsity pattern of the map `M`.
"""
function reduced_hess_x_log_pdf(M::HermiteMap, X; apply_rescaling::Bool = true)

        Nmax = maximum(i -> length(active_dim(M.C[i])), 1:M.Nx)
        NxX, Ne = size(X)
        Nx = M.Nx

        sparsity_pattern = spzeros(Nx, Nx)
        @inbounds for i=1:Nx
                dimi = active_dim(M.C[i])
                sparsity_pattern[dimi, dimi] .= 0.0
        end

        result = ntuple(x->copy(sparsity_pattern), Ne)

        d2cache = zeros(1, Nmax, Nmax)
        dcache  = zeros(1, Nmax)
        cache   = zeros(1)
        reduced_hess_x_log_pdf!(result, d2cache, dcache, cache, M, X; apply_rescaling = apply_rescaling)
end

## Optimization function

"""
$(TYPEDSIGNATURES)

Optimizes the `HermiteMap` `M` based on the ensemble matrix `X`. Several kind of optimization are implemented
depending on the desired robustness of the map

Several options for `optimkind` are available:
* `kfold` uses a k-fold cross validation procedure (the more robust choice)
* `split` splits the set of samples into a training and a testing
*  An `Int64` to determine the maximum number of features for each map component.
* `nothing` to simply optimize the existing coefficients in the basis expansion.

The following optionial settings can be tuned:
* `maxterms::Int64 = 1000`: a maximum number of terms for each map component
* `withconstant::Bool = false`: the option to remove the constant feature in the greedy optimization routine, if the zero feature is the constant function for the basis of interest.
* `withqr::Bool = false`: improve the conditioning of the optimization problem with a QR factorization of the feature basis (recommended option)
* `maxpatience::Int64 = 10^5`: for `optimkind = split`, the maximum number of extra terms that can be added without decreasing the validation error before the greedy optimmization get stopped.
* `verbose::Bool = false`: prints details of the optimization procedure, current component optimize, training and validation errors, number of features added.
* `apply_rescaling::Bool=true`: standardize the ensemble matrix `X` according to the `LinearTransform` `M.L`.
* `hessprecond::Bool=true`: use a preconditioner based on the Gauss-Newton of the Hessian of the loss function to accelerate the convergence.
* `start::Int64=1`: the first component of the map to optimize
* `P::Parallel = serial`: option to use multi-threading to optimize in parallel the different map components
* `ATMcriterion::String="gradient"`: sensitivty criterion used to select the feature to add to the expansion. The default uses the derivative of the cost function with respect to the coefficient
   of the features in the reduced margin of the current set of features.
"""
function optimize(M::HermiteMap, X::Array{Float64,2}, optimkind::Union{Nothing, Int64, String};
                  maxterms::Int64 = 100, withconstant::Bool = false, withqr::Bool = false, maxpatience::Int64 = 10^5,
                  verbose::Bool = false, apply_rescaling::Bool=true, hessprecond::Bool=true,
                  start::Int64=1, P::Parallel = serial, ATMcriterion::String="gradient")
        Nx = M.Nx

        @assert size(X,1) == Nx "Error dimension of the sample"

        # We can apply the rescaling to all the components once
        if apply_rescaling == true
         transform!(M.L, X)
        end

        if typeof(P) <: Serial
                # We can skip the evaluation of the map on the observation components
                for i=start:Nx
                        Xi = view(X,1:i,:)
                        M.C[i], _ = optimize(M.C[i], Xi, optimkind;
                                             maxterms = maxterms, withconstant = withconstant,
                                             withqr = withqr, maxpatience = maxpatience, verbose = verbose,
                                             hessprecond = hessprecond, ATMcriterion = ATMcriterion)
                end

        elseif typeof(P) <: Thread
                # We can skip the evaluation of the map on the observation components,
                # ThreadPools.@qthreads perform better than Threads.@threads for non-uniform tasks
                @inbounds ThreadPools.@qthreads for i=Nx:-1:start
                         Xi = view(X,1:i,:)
                         M.C[i], _ = optimize(M.C[i], Xi, optimkind;
                                              maxterms = maxterms, withconstant = withconstant,
                                              withqr = withqr, maxpatience = maxpatience, verbose = verbose,
                                              hessprecond = hessprecond, ATMcriterion = ATMcriterion)
                end
        end

        # We can apply the rescaling to all the components once
        if apply_rescaling == true
                itransform!(M.L, X)
        end

        return M
end

"""
$(TYPEDSIGNATURES)

Optimizes the `HermiteMap` `M` based on the ensemble matrix `X`.
In this routine, a greedy feature selection is performed for each map component with the entire ensemble matrix until optimkind[i] terms
are in the i-th map component.

The following optionial settings can be tuned:
* `maxterms::Int64 = 1000`: a maximum number of terms for each map component
* `withconstant::Bool = false`: the option to remove the constant feature in the greedy optimization routine, if the zero feature is the constant function for the basis of interest.
* `withqr::Bool = false`: improve the conditioning of the optimization problem with a QR factorization of the feature basis (recommended option)
* `maxpatience::Int64 = 10^5`: for `optimkind = split`, the maximum number of extra terms that can be added without decreasing the validation error before the greedy optimmization get stopped.
* `verbose::Bool = false`: prints details of the optimization procedure, current component optimize, training and validation errors, number of features added.
* `apply_rescaling::Bool=true`: standardize the ensemble matrix `X` according to the `LinearTransform` `M.L`.
* `hessprecond::Bool=true`: use a preconditioner based on the Gauss-Newton of the Hessian of the loss function to accelerate the convergence.
* `start::Int64=1`: the first component of the map to optimize
* `P::Parallel = serial`: option to use multi-threading to optimize in parallel the different map components
* `ATMcriterion::String="gradient"`: sensitivty criterion used to select the feature to add to the expansion. The default uses the derivative of the cost function with respect to the coefficient
   of the features in the reduced margin of the current set of features.
"""
function optimize(M::HermiteMap, X::Array{Float64,2}, optimkind::Array{Int64,1};
                  maxterms::Int64 = 100, withconstant::Bool = false, withqr::Bool = false,
                  verbose::Bool = false, apply_rescaling::Bool=true, hessprecond::Bool=true,
                  start::Int64=1, P::Parallel = serial, ATMcriterion = ATMcriterion)
        Nx = M.Nx

        @assert size(X,1) == Nx "Error dimension of the sample"
        @assert size(optimkind, 1) == Nx-start+1 "Error dimension of the components"

        # We can apply the rescaling to all the components once
        if apply_rescaling == true
                transform!(M.L, X)
        end

        if typeof(P) <: Serial
        # We can skip the evaluation of the map on the observation components
	        for i=start:Nx
		        Xi = view(X,1:i,:)
		        M.C[i], _ = optimize(M.C[i], Xi, optimkind[i-start+1];
                                             maxterms = maxterms, withconstant = withconstant, withqr = withqr,
                                             verbose = verbose, hessprecond = hessprecond, ATMcriterion = ATMcriterion)
	        end

        elseif typeof(P) <: Thread
	        # We can skip the evaluation of the map on the observation components,
	        # ThreadPools.@qthreads perform better than Threads.@threads for non-uniform tasks
	        @inbounds ThreadPools.@qthreads for i=Nx:-1:start
		         Xi = view(X,1:i,:)
		         M.C[i], _ = optimize(M.C[i], Xi, optimkind[i-start+1];
                                              maxterms = maxterms, withconstant = withconstant, withqr = withqr,
                                              verbose = verbose, hessprecond = hessprecond, ATMcriterion = ATMcriterion)
	        end
        end

        # We can apply the rescaling to all the components once
        if apply_rescaling == true
        	itransform!(M.L, X)
        end

        return M
end

"""
$(TYPEDSIGNATURES)

Partially-inverts the ensemble matrix `X` such that M_i([Ystar[1:Ny,j]; X[Ny+1:i,j]]) = F[Ny+i,j] for i = start:Nx
"""
function inverse!(X, F, M::HermiteMap, Ystar::AbstractMatrix{Float64}; apply_rescaling::Bool=true, start::Int64=1, P::Parallel = serial)

        Nx = M.Nx
        NxX, Ne = size(X)

        Ny, NeY = size(Ystar)
        @assert NxX == Nx
        @assert Ne == NeY
        @assert 1 <= Ny < Nx
        @assert size(F) == (Nx, Ne)

        @view(X[1:Ny,:]) .= Ystar

        # We can apply the rescaling to all the components once
        if apply_rescaling == true
            transform!(M.L, X)
        end
        # if P == serial
        # We can skip the evaluation of the map on the observations components
        @inbounds for k = start:Nx
            Fk = view(F,k,:)
            Xk = view(X,1:k,:)
            Sk = Storage(M[k].I.f, Xk)
            inverse!(Xk, Fk, M[k], Sk)

        end

        if apply_rescaling == true
            itransform!(M.L, X)
        end
        # else P == thread
        # There is a run-race problem, and the serial version is fast enough.
        #         nthread = Threads.nthreads()
        #         @time if nthread == 1
        #                 idx_folds = 1:Ne
        #         else
        #                 q = div(Ne, nthread)
        #                 r = rem(Ne, nthread)
        #                 @assert Ne == q*nthread + r
        #                 idx_folds = UnitRange{Int64}[i < nthread ? ((i-1)*q+1:i*q) : ((i-1)*q+1:i*q+r) for i in 1:nthread]
        #         end
        #
        #         @inbounds Threads.@threads for idx in idx_folds
        #                 for k = start:Nx
        #                 Fk = view(F,k,idx)
        #                 Xk = view(X,1:k,idx)
        #                 Sk = Storage(M[k].I.f, Xk)
        #                 inverse!(Xk, Fk, M[k], Sk)
        #                 end
        #         end
        # end
end

"""
$(TYPEDSIGNATURES)

Partially-inverts the ensemble matrix `X` such that M_i([ystar; X[Ny+1:i,j]]) = F[Ny+i,j] for i = start:Nx
"""
function inverse!(X, F, M::HermiteMap, ystar::AbstractVector{Float64}; apply_rescaling::Bool=true, start::Int64=1, P::Parallel = serial)

        Nx = M.Nx
        NxX, Ne = size(X)

        Ny = size(ystar,1)
        @assert NxX == Nx
        @assert 1 <= Ny < Nx
        @assert size(F) == (Nx, Ne)

        @view(X[1:Ny,:]) .= ystar

        # We can apply the rescaling to all the components once
        if apply_rescaling == true
            transform!(M.L, X)
        end

        # We can skip the evaluation of the map on the observations components
        @inbounds for k = start:Nx
            Fk = view(F,k,:)
            Xk = view(X,1:k,:)
            Sk = Storage(M[k].I.f, Xk)
            inverse!(Xk, Fk, M[k], Sk)
        end

        if apply_rescaling == true
            itransform!(M.L, X)
        end
end

"""
$(TYPEDSIGNATURES)

Inverts the ensemble matrix `X` such that `M(X) = F`
"""
function inverse!(X, F, M::HermiteMap; apply_rescaling::Bool=true, P::Parallel = serial)

        Nx = M.Nx
        NxX, Ne = size(X)

        @assert NxX == Nx
        @assert size(F) == (Nx, Ne)

        # We can apply the rescaling to all the components once
        if apply_rescaling == true
            transform!(M.L, X)
        end
        # We can skip the evaluation of the map on the observations components
        @inbounds for k = 1:Nx
            Fk = view(F,k,:)
            Xk = view(X,1:k,:)
            Sk = Storage(M[k].I.f, Xk)
            inverse!(Xk, Fk, M[k], Sk)
        end

        if apply_rescaling == true
            itransform!(M.L, X)
        end
end

"""
$(TYPEDSIGNATURES)

Partially-inverts the ensemble matrix `X` such that M_i([Ystar[1:Ny,j]; X[Ny+1:i,j]]) = F[Ny+i,j] for i = start:Nx
"""
function hybridinverse!(X, F, M::HermiteMap, Ystar::AbstractMatrix{Float64}; apply_rescaling::Bool=true, start::Int64=1, P::Parallel = serial)

        Nx = M.Nx
        NxX, Ne = size(X)

        Ny, NeY = size(Ystar)
        @assert NxX == Nx
        @assert Ne == NeY
        @assert 1 <= Ny < Nx
        @assert size(F) == (Nx, Ne)

        @view(X[1:Ny,:]) .= Ystar

        # We can apply the rescaling to all the components once
        if apply_rescaling == true
            transform!(M.L, X)
        end
        # if P == serial
        # We can skip the evaluation of the map on the observations components
        @inbounds for k = start:Nx
            Fk = view(F,k,:)
            Xk = view(X,1:k,:)
            Sk = Storage(M[k].I.f, Xk)
            hybridinverse!(Xk, Fk, M[k], Sk; P = P)
        end

        if apply_rescaling == true
            itransform!(M.L, X)
        end
end

"""
$(TYPEDSIGNATURES)

Partially-inverts the ensemble matrix `X` such that M_i([ystar; X[Ny+1:i,j]]) = F[Ny+i,j] for i = start:Nx
"""
function hybridinverse!(X, F, M::HermiteMap, ystar::AbstractVector{Float64}; apply_rescaling::Bool=true, start::Int64=1, P::Parallel = serial)

        Nx = M.Nx
        NxX, Ne = size(X)

        Ny = size(ystar,1)
        @assert NxX == Nx
        @assert 1 <= Ny < Nx
        @assert size(F) == (Nx, Ne)

        @view(X[1:Ny,:]) .= ystar

        # We can apply the rescaling to all the components once
        if apply_rescaling == true
            transform!(M.L, X)
        end

        # We can skip the evaluation of the map on the observations components
        @inbounds for k = start:Nx
            Fk = view(F,k,:)
            Xk = view(X,1:k,:)
            Sk = Storage(M[k].I.f, Xk)
            hybridinverse!(Xk, Fk, M[k], Sk; P = P)
        end

        if apply_rescaling == true
            itransform!(M.L, X)
        end
end

"""
$(TYPEDSIGNATURES)

Inverts the ensemble matrix `X` such that `M(X) = F`
"""
function hybridinverse!(X, F, M::HermiteMap; apply_rescaling::Bool=true, P::Parallel = serial)

        Nx = M.Nx
        NxX, Ne = size(X)

        @assert NxX == Nx
        @assert size(F) == (Nx, Ne)

        # We can apply the rescaling to all the components once
        if apply_rescaling == true
            transform!(M.L, X)
        end
        # We can skip the evaluation of the map on the observations components
        @inbounds for k = 1:Nx
            Fk = view(F,k,:)
            Xk = view(X,1:k,:)
            Sk = Storage(M[k].I.f, Xk)
            hybridinverse!(Xk, Fk, M[k], Sk; P = P)
        end

        if apply_rescaling == true
            itransform!(M.L, X)
        end
end


text:
struct FA{T1<:Real, F1<:FactorAnalysis, R<:AbstractMatrix}
    data
    mat::R
    fit::F1
    cor
    loadings::AbstractArray{T1}
    cov::AbstractArray{T1}
    projection::AbstractArray{T1}
    nfactors
    dimension
end

"""
    fa(M, nfactors::Int64 = 1; args...)
Factor analysis wrapper.

`M` is a dataframe to be fitted. `nfactors` is the number of the latent traits.
`cor_method` is method for calculating the correlation matrix.
`args` are arguments for FA model. See, `MultivariateStats.fit(FactorAnalysis, X; args...)`(https://multivariatestatsjl.readthedocs.io/en/latest/fa.html).
`maxoutdim` was preferentially fixed by `nfactors` and `mean` was fixed by 0. User can modified only `method`, `tol`, `tot` and `η`.
"""

function fa(M; cor_method = :Polychoric, nfactors = 1, args...)
    if cor_method == :Polychoric
        fa_polychoric(M, nfactors; args...)
    elseif cor_method == :Pearson
        fa_pearson(M, nfactors; args...)
    else
        throw(ArgumentError("`method` should be choosen from (:Pearson, :Polychoric)."))
    end
end

function fa_polychoric(M, nfactors = 1; args...)
    S = convert(Matrix{Float64}, polycor(M))
    n = size(S, 1)
    ft = fit(FactorAnalysis, S; mean = fill(0.0, n), maxoutdim = nfactors, args...)
    FA(M, S, ft, :Polychoric, loadings(ft), cov(ft), projection(ft), nfactors, n)
end
function fa_pearson(M, nfactors = 1; args...)
    S = convert(Matrix{Float64}, cor(Matrix(M)))
    n = size(S, 1)
    ft = fit(FactorAnalysis, S; mean = fill(0.0, n), maxoutdim = nfactors, args...)
    FA(M, S, ft, :Pearson, loadings(ft), cov(ft), projection(ft), nfactors, n)
end

"""
    cov(x::FA)
Extract covariance matrix from FA models.
"""
function cov(x::FA)
    x.cov
end

"""
    loadings(x::FA)
Extract factor loadings from FA models.
"""
function loadings(x::FA)
    x.loadings
end

function show(io::IO, x::FA)
    println(io, "Factor Analysis $(x.cor)")
end

function communalities(x::FA)
    return loadings(x) .^ 2
end


text:
export adjointTest

function adjointTest(m0,pFor,out::Bool=true)

end

import jInv.Utils.checkDerivative
function checkDerivative(m0,pFor::ForwardProbType;kwargs...)
    function testFun(x,v=[])
        fc = vec(getData(x,pFor)[1])
		if !(isempty(v))
			return fc, getSensMatVec(v,x,pFor)
		else
			return fc
		end
	end
    return checkDerivative(testFun,m0;kwargs...)
end

"""
function jInv.adjointTest

automatic adjoint test for forward problems. If sensitivity matrix is m x n then
this method generates vector v,m of lengths n and m, respectively and computes

err = abs(dot(v,JTw) - dot(w,Jv))  / abs(dot(w,Jv))

Input:

	sig::Vector           - current model
	pFor::ForwardProbType - forward problem

Optional Inputs:

	out::Bool            - controls verbosity (default=false)
	tol::Real            - tol on relative error (default=1e-10)

Output:

	passed              - true/false depending whether test passed or not
	err                 - absolute error
"""
function adjointTest(sig,pFor::ForwardProbType;out::Bool=false,tol::Real=1e-10)

    (out) && println("calling forward problem")
    dobs, pFor = getData(sig,pFor)

    (out) && println("calling getSensMatVec")
    v = getRandomTestDirection(sig)
    Jv = getSensMatVec(v,sig,pFor)
    w = getRandomTestDirection(Jv)
    t1 = dot(w,Jv)
    (out) && println("t1 = $(t1)")

    (out) && println("calling getSensTMatVec")
    JTw = getSensTMatVec(w,sig,pFor)
    t2 = dot(v,JTw)
    (out) && println("t2 = $(t2)")

    err = abs(t1-t2)/abs(t1)
    return err < tol, err
end