text:
println("\n\tEstimate autocorrelations");
println("\n\t...for N1...");
xcorr = [0:10:400];
ycorr1 = zeros( length( xcorr ), 5 );
for a in [1:5]
  # Calculate autocorrelations
  ycorr1[:,a] = autocorrelation(N1, v1, a, xcorr);
  # Normalize by dividing through by Ĉ[1,a]
  ycorr1[:,a] /= ycorr1[1,a];
end
println("\n\t...for N2 (See plots below)");
ycorr2 = zeros( length( xcorr ), 5 );
for a in [1:5]
  # Calculate autocorrelations
  ycorr2[:,a] = autocorrelation(N2, v2, a, xcorr);
  # Normalize by dividing through by Ĉ[1,a]
  ycorr2[:,a] /= ycorr2[1,a];
end

# Make plots for each alpha
println("\n\tPlacing estimated autocorrelation plots 'plots1' and 'plots2':");
using Gadfly
plots1 = Plot[];
plots2 = Plot[];
for i in [1:5]
  # Plot
  push!(plots1, Gadfly.plot ( x = xcorr, y = ycorr1[:,i], Geom.line, Geom.point,
      Guide.xlabel("Separation"), Guide.ylabel("Normalized Autocorrelation"),
      Guide.title(join(["Estimated Autocorrelation for Dataset",string(i),
                        "with Sample Size N1 = $N1"]," "," "))))
  push!(plots2, Gadfly.plot ( x = xcorr, y = ycorr2[:,i], Geom.line, Geom.point,
      Guide.xlabel("Separation"), Guide.ylabel("Normalized Autocorrelation"),
      Guide.title(join(["Estimated Autocorrelation for Dataset",string(i),
                        "with Sample Size N2 = $N2"]," "," "))))
end

# Estimate autocorrelation times
println("\n\tEstimating autocorrelation times τ1 and τ2");
τ1 = Float64[];
τ2 = Float64[];
for a in [1:5]
  corra = autocorrelation(N1, v1, a, ns);
  corra /= corra[1];                          # Normalize
  push!(τ1,sum(corra) - 0.5);                 # Get autocorrelation time
end
for a in [1:5]
  corra = autocorrelation(N2, v2, a, ns);
  corra /= corra[1];                          # Normalize
  push!(τ2,sum(corra) - 0.5);                 # Get autocorrelation time
end

println("\n\tEstimated Autocorrelation Times for N1 = $N1:\n");
for i in [1:5]
  println("\t\tτ̂1v$i:\t$(τ1[i])");
end
println("\n\tEstimated Autocorrelation Times for N2 = $N2:\n");
for i in [1:5]
  println("\t\tτ̂2v$i:\t$(τ2[i])");
end
println("\n\tThe estimated autocorrelation times are clearly far better for the larger sample.");


text:
using Test
using Documenter: doctest
using EponymKeywordSyntax
const E = EponymKeywordSyntax

doctest(EponymKeywordSyntax; manual=false)

@testset "@eponym" begin
    a = 1; b = 2
    @test (A=:A, b=b) == @eponym((;A=:A, b))
    @test (a=a, b=b) == @eponym((;a, b))
    @test (b=b, a=a) == @eponym((;b, a))
    kw = (c=3, d=4)
    @test (;a=a, kw...) == @eponym((;a, kw...))

    f(args...; kw...) = (args=args, kw=kw)

    @test f(1,2) == @eponym f(1,2)
    @test f(1,2; a=a, b=b) == @eponym f(1,2;a,b)
    @test f(1,2; a=a, b=b) == @eponym f(1,2;a,b=b)
    @test f(1,2; a=b, b=b)  == @eponym f(1,2;a=b,b)
    kw = (c=3, d=4)
    @test f(;kw...) == @eponym f(;kw...)
    @test f(;a=a, hello="world", kw...) == @eponym f(;a, hello="world", kw...)

    @test_throws ArgumentError E.eponymmacro(:((a)))
    @test_throws ArgumentError E.eponymmacro(:((a=1; b)))
end

ext: jl


text:
using Pkg: Pkg

Pkg.add([
    Pkg.PackageSpec(;
        name="KernelFunctions", uuid="ec8451be-7e33-11e9-00cf-bbf324bd1392", version="0.10"
    ),
    Pkg.PackageSpec(;
        name="CalibrationErrorsDistributions",
        uuid="20087e1a-bb94-462b-b900-33d17a750383",
        version="0.2",
    ),
    Pkg.PackageSpec(;
        name="CalibrationTests", uuid="2818745e-0823-50c7-bc2d-405ac343d48b", version="0.5"
    ),
])

max_stars_repo_name: devmotion/pycalibration

max_issues_repo_name: devmotion/pycalibration

max_forks_repo_name: devmotion/pycalibration


text:
function read_data(name, dims, side)
    s = Tuple(side for _ in 1:dims)
    array = zeros(UInt8, s)

    open(name) do io
        read!(io, array)
    end

    return BitArray(array)
end

const data = read_data("ceramic200", 3, 200)

@testset "Test filters (x3 .. x6)" begin
    for n in 3:6
        filter = filter_coeffs(n)
        @test 0.98 < AutoCorrelationResampling.csum(filter) < 1.02
    end
end

autocorrelation(array) = (array |> fft .|> abs2 |> ifft |> real) / length(array)

@testset "Test resampling" begin
    for n in [2, 3, 2//3, 3//2]
        ac = autocorrelation(data)
        ac_resampled = ac_resample(ac, n)
        ac_resampled_ft = fft(ac_resampled)

        @test (ac_resampled |> imag |> maximum) < +1e-8
        @test (ac_resampled |> real |> minimum) > -1e-8
    end
end


text:
"""
    ssrand(T::Type, ny::Int, nu::Int, nstates::Int; proper=false, stable=true, Ts=nothing)

Returns a random `StateSpace` model with `ny` outputs, `nu` inputs, and `nstates` states,
whose matrix elements are normally distributed.

It is possible to specify if the system should be `proper` or `stable`.

Specify a sample time `Ts` to obtain a discrete-time system.
"""
function ssrand(T::Type, ny::Int, nu::Int, nstates::Int; proper=false, stable=true, Ts=nothing)
    A = randn(T, nstates, nstates)
    if stable
        Λ = eigvals(A)
        A = isnothing(Ts) ? A - 1.1*maximum(real(Λ))*I : A*0.9/maximum(abs.(Λ))
    end

    B = randn(T, nstates, nu)
    C = randn(T, ny, nstates)
    D = (proper ? zeros(T, ny, nu) : randn(T, ny, nu))
    return isnothing(Ts) ? StateSpace(A, B, C, D) : StateSpace(A, B, C, D, Ts)
end
ssrand(ny::Int, nu::Int, nstates::Int; kwargs...) = ssrand(Float64, ny::Int, nu::Int, nstates::Int; kwargs...)
ssrand(n; kwargs...) = ssrand(n, n, 2*n; kwargs...)
ssrand(T::Type, n; kwargs...) = ssrand(T, n, n, 2*n; kwargs...)

"""
    DemoSystems

A module with standard test systems from the control literature.

The returned systems are of type `StateSpace` or `DelayLTISystem`.

Default parameter can be changed using keyword arguments.

SISO systems
============
* `DemoSystems.lag(;T=1)`     First-order system  `1/(sT+1)`
* `DemoSystems.fotd(;T=1, τ=1)`   First-order system with time delay `exp(-sτ)/(sT+1)`
* `DemoSystems.sotd(;T=1, T2=10, τ=1)`   Second-order non-resonant system with time delay `exp(-sτ)/(sT+1)/(sT2 + 1)`
* `DemoSystems.resonant(;ω0=1, ζ=0.25)`   Second-order resonant systems `ω0^2/(s^2 + 2ζ*ω0*s + ω0^2)`

MIMO systems
============
* `DemoSystems.woodberry()`     Wood--Berry distillation column
* `DemoSystems.doylesat(;a=10)`    The spinning body example by Doyle

Fore more information, see the help for the specific systems.
"""
module DemoSystems

using ControlSystems
using LinearAlgebra

"""
    lag(;T=1, τ=1)

Returns a first-order system `StateSpace` with transfer function `1/(sT+1)`.
"""
lag(;T=1) = ss(-1/T, 1, 1/T, 0)

"""
    fotd(;T=1, τ=1)

Returns a first-order system + a time delay `exp(-sτ)/(sT+1)` (`DelayLTISystem`).
"""
fotd(;T=1, τ=1) = ss(-1/T, 1, 1/T, 0) * delay(τ)

"""
    sotd(;T=1, T2=10, τ=1)

Returns a second-order system + a time delay `exp(-sτ)/(sT+1)/(sT2 + 1)` (`DelayLTISystem`).
"""
sotd(;T=1, T2=10, τ=1) = ss(-1/T, 1, 1/T, 0)*ss(-1/T2, 1, 1/T2, 0)*delay(τ)

"""
    resonant(;ω0=1, ζ=0.25)

Returns a second-order system + a time delay `exp(-sτ)/(sT+1)/(sT2 + 1)` (`DelayLTISystem`).
"""
resonant(;ω0=1, ζ=0.25) = ss([-ζ -ω0; ω0 -ζ], [ω0; 0], [0 ω0], 0) # QUESTION: Is this the form that we like? Perhhaps not.

"""
    `doylesat(;a=10) = ss([0 a; -a 0], I(2), [1 a; -a 1], 0)`

Returns a `StateSpace` model for the spinning body (satellite) example by Doyle,
which is a classic example that illustrates that robustness analysis of MIMO systems
is more involved than for SISO systems.

*References:*

**Zhou, K., J. C. Doyle, and K. Glover**, Robust and optimal control,
Prentice hall (NJ), 1996.

**Skogestad, S, and I. Postlethwaite**, Multivariable feedback control:
analysis and design, Wiley (NY), 2007.
"""
doylesat(;a=10) = ss([0 a; -a 0], I(2), [1 a; -a 1], 0)

"""
    `woodberry()`

Returns a `DelayLTISystem` model for the Wood--Berry distillation column,
which is a classic example from the literature on process control of MIMO process.

*References:*

**Wood, R. K., and M. W. Berry**, Terminal composition control of a binary
distillation column, Chemical Engineering Science, 28.9, 1973, pp. 1707-1717.

"""
function woodberry()
    s = tf("s")
    return [12.8/(1 + 16.7*s)*delay(1.0)   -18.9/(1 + 21*s) * delay(3.0)
            6.6/(1 + 10.9*s) * delay(7.0)  -19.4/(1 + 14.4*s) * delay(3.0)]
end

end


text:
module Utils

export mergedict, getiris
export skipmean,skipmedian,skipstd
export aggregatorclskipmissing

using Statistics
using DataFrames
using CSV

# closure to create aggregator closure with skipmissing features
function aggregatorclskipmissing(fn::Function)
  function skipagg(x::Union{AbstractArray,DataFrame})
    if length(collect(skipmissing(x))) == 0
      return missing
    else
      return fn(skipmissing(x))
    end
  end
  return skipagg
end

function skipmean(x::T) where {T<:Union{AbstractArray,DataFrame}}
  if length(collect(skipmissing(x))) == 0
    missing
  else
    mean(skipmissing(x))
  end
end

function skipmedian(x::T) where {T<:Union{AbstractArray,DataFrame}}
  if length(collect(skipmissing(x))) == 0
    missing
  else
    median(skipmissing(x))
  end
end

function skipstd(x::T) where {T<:Union{AbstractArray,DataFrame}}
  if length(collect(skipmissing(x))) == 0
    missing
  else
    std(skipmissing(x))
  end
end

function mergedict(first::Dict, second::Dict)
  target = copy(first)
  for (second_key, second_value) in second
    values_both_dict =
      typeof(second_value) <: Dict &&
      typeof(get(target, second_key, nothing)) <: Dict
    if values_both_dict
      target[second_key] = mergedict(target[second_key], second_value)
    else
      target[second_key] = second_value
    end
  end
  return target
end

function getiris()
  iris = CSV.read(joinpath(Base.@__DIR__,"../data","iris.csv"))
  return iris
end

end


text:
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020
# Copyright: Georgia Institute of Technology, 2020

module InvertibleNetworks

# Dependencies
using LinearAlgebra, Random
using Statistics, Wavelets
using JOLI
using NNlib, Flux, Zygote, ChainRulesCore

# Overloads and reexports
import Base.size, Base.length, Base.getindex, Base.reverse, Base.reverse!
import Base.+, Base.*, Base.-, Base./
import LinearAlgebra.dot, LinearAlgebra.norm, LinearAlgebra.adjoint
import Flux.glorot_uniform
import CUDA: CuArray

export clear_grad!, glorot_uniform

# Getters for DenseConvDims fields
# (need to redefine here as they are not public methods in NNlib)
input_size(c::DenseConvDims) = c.I
kernel_size(::DenseConvDims{N,K,C_in,C_out,S,P,D,F}) where {N,K,C_in,C_out,S,P,D,F} = K
channels_in(::DenseConvDims{N,K,C_in,C_out,S,P,D,F}) where {N,K,C_in,C_out,S,P,D,F} = C_in
channels_out(::DenseConvDims{N,K,C_in,C_out,S,P,D,F}) where {N,K,C_in,C_out,S,P,D,F} = C_out

function DCDims(X::AbstractArray{Float32, N}, W::AbstractArray{Float32, N}; stride=1, padding=1, nc=nothing) where N
    sw = size(W)
    isnothing(nc) && (nc = sw[N-1])
    sx = (size(X)[1:N-2]..., nc, size(X)[end])
    return DenseConvDims(sx, sw; stride=Tuple(stride for i=1:N-2), padding=Tuple(padding for i=1:N-2))
end

# Utils
include("utils/parameter.jl")
include("utils/objective_functions.jl")
include("utils/dimensionality_operations.jl")
include("utils/activation_functions.jl")
include("utils/test_distributions.jl")
include("utils/neuralnet.jl")
include("utils/invertible_network_sequential.jl")
# AD rules
include("utils/chainrules.jl")

# Single network layers (invertible and non-invertible)
include("conditional_layers/conditional_layer_residual_block.jl")
include("layers/layer_residual_block.jl")
include("layers/layer_flux_block.jl")
include("layers/layer_resnet.jl")
include("layers/layer_affine.jl")
include("layers/invertible_layer_actnorm.jl")
include("layers/invertible_layer_conv1x1.jl")
include("layers/invertible_layer_basic.jl")
include("layers/invertible_layer_irim.jl")
include("layers/invertible_layer_glow.jl")
include("layers/invertible_layer_hyperbolic.jl")
include("layers/invertible_layer_hint.jl")

# Invertible network architectures
include("networks/invertible_network_hint_multiscale.jl")
include("networks/invertible_network_irim.jl")  # i-RIM: Putzky and Welling (2019)
include("networks/invertible_network_glow.jl")  # Glow: Dinh et al. (2017), Kingma and Dhariwal (2018)
include("networks/invertible_network_hyperbolic.jl")    # Hyperbolic: Lensink et al. (2019)

# Conditional layers and nets
include("conditional_layers/conditional_layer_hint.jl")
include("networks/invertible_network_conditional_hint.jl")
include("networks/invertible_network_conditional_hint_multiscale.jl")

# Jacobians
include("utils/jacobian.jl")

# gpu
include("utils/cuda.jl")

end


text:
a = Dict( "finame"  => "test",
            "val1"  => 1,
            "val2"  => 100)

for key in keys(a)
    println(@sprintf("%10s = %10s", key, string(a["$key"])))
end

a = [1,2,3]
b = mapreduce(identity, +, a)
b = mapreduce(x->x^2, +, a)

c = zeros(5)
fill!(c, 3)

m = 10
v = collect(1:100)
TN=10
function _f_kurt1(v, m)
    return (v - m)^2
end

z2 = map(x -> _f_kurt1(x, m), v[1:TN])

a = [1,2,3]
map(x -> x^2, a)

v = collect(1:10)
TN = 3
n = length(v)
A = Array{Float64, 2}(undef, TN, n-TN+1)
for k = TN:n
    icount = 1
    println(k)
    A[:, k-TN+1] = view(v, k-TN+1:k)
end

m = mean!(ones(n-TN+1)', A)

varm(A, m.parent', dims=1, corrected = false)

A = v

using StatsBase
c = [false, false, false, true]
a = [1,2,3,4]
mean(a, weights(c))

using Dierckx

x = 1:10
y = collect(1:2:20)
spl = Spline1D(x, y; k=1, bc="nearest")

t = collect(1:0.1:10)
y1 = evaluate(spl, t)


text:
"""
Methods to help handling of identical atoms in molecules
"""
module IdenticalTools

export AbstractIdentical,
       areidentical,
       Identical

abstract type AbstractIdentical end

"""
    Identical <: AbstractIdentical

Struct to handle identical information for atoms or any object in general

# Fields
- `identical::Vector{Set}` :
"""
mutable struct Identical <: AbstractIdentical
    identical::Vector{Set}
    Identical() = new([])
    Identical(n::Integer) = new(Set.(1:n))
end

"""
    Base.push!(a::AbstractIdentical, x)

Adds identical information
"""
function Base.push!(a::AbstractIdentical, x)
    s = Set(x)
    if isempty(a.identical)
        push!(a.identical, s)
        return a
    end
    n = Set()
    d = []
    for i in eachindex(a.identical)
        if !isempty( intersect(s,a.identical[i]) )
            if isempty(n)
                n = a.identical[i]
                union!(n,s)
            else
                union!(n,a.identical[i])
                push!(d,i)
            end
        end
    end
    if isempty(d) && isempty(n)
        push!(a.identical,s)
    else
        deleteat!(a.identical, d)
    end
    return a
end

"""
    areidentical(a::AbstractIdentical, x)

Tests are objects x

# Arguments
- 'identical::AbstractIdentical' : object holding identical data
- 'x' : collectable holiding keys of object to be tested for identical information
"""
function areidentical(identical::AbstractIdentical, x)
    s = Set(x)
    for y in identical.identical
        if issubset(s,y)
            return true
        end
    end
    return false
end

end #module


text:
using Documenter, LiveServer
using LiveServer: SimpleWatcher, WatchedFile, set_callback!, file_changed_callback

const deploy = "deploy" in ARGS

make() = makedocs(
    format=Documenter.HTML(
        prettyurls = deploy,
        canonical = ("deploy" in ARGS) ? "https://rogerluo.me/Brochure.jl/dev/" : nothing,
    ),
    sitename="Julia编程指南",
    authors="罗秀哲 Xiuzhe (Roger) Luo",
    pages=[
        "简介" => "index.md",
        "在开始之前" => "before-we-start.md",
        "程序的编写流程" => "workflow.md",
        "快速入门" => "quick-start.md",
        "实现你自己的稀疏矩阵" => "define-your-own-matrix.md",
        "实现你自己的自动微分" => "automatic-differentiation.md",
    ],
)

function scan_files!(dw::SimpleWatcher)
    for (root, _, files) in walkdir("src"), file in files
        push!(dw.watchedfiles, WatchedFile(joinpath(root, file)))
    end
end

function update_book_callback(fp::AbstractString)
    if splitext(fp)[2] == ".md"
        make()
    end
    file_changed_callback(fp)
    return nothing
end

function serve_book(verbose=false)
    watcher = SimpleWatcher()
    scan_files!(watcher)
    set_callback!(watcher, update_book_callback)
    make()
    serve(watcher, dir="build", verbose=verbose)
    return nothing
end

if ("s" in ARGS) || ("serve" in ARGS) || isinteractive()
    serve_book("verbose" in ARGS)
end

if deploy
    make()
    deploydocs(
        repo = "github.com/Roger-luo/Brochure.jl.git",
        target = "build",
    )
end


text:
#=##############################################################################
# # DESCRIPTION
#     Examples of VTK formatting and usage of VTKtools
# # AUTHORSHIP
#   * Author    : Eduardo J Alvarez
#   * Email     : Edo.AlvarezR@gmail.com
#   * Created   : Nov 2017
#   * License   : MIT License
#   * Modified  : Kevin Moore moorekevin711@gmail.com
=###############################################################################

function save_vtk_opt(p_stress,p_strain,p_stresslayer,p_c_bucklingplot,p_xloc,p_yloc,p_zloc,numprops,w_stress,w_strain,w_stresslayer,w_c_bucklingplot,w_xloc,w_yloc,w_zloc;
    file_name="./VTK_output/temp_kevin00",
    p_radii = [0.0,0.2,0.4,0.6,0.8,1.0],
    chord = [1.619, 1.473, 1.456, 0.789, 0.789, 0.789],
    p_twist_d = [1.619, 1.473, 1.456, 0.789, 0.789, 0.789], #twist includes pitch
    p_dihedral_d = pi/180*vcat([0 for i in 1:5], [23.20]),
    p_sweep_d = pi/180*vcat([23.20 for i in 1:6], []),
    Rtip = 1.25,
    p_airfoilx = 0.0, #TODO, input default
    p_airfoily = 0.0,
    p_strainx = 0.0,
    p_strainy = 0.0,
    w_twist_d = [1.619, 1.473, 1.456, 0.789, 0.789, 0.789],
    w_dihedral_d = pi/180*vcat([0 for i in 1:5], [23.20]),
    w_sweep_d = pi/180*vcat([23.20 for i in 1:6], []),
    w_strainx = 0.0,
    w_strainy = 0.0,
    spanloc = [0.0,1.0],
    halfspan = 10.0,
    printiter = 1,
    proploc = 0.0,
    numblades = 2,
    plots = false )

    numblades = Int(round(numblades))

    # GEOMETRY DEFINITION
    # dihedral = p_dihedral_d * pi/180
    # sweep = p_sweep_d * pi/180

    forward_offset = Rtip/10.0

    y_pos = copy(p_yloc)#Rtip*p_radii # span position of each section
    x_pos = zeros(p_yloc)-forward_offset#[y*tan(sweep[i]) for (i,y) in enumerate(y_pos)] #Already translated
    z_pos = zeros(p_yloc)#[y*tan(dihedral[i]) for (i,y) in enumerate(y_pos)]
    # PARAMETERS
    r = fill(1.0,length(p_radii))          # Expansion ratio in both surfaces of each airfoil

    azimuth_deg_spacing = round(Int, 360/numblades)
    azimuth_deg = 0.0
    for i = 1#:numprops
        for j = 1:numblades

            generateprop!(p_strainx,p_strainy,p_strain,p_stresslayer,p_stress,p_c_bucklingplot,p_airfoilx,p_airfoily,x_pos,y_pos,z_pos,proploc[i,:],p_sweep_d,p_twist_d,i,j,p_radii,azimuth_deg,file_name;plots=plots)
            azimuth_deg = azimuth_deg_spacing + azimuth_deg

        end
    end

    # ############################################
    # ########    Generate Wing  #################
    # ############################################
    #
    # y_pos = copy(w_yloc)#Rtip*w_radii # span position of each section
    # x_pos = zeros(w_yloc)#[y*tan(sweep[i]) for (i,y) in enumerate(y_pos)] #Already translated
    # z_pos = zeros(w_yloc)#[y*tan(dihedral[i]) for (i,y) in enumerate(y_pos)]
    # # PARAMETERS
    # r = fill(1.0,length(spanloc))          # Expansion ratio in both surfaces of each airfoil
    #
    # sections = fill(Tuple{Float64,Int64,Float64,Bool}[(1.0, 10, 1.0, false)],(length(spanloc)-1))
    #
    # # Leading edge position of each airfoil
    # Os = [ [x_pos[i], y_pos[i], z_pos[i]] for i in 1:size(w_strainx)[1]]
    # # Orientation of chord of each airfoil (yaw, pitch, roll)
    # orien = zeros(length(w_twist_d),3) #twist already included
    # for i = 1:length(orien[:,1])
    #     orien[i,:] = [0.0,w_sweep_d[i],270.0]
    #     if i == 1
    #         orien[i,:] = [0.0,0.0,270.0]
    #     end
    # end
    #
    # crosssections = []        # It will store here the cross sections for lofting
    # strains = []
    # stress1 = []
    # stress2 = []
    # stressplot = []
    # buckle = []
    # point_datas = []
    #
    # # Processes each airfoil geometry - STRESS/STRAIN CALCS
    # for i = 1:length(w_strainx[:,1])
    #
    #     # Read airfoil file
    #     x = w_strainx[i,:]
    #     y = w_strainy[i,:]
    #
    #     # Scales the airfoil acording to its chord length
    #     new_x = x#new_x #already scaled
    #     new_y = y#new_y
    #
    #     # plot_airfoil(new_x, new_y; style=styles[i], label="airfoil")
    #     if plots
    #         figure("wing_airfoil")
    #         plot(new_x,new_y)
    #     end
    #
    #     # Reformats into points
    #     npoints = size(new_x)[1]
    #     airfoil = Array{Float64, 1}[[new_x[j], new_y[j], 0] for j in 1:npoints]
    #
    #     # Positions the airfoil along the blade in the right p_orientation
    #     Oaxis = VTKtools.rotation_matrix(orien[i,1], orien[i,2], orien[i,3])
    #     invOaxis = inv(Oaxis)
    #     airfoil = VTKtools.countertransform(airfoil, invOaxis, Os[i])
    #
    #     push!(crosssections, airfoil)
    #
    #     push!(strains, w_strain[i,:])
    #     push!(stress1, w_stresslayer[1,i,:])
    #     push!(stress2, w_stresslayer[2,i,:])
    #     push!(stressplot, w_stress[i,:])
    #     push!(buckle, w_c_bucklingplot[i,:])
    #     push!(point_datas, [j for j in npoints*(i-1)+1:npoints*i])
    # end
    #
    # # Generates cells in VTK Legacy format
    # points, vtk_cells, strain = VTKtools.multilines2vtkmulticells(crosssections, sections;point_datas=strains)
    # _, _, stresses1 = VTKtools.multilines2vtkmulticells(crosssections, sections;point_datas=stress1)
    # _, _, stresses2 = VTKtools.multilines2vtkmulticells(crosssections, sections;point_datas=stress2)
    # _, _, stresses6 = VTKtools.multilines2vtkmulticells(crosssections, sections;point_datas=stressplot)
    # _, _, bucklecon = VTKtools.multilines2vtkmulticells(crosssections, sections;point_datas=buckle)
    # points, vtk_cells, point_datas = VTKtools.multilines2vtkmulticells(crosssections, sections;point_datas=point_datas)
    #
    # # Formats the point data for generateVTK
    # data = []
    # push!(data, Dict(
    # "field_name" => "Strain",
    # "field_type" => "scalar",
    # "field_data" => strain))
    #
    # push!(data, Dict(
    # "field_name" => "uni_stress",
    # "field_type" => "scalar",
    # "field_data" => stresses1))
    #
    # push!(data, Dict(
    # "field_name" => "weave_stress",
    # "field_type" => "scalar",
    # "field_data" => stresses2))
    #
    # push!(data, Dict(
    # "field_name" => "stress",
    # "field_type" => "scalar",
    # "field_data" => stresses6))
    #
    # push!(data, Dict(
    # "field_name" => "Buckling_Constraint",
    # "field_type" => "scalar",
    # "field_data" => bucklecon))
    #
    # push!(data, Dict(
    # "field_name" => "PointIndex",
    # "field_type" => "scalar",
    # "field_data" => point_datas))
    #
    # # Generates the vtk file
    # VTKtools.generateVTK("$(file_name)wingR.$printiter", points; cells=vtk_cells, point_data=data, keep_points=false)
    #
    #
    # #generate other wing
    # Os = [ [x_pos[i], -y_pos[i], z_pos[i]] for i in 1:size(w_strainx)[1]]
    # # Orientation of chord of each airfoil (yaw, pitch, roll)
    # orien = zeros(length(w_twist_d),3) #twist already included
    # for i = 1:length(orien[:,1])
    #     orien[i,:] = [0.0,-w_sweep_d[i],270.0]
    #     if i == 1
    #         orien[i,:] = [0.0,0.0,270.0]
    #     end
    # end
    #
    # crosssections = []        # It will store here the cross sections for lofting
    # strains = []
    # stress1 = []
    # stress2 = []
    # stressplot = []
    # buckle = []
    # point_datas = []
    # for i = length(w_strainx[:,1]):-1:1
    #
    #     # Read airfoil file
    #     x = w_strainx[i,:]
    #     y = w_strainy[i,:]
    #
    #     # Scales the airfoil acording to its chord length
    #     new_x = x#new_x #already scaled
    #     new_y = y#new_y
    #
    #     # plot_airfoil(new_x, new_y; style=styles[i], label="airfoil")
    #     if plots
    #         figure("wing_airfoil")
    #         plot(new_x,new_y)
    #     end
    #
    #     # Reformats into points
    #     npoints = size(new_x)[1]
    #     airfoil = Array{Float64, 1}[[new_x[j], new_y[j], 0] for j in 1:npoints]
    #
    #     # Positions the airfoil along the blade in the right p_orientation
    #     Oaxis = VTKtools.rotation_matrix(orien[i,1], orien[i,2], orien[i,3])
    #     invOaxis = inv(Oaxis)
    #     airfoil = VTKtools.countertransform(airfoil, invOaxis, Os[i])
    #
    #     push!(crosssections, airfoil)
    #
    #     push!(strains, w_strain[i,:])
    #     push!(stress1, w_stresslayer[1,i,:])
    #     push!(stress2, w_stresslayer[2,i,:])
    #     push!(stressplot, w_stress[i,:])
    #     push!(buckle, w_c_bucklingplot[i,:])
    #     push!(point_datas, [j for j in npoints*(i-1)+1:npoints*i])
    # end
    #
    # # Generates cells in VTK Legacy format
    # points, vtk_cells, strain = VTKtools.multilines2vtkmulticells(crosssections, sections;point_datas=strains)
    # _, _, stresses1 = VTKtools.multilines2vtkmulticells(crosssections, sections;point_datas=stress1)
    # _, _, stresses2 = VTKtools.multilines2vtkmulticells(crosssections, sections;point_datas=stress2)
    # _, _, stresses6 = VTKtools.multilines2vtkmulticells(crosssections, sections;point_datas=stressplot)
    # _, _, bucklecon = VTKtools.multilines2vtkmulticells(crosssections, sections;point_datas=buckle)
    # points, vtk_cells, point_datas = VTKtools.multilines2vtkmulticells(crosssections, sections;point_datas=point_datas)
    #
    # # Formats the point data for generateVTK
    # data = []
    # push!(data, Dict(
    # "field_name" => "Strain",
    # "field_type" => "scalar",
    # "field_data" => strain))
    #
    # push!(data, Dict(
    # "field_name" => "uni_stress",
    # "field_type" => "scalar",
    # "field_data" => stresses1))
    #
    # push!(data, Dict(
    # "field_name" => "weave_stress",
    # "field_type" => "scalar",
    # "field_data" => stresses2))
    #
    # push!(data, Dict(
    # "field_name" => "stress",
    # "field_type" => "scalar",
    # "field_data" => stresses6))
    #
    # push!(data, Dict(
    # "field_name" => "Buckling_Constraint",
    # "field_type" => "scalar",
    # "field_data" => bucklecon))
    #
    # push!(data, Dict(
    # "field_name" => "PointIndex",
    # "field_type" => "scalar",
    # "field_data" => point_datas))
    #
    # # Generates the vtk file
    # VTKtools.generateVTK("$(file_name)wingL.$printiter", points; cells=vtk_cells, point_data=data, keep_points=false)
    #
    #
    # # #TEST vector output
    # # p0 = [0,0,0.0]
    # # p1 = [0,1,0.0]
    # # c0 = [p0,p1,p1,p0]
    # # d0 = [1.0,1.0,1.0,1.0]
    # # points, vtk_cells, lift = VTKtools.lines2vtk([c0];
    # #                                        values=[d0])
    # # data = []
    # # push!(data, Dict(
    # # "field_name" => "Strain",
    # # "field_type" => "vector",
    # # "field_data" => lift))
    # #
    # # VTKtools.generateVTK("$(file_name)lift.$printiter", points;
    # #          cells=vtk_cells, point_data=data)

    println("VTK Saved")

end

function generateprop!(p_strainx,p_strainy,p_strain,stresslayer,stress,p_c_bucklingplot,p_airfoilx,p_airfoily,x_pos,y_pos,z_pos,proploc,p_sweep_d,p_twist_d,propnum,bladenum,p_radii,azimuth_deg,file_name;plots = false)

    sections = fill(Tuple{Float64,Int64,Float64,Bool}[(1.0, 10, 1.0, false)],((length(p_radii)-1)))#*2+1))

    crosssections = []        # It will store here the cross sections for lofting
    strains = []
    stress1 = []
    stress2 = []
    stressplot = []
    buckle = []
    point_datas = []

    # Leading edge position of each airfoil
    Os = [ [x_pos[i]+proploc[1], y_pos[i]+proploc[2], z_pos[i]]+proploc[3] for i in 1:size(p_airfoilx)[1]]
    # Orientation of chord of each airfoil (yaw, pitch, roll)
    orien = zeros(length(p_twist_d),3) #twist already included
    for i = 1:length(orien[:,1])
        orien[i,:] = [-90.0,0.0,270.0]
    end

    # Processes each airfoil geometry - STRESS/STRAIN CALCS
    for i = 1:length(p_strainx[:,1])

        # Read airfoil file
        x = p_strainx[i,:]
        y = p_strainy[i,:]

        # Scales the airfoil acording to its chord length
        new_x = x#new_x #already scaled
        new_y = y#new_y

        # plot_airfoil(new_x, new_y; style=styles[i], label="airfoil")
        if plots
            figure("prop_airfoil")
            plot(new_x,new_y)
        end

        # Reformats into points
        npoints = size(new_x)[1]
        airfoil = Array{Float64, 1}[[new_x[j], new_y[j], 0] for j in 1:npoints]

        # Positions the airfoil along the blade in the right p_orientation
        Oaxis = VTKtools.rotation_matrix(orien[i,1], orien[i,2], orien[i,3])
        invOaxis = inv(Oaxis)
        airfoil = VTKtools.countertransform(airfoil, invOaxis, Os[i])

        spinaxis = VTKtools.rotation_matrix(0.0, 0.0, azimuth_deg)
        # invspinaxis = inv(spinaxis)
        airfoil = transform(airfoil, spinaxis, [proploc[1],proploc[2],proploc[3]])

        push!(crosssections, airfoil)

        push!(strains, p_strain[i,:])
        push!(stress1, stresslayer[1,i,:])
        push!(stress2, stresslayer[2,i,:])
        push!(stressplot, stress[i,:])
        push!(buckle, p_c_bucklingplot[i,:])
        push!(point_datas, [j for j in npoints*(i-1)+1:npoints*i])
    end

    # #generate other blade
    # Os = [ [x_pos[i]+proploc[1], -y_pos[i]+proploc[2], z_pos[i]]+proploc[3] for i in 1:size(p_airfoilx)[1]]
    # for i = length(p_strainx[:,1]):-1:1
    #
    #     # Read airfoil file
    #     x = p_strainx[i,:]
    #     y = p_strainy[i,:]
    #
    #     # Scales the airfoil acording to its chord length
    #     new_x = x#new_x #already scaled
    #     new_y = y#new_y
    #
    #     # plot_airfoil(new_x, new_y; style=styles[i], label="airfoil")
    #     if plots
    #         figure("prop_airfoil")
    #         plot(new_x,new_y)
    #     end
    #
    #     # Reformats into points
    #     npoints = size(new_x)[1]
    #     airfoil = Array{Float64, 1}[[-new_x[j], new_y[j], 0] for j in 1:npoints]
    #
    #     # Positions the airfoil along the blade in the right p_orientation
    #     Oaxis = VTKtools.rotation_matrix(orien[i,1], orien[i,2], orien[i,3])
    #     invOaxis = inv(Oaxis)
    #     airfoil = VTKtools.countertransform(airfoil, invOaxis, Os[i])
    #
    #     spinaxis = VTKtools.rotation_matrix(0.0, 0.0, 90.0)
    #     # invspinaxis = inv(spinaxis)
    #     airfoil = transform(airfoil, spinaxis, [proploc[1],proploc[2],proploc[3]])
    #
    #
    #     push!(crosssections, airfoil)
    #
    #     push!(strains, p_strain[i,:])
    #     push!(stress1, stresslayer[1,i,:])
    #     push!(stress2, stresslayer[2,i,:])
    #     push!(stressplot, stress[i,:])
    #     push!(buckle, p_c_bucklingplot[i,:])
    #     push!(point_datas, [j for j in npoints*(i-1)+1:npoints*i])
    # end

    # Generates cells in VTK Legacy format
    points, vtk_cells, strain = VTKtools.multilines2vtkmulticells(crosssections, sections;point_datas=strains)
    _, _, stresses1 = VTKtools.multilines2vtkmulticells(crosssections, sections;point_datas=stress1)
    _, _, stresses2 = VTKtools.multilines2vtkmulticells(crosssections, sections;point_datas=stress2)
    _, _, stresses6 = VTKtools.multilines2vtkmulticells(crosssections, sections;point_datas=stressplot)
    _, _, bucklecon = VTKtools.multilines2vtkmulticells(crosssections, sections;point_datas=buckle)
    _, _, point_datas = VTKtools.multilines2vtkmulticells(crosssections, sections;point_datas=point_datas)
    # Formats the point data for generateVTK
    data = []
    push!(data, Dict(
    "field_name" => "Strain",
    "field_type" => "scalar",
    "field_data" => strain))

    push!(data, Dict(
    "field_name" => "uni_stress",
    "field_type" => "scalar",
    "field_data" => stresses1))

    push!(data, Dict(
    "field_name" => "weave_stress",
    "field_type" => "scalar",
    "field_data" => stresses2))

    push!(data, Dict(
    "field_name" => "stress",
    "field_type" => "scalar",
    "field_data" => stresses6))

    push!(data, Dict(
    "field_name" => "Buckling_Constraint",
    "field_type" => "scalar",
    "field_data" => bucklecon))

    push!(data, Dict(
    "field_name" => "PointIndex",
    "field_type" => "scalar",
    "field_data" => point_datas))

    # Generates the vtk file
    VTKtools.generateVTK("$(file_name)prop$(propnum)bladenum$(bladenum).$printiter", points; cells=vtk_cells, point_data=data, keep_points=false)

end

function transform(V::Array{Float64,1},
    M::Array{Float64,2}, T::Array{Float64,1})
    return M*(V-T)+T
end

function transform(Vs::Array{Array{Float64,1},1},
    M::Array{Float64,2}, T::Array{Float64,1})
    out = Array{Float64,1}[]
    for V in Vs
        push!(out, transform(V, M, T))
    end
    return out
end

#Generate prop with all airfoil points
# y_pos = p_yloc#Rtip*p_radii # span position of each section
# x_pos = p_xloc#[y*tan(sweep[i]) for (i,y) in enumerate(y_pos)]
# z_pos = p_zloc#[y*tan(dihedral[i]) for (i,y) in enumerate(y_pos)]
# # PARAMETERS
# r = fill(1.0,length(p_radii))          # Expansion ratio in both surfaces of each airfoil
#
# sections = fill(Tuple{Float64,Int64,Float64,Bool}[(1.0, 10, 1.0, false)],length(p_radii)-1)
#
# # Leading edge position of each airfoil
# Os = [ [x_pos[i], y_pos[i], z_pos[i]] for i in 1:size(p_airfoilx)[1]]
# # Orientation of chord of each airfoil (yaw, pitch, roll)
# orien = zeros(length(p_twist_d),3) #twist already included
# for i = 1:length(orien[:,1])
#     orien[i,:] = [p_twist_d[i],p_sweep_d[i],270.0]
# end
#
#
# crosssections = []        # It will store here the cross sections for lofting
# point_datas = []
#
# # Processes each airfoil geometry
# for i = 1:length(p_airfoilx[:,1])
#
#     # Read airfoil file
#     x = p_airfoilx[i,:]
#     y = p_airfoily[i,:]
#
#     # Scales the airfoil acording to its chord length
#     new_x = chord[i]*x#new_x
#     new_y = chord[i]*y#new_y
#
#     # plot_airfoil(new_x, new_y; style=styles[i], label="airfoil")
#     if plots
#         figure("wing_airfoil")
#         plot(new_x,new_y)
#     end
#
#     # Reformats into points
#     npoints = size(new_x)[1]
#     airfoil = Array{Float64, 1}[[new_x[j], new_y[j], 0] for j in 1:npoints]
#
#     # Positions the airfoil along the blade in the right p_orientation
#     Oaxis = VTKtools.rotation_matrix(orien[i,1], orien[i,2], orien[i,3])
#     invOaxis = inv(Oaxis)
#     airfoil = VTKtools.countertransform(airfoil, invOaxis, Os[i])
#
#     push!(crosssections, airfoil)
#     push!(point_datas, [j for j in npoints*(i-1)+1:npoints*i])
#
# end
#
# # Generates cells in VTK Legacy format
# points, vtk_cells, pt_idx = VTKtools.multilines2vtkmulticells(crosssections, sections;
# point_datas=point_datas)
#
# data = []
# push!(data, Dict(
# "field_name" => "Point_index",
# "field_type" => "scalar",
# "field_data" => pt_idx
# )
# )
#
# # Generates the vtk file for geometry
# VTKtools.generateVTK("$(file_name)_geom$printiter", points; cells=vtk_cells, point_data=data, keep_points=false)


text:
#
# Squared distance function, considering Periodic Boundary Condition
#

function r2pbc(p1, p2, bdim)
    r2 = 0.
    for i = 1:2
        Δx = abs(p1[i] - p2[i]) % bdim
        if Δx > bdim/2
            Δx -= bdim
        end
        r2 += Δx^2
    end
    return r2
end


text:
mutable struct SystemBuildStats <: PowerSystemCaseBuilderType
    count::Int
    initial_construct_time::Float64
    serialize_time::Float64
    min_deserialize_time::Float64
    max_deserialize_time::Float64
    total_deserialize_time::Float64
end

function SystemBuildStats(;
    count,
    initial_construct_time,
    serialize_time,
    min_deserialize_time,
    max_deserialize_time,
    total_deserialize_time,
)
    return SystemBuildStats(
        count,
        initial_construct_time,
        serialize_time,
        min_deserialize_time,
        max_deserialize_time,
        total_deserialize_time,
    )
end

function SystemBuildStats(initial_construct_time::Float64, serialize_time::Float64)
    return SystemBuildStats(1, initial_construct_time, serialize_time, 0.0, 0.0, 0.0)
end

function update_stats!(stats::SystemBuildStats, deserialize_time::Float64)
    stats.count += 1
    if stats.min_deserialize_time == 0 || deserialize_time < stats.min_deserialize_time
        stats.min_deserialize_time = deserialize_time
    end
    if deserialize_time > stats.max_deserialize_time
        stats.max_deserialize_time = deserialize_time
    end
    stats.total_deserialize_time += deserialize_time
end

avg_deserialize_time(stats::SystemBuildStats) = stats.total_deserialize_time / stats.count

max_stars_repo_name: NREL-SIIP/PowerSystemBuilder

max_issues_repo_name: NREL-SIIP/PowerSystemBuilder

max_forks_repo_name: NREL-SIIP/PowerSystemBuilder


text:
import LIBSVM

function LIBSVMModel(
        ;
        single_label_levels::AbstractVector = [],
        name::AbstractString = "",
        isclassificationmodel::Bool = false,
        isregressionmodel::Bool = false,
        svmtype::Type=LIBSVM.SVC,
        kernel::LIBSVM.Kernel.KERNEL = LIBSVM.Kernel.RadialBasis,
        degree::Integer = 3,
        gamma::AbstractFloat = 0.1,
        coef0::AbstractFloat = 0.0,
        cost::AbstractFloat = 1.0,
        nu::AbstractFloat = 0.5,
        epsilon::AbstractFloat = 0.1,
        tolerance::AbstractFloat = 0.001,
        shrinking::Bool = true,
        weights::Union{Dict, Nothing} = nothing,
        cachesize::AbstractFloat = 100.0,
        verbose::Bool = true,
        )
    hyperparameters = Dict()
    hyperparameters[:svmtype] = svmtype
    hyperparameters[:kernel] = kernel
    hyperparameters[:degree] = degree
    hyperparameters[:gamma] = gamma
    hyperparameters[:coef0] = coef0
    hyperparameters[:cost] = cost
    hyperparameters[:nu] = nu
    hyperparameters[:epsilon] = epsilon
    hyperparameters[:tolerance] = tolerance
    hyperparameters[:shrinking] = shrinking
    weights = fix_type(weights)
    hyperparameters[:weights] = weights
    hyperparameters[:cachesize] = cachesize
    hyperparameters[:verbose] = verbose
    hyperparameters = fix_type(hyperparameters)
    underlyingsvm = FitNotYetRunUnderlyingObject()
    result = LIBSVMModel(
        name,
        isclassificationmodel,
        isregressionmodel,
        single_label_levels,
        hyperparameters,
        underlyingsvm,
        )
    return result
end

"""
"""
function get_underlying(
        x::LIBSVMModel;
        saving::Bool = false,
        loading::Bool = false,
        )
    result = x.underlyingsvm
    return result
end

"""
"""
function fit!(
        estimator::LIBSVMModel,
        featuresarray::AbstractArray,
        labelsarray::AbstractArray,
        )
    if estimator.isclassificationmodel && !estimator.isregressionmodel
        probability = true
    elseif !estimator.isclassificationmodel && estimator.isregressionmodel
        probability = false
    else
        error(
            "Could not figure out if model is classification or regression"
            )
    end
    @info(string("Starting to train LIBSVM model."))
    svm = try
        LIBSVM.svmtrain(
            featuresarray,
            labelsarray;
            probability = probability,
            estimator.hyperparameters...
            )
    catch e
        @warn(
            string(
                "While training LIBSVM model, ignored error: ",
                e,
                )
            )
        FitFailedUnderlyingObject()
    end
    # svm =
    @info(string("Finished training LIBSVM model."))
    estimator.underlyingsvm = svm
    estimator.levels = estimator.underlyingsvm.labels
    return estimator
end

"""
"""
function predict(
        estimator::LIBSVMModel,
        featuresarray::AbstractArray,
        )
    if estimator.isclassificationmodel && !estimator.isregressionmodel
        probabilitiesassoc = predict_proba(
            estimator,
            featuresarray,
            )
        predictionsvector = single_labelprobabilitiestopredictions(
            probabilitiesassoc
            )
        return predictionsvector
    elseif !estimator.isclassificationmodel && estimator.isregressionmodel
        if is_nothing(estimator.underlyingsvm)
            predicted_values = fill(Float64(0), size(featuresarray, 2))
        else
            predicted_values, decision_values = LIBSVM.svmpredict(
                estimator.underlyingsvm,
                featuresarray,
                )
            if !(typeof(predicted_values) <: AbstractVector)
                error("!(typeof(predicted_values) <: AbstractVector)")
            end
        end
        return predicted_values
    else
        error(
            "Could not figure out if model is classification or regression"
            )
    end
end

"""
"""
function predict(
        estimator::LIBSVMModel,
        featuresarray::AbstractArray,
        positive_class::AbstractString,
        threshold::AbstractFloat,
        )
    if estimator.isclassificationmodel && !estimator.isregressionmodel
        probabilitiesassoc = predict_proba(
            estimator,
            featuresarray,
            )
        predictionsvector = single_labelprobabilitiestopredictions(
            probabilitiesassoc,
            positive_class,
            threshold,
            )
        return predictionsvector
    else
        error(
            "Can only use the `threshold argument` with classification models"
            )
    end
end

"""
"""
function predict_proba(
        estimator::LIBSVMModel,
        featuresarray::AbstractArray,
        )
    if estimator.isclassificationmodel && !estimator.isregressionmodel
        if is_nothing(estimator.underlyingsvm)
            decision_values = fill(
                Int(0),
                size(featuresarray, 2),
                length(estimator.underlyingsvm.labels),
                )
            decision_values[:, 1] = 1
        else
            predicted_labels, decision_values =
                LIBSVM.svmpredict(estimator.underlyingsvm,featuresarray,)
            decision_values = transpose(decision_values)
        end
        result = Dict()
        for i = 1:length(estimator.underlyingsvm.labels)
            result[estimator.underlyingsvm.labels[i]] =
                decision_values[:, i]
        end
        result = fix_type(result)
        return result
    elseif !estimator.isclassificationmodel && estimator.isregressionmodel
        error("predict_proba is not defined for regression models")
    else
        error(
            "Could not figure out if model is classification or regression"
            )
    end
end

"""
"""
function single_labelmulticlassdataframesvmclassifier_LIBSVM(
        feature_names::AbstractVector,
        single_label_name::Symbol,
        single_label_levels::AbstractVector;
        name::AbstractString = "",
        svmtype::Type=LIBSVM.SVC,
        kernel::LIBSVM.Kernel.KERNEL = LIBSVM.Kernel.RadialBasis,
        degree::Integer = 3,
        gamma::AbstractFloat = 0.1,
        coef0::AbstractFloat = 0.0,
        cost::AbstractFloat = 1.0,
        nu::AbstractFloat = 0.5,
        epsilon::AbstractFloat = 0.1,
        tolerance::AbstractFloat = 0.001,
        shrinking::Bool = true,
        weights::Union{Dict, Nothing} = nothing,
        cachesize::AbstractFloat = 100.0,
        verbose::Bool = true,
        feature_contrasts::Union{Nothing, AbstractFeatureContrasts} =
            nothing,
        )
    dftransformer = DataFrame2LIBSVMTransformer(
        feature_names,
        single_label_name;
        levels = single_label_levels,
        )
    svmestimator = LIBSVMModel(
        ;
        name = name,
        single_label_levels = single_label_levels,
        isclassificationmodel = true,
        isregressionmodel = false,
        svmtype=svmtype,
        kernel = kernel,
        degree = degree,
        gamma = gamma,
        coef0 = coef0,
        cost = cost,
        nu = nu,
        epsilon = epsilon,
        tolerance = tolerance,
        shrinking = shrinking,
        weights = weights,
        cachesize = cachesize,
        verbose = verbose,
        )
    probapackager = ImmutablePackageSingleLabelPredictProbaTransformer(
        single_label_name,
        )
    predpackager = ImmutablePackageSingleLabelPredictionTransformer(
        single_label_name,
        )
    finalpipeline = dftransformer |>
                    svmestimator |>
                    probapackager |>
                    predpackager
    finalpipeline.name = name
    if !is_nothing(feature_contrasts)
        set_feature_contrasts!(finalpipeline, feature_contrasts)
    end
    return finalpipeline
end

"""
"""
function single_labelmulticlassdataframesvmclassifier(
        feature_names::AbstractVector,
        single_label_name::Symbol,
        single_label_levels::AbstractVector;
        package::Symbol = :none,
        name::AbstractString = "",
        svmtype::Type=LIBSVM.SVC,
        kernel::LIBSVM.Kernel.KERNEL = LIBSVM.Kernel.RadialBasis,
        degree::Integer = 3,
        gamma::AbstractFloat = 0.1,
        coef0::AbstractFloat = 0.0,
        cost::AbstractFloat = 1.0,
        nu::AbstractFloat = 0.5,
        epsilon::AbstractFloat = 0.1,
        tolerance::AbstractFloat = 0.001,
        shrinking::Bool = true,
        weights::Union{Dict, Nothing} = nothing,
        cachesize::AbstractFloat = 100.0,
        verbose::Bool = true,
        feature_contrasts::Union{Nothing, AbstractFeatureContrasts} = nothing,
        )
    if package == :LIBSVM
        result = single_labelmulticlassdataframesvmclassifier_LIBSVM(
            feature_names,
            single_label_name,
            single_label_levels;
            name = name,
            svmtype=svmtype,
            kernel = kernel,
            degree = degree,
            gamma = gamma,
            coef0 = coef0,
            cost = cost,
            nu = nu,
            epsilon = epsilon,
            tolerance = tolerance,
            shrinking = shrinking,
            weights = weights,
            cachesize = cachesize,
            verbose = verbose,
            feature_contrasts = feature_contrasts
        )
        return result
    else
        error("$(package) is not a valid value for package")
    end
end

"""
"""
function single_labeldataframesvmregression_LIBSVM(
        feature_names::AbstractVector,
        single_label_name::Symbol;
        name::AbstractString = "",
        svmtype::Type=LIBSVM.EpsilonSVR,
        kernel::LIBSVM.Kernel.KERNEL = LIBSVM.Kernel.RadialBasis,
        degree::Integer = 3,
        gamma::AbstractFloat = 0.1,
        coef0::AbstractFloat = 0.0,
        cost::AbstractFloat = 1.0,
        nu::AbstractFloat = 0.5,
        epsilon::AbstractFloat = 0.1,
        tolerance::AbstractFloat = 0.001,
        shrinking::Bool = true,
        weights::Union{Dict, Nothing} = nothing,
        cachesize::AbstractFloat = 100.0,
        verbose::Bool = true,
        feature_contrasts::Union{Nothing, AbstractFeatureContrasts} =
                nothing,
        )
    dftransformer = DataFrame2LIBSVMTransformer(
        feature_names,
        single_label_name,
        )
    svmestimator = LIBSVMModel(
        ;
        name = name,
        isclassificationmodel = false,
        isregressionmodel = true,
        svmtype=svmtype,
        kernel = kernel,
        degree = degree,
        gamma = gamma,
        coef0 = coef0,
        cost = cost,
        nu = nu,
        epsilon = epsilon,
        tolerance = tolerance,
        shrinking = shrinking,
        weights = weights,
        cachesize = cachesize,
        verbose = verbose,
        )
    predpackager = ImmutablePackageSingleLabelPredictionTransformer(
        single_label_name,
        )
    finalpipeline = dftransformer |>
                    svmestimator |>
                    predpackager
    finalpipeline.name = name
    if !is_nothing(feature_contrasts)
        set_feature_contrasts!(finalpipeline, feature_contrasts)
    end
    return finalpipeline
end

"""
"""
function single_labeldataframesvmregression(
        feature_names::AbstractVector,
        single_label_name::Symbol,
        ;
        package::Symbol = :none,
        name::AbstractString = "",
        svmtype::Type=LIBSVM.EpsilonSVR,
        kernel::LIBSVM.Kernel.KERNEL = LIBSVM.Kernel.RadialBasis,
        degree::Integer = 3,
        gamma::AbstractFloat = 0.1,
        coef0::AbstractFloat = 0.0,
        cost::AbstractFloat = 1.0,
        nu::AbstractFloat = 0.5,
        epsilon::AbstractFloat = 0.1,
        tolerance::AbstractFloat = 0.001,
        shrinking::Bool = true,
        weights::Union{Dict, Nothing} = nothing,
        cachesize::AbstractFloat = 100.0,
        verbose::Bool = true,
        feature_contrasts::Union{Nothing, AbstractFeatureContrasts} = nothing,
        )
    if package == :LIBSVM
        result = single_labeldataframesvmregression_LIBSVM(
            feature_names,
            single_label_name;
            name = name,
            svmtype=svmtype,
            kernel = kernel,
            degree = degree,
            gamma = gamma,
            coef0 = coef0,
            cost = cost,
            nu = nu,
            epsilon = epsilon,
            tolerance = tolerance,
            shrinking = shrinking,
            weights = weights,
            cachesize = cachesize,
            verbose = verbose,
            feature_contrasts = feature_contrasts
        )
        return result
    else
        error("$(package) is not a valid value for package")
    end
end


text:
using ..Tagging, ..IDS, ..Traits

export Member, nth

"The `id`th member of the family `f`"
struct Member{F, ID} <: AbstractVariable
  id::ID
  class::F
end

Base.:(==)(m1::Member, m2::Member) =
  m1.id == m2.id && m1.class == m2.class

Base.hash(x::Member) = hash(x.id)

@inline Var.recurse(x::Member, ω) = x.class(x.id, ω)

"""
`nth(f, id)`

`id`th member of sequence `(f_1, f_2, ..., f_n)`

Each element `f_i` is Conditionally independent of all other `f_j` given parents
"""
@inline nth(f, id) = Member(id, f)

"`id ~ f` is an alias for `nth(f, i)`"
@inline Base.:~(id, f) = nth(f, id)

@inline Base.:~(f) = Var.Class(f)

## Display
Base.show(io::IO, x::Member) = print(io, x.id,"@",x.class)

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
"""
    CharXRay

Represents a specific known characteristic X-ray in a specific element.
"""
struct CharXRay
    z::Int
    transition::Transition
    function CharXRay(z::Int, transition::Transition)
        @assert charactericXRayAvailable(
            z,
            transition.innershell.index,
            transition.outershell.index
        ) "$z $transition unknown!"
        "$(symbol(elements[z])) $(transition) does not occur."
        return new(z, transition)
    end
end

innerindex(cxr::CharXRay) = cxr.transition.innershell.index
outerindex(cxr::CharXRay) = cxr.transition.outershell.index

Base.hash(cxr::CharXRay, h::UInt)::UInt = hash(cxr.z, hash(cxr.transition, h))

Base.isequal(cxr1::CharXRay, cxr2::CharXRay) =
    isequal(cxr1.z, cxr2.z) && isequal(cxr1.transition, cxr2.transition)

Base.isless(cxr1::CharXRay, cxr2::CharXRay) =
    isequal(cxr1.z, cxr2.z) ? isless(cxr1.transition, cxr2.transition) : cxr1.z < cxr2.z

characteristic(elm::Element, tr::Transition) = CharXRay(z(elm), tr)

Base.show(io::IO, cxr::CharXRay) = print(io, elements[cxr.z].symbol, " ", cxr.transition)

"""
    inner(cxr::CharXRay)

Return the inner AtomicSubShell associated with the specified CharXRay.
"""
inner(cxr::CharXRay) = AtomicSubShell(cxr.z, cxr.transition.innershell)

"""
    outer(cxr::CharXRay)

Return the outer AtomicSubShell associated with the specified CharXRay.
"""
outer(cxr::CharXRay) = AtomicSubShell(cxr.z, cxr.transition.outershell)

"""
    transition(cxr::CharXRay)

Return the Transition structure associated with the specified CharXRay.
"""
transition(cxr::CharXRay) = cxr.transition

"""
    shell(cxr::CharXRay)

Return the shell, Shell(1), Shell(2), ... associated with the inner shell.
"""
shell(cxr::CharXRay) = shell(cxr.transition)

"""
    element(cxr::CharXRay)

Return the element for this CharXRay.
"""
element(cxr::CharXRay) = PeriodicTable.elements[cxr.z]

ν(cxr::CharXRay) = energy(cxr) / plancksConstant
ω(cxr::CharXRay) = 2π * ν(cxr)

"""
    λ(cxr::CharXRay)

X-ray wavelength in cm.
"""
λ(cxr::CharXRay) = hc / energy(cxr)
λ(energy::Real) = hc / energy
"""
    wavenumber(cxr::CharXRay)

X-ray wavenumber in cm¯¹.
"""
wavenumber(cxr::CharXRay) = 1.0 / λ(cxr)
wavenumber(energy::Real) = 1.0 / λ(energy)

abstract type WeightNormalization end
"""
`NormalizeByShell` normalizes the sum of all the weights associated with a shell to unity.
Example:

    sum(cxr=>weight(NormalizeByShell, cxr), characteristic(n"Fe", ltransitions))==1.0
"""
struct NormalizeByShell <: WeightNormalization end
"""
`NormalizeBySubShell` normalizes the sum of all the weights associated with a sub-shell to unity.

Example:

    sum(cxr=>weight(NormalizeBySubShell, cxr), characteristic(n"Fe", ltransitions))==1.0+1.0+1.0
"""
struct NormalizeBySubShell <: WeightNormalization end
"""
`NormalizeToUnity` normalizes intensities such that the most intense line in each shell to 1.0.

Example:

    sum(cxr=>weight(NormalizeBySubShell, cxr), n"Fe K-L3")==1.0
"""
struct NormalizeToUnity <: WeightNormalization end

"""
    weight(::Type{<:WeightNormalization}, cxr::CharXRay)

where `WeightNormalization` is one of the following:

  * `NormalizeByShell` normalizes the sum of all the weights associated with a shell to unity.
  * `NormalizeBySubShell` normalizes the sum of all the weights associated with each sub-shell to unity.
  * `NormalizeToUnity` normalizes intensities such that the most intense line in each shell to 1.0.

Computes a rough estimate of the relative intensity of `cxr` relative to the other characteristic
X-rays in its shell, sub-shell etc. The different `WeightNormalization` modes reflect different ways
that the `weight(...)` function is used.

The difference between `fluorescenceyield(...)` and `weight(...)` is that

  * fluorescenceyield assumes that a sub-shell in the atom is already ionized
  * weight also considers the relative likelihood of ionization of each sub-shell assuming an overvoltage of 4.0.
"""
function weight(::Type{NormalizeBySubShell}, cxr::CharXRay)
    inn, out = cxr.transition.innershell.index, cxr.transition.outershell.index
    return get(xrayweights()[cxr.z], (inn, inn, out), 0.0)[2]
end
function weight(::Type{NormalizeByShell}, cxr::CharXRay)
    inn, out = cxr.transition.innershell.index, cxr.transition.outershell.index
    return get(xrayweights()[cxr.z], (inn, inn, out), 0.0)[3]
end
function weight(::Type{NormalizeToUnity}, cxr::CharXRay)
    inn, out = cxr.transition.innershell.index, cxr.transition.outershell.index
    return get(xrayweights()[cxr.z], (inn, inn, out), 0.0)[4]
end

"""
    brightest(elm::Element, shell)

Return the brightest transition among the shell of transitions for the
specified element.  (group="K"|"Ka"|"Kb"|"L" etc)
"""
brightest(elm::Element, transitions) = brightest(characteristic(elm, transitions))

"""
    characteristic(elm::Element, iter::Tuple{Vararg{Transition}}, minweight=0.0, maxE=1.0e6)

"""
characteristic(
    elm::Element,
    iter::Tuple{Vararg{Transition}},
    minweight = 0.0,
    maxE = 1.0e6,
) = characteristic(
    elm,
    iter,
    cxr -> (weight(NormalizeToUnity, cxr) > minweight) && (energy(inner(cxr)) <= maxE),
)
characteristic(
    elm::Element,
    iter::AbstractVector{Transition},
    minweight = 0.0,
    maxE = 1.0e6,
) = characteristic(
    elm,
    iter,
    cxr -> (weight(NormalizeToUnity, cxr) > minweight) && (energy(inner(cxr)) <= maxE),
)

"""
    characteristic(elm::Element, iter::Tuple{Vararg{Transition}}, filterfunc::Function)
    characteristic(elm::Element, iter::AbstractVector{Transition}, filterfunc::Function)
    characteristic(elm::Element, iter::AbstractVector{Transition}, minweight = 0.0, maxE = 1.0e6)
    characteristic(elm::Element, iter::Tuple{Vararg{Transition}}, minweight = 0.0, maxE = 1.0e6)
    characteristic(ass::AtomicSubShell)

The collection of available CharXRay for the specified `Element` or `AtomicSubShell`.  `filterfunc(...)`
    * ` maxE` is compared to the edge energy.
    * `minWeight` is compared to the weight

Example:

    characteristic(n"Fe",ltransitions,0.01)
    characteristic(n"Fe",ltransitions,cxr->energy(cxr)>700.0)
    characteristic(n"Fe L3")
"""
function characteristic(
    elm::Element,
    iter::Tuple{Vararg{Transition}},
    filterfunc::Function,
)::Vector{CharXRay}
    res = CharXRay[]
    for tr in filter(t->has(elm, t), iter)
        cxr=characteristic(elm, tr)
        filterfunc(cxr) && push!(res, cxr)
    end
    res
end

characteristic(
    elm::Element,
    iter::AbstractVector{Transition},
    filterfunc::Function,
)::Vector{CharXRay} = map(
    tr -> characteristic(elm, tr),
    filter(tr -> has(elm, tr) && filterfunc(characteristic(elm, tr)), collect(iter)),
)

characteristic(ass::AtomicSubShell) =
    characteristic(element(ass), filter(tr -> inner(tr) == ass.subshell, alltransitions))

"""
    splitbyshell(cxrs)

Splits a collection of CharXRay into a dictionary where the key is the inner
AtomicSubShell and the values are a vector of CharXRay for that inner shell.
"""
function splitbyshell(cxrs)
    res = Dict{AtomicSubShell,Vector{CharXRay}}()
    for cxr in cxrs
        shell = inner(cxr)
        if haskey(res, shell)
            push!(res[shell], cxr)
        else
            res[shell] = [cxr]
        end
    end
    return res
end

brightest(cxrs::Vector{CharXRay})::CharXRay = cxrs[findmax(weight.(NormalizeToUnity, cxrs))[2]]

"""
    name(cxrs::AbstractVector{CharXRay}, byfamily=false)

An abbeviated name for a collection of CharXRay.
"""
function name(cxrs::AbstractVector{CharXRay}, byfamily::Bool=false)::String
    if byfamily
        if all(transition(cxr) in kalpha for cxr in cxrs)
            "Kα"
        elseif all(transition(cxr) in kbeta for cxr in cxrs)
            "Kβ"
        elseif all(transition(cxr) in ktransitions for cxr in cxrs)
            "K"
        elseif all(transition(cxr) in lalpha for cxr in cxrs)
            "Lα"
        elseif all(transition(cxr) in lbeta for cxr in cxrs)
            "Lβ"
        elseif all(transition(cxr) in ltransitions for cxr in cxrs)
            "L"
        elseif all(transition(cxr) in malpha for cxr in cxrs)
            "Mα"
        elseif all(transition(cxr) in mbeta for cxr in cxrs)
            "Mβ"
        elseif all(transition(cxr) in mtransitions for cxr in cxrs)
            "M"
        elseif all(transition(cxr) in ntransitions for cxr in cxrs)
            "N"
        else
            "Unknown"
        end
    else
        res = []
        elms = Set{Element}(element.(cxrs))
        for elm in elms
            for sh in Set(shell.(filter(cxr -> element(cxr) == elm, cxrs)))
                fc = filter(cxr -> (element(cxr) == elm) && (shell(cxr) == sh), cxrs)
                br, cx = brightest(fc), length(fc)
                other = cx > 2 ? "others" : "other"
                push!(res, cx > 1 ? "$(br) + $(cx-1) $(other)" : "$(br)")
            end
        end
        return join(res, ", ")
    end
end

function Base.show(io::IO, cxrs::AbstractVector{CharXRay})
    print(io, name(cxrs))
end

function NeXLUncertainties.asa(::Type{DataFrame}, cxrs::AbstractVector{CharXRay})
    cc = sort(cxrs)
    return DataFrame(
        XRay = cc,
        Inner = inner.(cc),
        Outer = outer.(cc),
        Energy = energy.(cc),
        WgtByShell = map(cxr->weight(NormalizeByShell, cxr), cc),
        WgtBySubShell = map(cxr->weight(NormalizeBySubShell, cxr), cc),
        Weight = map(cxr->weight(NormalizeToUnity, cxr), cc),
    )
end


text:
# data generation functions for motif finding

DNA_arr = ["A", "T", "C", "G"]

# generate a random DNA sequence of length n
genseq(n) = join([DNA_arr[i] for i in rand(1:4, 10)])

"""
Generate *count* sequences of length *len*, each with a motif of length
*motiflen* is inserted into a random index. *num_mutations* mutations are
added randomly to each motif in each sequence.

Return an array of sequences, the motif, and the indices of the inserted
motifs in each sequence.
"""
function gendata(count, len, motiflen, num_mutations)

    motif_arr = [DNA_arr[i] for i in rand(1:4, motiflen)]
    seqs_arr = [DNA_arr[i] for i in rand(1:4, count, len)]

    motif_ind = rand(1:len - motiflen + 1, count)

    # insert (possibly mutated) motifs into sequences
    for i = 1:count
        if num_mutations == 0
            seqs_arr[i, motif_ind[i]:motif_ind[i]+motiflen-1] = motif_arr
        else
            motif_arrᵢ = copy(motif_arr)
            for d = 1:num_mutations
                # add a random mutation to motif
                motif_arrᵢ[rand(1:motiflen)] = DNA_arr[rand(1:4)]

            end
            seqs_arr[i, motif_ind[i]:motif_ind[i]+motiflen-1] = motif_arrᵢ
        end
    end

    return [join(seqs_arr[i, :]) for i in 1:count], join(motif_arr), motif_ind
end

max_stars_repo_name: dillondaudert/JuliaMotifs

max_issues_repo_name: dillondaudert/JuliaMotifs

max_forks_repo_name: dillondaudert/JuliaMotifs


text:
### GROUP MEMBERS : Ashish Kumar, Pallavi Wats and Firat Melih Yilmaz

using Random
using LinearAlgebra
using BenchmarkTools
using Distributions
using JLD
using DataFrames
using CSV
using DelimitedFiles
using TypedTables
using Latexify
using LaTeXStrings
using FreqTables

############## Q1 ################
#1.(a)
Random.seed!(1234)
A = rand(Uniform(-1,2), 10,7)
B = rand(Normal(-2,15),10,7)
C = hcat(A[1:5,1:5],B[1:5,6:7])
D = similar(A)
for i in 1:7, j in 1:7
    if A[i,j] <= 0
        D[i,j]=A[i,j]
    else D[i,j]=0
    end
end

#1.(b)
#list the number of elements of the matrix
println(length(A))

#describe the size of the matrix
size(A)

#list of elements in A
for i in 1:10, j in 1:7
println(A[i,j])
end

#1.(c)
#number of unique elements in the matrix D
length(unique(D))

#list of unique elements of matrix D
vector_D = D[:];
println(unique(vector_D))

#1.(d)
#convert a matrix to a vector
#using reshape function
E = reshape(B,70)

#simpler way
E  = B[:]

#1.(e)
#creating a three-dimensional vector using A in the first column of the 3rd dimension and B in the second column of the 3rd dimension
F = cat(A,B,dims=3)

#1.(f)
#using permutedims to transform F from (10*7*2) to (2*10*7)
#F = permutedims(F,(2,10,7)) # not working

#using reshape
F = reshape(F,(2,10,7))

#1.(g)
#kronecker product of B and C
G = kron(A,B)

#kron(C,F)   ## kron does not work well for 3-D matrices

#1.(h)
#saving matrix A,B,C,D,E,F,G as a .jld file called matrixpractice

#saving all the matrices with their respective unique identifiers
save("/Users/prachi/Desktop/Econometrics6343/PS1/matrixpractice.jld","A",A,"B",B,"C",C,"D",D,"E",E,"F",F,"G",G)

#if I have to fetch only Matrix A
#load("/Users/prachi/Desktop/Econometrics6343/PS1/matrixpractice.jld","A")

#1.(i)
#save only the matrices A, B, C, and D as a .jld file called firstmatrix.
save("/Users/prachi/Desktop/Econometrics6343/PS1/firstmatrix.jld","A",A,"B",B,"C",C,"D",D)

#1.(j)
#export C as a .csv file called Cmatrix -- you will first need to transform C into a Dataframe
#converting to dataframe
C = convert(DataFrame,C)

#write a csv file
CSV.write("/Users/prachi/Desktop/Econometrics6343/PS1/Cmatrix.csv",C)

#1.(h)
#export D as a tab-delimited .dat file called Dmatrix -- you will first need to transform D into a DataFrame
D = convert(DataFrame,D)

#write a tab delimited .dat file
CSV.write("/Users/prachi/Desktop/Econometrics6343/PS1/Dmatrix.dat", D)

#1.(k)
#wrap it in function
function q1()
    Random.seed!(1234)
    A = rand(10,7) .* (10-(-5)).-5
    B = randn(10,7) .* (225-(-2)).-2
    C = hcat(A[1:5,1:5],B[1:5,6:7])
    D = similar(A)
    for i in 1:10, j = 1:7
        if A[i,j] <= 0
            D[i,j] = A[i,j]
        else
            D[i,j] = 0
        end
    end

    return A,B,C,D
end

A,B,C,D = q1()

############# Q2 #################
#2.(a)
#loop that computes element-by-element product of Matrix A and B, called AB
AB = zeros(10,7)
for i in 1:length(A)
    AB[i]=A[i]*B[i]
end

#without a loop element-by-element product of A and B
AB2 = A.*B

#2.(b)
#write a loop to create a column vector called Cprime which contains only elements of C that are between [-5,5]
Cprime = C[:]
for item in eachrow(C), i in item
    if -5<= i <= 5
        append!(Cprime,i)
    end
end

#do above without a loop
Cprime2 = C[-5 .≤ C .≤ 5];

#2.(c)
#use a loop to craete a 3D array X of dimension 15*169*5
X = zeros((15169,6,5))
for i in 1:size(X)[3], j in 1:size(X)[2]
    if j == 1
        X[:,j,i] .= 1
    elseif  j==2
        a2 = rand(Binomial(1,(0.75*(6-i))/5),15169)
        X[:,j,i] = a2
    elseif j==3
        a3 = rand(Normal(15+i-1,5*(i-1)),15169)
        X[:,j,i] = a3
    elseif j == 4
        a4 = rand(Normal((pi*(6-i))/3,1/ℯ),15169)
        X[:,j,i] = a4
    elseif  j == 5
        a5 = rand(Normal(12,2.19),15169)
        X[:,j,i] = a5
    elseif j == 6
        a6 = rand(Binomial(20,0.5),15169)
        X[:,j,i]  = a6
    end
end

#2.(d)
#use comprehension to create a matrix β (6*5) whose elements eveolve across time
K = 6
T = 5
β = zeros(K,T)
β[1,:] = [1:0.25:2;]
β[2,:] = [ log(x) for x in 1:T ]
β[3,:] = [-sqrt(x) for x in 1:T ]
β[4,:] = [ exp(x)-exp(x+1) for x in 1:T ]
β[5,:] = [x for x in 1:T]
β[6,:] = [x/3 for x in 1:T]

#2.(e)
#use comprehensions to create matrix Y (N*T) defined by ...

Y = zeros((15169,5,5))
for i in 1:5
    ϵ = rand(Normal(0,0.36),(15169,5))
    Y[:,:,i] = X[:,:,i]*β + ϵ
end

################ Q3 ######################
#3.(a)
#import nlsw88 as DataFrames
## first import nlsw as CSV file
nlsw = CSV.read("/Users/prachi/Desktop/Econometrics6343/PS1/nlsw88.csv")

#check if there are missing values
ismissing(nlsw)    #false -- no missing data

#convert nlsw88 into a DataFrame
df = convert(DataFrame, nlsw)

save("/Users/prachi/Desktop/Econometrics6343/PS1/nlsw88.jld", "df", df)

#make nullable
allowmissing!(df)

#describe data
describe(df)

#3.(b)
#what % of the sample have never been married?
m = freqtable(df, :never_married)

never_married = (sum(df["never_married"])/(sum(df["never_married"]) + sum(df["married"]))) * 100;
println("Percentage of the sample that has never been married: $never_married %")

#3.(c)
#use tabulate() to report what percentage of sample is in each race category
#using freqtable instead
race = prop(freqtable(df["race"]));

#3.(d)
#use the describe() function to create a matrix called summarystats which lists the mean, median, standard deviation, min, max, number of unique elements, and IQR
summarystats = describe(df);

#3.(e)
#show the joint distribution of industry and occupation using a cross-tabulation

joint_table = freqtable(df, :industry, :occupation);

#3.(g)
#wrap a function definition around all of the code for question 3
#call the function q3(). The function should have no inputs and no outputs
function q3()
    df = CSV.read("/Users/prachi/Desktop/Econometrics6343/PS1/nlsw88.csv");
    @save("nlsw88.jld", df)
    never_married = (sum(df["never_married"])/(sum(data["never_married"]) + sum(data["married"]))) * 100;
    race = prop(freqtable(df["race"]));
    summarystats = describe(df);
end

############## Q4 ###############
# Practice with functions

#4.(a)
#load firstmatrix.jld
load("/Users/prachi/Desktop/Econometrics6343/PS1/firstmatrix.jld")

#write a function called matrixops that takes as inputs the matrices A and B from question (a) of problem 1 and has three outputs: (i) the element-by-element product of the inputs,
# (ii) the product A′B, and (iii) the sum of all the elements of A + B.
#4.(d)
function matrixops(A, B)
    i = A .+ B
    ii = A' * B
    iii = sum(A + B)

    return i, ii, iii
end

matrixops(A,B);

#4.(c)
function matrixops(A, B)
    """
    Takes  matrix A and B as inputs and returns following outputs:

    i = A .+ B
    ii = A' * B
    iii = sum(A + B)
    """
    i = A .+ B
    ii = A' * B
    iii = sum(A + B)

    return i, ii, iii
end

#4.(e)
function matrixops(A, B)
    """
    Takes  matrix A and B as inputs and returns following outputs:

    i = A .+ B
    ii = A' * B
    iii = sum(A + B)
    """
    if size(A) != size(B)
        throw(DimensionMismatch("inputs must have the same size."))
    i = A .+ B
    ii = A' * B
    iii = sum(A + B)

    return i, ii, iii
end


text:
using FieldTraits, Base.Test

"`zap(x)` sends `x` to oblivion" zap(x) = nothing

@testset "UsageError" begin
    ue = FieldTraits.UsageError(zap, 2)
    @test ue.message == "`zap(x)` sends `x` to oblivion\n"
    @test ue.value == 2
    ue = FieldTraits.UsageError("a custom error message", 2)
    @test ue.message == "a custom error message"
end

@field FIscale
@field FIposition
@field FInotfield
@composed type FIcomposed
    FIscale
    FIposition
end
@testset "fieldindex" begin
    @test FieldTraits.fieldindex(FIcomposed, FIposition) === (Val{2}(),)
    @test FieldTraits.fieldindex(FIcomposed, FInotfield) === (Val{0}(),)
end


text:
type SobolDiagonalOrdering <: SobolOrdering end
type SobolStepsOrdering <: SobolOrdering end
type SobolFactorsOrdering <: SobolOrdering end

type SobolBrownianGenerator{O <: SobolOrdering} <: BrownianGenerator
  factors::Int
  steps::Int
  ordering::O
  generator::SobolInverseCumulativeRSG
  bridge::BrownianBridge
  lastStep::Int
  orderedIndices::Vector{Vector{Int}}
  bridgedVariates::Vector{Vector{Float64}}
end

function SobolBrownianGenerator(factors::Int, steps::Int, ordering::SobolOrdering, ::Int)
  generator = SobolInverseCumulativeRSG(factors * steps)
  bridge = BrownianBridge(steps)
  orderedIndices = Vector{Int}[Vector{Int}(steps) for _ = 1:factors]
  bridgedVariates = Vector{Float64}[Vector{Float64}(steps) for _ = 1:factors]

  fill_by_ordering!(ordering, orderedIndices, factors, steps)

  return SobolBrownianGenerator(factors, steps, ordering, generator, bridge, 0, orderedIndices, bridgedVariates)
end

function fill_by_ordering!(::SobolDiagonalOrdering, M::Vector{Vector{Int}}, factors::Int, steps::Int)
  # starting position of current diagonal
  i1 = j1 = 1

  # current position
  i = j = 1
  counter = 1
  while counter <= factors * steps
    M[i][j] = counter
    counter += 1
    if i == 1 || j == steps
      # we have completed a diagonal and have to start a new one
      if i1 < factors
        # we have to start path of the next factor
        i1 += 1
        j1 = 1
      else
        # we have to move along the path of the last factor
        i1 = factors
        j1 += 1
      end
      i = i1
      j = j1
    else
      # we move along the diagonal
      i = i - 1
      j = j + 1
    end
  end

  return M
end

function next_path!(sob::SobolBrownianGenerator)
  vals, wgt = next_sequence!(sob.generator)
  for i = 1:sob.factors
    transform!(sob.bridge, vals[sob.orderedIndices[i]], sob.bridgedVariates[i])
  end

  sob.lastStep = 1

  return wgt
end

function next_step!(sob::SobolBrownianGenerator, output::Vector{Float64})
  for i = 1:sob.factors
    output[i] = sob.bridgedVariates[i][sob.lastStep]
  end
  sob.lastStep += 1

  return 1.0
end

type SobolBrownianGeneratorFactory{O <: SobolOrdering} <: BrownianGeneratorFactory
  ordering::O
  seed::Int
end

create(sob::SobolBrownianGeneratorFactory, factors::Int, steps::Int) = SobolBrownianGenerator(factors, steps, sob.ordering, sob.seed)

# Clone #
clone(sob::SobolBrownianGenerator) = SobolBrownianGenerator(sob.factors, sob.steps, sob.ordering, sob.generator,
                                                          sob.bridge, sob.lastStep, deepcopy(sob.orderedIndices), deepcopy(sob.bridgedVariates))


text:
"""
    initial_evaluations(dd_function::DataFrame, limits::Tuple{T,T}, step_right_to_left::Bool, x_name::Symbol, output_for_envelope::Symbol; do_all_at_ends::Bool = false, additional_points_for_all_functions::AbstractArray) where T<:Real
This function does evaluations on the from edge of the limits. It will do any required evaluations on the to edge (or all if the optional parameter specified it).
It will additionally evaluate all functions at each of the points specified in the do_all_points optional parameter.
"""
function initial_evaluations(dd_function::DataFrame, limits::Tuple{T,T}, step_right_to_left::Bool, x_name::Symbol, output_for_envelope::Symbol; do_all_at_ends::Bool = false, additional_points_for_all_functions::AbstractArray = Array{T,1}([])) where T<:Real
    from_side = limits[1+convert(Int,step_right_to_left)]
    to_side = limits[2-convert(Int,step_right_to_left)]
    dd_eval = DataFrame()
    for i in 1:size(dd_function)[1]
        func_name = dd_function[i,:name]
        type = dd_function[i,:type]
        # Everything gets done on the from side
        # On the to side we need decreasings if we step right to left and increasings if we step left to right.
        output, dd_eval = get_function_value(dd_eval, dd_function, from_side, func_name, x_name, output_for_envelope)
        if ((type == :Decreasing) & step_right_to_left) | ((type == :Increasing) & (!step_right_to_left)) | do_all_at_ends
            output, dd_eval = get_function_value(dd_eval, dd_function, to_side, func_name, x_name, output_for_envelope)
        end
    end
    # Now doing additional optional points.
    for x in additional_points_for_all_functions
        for i in 1:size(dd_function)[1]
            func_name = dd_function[i,:name]
            output, dd_eval = get_function_value(dd_eval, dd_function, x, func_name, x_name, output_for_envelope)
        end
    end
    return dd_eval
end

function maximum_possible_value(target::Real, dd_function::DataFrame, dd_eval::DataFrame,
                                output_for_envelope::Symbol, x_name::Symbol)
    dic = Dict{Symbol,typeof(dd_eval[!,output_for_envelope][1])}()
    for i in 1:size(dd_function)[1]
        name = dd_function[i,:name]
        type = dd_function[i,:type]
        reduced_evals = dd_eval[dd_eval[!,:name] .== name,:]
        if type == :Constant
            dic[name] = dd_function[i,:func][String(output_for_envelope)]
        elseif type == :Increasing
            index     = searchsortedlast(reduced_evals[!,x_name], target)+1
            if index > length(reduced_evals[!,x_name])
                dic[name] = Inf
            else
                dic[name] = reduced_evals[index,output_for_envelope]
            end
            dic[name] = reduced_evals[index,output_for_envelope]
        elseif type == :Decreasing
            index     = searchsortedfirst(reduced_evals[!,x_name], target)-1
            if index < 1
                dic[name] = Inf
            else
                dic[name] = reduced_evals[index,output_for_envelope]
            end
        else
            dic[name] = Inf
        end
    end
    return dic
end

function get_func(dd_function::DataFrame, func_name::Symbol)
    return dd_function[dd_function[:name] .== func_name,:func][1]
end

function functions_potentially_above_mark(marc::Real, max_vals::Dict{Symbol,T}) where T<:Real
    result = Array{Symbol,1}()
    for k in keys(max_vals)
        if max_vals[k] > marc
            push!(result, k)
        end
    end
    return result
end

function top_function_at_point(dd_eval::DataFrame, x_name::Symbol, output_for_envelope::Symbol, point::Real, tol = 1e-14)
    dd_at_point = dd_eval[abs.(dd_eval[!,x_name] .- point) .< tol,:]
    len = size(dd_at_point)[1]
    if len == 0 error("No functions evaluated at this point") end
    return sort(dd_at_point, output_for_envelope)[len, :name]
end

"""
This crawls from one side to the other to get the top function value at each point in a grid.
"""
function crawler(dd_function::DataFrame, limits::Tuple{T,T}, stepp::R, x_name::Symbol, output_for_envelope::Symbol;
                 step_right_to_left::Bool = true, do_all_at_ends::Bool = false,additional_points_for_all_functions::AbstractArray = Array{promote_type(T,R),1}([])) where T<:Real where R<:Real
    grid = limits[1]:step:limits[2]
    if step_right_to_left grid = reverse(grid) end
    return crawler(dd_function, grid, x_name, output_for_envelope; do_all_at_ends = do_all_at_ends)
end
function crawler(dd_function::DataFrame, grid::Union{Array,StepRangeLen,UnitRange}, x_name::Symbol, output_for_envelope::Symbol; do_all_at_ends::Bool = false, additional_points_for_all_functions::AbstractArray = Array{typeof(grid[1]),1}([]))
    limits = (min(grid[1], grid[length(grid)]), max(grid[1], grid[length(grid)]))
    step_right_to_left = issorted(reverse(grid))
    if !step_right_to_left & !issorted(grid) error("The input grid must be strictly increasing or strictly decreasing") end
    dd_eval = initial_evaluations(dd_function, limits, step_right_to_left, x_name, output_for_envelope; do_all_at_ends = do_all_at_ends, additional_points_for_all_functions = additional_points_for_all_functions)
    top_function_at_last_point = top_function_at_point(dd_eval, x_name, output_for_envelope, grid[1])
    for i in 2:length(grid)
        target = grid[i]
        max_vals =  maximum_possible_value(target, dd_function, dd_eval, output_for_envelope, x_name)
        result, dd_eval = get_function_value(dd_eval, dd_function, target, top_function_at_last_point, x_name, output_for_envelope)
        #result, dd_eval = evaluate_function_and_add_to_dd(dd_eval, get_func(dd_function, top_function_at_last_point), target, top_function_at_last_point, x_name, output_for_envelope)
        other_functions_to_evaluate = functions_potentially_above_mark(result, max_vals)
        if length(other_functions_to_evaluate) > 0
            for f in other_functions_to_evaluate
                _, dd_eval = get_function_value(dd_eval, dd_function, target, f, x_name, output_for_envelope)
                #_, dd_eval = evaluate_function_and_add_to_dd(dd_eval, get_func(dd_function, f), target, f, x_name, output_for_envelope)
            end
        end
        top_function = top_function_at_point(dd_eval, x_name, output_for_envelope, grid[i])
        top_function_at_last_point = top_function
    end
    return dd_eval
end

# Intercept finders
"""
get_intercept_with_constant(dd_eval::DataFrame, func_name::Symbol, type::Symbol, constant::Real, dd_function::DataFrame, interval::Tuple{R,R}, x_name::Symbol, output_for_envelope::Symbol,
                                     bracketing_parameter::Real, max_interval_width::Real; recursion_count = 1, recursion_limit = 100) where R <:Real
Get an intercept between a function and a constant.
"""
function get_intercept_with_constant(dd_eval::DataFrame, func_name::Symbol, type::Symbol, constant::Real, dd_function::DataFrame, interval::Tuple{R,R}, x_name::Symbol, output_for_envelope::Symbol,
                                     bracketing_parameter::Real, max_interval_width::Real; recursion_count = 1, recursion_limit = 1000) where R <:Real
    if recursion_count > recursion_limit
        error("We have gone over the recursion limit in the function get_intercept_with_constant. This can happen if bracketing_parameter is too close to 1.")
    end
    @assert (type == :Increasing) | (type == :Decreasing) "No constant or unspecified functions can be used in the get intercept with constant function."
    # Making the first function
    dd_1 = dd_eval[dd_eval[!,:name] .== func_name , :]
    s1 =  Schumaker(dd_1[:,x_name], dd_1[:,output_for_envelope])
    # Finding roots
    roots = find_roots(s1; root_value = constant, interval = interval)[:roots]
    num_roots = length(roots)
    if num_roots == 0
        bson(string("C:\\Dropbox\\Stuart\\Julia_Library\\debug\\FunctionEnvelope_NoRootsInInterval.bson"),
             Dict(["dd_eval", "func_name", "type", "constant", "dd_function", "interval", "x_name", "output_for_envelope", "bracketing_parameter", "max_interval_width", "recursion_count", "recursion_limit"] .=>
                    [dd_eval, func_name, type, constant, dd_function, interval, x_name, output_for_envelope, bracketing_parameter, max_interval_width, recursion_count, recursion_limit]))
        # This is probably happening when we are rootfinding. Unfortunately we have got into an interval with no intercept. This can happen when our most recent evaluation changes the spline so a
        # root is narrowly excluded. We will just pick an edge.
        roots = [abs(s1(interval[1])) < abs(s1(interval[2])) ? interval[1] : interval[2]]
        return roots[1], dd_eval
        #error("No roots were found. Potentially more evaluations of this function are needed.")
    elseif num_roots > 1 # This can happen when a function is nearly flat which results in the schumaker spline adding some turns to it which makes it cross the constant more than once.
        # We will jusst select the first which should be fine.
        roots = [roots[1]]
        #bson(string("C:\\Dropbox\\Stuart\\Julia_Library\\debug\\FunctionEnvelope_MultipleRootsInInterval.bson"),
        #     Dict(["dd_eval", "func_name", "type", "constant", "dd_function", "interval", "x_name", "output_for_envelope", "bracketing_parameter", "max_interval_width", "recursion_count", "recursion_limit"] .=>
        #            [dd_eval, func_name, type, constant, dd_function, interval, x_name, output_for_envelope, bracketing_parameter, max_interval_width, recursion_count, recursion_limit]))
        #println("func_name = ", func_name)
        #println("type = ", type)
        #println("constant = ", constant)
        #error("This package does not support multiple roots in an interval.")
    end
    if abs(interval[2] - interval[1]) < max_interval_width
        return roots[1], dd_eval
    end
    # Finding point to try.
    point_to_try =  (1-bracketing_parameter) * (sum(interval)/2) + bracketing_parameter * (roots[1])
    # Now trying this point and saving it in dd_eval.
    result, dd_eval = get_function_value(dd_eval, dd_function, point_to_try, func_name, x_name, output_for_envelope)
    # Now using the result to make a smaller interval.
    if ((result > constant) & (type == :Increasing)) | ((result <= constant) & (type == :Decreasing))  # Intercept is left
        new_interval = Tuple{R,R}([interval[1], point_to_try])
        return get_intercept_with_constant(dd_eval, func_name, type, constant, dd_function, new_interval, x_name, output_for_envelope, bracketing_parameter, max_interval_width; recursion_count = recursion_count + 1, recursion_limit = recursion_limit)
    else
        new_interval = Tuple{R,R}([point_to_try, interval[2]])
        return get_intercept_with_constant(dd_eval, func_name, type, constant, dd_function, new_interval, x_name, output_for_envelope, bracketing_parameter, max_interval_width; recursion_count = recursion_count + 1, recursion_limit = recursion_limit)
    end
end

"""
get_intercept(dd_eval::DataFrame, increasing_function::Symbol, decreasing_function::Symbol, dd_function::DataFrame, interval::Tuple{R,R}, x_name::Symbol,
              output_for_envelope::Symbol, bracketing_parameter::Real, max_interval_width::Real; recursion_count = 1, recursion_limit = 100) where R <:Real
Get an intercept between two functions.
"""
function get_intercept(dd_eval::DataFrame, increasing_function::Symbol, decreasing_function::Symbol, dd_function::DataFrame, interval::Tuple{R,R}, x_name::Symbol,
                       output_for_envelope::Symbol, bracketing_parameter::Real, max_interval_width::Real; recursion_count = 1, recursion_limit = 100) where R <:Real
    if recursion_count > recursion_limit
        error("We have gone over the recursion limit in the function get_intercept. This can happen if bracketing_parameter is too close to 1.")
    end
    # Making the first function
    dd_1 = dd_eval[dd_eval[:,:name] .== increasing_function , :]
    s1 =  Schumaker(dd_1[:,x_name], dd_1[:,output_for_envelope])
    # Making the second function
    dd_2 = dd_eval[dd_eval[:,:name] .== decreasing_function , :]
    s2 =  Schumaker(dd_2[:,x_name], dd_2[:,output_for_envelope])
    # Finding roots
    roots = get_intersection_points(s1, s2)
    root_in_interval = roots[(roots .>= interval[1]-10*eps()) .& (roots .<= interval[2]+10*eps())]
    point_to_try = (sum(interval)/2)
    if length(root_in_interval) == 1 # 2 or 0 roots can happen sometimes if the spline is a bit inaccurate.
        point_to_try =  (1-bracketing_parameter) * (sum(interval)/2) + bracketing_parameter * (root_in_interval[1])
        if abs(interval[2] - interval[1]) < max_interval_width
            return root_in_interval[1], dd_eval
        end
    else
        if abs(interval[2] - interval[1]) < max_interval_width
            return (sum(interval)/2), dd_eval
        end
    end
    # Evaluation of function.
    result_inc, dd_eval = get_function_value(dd_eval, dd_function, point_to_try, increasing_function, x_name, output_for_envelope)
    #result1, dd_eval = evaluate_function_and_add_to_dd(dd_eval, get_func(dd_function, increasing_function), point_to_try, increasing_function, x_name, output_for_envelope)
    result_dec, dd_eval = get_function_value(dd_eval, dd_function, point_to_try, decreasing_function, x_name, output_for_envelope)
    #result2, dd_eval = evaluate_function_and_add_to_dd(dd_eval, get_func(dd_function, decreasing_function), point_to_try, decreasing_function, x_name, output_for_envelope)
    if result_inc > result_dec # We are to the right of the intercept as the increasing function is higher.
        new_interval = Tuple{R,R}([interval[1], point_to_try])
        return get_intercept(dd_eval, increasing_function, decreasing_function, dd_function, new_interval, x_name, output_for_envelope, bracketing_parameter, max_interval_width; recursion_count = recursion_count + 1, recursion_limit = recursion_limit)
    else
        new_interval = Tuple{R,R}([point_to_try, interval[2]])
        return get_intercept(dd_eval, increasing_function, decreasing_function, dd_function, new_interval, x_name, output_for_envelope, bracketing_parameter, max_interval_width; recursion_count = recursion_count + 1, recursion_limit = recursion_limit)
    end
end

"""
    get_upper_envelope(dd_eval::DataFrame, dd_function::DataFrame, x_name::Symbol, output_for_envelope::Symbol, bracketing_parameter::Real, max_interval_width::Real;
                             grid::Array{R,1} = sort(unique(dd_eval[x_name]))) where R<:Real
"""
function get_upper_envelope(dd_eval::DataFrame, dd_function::DataFrame, x_name::Symbol, output_for_envelope::Symbol, bracketing_parameter::Real, max_interval_width::Real;
                                      grid::Union{StepRangeLen,Array{R,1}} = sort(unique(dd_eval[x_name]))) where R<:Real
    interval_dd = DataFrame(interval_start = grid[1], interval_end = NaN, func = top_function_at_point(dd_eval, x_name, output_for_envelope, grid[1]))
    i = 1 # This is the row index for interval_dd
    for j in 2:length(grid)
        # We need to ensure that the previous function has been evaluated at this point. Otherwise we might get a dodgy "change"
        _, dd_eval = get_function_value(dd_eval, dd_function, grid[j], interval_dd[i, :func], x_name, output_for_envelope)
        top_func = top_function_at_point(dd_eval, x_name, output_for_envelope, grid[j])
        if top_func != interval_dd[i, :func]
            old_top_func = interval_dd[i, :func]
            old_func_type = dd_function[dd_function[:,:name] .== old_top_func,:type][1]
            top_func_type = dd_function[dd_function[:,:name] .== top_func    ,:type][1]
            intercept = NaN
            interval = Tuple([grid[j-1], grid[j]])
            if (old_func_type == :Constant) && (top_func_type == :Constant)
                error("This does not make sense. How did they cross?")
            elseif (old_func_type == :Constant)
                old_func_const = dd_function[dd_function[:,:name] .== old_top_func,:func][1][String(output_for_envelope)]
                intercept, dd_eval = get_intercept_with_constant(dd_eval, top_func, top_func_type, old_func_const, dd_function, interval, x_name, output_for_envelope, bracketing_parameter, max_interval_width)
            elseif (top_func_type == :Constant)
                top_func_const = dd_function[dd_function[:,:name] .== top_func,:func][1][String(output_for_envelope)]
                intercept, dd_eval = get_intercept_with_constant(dd_eval, old_top_func, old_func_type, top_func_const, dd_function, interval, x_name, output_for_envelope, bracketing_parameter, max_interval_width)
            else
                intercept, dd_eval = get_intercept(dd_eval, top_func, old_top_func, dd_function, interval, x_name, output_for_envelope, bracketing_parameter, max_interval_width)
            end
            interval_dd[i, :interval_end] = intercept
            i+=1
            new_interval = DataFrame(interval_start = intercept, interval_end = NaN, func = top_func)
            interval_dd = vcat(interval_dd, new_interval)
        end
    end
    interval_dd[i,:interval_end] = grid[length(grid)]
    interval_dd[!,:length] .= interval_dd[!,:interval_end] .- interval_dd[!,:interval_start]
    return interval_dd, dd_eval
end


text:
"""
run benchmark measurements
"""
function runbench(f::Function,s,t, b::Benchmark, stats::Bool)
    #print(string(f),": ")
    sf = zeros(Int,Int(typemax(StatsFields)))
    starttime =  time_ns()
    pattern = f(t)
    inittime =  time_ns()
    pos = 0
    if stats
        while (pos = f(s,pattern,pos+1,sf))>0 end
    else
        while (pos = f(s,pattern,pos+1,nothing))>0 end
    end
    kerneltime = time_ns()
    if sizeof(t)>3 # do not report pattern size up to 3
        sf[Int(SFinit)] = (inittime - starttime)%Int
        sf[Int(SFkernel)] = (kerneltime-inittime)%Int
        # patternsize 2 is 1st run, ignore it (may include compilation)
        record(b,f,sizeof(t),sf)
    end
    nothing
end

"""
    benchmark(patternsize::Int, textsize::Int, source::String, functions::Vector{Function}, result::String, stats::Bool)

patternsize: largest size of search pattern in benchmark loop.
Loop runs from 2 (not measured) up to patternsize

textsize: length of the byte sequence to search in. Source must supply al least
so many bytes

source: either an UInt8 and an Int literal, separated by comma, or a file path.
If source starts with "0x" benchmark expects an UInt8 and Int litaral. The UInt8
specifies the highest value in the alphabet for a randomly generated search
sequence, the Int gives a random seed.

functions: the list of search functions to include in the benchmark

result: a file path. benchmarks writes files <result>.png and <result>.csv.
The former file shows the search kernel timings, the latter saves all collected
statistics in CSV format, including search timings.

stats: if true, additional measures are talen, like the number of main loop
cycles or the number of successful bloom skips. Thanks to multiple dispatch,
code collecting statistics is optimized away if stats==false. Because we are
benchmarking very lowlevel ans simple operations, code collecting statistics
can add overhead which results in senseless benchmark timings. Use stats==true
mostly in the design and anaysis phase of a search function, and switch stats
off for final benchmarking.

"""
Base.@propagate_inbounds function benchmark(patternsize::Int, textsize::Int, source::String, functions::Vector{Function}, result::String, stats::Bool)
    b = Benchmark(patternsize,source,result)
    s = Vector{UInt8}(undef,textsize)
    if startswith(source,"0x")
        # generate random sequence
        args = split(source,',')
        alphabetsize = parse(UInt8,args[1])
        seed = parse(Int,args[2])
        Random.seed!(seed)
        Random.rand!(s,0x00:alphabetsize) # array to search in
    else
        # read file
        f = open(source)
        read!(f,s)
        close(f)
    end
    #ss = String(copy(s))
    for j in 2:patternsize
        tj = s[textsize-patternsize+1:textsize-patternsize+j]
        #println()
        #println("benchmark with pattern size $j")
        for f in functions
            runbench(f,s,tj,b,stats)
        end
    end
    io = open(result * ".csv", "w")
    write(io,b)
    close(io)

    #= initializing timings need more repititions to be statistically valid
    plotinit =plot()
    title!("$result: initialization timing")
    for f in functions
        plot!(line(b,f,SFinit), label=string(f))
    end
    png(result*"_init")
    =#
    plotkernel =plot()
    # remove (windows) path part from result for title to keep plot title intact
    title = result
    endpath = findlast("\\",result)
    if endpath != nothing
        title = result[endpath+1:end]
    end
    title!(title*": kernel timing")
    for f in functions
        plot!(line(b,f,SFkernel), label=string(f))
    end
    png(result)
    display(plotkernel)
    nothing
end

#benchmark(128,3,100,"r:\\testfile.txt")


text:
module BFZ

using CompEcon: qnwnorm
using Interpolations

type BFZModel
    A::Float64
    B::Float64
    vbar::Float64
    lgbar::Float64
    ρ::Float64
    α::Float64
    β::Float64
    η::Float64
    ν::Float64
    δ::Float64
    ϵ::Vector{Float64}
    Π::Vector{Float64}
    x::AbstractVector
    k::AbstractVector
    xp::Matrix{Float64}
    gp::Matrix{Float64}
end

function BFZModel(;ρ=-1.0, α=-9.0, β=0.99, η=1/3, ν=0.1, δ=0.025, lgbar=0.004,
                   A=0.0, B=1.0, vbar=0.015^2, nk=35, nx=15, nϵ1=5)

    # quadrature weights and nodes
    ϵ, Π = qnwnorm(nϵ1, 0.0, 1.0)

    # grids for x and k
    x = linspace(-0.05, 0.05, nx)
    k = linspace(15.0, 30.0, nk)

    # grids for x_{t+1} and g_{t+1} -- apply law of motion
    xp = Array(Float64, nϵ1, nx)
    gp = Array(Float64, nϵ1, nx)

    for i=1:nx
        _x = x[i]
        for j=1:nϵ1
            xp[j, i] = A*_x + B*sqrt(vbar)*ϵ[j]
            gp[j, i] = exp(lgbar + xp[j, i])
        end

    end

    BFZModel(A, B, vbar, lgbar, ρ, α, β, η, ν, δ, ϵ, Π, x, k, xp, gp)
end

## Helpful functions
_Y(m::BFZModel, k) = (m.η .* k.^m.ν + (1 - m.η)).^(1/m.ν)

"""
Initializes the guess of the value function on the grid to be a small fraction
of consuming production each period. Use the same value for all x
"""
function initialize(m::BFZModel)
    JT = (((1-m.β) * (_Y(m, m.k) + (1-m.δ)*m.k).^(m.ρ)).^(1/m.ρ)) ./ 500
    Jvals = repeat(JT, inner=[1, length(m.x)])
    itp = scale(interpolate(Jvals, BSpline(Quadratic(Line())), OnGrid()),
                m.k, m.x)
    Jvals, itp
end

function envelope_method(m::BFZModel;maxiter=500, tol=1e-4)
    # simplify notation
    β, δ, ρ, η, ν, α = m.β, m.δ, m.ρ, m.η, m.ν, m.α
    nshocks = length(m.ϵ)
    nk, nx = length(m.k), length(m.x)

    # initialize interpolant for value function
    Jvals, itp = initialize(m)
    χpolicy = Array(Float64, nk, nx)
    Jupd = similar(χpolicy)

    # intermediate helper variables
    err = 100.0
    iter = 0
    start_time = time()

    # Solve this bad boy
    while err>tol && iter<maxiter

        iter += 1
        # Given the current vf, find the χ that solves the env condition
        for i_x=1:nx
            # Pull out current x
            xt = m.x[i_x]

            for i_k=1:nk
                # Pull out current k
                kt = m.k[i_k]

                yt = _Y(m, kt)  # production

                dJ_k = gradient(itp, kt, xt)[1]  # evaluate dJ/dk
                dJ_k = dJ_k < 1e-12 ? 0. : dJ_k

                Jt = itp[kt, xt] # evaluate J

                # use envelope and budget constraint to get c and χ, then save χ
                # in policy matrix
                ct = (dJ_k/((1-β)*Jt^(1-ρ)*(yt^(1-ν)*η*kt^(ν-1)+1-δ)))^(1./(ρ-1))
                χt = yt + (1 - δ)*kt - ct
                χpolicy[i_k, i_x] = χt

                # calculate the certainty equivalent
                expterm = 0.
                for i_eps=1:nshocks
                    xp = m.xp[i_eps, i_x]
                    gp = m.gp[i_eps, i_x]
                    jp = itp[χt/gp, xp]
                    expterm += m.Π[i_eps] * (jp * gp)^α
                end
                μ = expterm^(1/α)

                # evaluate value function
                Jupd[i_k, i_x] = ((1 - β)*ct^ρ + β*μ^ρ)^(1/ρ)

            end
        end

        err = maxabs(Jvals - Jupd)
        copy!(Jvals, Jupd)
        itp = scale(interpolate(Jupd, BSpline(Quadratic(Line())), OnGrid()),
                    m.k, m.x)
        @show iter, err, time() - start_time

    end

    return itp
end

end  # module

max_stars_repo_name: NYUEcon/GrowthModels

max_issues_repo_name: NYUEcon/GrowthModels

max_forks_repo_name: NYUEcon/GrowthModels


text:
# Path planning simulation by reinforcement learning
# learning method: q learning

module AnimeQLearning
  using Plots
  pyplot()

  include(joinpath(split(@__FILE__, "src")[1], "src/model/world/puddle_world.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/model/map/map.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/model/robot/warp_robot/warp_robot.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/localization/extended_kalman_filter/extended_kalman_filter.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/model/goal/goal.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/model/puddle/puddle.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/model/agent/q_agent.jl"))

  function main(delta_time=0.1, end_time=20000; is_test=false)
    # save path of gif file
    path = "src/decision_making/reinforcement_learning/q_learning/anime_q_learning.gif"

    # simulation world
    world = PuddleWorld(-5.0, 5.0, -5.0, 5.0,
                        delta_time=delta_time,
                        end_time=end_time,
                        is_test=is_test,
                        save_path=path)

    # map including landmarks
    map = Map()
    landmark_positions = [[-4.0, 2.0],[2.0, -3.0],[4.0, 4.0],[-4.0, -4.0]]
    for (idx, pos) in enumerate(landmark_positions)
      add_object(map, Object(pos[1], pos[2], id=idx))
    end
    append(world, map)

    # goal
    goal = Goal(-3.0, -3.0)
    append(world, goal)

    # puddles
    append(world, Puddle([-2.0, 0.0], [0.0, 2.0], 0.1))
    append(world, Puddle([-0.5, -2.0], [2.5, 1.0], 0.1))

    # robot including sensor, agent, estimator
    init_pose = [3.0, 3.0, 0.0]
    sensor = Sensor(map, dist_noise_rate=0.1, dir_noise=pi/90)
    ekf = ExtendedKalmanFilter(init_pose, env_map=map,
                               dist_dev_rate=0.14,
                               dir_dev=0.05)
    agent = QAgent(delta_time=delta_time, estimator=ekf)
    robot = WarpRobot(init_pose, 0.2, "red",
                      agent, delta_time,
                      noise_per_meter=5,
                      noise_std=pi/30,
                      bias_rate_stds=[0.1, 0.1],
                      sensor=sensor)
    append(world, robot)

    draw(world)
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy