text:
@with_kw mutable struct CoVParameters
    β::Float64 = 2.5/3.6
    γ::Float64 = 1/3.6
    σ::Float64 = 1/2.
    δ::Float64 = 0.05#Proportion of symptomatic/diseased vs non-symptomatic cases
    τ::Float64 = 1/15. #treatment/isolation rate for symptomatics
    μ₁::Float64 = 0.01#Excess mortality due to disease
    ϵ::Float64 = 0.1 #Relative infectiousness of undetectable infecteds
    ρ::Vector{Float64} = [0.01 for i in 1:n] #spatial coupling
    T::Matrix{Float64} = zeros(n,n)#transmission matrix
    ext_inf_rate::Float64 = 0. #Scales the contact rate with the infecteds arriving via air
    into_mom::Vector{Int}#Number of people flying into Mombassa each day
    into_nai::Vector{Int}#Number of people flying into Nairobi each day
    global_prev::Vector{Float64}
    dt::Float64 = 1. #Useful for the non-negative method
    Î::Vector{Float64} = zeros(n) #For inplace calculations
    N̂::Vector{Float64} = zeros(n)#For inplace calculations
    λ_urb::Vector{Float64} = zeros(n)#For inplace calculations
    λ_rur::Vector{Float64} = zeros(n)#For inplace calculations
    dN::Vector{Int64} = zeros(Int64,n*n_t)#For inplace calculations
    poi_rates::Vector{Float64} = zeros(n*n_t)#For inplace calculations
    dc::SparseMatrixCSC{Int64,Int64} = sparse(zeros(Int64,n*n_s*2,n*n_t))#For inplace calculations
end

"""
    mutable struct CoVParameters_AS

Struct for containing relevant epidemilogical parameters for the age-structured version of KenyaCoV
"""
@with_kw mutable struct CoVParameters_AS
    #Epidemiological parameters and social contact/mixing rates
    β::Float64 = 1. #Basic transmission probability per contact OR infectious contact rate (can be greater than one depending on scaling)
    c_t::Function = t -> 1. #Time varying basic contact rate
    γ::Float64 = 1/2. #recovery rate for mild and asymptomatic cases (Consensus estimate mean 4 days infectious, default is 1/σ₂ + 1/γ = 4 days)
    σ₁::Float64 = 1/5. #end of latency rate (Consensus estimate mean 5 days)
    σ₂::Float64 = 1/2. #end of pre-symptomatic rate (Consensus estimate mean 4 days infectious, default is 1/σ₂ + 1/γ = 4 days)
    δ::Float64 = 0.9#Proportion of over 80s who get identified
    τ::Float64 = 1/5. #Rate of hospitalisation treatment, conditional on eventually needing it (V category)
    τ_initial::Float64 = 0. # isolation rate for symptomatics at beginning of epidemic
    clear_quarantine = 0. # Average time to end isolation NOT USED IN THIS VERSION
    μ₁::Float64 = 0.0#Excess mortality due to disease NOT USED IN THIS VERSION
    ϵ::Float64 = 0.1 #Relative infectiousness of undetectable/undetected infecteds both pre-symptomatic and asymptomatic
    ϵ_D::Float64 = 1.#Relative infectiousness of mild infecteds due to social avoidance
    ϵ_V::Float64 = 0.4#Relative infectiousness of mild, then severe, infecteds due to social avoidance
    rel_detection_rate::Vector{Float64} = ones(n_a) #relative symptomatic rate
    hₐ::Vector{Float64} = zeros(n_a)   #proportion of severe cases if symptomatic
    ICUₐ::Vector{Float64} = zeros(n_a) #proportion of hospitalised cases that become critical
    χ::Vector{Float64} = ones(n_a) #relative susceptibility
    ρ::Vector{Float64} = [0.01 for i in 1:n] #Time spent outside of area
    T::Matrix{Float64} = zeros(n,n)#Probability distributions of location for mobile individuals
    M::Matrix{Float64} = zeros(n_a,n_a) #Age mixing matrix
    M_ho::Matrix{Float64} = zeros(n_a,n_a) #Age mixing matrix --- if only home contacts are made
    ext_inf_rate::Float64 = 0. #Scales the contact rate with the infecteds arriving via air
    into_mom::Vector{Int}#Number of people flying into Mombassa each day
    into_nai::Vector{Int}#Number of people flying into Nairobi each day
    global_prev::Vector{Float64}#Chance that each person arriving was infected on that day
    #Control variables
    isolating_detecteds::Bool = false #This determines if people are still being isolated
    lockdown::Bool = false  #This determines if social distancing and travel restrictions are still in force
    #Calculation variables
    dt::Float64 = 1. #Useful for the non-negative method
    Î::Matrix{Float64} = zeros(n_wa,n_a) #For inplace calculations
    N̂::Matrix{Float64} = zeros(n_wa,n_a)#For inplace calculations
    λ::Matrix{Float64} = zeros(n_wa,n_a)#For inplace calculations
    λ_loc::Matrix{Float64} = zeros(n_wa,n_a)#For inplace calculations
    dN::Vector{Int64} = zeros(Int64,n_wa*n_a*n_ta)#For inplace calculations
    poi_rates::Vector{Float64} = zeros(n_wa*n_a*n_ta)#For inplace calculations
    dc::SparseMatrixCSC{Int64,Int64} = sparse(zeros(Int64,n_wa*n_a*n_s,n_wa*n_a*n_ta))#For inplace calculations
    du_linear::Vector{Int64} = zeros(Int64,n_wa*n_a*n_s) #for inplace calculations
end


text:
const minfloat = realmin(Float64)
const oneMeps  = 1. - eps()
const llmaxabs = log(-log(minfloat))
const logeps   = log(eps())
abstract Link                           # Link types define linkfun, linkinv, mueta,
                                        # valideta and validmu.

chkpositive(x::Real) = isfinite(x) && 0. < x ? x : error("argument must be positive")
chkfinite(x::Real) = isfinite(x) ? x : error("argument must be finite")
clamp01(x::Real) = clamp(x, minfloat, oneMeps)
chk01(x::Real) = 0. < x < 1. ? x : error("argument must be in (0,1)")

type CauchitLink  <: Link end
linkfun (l::CauchitLink,   mu::Real) = tan(pi * (mu - 0.5))
linkinv (l::CauchitLink,  eta::Real) = 0.5 + atan(eta) / pi
mueta   (l::CauchitLink,  eta::Real) = 1. /(pi * (1 + eta * eta))
valideta(l::CauchitLink,  eta::Real) = chkfinite(eta)
validmu (l::CauchitLink,   mu::Real) = chk01(mu)

type CloglogLink  <: Link end
linkfun (l::CloglogLink,   mu::Real) = log(-log(1. - mu))
linkinv (l::CloglogLink,  eta::Real) = -expm1(-exp(eta))
mueta   (l::CloglogLink,  eta::Real) = exp(eta) * exp(-exp(eta))
valideta(l::CloglogLink,  eta::Real) = abs(eta) < llmaxabs? eta: error("require abs(eta) < $llmaxabs")
validmu (l::CloglogLink,   mu::Real) = chk01(mu)

type IdentityLink <: Link end
linkfun (l::IdentityLink,  mu::Real) = mu
linkinv (l::IdentityLink, eta::Real) = eta
mueta   (l::IdentityLink, eta::Real) = 1.
valideta(l::IdentityLink, eta::Real) = chkfinite(eta)
validmu (l::IdentityLink,  mu::Real) = chkfinite(mu)

type InverseLink  <: Link end
linkfun (l::InverseLink,   mu::Real) =  1. / mu
linkinv (l::InverseLink,  eta::Real) =  1. / eta
mueta   (l::InverseLink,  eta::Real) = -1. / (eta * eta)
valideta(l::InverseLink,  eta::Real) = chkpositive(eta)
validmu (l::InverseLink,  eta::Real) = chkpositive(mu)

type LogitLink    <: Link end
linkfun (l::LogitLink,     mu::Real) = log(mu / (1 - mu))
linkinv (l::LogitLink,    eta::Real) = 1. / (1. + exp(-eta))
mueta   (l::LogitLink,    eta::Real) = (e = exp(-abs(eta)); f = 1. + e; e / (f * f))
valideta(l::LogitLink,    eta::Real) = chkfinite(eta)
validmu (l::LogitLink,     mu::Real) = chk01(mu)

type LogLink      <: Link end
linkfun (l::LogLink,       mu::Real) = log(mu)
linkinv (l::LogLink,      eta::Real) = exp(eta)
mueta   (l::LogLink,      eta::Real) = eta < logeps ? eps() : exp(eta)
valideta(l::LogLink,      eta::Real) = chkfinite(eta)
validmu (l::LogLink,       mu::Real) = chkpositive(mu)

type ProbitLink   <: Link end
linkfun (l::ProbitLink,    mu::Real) = ccall((:qnorm5, Rmath), Float64,
                                             (Float64,Float64,Float64,Int32,Int32),
                                             mu, 0., 1., 1, 0)
linkinv (l::ProbitLink,   eta::Real) = (1. + erf(eta/sqrt(2.))) / 2.
mueta   (l::ProbitLink,   eta::Real) = exp(-0.5eta^2) / sqrt(2.pi)
valideta(l::ProbitLink,   eta::Real) = chkfinite(eta)
validmu (l::ProbitLink,    mu::Real) = chk01(mu)
                                        # Vectorized methods, including validity checks
function linkfun{T<:Real}(l::Link, mu::AbstractArray{T,1})
    [linkfun(l, validmu(l, m)) for m in mu]
end

function linkinv{T<:Real}(l::Link, eta::AbstractArray{T,1})
    [linkinv(l, valideta(l, et)) for et in eta]
end

function mueta{T<:Real}(l::Link, eta::AbstractArray{T,1})
    [mueta(l, valideta(l, et)) for et in eta]
end

canonicallink(d::Gamma)     = InverseLink()
canonicallink(d::Normal)    = IdentityLink()
canonicallink(d::Bernoulli) = LogitLink()
canonicallink(d::Poisson)   = LogLink()

ext: jl


text:
# Working with Tasks

*Warning.* The formal task API described here may be depreciated in
the future. It is now recommended that the user use the `matching`
method desribed in [Model Search](model_search.md) to match models to
machine learning tasks (which can be understood there as simply
`Bool`-valued functions on models).

---

In MLJ a *task* is a synthesis of three elements: *data*, an
*interpretation* of that data, and a *learning objective*. Once one has a
task one is ready to choose learning models.

### Scientific types and the interpretation of data

Generally the columns of a table, such as a DataFrame, represents real
quantities. However, the nature of a quantity is not always clear from
the representation. For example, we might count phone calls using the
`UInt32` type but also use `UInt32` to represent a categorical
feature, such as the species of conifers. MLJ mitigates such ambiguity
with the use of scientific types. See [Getting Started](index.md)) for details

Explicitly specifying scientific types during the construction of a
MLJ task is the user's opportunity to articulate how the supplied data
should be interpreted.

### Learning objectives

In MLJ specifying a learning objective means specifying: (i) whether
learning is supervised or not; (ii) whether, in the supervised case,
predictions are to be probabilistic or deterministic; and (iii) what
part of the data is relevant and what role is each part to play.

### Sample usage

Load a built-in task:

```@setup 1
import Base.eval
using MLJ
MLJ.color_off()
```

```@example 1
using MLJ
using CSV
task = load_iris()
```

Extract input and target:

```@example 1
X, y = task()
X[1:3, :]
```

Now, starting with some tabular data...

```@example 1
using RDatasets
df = dataset("boot", "channing");
first(df, 4)
```

...we can check MLJ's interpretation of that data:

```@example 1
schema(df)
```

And construct a task by wrapping the data in a learning objective, and
coercing the data into a form MLJ will correctly interpret. (The middle three
fields of `df` refer to ages, in months, the last is a flag.):

```@example 1
task = supervised(data=df,
                  target=:Exit,
                  ignore=:Time,
                  is_probabilistic=true,
                  types=Dict(:Entry=>Continuous,
                             :Exit=>Continuous,
                             :Cens=>Multiclass))
schema(task.X)
```

Shuffle the rows of a task:

```@example 1
task = load_iris()
using Random
rng = MersenneTwister(1234)
shuffle!(rng, task) # rng is optional
task[1:4].y
```

Counting and selecting rows of a task:

```@example 1
nrows(task)
```

```@example 1
task[1:2].y
```

Listing the models available to complete a task:

```@example 1
models(task)
```

Binding a model to a task and evaluating performance:

```@repl 1
tree = @load DecisionTreeClassifier verbosity=1
mach = machine(tree, task)
evaluate!(mach, operation=predict_mode, resampling=Holdout(), measure=misclassification_rate, verbosity=0)
```

### API Reference

```@docs
supervised
```

```@docs
unsupervised
```


text:
# This script includes functions that are used to enumerate paths explictly, not using sampling.

type PathEnum
    path::Array
    length::Real
end

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

function path_enumeration(origin, destination, adj_mtx)
    link_length_dict = getLinkLengthDict(adj_mtx)
    path_enums = Array{PathEnum,1}(0)
    return path_enumeration(origin, destination, adj_mtx, link_length_dict)
end

function path_enumeration(origin, destination, start_node::Array, end_node::Array, link_length::Array)
    adj_mtx = getAdjacency(start_node, end_node)
    link_length_dict = getLinkLengthDict(start_node, end_node, link_length)
    return path_enumeration(origin, destination, adj_mtx, link_length_dict)
end

function path_enumeration(origin, destination, adj_mtx, link_length_dict::Dict)
    path_enums = Array{PathEnum,1}(0)
    current_list = [origin]
    find_neighbor!(origin, destination, adj_mtx, current_list, path_enums, link_length_dict)
    return path_enums
end

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

function find_neighbor!(current, destination, adj_mtx, current_list, path_enums, link_length_dict)

    for i=1:size(adj_mtx, 1)
        visited = length(find(x-> x==i, current_list)) > 0

        if ! visited && adj_mtx[current, i] == 1
            # Create a new object for the new list
            new_list = [current_list; i]

            if i == destination
                # Obtained a new path
                # println(fid, new_list)
                push!(path_enums, PathEnum(new_list, getPathLength(new_list, link_length_dict)) )
            else
                # Branch another neighbor search
                find_neighbor!(i, destination, adj_mtx, new_list, path_enums, link_length_dict)
            end
        end

    end

end

function actual_cumulative_count(path_enums::Array{PathEnum,1}, option=:unique)
    path_lengths = Array{Float64,1}(0)
    for enum in path_enums
        push!(path_lengths, enum.length)
    end

    N_data = 100
    # option == :uniform
    x_data = collect(linspace(minimum(path_lengths), maximum(path_lengths), N_data))
    if option == :unique
        x_data = sort(unique(path_lengths))
    elseif option == :first_quarter
        x_q1 = linspace(minimum(path_lengths), 0.25*maximum(path_lengths), N_data/2)
        x_q234 = linspace(0.25*maximum(path_lengths), maximum(path_lengths), N_data/2)
        x_data = append!(collect(x_q1), collect(x_q234))
    end

    y_data = similar(x_data)

    for i=1:length(x_data)
        y_data[i] = count(x-> x<=x_data[i], path_lengths)
    end

    return x_data, y_data
end


text:
using TOML

# Parse config files
global_cfg = TOML.parsefile(joinpath("main", "config.toml"))

# Track config
# Check for existence of track config file
track_cfg_file = joinpath("track-repo", ".github", "org-wide-files-config.toml")
if isfile(track_cfg_file)
    track_cfg = TOML.parsefile(track_cfg_file)

    # Configlet
    if haskey(track_cfg, "configlet")
        # Iterate through allowed configlet settings
        for key in global_cfg["configlet_configurable_keys"]
            haskey(track_cfg["configlet"], key) || continue

            if track_cfg["configlet"]["fmt"]
                @info "Track wants to run configlet fmt; appending to configlet workflow"

                open(joinpath("track-repo", ".github", "workflows", "configlet.yml"), "a") do io
                    write(io, "    with:\n      fmt: true\n")
                end
            end
        end
    end
else
    @info "$(ENV["TRACK"]) does not have an org-wide-files-config.toml"
end

# Appends
# Iterate through files listed in global config file
for file in global_cfg["appendable_files"]
    append_path = joinpath("track-repo", ".appends", file)

    # Check for existance of append file
    isfile(append_path) || (@info "$(ENV["TRACK"]) does not have an append for $file"; continue)

    # Read content of append file
    append_content = read(append_path, String)

    # Read org-wide file
    org_wide_content = readchomp(joinpath("track-repo", file))

    # Append file
    open(joinpath("track-repo", file), "w") do io
        # Write org-wide content
        write(io, org_wide_content)

        # Add or remove trailing newlines depending on global config
        # The default assumption is one trailing newline
        newlines = haskey(global_cfg["files"], file) ? get(global_cfg["files"][file], "trailing_newlines", 1) : 1
        write(io, "\n" ^ newlines)

        # Write append content
        write(io, append_content)

        @info "Appended $file" append_content
    end
end

max_stars_repo_name: ErikSchierboom/org-wide-files

max_issues_repo_name: ErikSchierboom/org-wide-files

max_forks_repo_name: ErikSchierboom/org-wide-files


text:
function MORFE_integrate_rdyn_frc(analysis_name,zero_amplitude,harmonics_init,param_init,cont_param,
                                  time_integration_length,forward=true,MaxNumPoints=100,minstep=1e-8,maxstep=20.0,ncol=4.0,ntst=40.0,analysis_number=1)
    # check for continuation direction
    if (forward)
      var = 0
    else
      var = 1
    end
    # output folder
    analysis_folder = "./"*analysis_name*"/matcont_automatic"
    analysis_output = "./"*analysis_name*"/frc_"*string(analysis_number)*".mat"
    #
    nΩ = Int(size(harmonics_init)[1]/2)
    control_parameters = "mu"
    for i = 1:nΩ
      control_parameters *= ",beta"*string(i)
    end
    for i = 1:nΩ
      control_parameters *= ",w"*string(i)
    end
    # initialise normal coordinates and harmonics amplitudes
    X0 = zeros(Float64,size(zero_amplitude)[1]+size(harmonics_init)[1])
    X0[1:size(zero_amplitude)[1]] = zero_amplitude
    X0[size(zero_amplitude)[1]+1:end] = harmonics_init
    # inizialise control parameters
    param = ""
    param *= "mu = "*string(param_init[1])
    for i = 1:nΩ
      param *= ";beta"*string(i) * " = " * string(param_init[i+1])
    end
    for i = 1:nΩ
      param *= ";w"*string(i) * " = " * string(param_init[i+1+nΩ])
    end
    param *= ";"
    #
    mat"""
    warning off
    addpath(genpath('./MatCont7p3'))
    addpath($analysis_folder)
    ndofs = size($X0);
    ndofs = ndofs(1);
    init
    X0   = $X0;
    tfin = $time_integration_length;
    """
    eval_string(param)
    mat"""
    ndofs = size(X0);
    ndofs = ndofs(1);
    hls=feval(@MORFEsystem);
    options=odeset('RelTol',1e-8);
    """

    param = "[t,y]=ode45(hls{2},[0,tfin],...s
                X0,options,"*control_parameters*");"
    eval_string(param)

    mat"""
    x1 = y(end,:);
    [vl,pk]=findpeaks(y(:,1));
    period=t(pk(end))-t(pk(end-1));
    """

    param = "[t,y] = ode45(hls{2},[0 period],x1,options,"*control_parameters*");"
    eval_string(param)

    mat"""
    active_pars=$cont_param;
    ncol=$ncol;
    ntst=$ntst;
    tolerance=1e-4;
    """
    param = "[x0,v0]=initOrbLC(@MORFEsystem,...
             t,y,...
             ["*control_parameters*"],[active_pars],ntst,ncol,...
             tolerance);"
    eval_string(param)

    mat"""
    opt=contset;
    opt=contset(opt,'MaxNumPoints',$MaxNumPoints);
    opt=contset(opt, 'InitStepsize' , 0.1);
    opt=contset(opt,'MaxStepsize'  , $maxstep);
    opt=contset(opt,'MinStepsize'  , $minstep);
    opt=contset(opt,'Backward',0);
    opt=contset(opt,'FunTolerance', 1e-6);
    opt=contset(opt,'VarTolerance', 1e-6);
    opt=contset(opt,'ActiveParams', active_pars);
    [xlcc,vlcc,slcc,hlcc,flcc]=cont(@limitcycle,x0,v0,opt);

    %FRF=[];
    %omega=[];
    %for i=1:size(xlcc,2)
    %   x1=xlcc(1:4:end-2,i);
    %   y1=xlcc(2:4:end-2,i);
    %   omega(i)=(xlcc(end-1,i).^-1)*2*pi;
    %   FRF(i)=max(x1);
    %   FRFv(i)=max(y1);
    %end

    %figure(1)
    %hold on
    %plot(omega,FRF)

    save($analysis_output,'xlcc','ncol','ntst')
    """

    #
    return nothing
    #
  end

function MORFE_integrate_rdyn_backbone(analysis_name,zero_amplitude,time_integration_length,forward=true,MaxNumPoints=100,minstep=1e-8,maxstep=20.0,ncol=4.0,ntst=40.0,analysis_number=1)

  if (forward)
    var = 0
  else
    var = 1
  end

  analysis_folder = "./"*analysis_name*"/matcont_automatic"
  analysis_output = "./"*analysis_name*"/backbone_"*string(analysis_number)*".mat"

  mat"""
  warning off
  addpath(genpath('./MatCont7p3'))
  addpath($analysis_folder)

  init

  mu = 0.0;
  ndofs = size($zero_amplitude);
  ndofs = ndofs(1);

  X0 = $zero_amplitude;

  tfin=$time_integration_length;

  hls=feval(@MORFEsystem);
  options=odeset('RelTol',1e-8);

  [t,y]=ode45(hls{2},[0,tfin],...s
  X0,options,mu);

  x1 = y(end,:);

  [vl,pk]=findpeaks(y(:,1));
  period=t(pk(end))-t(pk(end-1));
  [t,y] = ode45(hls{2},[0 period],x1,options,mu);

  active_pars=[1];
  ncol=$ncol;
  ntst=$ntst;
  tolerance=1e-4;
  [x0,v0]=initOrbLC(@MORFEsystem,...
  t,y,...
  [mu],active_pars,ntst,ncol,...
  tolerance);

  opt=contset;
  opt=contset(opt,'MaxNumPoints',$MaxNumPoints);
  opt=contset(opt, 'InitStepsize' , 0.1);
  opt=contset(opt,'MaxStepsize'  , $maxstep);
  opt=contset(opt,'MinStepsize'  , $minstep);
  opt=contset(opt,'Backward',$var);
  opt=contset(opt,'FunTolerance', 1e-6);
  opt=contset(opt,'VarTolerance', 1e-6);
  [xlcc,vlcc,slcc,hlcc,flcc]=cont(@limitcycle,x0,v0,opt);

  %FRF=[];
  %omega=[];
  %for i=1:size(xlcc,2)
  %    x1=xlcc(1:2:end-2,i);
  %    y1=xlcc(2:2:end-2,i);
  %    omega(i)=(xlcc(end-1,i).^-1)*2*pi;
  %    FRF(i)=max(x1);
  %    FRFv(i)=max(y1);
  %end

  %figure(1)
  %hold on
  %plot(omega,FRF)

  save($analysis_output,'xlcc','ncol','ntst')

  """

  #
  return nothing
  #

end

function MORFE_compute_backbone_modal(analysis_name,analysis_number=1)
  #
  comb = readdlm(analysis_name*"/comb_a_0/monomials.txt")
  maps = readdlm(analysis_name*"/manifold_a_0/psi.txt")
  #
  po_dir = "./"*analysis_name*"/backbone_"*string(analysis_number)*".mat"
  #
  nc = size(comb)[1]
  ndofs = size(comb)[2]
  neig = size(maps)[2]
  #
  data = matread(po_dir)
  #
  xlcc = data["xlcc"]
  ncol = data["ncol"]
  ntst = data["ntst"]
  #
  npoints = Int(ncol*ntst)
  npo = size(xlcc)[2]
  vars = zeros(Float64,ndofs)
  field = zeros(Float64,neig)
  frf = zeros(npo,neig+1)
  #
  for po = 1:npo
  #
    ω = 2*π/xlcc[end-1,po]
    frf[po,1] = ω
    for p = 1:npoints
      # extract variables
      fill!(field,0.0)
      #
      for d = 1:ndofs
        vars[d] = xlcc[d+(p-1)*ndofs,po]
      end
      #
      for c = 1:nc
        #
        Πi = 1.0
        for d = 1:ndofs
          Πi *= vars[d]^comb[c,d]
        end
        #
        for d = 1:neig
          field[d] += maps[c,d]*Πi
        end
        #
      end
      #
      for d = 1:neig
        if (abs(field[d])>frf[po,1+d])
          frf[po,1+d] = abs(field[d])
        end
      end
      #
    end
    #
  end
  #
  return frf
  #
end

# reconstruction with only autonomous vector till the paper on the forcing is
# becomes public
function MORFE_compute_frc_modal(analysis_name,Ω_list,analysis_number=1)
  #
  comb = readdlm(analysis_name*"/comb_a_0/monomials.txt")
  maps = readdlm(analysis_name*"/manifold_a_0/psi.txt")
  #
  po_dir = "./"*analysis_name*"/frc_"*string(analysis_number)*".mat"
  #
  nc = size(comb)[1]
  ndofs = size(comb)[2]
  nΩ = 2*size(Ω_list)[1]
  neig = size(maps)[2]
  #
  data = matread(po_dir)
  #
  xlcc = data["xlcc"]
  ncol = data["ncol"]
  ntst = data["ntst"]
  #
  npoints = Int(ncol*ntst)
  npo = size(xlcc)[2]
  vars = zeros(Float64,ndofs)
  field = zeros(Float64,neig)
  frf = zeros(npo,neig+1)
  #
  for po = 1:npo
  #
    ω = xlcc[end,po]
    frf[po,1] = ω
    for p = 1:npoints
      # extract variables
      fill!(field,0.0)
      #
      for d = 1:ndofs
        vars[d] = xlcc[d+(p-1)*(ndofs+nΩ),po]
      end
      #
      for c = 1:nc
        #
        Πi = 1.0
        for d = 1:ndofs
          Πi *= vars[d]^comb[c,d]
        end
        #
        for d = 1:neig
          field[d] += maps[c,d]*Πi
        end
        #
      end
      #
      for d = 1:neig
        if (abs(field[d])>frf[po,1+d])
          frf[po,1+d] = abs(field[d])
        end
      end
      #
    end
    #
  end
  #
  return frf
  #
end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

"""
    AbstractCallback

Abstract type of callback functions used in training.
"""
abstract type AbstractCallback end

"""
    AbstractBatchCallback

Abstract type of callbacks to be called every mini-batch.
"""
abstract type AbstractBatchCallback <: AbstractCallback end

"""
    AbstractEpochCallback

Abstract type of callbacks to be called every epoch.
"""
abstract type AbstractEpochCallback <: AbstractCallback end

mutable struct BatchCallback <: AbstractBatchCallback
  frequency :: Int
  call_on_0 :: Bool
  callback  :: Function
end

"""
    every_n_batch(callback :: Function, n :: Int; call_on_0 = false)

A convenient function to construct a callback that runs every `n` mini-batches.

# Arguments
* `call_on_0::Bool`: keyword argument, default false. Unless set, the callback
          will *not* be run on batch 0.

For example, the [`speedometer`](@ref) callback is defined as

```julia
every_n_batch(frequency, call_on_0=true) do state :: OptimizationState
  if state.curr_batch == 0
    # reset timer
  else
    # compute and print speed
  end
end
```

See also [`every_n_epoch`](@ref) and [`speedometer`](@ref).
"""
function every_n_batch(callback::Function, n::Int; call_on_0::Bool = false)
  BatchCallback(n, call_on_0, callback)
end
function (cb :: BatchCallback)(state :: OptimizationState)
  if state.curr_batch == 0
    if cb.call_on_0
      cb.callback(state)
    end
  elseif state.curr_batch % cb.frequency == 0
    cb.callback(state)
  end
end

"""
    speedometer(;frequency=50)

Create an `AbstractBatchCallback` that measure the training speed
   (number of samples processed per second) every k mini-batches.

# Arguments
* `frequency::Int`: keyword argument, default 50. The frequency (number of
          min-batches) to measure and report the speed.
"""
function speedometer(;frequency::Int = 50)
  cl_tic = 0
  every_n_batch(frequency, call_on_0 = true) do state::OptimizationState
    if state.curr_batch == 0
      # reset timer
      cl_tic = time()
    else
      speed = frequency * state.batch_size / (time() - cl_tic)
      @info(format("Speed: {1:>6.2f} samples/sec", speed))
      cl_tic = time()
    end
  end
end

mutable struct EpochCallback <: AbstractEpochCallback
  frequency :: Int
  call_on_0 :: Bool
  callback  :: Function
end

"""
    every_n_epoch(callback :: Function, n :: Int; call_on_0 = false)

A convenient function to construct a callback that runs every `n` full data-passes.

* `call_on_0::Bool`: keyword argument, default false. Unless set, the callback
          will *not* be run on epoch 0. Epoch 0 means no training has been performed
          yet. This is useful if you want to inspect the randomly initialized model
          that has not seen any data yet.

See also [`every_n_batch`](@ref).
"""
every_n_epoch(callback::Function, n::Int; call_on_0::Bool = false) =
  EpochCallback(n, call_on_0, callback)

function (cb::EpochCallback)(model::Any, state::OptimizationState,
                             metric::Vector{Tuple{Symbol, T}}) where T<:Real
  if state.curr_epoch == 0
    if cb.call_on_0
      cb.callback(model, state, metric)
    end
  elseif state.curr_epoch % cb.frequency == 0
    cb.callback(model, state, metric)
  end
end

"""
    do_checkpoint(prefix; frequency=1, save_epoch_0=false)

Create an `AbstractEpochCallback` that save checkpoints of the model to disk.
The checkpoints can be loaded back later on.

# Arguments
* `prefix::AbstractString`: the prefix of the filenames to save the model.
  The model architecture will be saved to prefix-symbol.json,
  while the weights will be saved to prefix-0012.params,
  for example, for the 12-th epoch.
* `frequency::Int`: keyword argument, default is 1.
  The frequency (measured in epochs) to save checkpoints.
* `save_epoch_0::Bool`: keyword argument, default false. Whether we should save a
  checkpoint for epoch 0 (model initialized but not seen any data yet).
"""
function do_checkpoint(prefix::AbstractString;
                       frequency::Int = 1, save_epoch_0::Bool = false)
  mkpath(dirname(prefix))
  every_n_epoch(frequency, call_on_0=save_epoch_0) do model, state, metric
    save_checkpoint(model, prefix, state)
  end
end


text:
# Before running this, please make sure to activate and instantiate the environment
# corresponding to [this `Project.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Project.toml) and [this `Manifest.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Manifest.toml)
# so that you get an environment which matches the one used to generate the tutorials:
#
# ```julia
# cd("DataScienceTutorials") # cd to folder with the *.toml
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

using MLJ
using PrettyPrinting

@load KNNRegressor
# input
X = (age    = [23, 45, 34, 25, 67],
     gender = categorical(['m', 'm', 'f', 'm', 'f']))
# target
height = [178, 194, 165, 173, 168];

scitype(X.age)

pipe = @pipeline(X -> coerce(X, :age=>Continuous),
                OneHotEncoder(),
                KNNRegressor(K=3),
                target = UnivariateStandardizer());

pipe.knn_regressor.K = 2
pipe.one_hot_encoder.drop_last = true;

evaluate(pipe, X, height, resampling=Holdout(),
         measure=rms) |> pprint

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
# 这是一个风洞配平程序

function uisol(uitmp,θ7,θlat,θlon)  # 均匀入流求解力和力矩的流程
  # 用于雅可比矩阵中求平均周期力
  ψxx = 0.0
  iter = 1
  sumf = 0.0
  blat = 0.0
  blon = 0.0
  θ0 = Array{Float64}(NR,Nbe)
  θcp = zeros(Float64,NR)
  for k in 1:NR
    θcp[k] = θ7
    for i in 1:Nbe
      θ0[k,i] = θcp[k]+θtw*(rb[k,i]/R-0.7)
    end
  end

  while iter<=npsi
    vind_r = uitmp[1]
    vall_r = uitmp[2]
    beta = uitmp[4]
    dbeta = uitmp[5]
    blat = blat+uitmp[4]*sin(ψxx)
    blon = blon+uitmp[4]*cos(ψxx)

    # This file is used for get angle-of-attack
    aoatmp = aoaget(vall_r,ψxx,θ0,θlat,θlon,beta)
    α_aero = aoatmp[1]
    θ = aoatmp[2]

    # This file is to get the Cl and Cd of Blade Element
    clcdtmp = clcdget(α_aero,vall_r)
    Cl = clcdtmp[1]
    Cd = clcdtmp[2]

    # This file calculate the whole aerodynamic force of the rotor in the hub
    rftmp = rotorforce(ψxx,vall_r,α_aero,θ,Cl,Cd)
    sumf = sumf+rftmp[3]
    Mbeta_aero = rftmp[5][1]

    # The file calculate the uniform induced velocity
    uitmp = uniforminflow(ψxx,θcp,θlat,θlon,beta,dbeta,Mbeta_aero)

    ψxx = ψxx+dψ
    iter = iter+1
  end
  sumf = sumf/npsi
  blat = blat/npsi*2
  blon = blon/npsi*2

  return sumf,blat,blon
end

function yagb(uitmp,θcp,θlat,θlon,epsi=1/180*pi)
  # 前向差分雅可比矩阵

  # 后向差分的雅可比矩阵

  # 中心差分雅可比矩阵
  thecp = θcp+epsi/2
  uisoltmp = uisol(uitmp,thecp,θlat,θlon)
  rotfor = uisoltmp[1]
  blatfor = uisoltmp[2]
  blonfor = uisoltmp[3]

  thecp = θcp-epsi/2
  uisoltmp = uisol(uitmp,thecp,θlat,θlon)
  rotbac = uisoltmp[1]
  blatbac = uisoltmp[2]
  blonbac = uisoltmp[3]

  dthecp = (rotfor-rotbac)/epsi
  dtheblat = (blatfor-blatbac)/epsi
  dtheblon = (blonfor-blonbac)/epsi

  lat = θlat+epsi/2
  uisoltmp = uisol(uitmp,θcp,lat,θlon)
  rotfor = uisoltmp[1]
  blatfor = uisoltmp[2]
  blonfor = uisoltmp[3]

  lat = θlat-epsi/2
  uisoltmp = uisol(uitmp,θcp,lat,θlon)
  rotbac = uisoltmp[1]
  blatbac = uisoltmp[2]
  blonbac = uisoltmp[3]

  dlat = (rotfor-rotbac)/epsi
  dlatblat = (blatfor-blatbac)/epsi
  dlatblon = (blonfor-blonbac)/epsi

  lon = θlon+epsi/2
  uisoltmp = uisol(uitmp,θcp,θlat,lon)
  rotfor = uisoltmp[1]
  blatfor = uisoltmp[2]
  blonfor = uisoltmp[3]

  lon = θlon-epsi/2
  uisoltmp = uisol(uitmp,θcp,θlat,lon)
  rotbac = uisoltmp[1]
  blatbac = uisoltmp[2]
  blonbac = uisoltmp[3]

  dlon = (rotfor-rotbac)/epsi
  dlonblat = (blatfor-blatbac)/epsi
  dlonblon = (blonfor-blonbac)/epsi

  Myg = [dthecp  dlat  dlon;
          dtheblat  dlatblat  dlonblat;
          dtheblon  dlatblon  dlonblon]

  return Myg
end

function trimwt(uitmp,rot,beta_lat,beta_lon,θcp,θlat,θlon,epsi=1/180*π) # 纵横向挥舞配平法
  fzt = abs(rot*cos(αs)-T) # 力的配平条件
  epsi = 1/180*π # 雅可比求解的变化小量
  # 求解雅可比矩阵
  # 输入量：θ0,θ_lat,θ_lon
  # 输出量：rot,beta_lat,beta_lon

  if fzt<10&&abs(beta_lat)<1e-3&&abs(beta_lon)<1e-3 # 配平力和挥舞
  # if abs(beta_lat)<1e-3&abs(beta_lon)<1e-3 # 仅配平挥舞
    return true,θ0,θlat,θlon,θcp
  else
    Myg = yagb(uitmp,θcp[1],θlat,θlat,epsi) # 时变雅可比矩阵
    dc = inv(Myg)*[T-rot,-beta_lat,-beta_lon]
    for k = 1:NR
      θcp[k] = θcp[k]+dc[1]
      θlat[k] = θlat[k]+dc[2]
      θlon[k] = θlon[k]+dc[3]
      for i in 1:Nbe
        θ0[k,i] = θcp[k]+θtw*(rb[k,i]/R-0.7)
      end
    end
  end

  return false,θ0,θlat,θlon,θcp
end

max_stars_repo_name: dingzhw/RotorTrim

max_issues_repo_name: helicopterliu/RotorTrim

max_forks_repo_name: helicopterliu/RotorTrim


text:
# List of vector to vector conversion functions
# should accept a function call like: f(x) where x is a input vector
const vec2vec = [
                 :mcep,
                 :gcep,
                 :mgcep,
                 :uels,
                 :fftcep,
                 :lpc,

                 :mfcc,

                 :lpc2c,
                 :lpc2lsp,
                 :lpc2par,
                 :par2lpc,
                 :lsp2sp,

                 :mc2b,
                 :b2mc,
                 :b2c,
                 :c2acr,
                 :c2ir,
                 :ic2ir,
                 :c2ndps,
                 :ndps2c,
                 :gc2gc,
                 :gnorm,
                 :ignorm,
                 :freqt,
                 :frqtr,
                 :mgc2mgc,
                 :mgc2sp,
                 :mgclsp2sp
                 ]

# extend vector to vector conversion for matrix input (col-wise)
# should accept a function call like: f(x) where x is a input matrix
for f in vec2vec
    @eval begin
        function $f(x::StridedMatrix{Cdouble}, args...; kargs...)
            outbuf = $f(view(x, :, 1), args...; kargs...)
            ret = Array{eltype(outbuf)}(length(outbuf), size(x, 2))
            copy!(ret, 1, outbuf, 1, length(outbuf))
            for i = 2:size(x, 2)
                @inbounds ret[:, i] = $f(view(x, :, i), args...; kargs...)
            end
            ret
        end
    end
end


text:
using Widgets, Observables, OrderedCollections
using Widgets: div, @nodeps, Widget, dropdown, button, textbox
using Observables: @map, AbstractObservable
import FileIO

_empty(s::AbstractString) = s == ""
_empty(s::Symbol) = s == Symbol("")
_empty(p::Pair) = _empty(last(p))
exclude_empty(v) = Iterators.filter(!_empty, v)

function gui(df)
    df isa AbstractObservable || (df = Observable{Any}(df))
    t = map(columntable, df)
    names = @map collect(propertynames(&t))
    maybe_names = @map vcat(Symbol(""), &names)
    x = dropdown(names, label = "First axis")
    y = dropdown(maybe_names, label = "Second axis")
    z = dropdown(maybe_names, label = "Third axis")
    plot_func = dropdown([plot, scatter, lines, barplot, heatmap, surface, wireframe,
        volume, contour, boxplot, violin], label = "Plot function")
    analysis_options = OrderedDict(
        "" => tuple,
        "density" => density,
        "histogram" => histogram,
        "linear" => linear,
        "smooth" => smooth
    )
    analysis = dropdown(analysis_options, label = "Analysis")
    group_attr = [:color, :marker, :linestyle]
    style_attr = [:color, :markersize]
    groups = [dropdown(maybe_names, label = string(l)) for l in group_attr]
    styles = [dropdown(maybe_names, label = string(l)) for l in style_attr]

    output = Observable{Any}(text("Welcome to the StatsMakie GUI", align = (:center, :center)))

    plot_button = button("Plot")
    save_button = button("Save")
    save_name = textbox(placeholder = "Save as...")

    ui = Widget{:gui}(
        OrderedDict(
            "plot_button" => plot_button,
            "save_button" => save_button,
            "save_name" => save_name,
            "table" => t,
            "plot_func" => plot_func,
            "analysis" => analysis,
            "x" => x,
            "y" => y,
            "z" => z,
            "group" => groups,
            "style" => styles
        ),
        output = output
    )
    map!(output, plot_button) do _
        vars = exclude_empty([x[], y[], z[]])
        grps = exclude_empty([(key => val[]) for (key, val) in zip(group_attr, groups)])
        stls = exclude_empty((key => val[]) for (key, val) in zip(style_attr, styles))
        g = Group(; grps...)
        s = Style(; stls...)
        plot_func[](analysis[], Data(t[]), g, s, vars...)
    end
    on(save_button) do _
        save_plot(save_name[], output[])
    end
    row_style = Dict("display" => "flex", "direction" => "row")
    column_style = Dict("display" => "flex", "direction" => "column")
    title = t -> Widgets.node("p", t, style = Dict("font-weight" => "bold"))
    hspace = div(style = Dict("width" => "1em"))
    vspace = div(style = Dict("height" => "1em"))
    @layout! ui begin
        top_row = div(
            :plot_func,
            hspace,
            :analysis,
            hspace,
            :x,
            hspace,
            :y,
            hspace,
            :z,
            style = row_style
        )
        button_row = div(
            :plot_button,
            hspace,
            :save_button,
            hspace,
            :save_name,
            style = row_style
        )
        group_column = div(title("Group"), vspace, :group...)
        style_column = div(title("Style"), vspace, :style...)
        bottom_row = div(_.output, hspace, group_column, hspace, style_column, style = row_style)
        div(top_row, vspace, button_row, vspace, bottom_row)
    end
end

function save_plot(name, scene)
    fp = joinpath(homedir(), ".StatsMakie", "plots")
    ispath(fp) || mkpath(fp)
    FileIO.save(joinpath(fp, name), scene)
end


text:
using BinaryBuilder

# Collection of sources required to build PALISADE
name = "PALISADE"
version = v"1.4.1"
sources = [
    "https://git.njit.edu/palisade/PALISADE/-/archive/PALISADE-v1.4.1/PALISADE-PALISADE-v1.4.1.tar.gz" =>
    "cdeddfc6d6f059047d626ff0e221b3455058bd5bcabdef8602ba53b5586e3377",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/PALISADE*

# We need coreutils
apk add coreutils

# Don't build GMP, we've already got it
mkdir -p third-party/lib
touch third-party/lib/lib{gmp,ntl}.a

export CFLAGS="${CFLAGS} -I${prefix}/include"
export LDFLAGS="${LDFLAGS} -L${prefix}/lib"

FLAGS=(PREFIX="${prefix}" GMP_UNPACK_NEEDED= NTL_UNPACK_NEEDED= THIRDPARTYINCLUDE="-I${prefix}/include -I$(pwd)/third-party/include -I$(pwd)/third-party/include/rapidjson" NTLLIB=${prefix}/lib/libntl.so GMPLIB=${prefix}/lib/libgmp.so CXX="ccache ${CXX}")

make -j${nproc} "${FLAGS[@]}"
make installcore installpke installabe installsignature "${FLAGS[@]}"
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64; compiler_abi=CompilerABI(:gcc6)),
    Linux(:x86_64; libc=:musl, compiler_abi=CompilerABI(:gcc6)),
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libsnark", :libsnark),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/GMP-v6.1.2-1/build_GMP.v6.1.2.jl",
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/ntl-v10.5.0%2B0/build_ntl.v10.5.0.jl",
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
# N-body integrator with stupid name

module ImprovedGarbanzo
include("common.jl")

# using .common
using Plots
using LinearAlgebra
using ProgressMeter
# pyplot()

mode = :runge_kutta
ε = 0.0001

function rk4!(p::Body, a, h)
    k1_v = a
    k1_r = p.v

    k2_v = a .* (k1_r*h/2)
    k2_r = p.v .* k1_v*h/2

    k3_v = a .* (k2_r*h/2)
    k3_r = p.v .* k2_v*h/2

    k4_v = a .* (k3_r*h)
    k4_r = p.v .* k3_v*h

    p.r += h/6 * (k1_r + 2*k2_r + 2*k3_r + k4_r)
    p.v += h/6 * (k1_v + 2*k2_v + 2*k3_v + k4_v)
end

function euler!(p::Body, a, h)
    p.r += h * p.v
    p.v += h * a
end

function leapfrog!(p::Body, a, h)
    p.r += 0.5h * p.v
    p.v += h * a
    p.r += 0.5h * p.v
end

function step!(ensemble, h)
    global mode
    for p1 in ensemble
        a = [0.,0.,0.]
        for p2 in ensemble
            if p1 == p2
                continue
            end
            dr = p1.r - p2.r
            r2 = dot(dr, dr)
            r3 = r2 * √r2 + ε
            a -= (dr / r3)
        end
        if mode == :fwd_euler
            euler!(p1, a, h)
        else
            rk4!(p1, a, h)
        end
    end
end

"""
Performs full simulation during `n` steps.
it calls `output` function every  `output_n` steps,
where `output_n` = `n` / `1000`
"""
function perform!(n, ensemble::Vector{Body}, output, h=0.0001)
    hidden_n = 100
    output_n = Integer(n / hidden_n)

    println(stderr, "Performing simulation")
    prog = Progress(output_n, 1, "Computing...")
    @gif for i=1:output_n
        for j=1:hidden_n
            step!(ensemble, h)
        end
        output(ensemble)
        next!(prog)
    end
end

function prep_plot(ensemble)
    lim = (-2, 2)
    marker = (:cross, 3, stroke(0))
    global plt = scatter3d(1, m = marker,
        ylim=lim, xlim=lim, zlim=lim)

    for i=2:length(ensemble)
        scatter3d!(1, m = marker)
    end
end

function print_out(ensemble::Vector{Body})
    println(stderr, ensemble[1])
end

"""
Adds subsequent postitions of particles in the system
to draw trajectory trails.
"""
function plot_trajectory_out(ensemble::Vector{Body})
    for (i, p) in enumerate(ensemble)
        push!(plt.series_list[i], p.r[1], p.r[2], p.r[3])
    end
end

"""
Every time it creates new plot to present
only current position of particles in the system.
"""
function plot_distribution_out(ensemble::Vector{Body})
    lim = (-2, 2)
    marker = (:cross, 3, stroke(0))
    global plt = scatter3d(1, m = marker,
            ylim=lim, xlim=lim, zlim=lim)
    for (i, p) in enumerate(ensemble)
        push!(plt, p.r[1], p.r[2], p.r[3])
    end
end

"""
_Example_
Solves system for `N` random bodies
"""
function ex_random(N, new_mode)
    global mode = new_mode
    e = [Body(
            2rand(Float32, 3) .- 1.,
            (2rand(Float32, 3) .- 1.)*4,
            1.)
        for i=1:N]

    perform!(100000, e, plot_distribution_out, 0.0001)
end

# perform!(150000, e, plot_trajectory_out)
# ex_random(100)
# plot(plt, size = (500, 500))

end

doc = """ImprovedGarbanzo
Simple N-body simulator

Usage:
  ImprovedGarbanzo.jl perform <N>
  ImprovedGarbanzo.jl perform <N>  [--euler|--runge-kutta]
  ImprovedGarbanzo.jl -h | --help
  ImprovedGarbanzo.jl --version

Options:
  -h --help     Show this screen.
  --version     Show version.
  --euler
  --runge-kutta

"""
using DocOpt  # import docopt function

function cli()

  arguments = docopt(doc, version=v"1.0.0")

  # aliases = Dict(
  #   "p" => "perform",
  # )
  #
  # for (key, value) in arguments
  #   if value != true ; continue ; end
  #
  #   if haskey(aliases, key)
  #     key = aliases[key]
  #   end
  # end

  if haskey(arguments, "perform")
      println(arguments)
      N = parse(Int ,arguments["<N>"])
      if arguments["--euler"]
          mode = :fwd_euler
      elseif arguments["--runge-kutta"]
          mode = :runge_kutta
      end
      ImprovedGarbanzo.ex_random(N, mode)
      return
  end

  error("Invalid cli input.")

end

# cli()


text:
using .Gadfly

function plotMACs(z::Integer)
    colors, name = [ "red", "green", "blue", "orange" ], [ "[μ/ρ]ₚₑ", "[μ/ρ]ci", "[μ/ρ]ₜₒₜ", "[μ/ρ]ₖ" ]
    lpe = layer(loge->log10(mac(PhotoElectricMAC, z, 10.0^loge)), log10(50.0), log10(4.0e5), Theme(default_color = colors[1]))
    lci = layer(loge->log10(mac(CoherentIncoherentMAC, z, 10.0^loge)), log10(50.0), log10(4.0e5), Theme(default_color = colors[2]))
    ltot = layer(loge->log10(mac(TotalMAC, z, 10.0^loge)), log10(50.0), log10(4.0e5), Theme(default_color = colors[3]))
    lK = layer(loge->log10(mac(KMAC,z, 10.0^loge)), log10(50.0), log10(4.0e5), Theme(default_color = colors[4]))
    plot(lpe, lci, ltot, lK,
        Guide.title("Z = $z"),
        Guide.manual_color_key("Type", name, colors),
        Guide.xlabel("log₁₀(E) (eV)"), Guide.ylabel("log₁₀([μ/ρ]) (g/cm²)"),
        Coord.cartesian(xmin = log10(50.0), xmax = log10(4.0e5), ymin = -2.0))
end

function plotFormFactors(z::Integer)
    colors, name = [ "red", "blue" ], [ "f₁(E)", "f₂(E)" ]
    lf1 = layer(loge->formfactors(z, 10.0^loge)[1], log10(50.0), log10(4.0e5), Theme(default_color = colors[1]))
    lf2 = layer(loge->formfactors(z, 10.0^loge)[2], log10(50.0), log10(4.0e5), Theme(default_color = colors[2]))
    plot(lf1, lf2,
        Guide.title("Z = $z"),
        Guide.manual_color_key("Type", name, colors),
        Guide.xlabel("log₁₀(E) (eV)"), Guide.ylabel("log₁₀(f(E)) (e/atom)"),
        Coord.cartesian(xmin = log10(50.0), xmax = log10(4.0e5), ymin = 0.0))
end

function plotEdgeCount()
    colors, name = [ "red", "blue" ], [ "Count", "Last" ]
    l1 = layer(x = eachelement(), y = length.(eachedge.(eachelement())), Theme(default_color = colors[1]))
    l2 = layer(x = eachelement(), y = maximum.(eachedge.(eachelement())), Theme(default_color = colors[2]))
    plot(l1, l2, Guide.title("Edges"), Guide.manual_color_key("Quantity", name, colors), Guide.xlabel("Z"), Guide.ylabel("Quantity"))
end

plotAtomicWeight() =
    plot(x = eachelement(), y = atomicweight.(eachelement()),
    Guide.title("Nominal Atomic Weight"),
    Guide.xlabel("Z"), Guide.ylabel("g/mole"))

plotCrossSectionFactor() =
    plot(x = eachelement(), y = 1.0e24 * crosssectionfactor.(eachelement()),
        Guide.xlabel("Z"), Guide.ylabel("barns/atom"),
        Guide.title("Cross Section Factor"),
        Coord.cartesian(xmin = 0, xmax = 92, ymin = 0.0, ymax = 400.0))

plotDensity() =
    plot(x = eachelement(), y = density.(eachelement()), #
    Guide.title("Nominal Density"), #
    Guide.xlabel("Z"), Guide.ylabel("g/cm³"))

function plotRelativisticCorrection()
    l1 = layer(x = eachelement(), y = map(z->relativisticcorrection(z)[1], eachelement()), Theme(default_color = "red"))
    l2 = layer(x = eachelement(), y = map(z->relativisticcorrection(z)[2], eachelement()), Theme(default_color = "blue"))
    plot(l1, l2, Guide.xlabel("Z"), Guide.ylabel("e/atom"), #
        Guide.manual_color_key("Type", ["1", "2" ], [ "red", "blue"]), #
        Guide.title("fᵣₑₗ(H82,3/5CL)"), #
        Coord.cartesian(xmin = 0, xmax = 92, ymin = -2.6, ymax = 0.0))
end

plotNuclearThompsonCorrection() =
    plot(x = eachelement(), y = nuclearthompsoncorrection.(eachelement()),
        Guide.title("Nuclear Thompson Correction"),
        Guide.xlabel("Z"), Guide.ylabel("e/atom"))

plotJumpRatios(shell::Int) =
    plot(x = eachelement(), y = jumpratio.(eachelement(), shell),
        Guide.title("Jump Ratio"),
        Guide.xlabel("Z"), Guide.ylabel("Ratio"))


text:
using FlatBuffers
using Base.Test

if !isdefined(Core, :String)
    typealias String UTF8String
end

include("internals.jl")
CheckByteLayout()
CheckManualBuild()
CheckVtableDeduplication()
CheckNotInObjectError()
CheckStringIsNestedError()
CheckByteStringIsNestedError()
CheckStructIsNotInlineError()
CheckFinishedBytesError()
CheckCreateByteVector()
checkFuzz(100, 100, true)

include("monster.jl")

vec3 = Example.Vec3(1.0, 2.0, 3.0, 3.0, Example.Color(1), Example.Test(5, 6))
test4 = Example.Test[Example.Test(10, 20), Example.Test(30, 40)]
testArrayOfString = ["test1","test2"]

mon = Example.Monster(vec3, 150, 80, "MyMonster", false, collect(0x00:0x04),
        Example.Blue, test4, testArrayOfString, Example.Monster[],
        UInt8[], Example.Stat("",0,0), false, 0, 0, 0, 0, 0, 0, 0, 0,
        Bool[], 0, 0, 0)
b = FlatBuffers.build!(mon)
monst = FlatBuffers.read(b)

@test mon.pos == monst.pos

# create test types
# types (Scalar, Enum, immutable, T, String, Vector{UInt8})
type TestInt8T
    x::Int8
end

inst1 = TestInt8T(1)

b = FlatBuffers.Builder(TestInt8T)
FlatBuffers.build!(b, inst1)
t = FlatBuffers.Table(b)
inst1_2 = FlatBuffers.read(t)

@test inst1.x === inst1_2.x

immutable TestInt8I
    x::Int8
end

inst2 = TestInt8I(2)

type TestInt8A
    x::Vector{Int8}
end

inst3 = TestInt8A([1,2,3])

b = FlatBuffers.Builder(TestInt8A)
FlatBuffers.build!(b, inst3)
t = FlatBuffers.Table(b)
inst3_2 = FlatBuffers.read(t)

@test inst3.x == inst3_2.x

type TestMixT
    x::Int8
    y::String
    z::Vector{Int8}
end

inst4 = TestMixT(10,"hey there sailor",[1,2,3])

b = FlatBuffers.Builder(TestMixT)
FlatBuffers.build!(b, inst4)
t = FlatBuffers.Table(b)
inst4_2 = FlatBuffers.read(t)

@test inst4.x == inst4_2.x && inst4.y == inst4_2.y && inst4.z == inst4_2.z

# simple sub-table/type (Stat)
type TestSubT
    x::TestInt8T
    y::TestInt8I
    z::TestInt8A
end

inst5 = TestSubT(inst1, inst2, inst3)

b = FlatBuffers.Builder(TestSubT)
FlatBuffers.build!(b, inst5)
t = FlatBuffers.Table(b)
inst5_2 = FlatBuffers.read(t)

@test inst5.x.x == inst5_2.x.x && inst5.y.x == inst5_2.y.x && inst5.z.x == inst5_2.z.x

# vtable duplicates
type TestDupT
    x::TestInt8T
    y::TestInt8I
    z::TestInt8T
end

inst6 = TestDupT(inst1, inst2, TestInt8T(2))

b = FlatBuffers.Builder(TestDupT)
FlatBuffers.build!(b, inst6)
t = FlatBuffers.Table(b)
inst6_2 = FlatBuffers.read(t)

@test inst6.x.x == inst6_2.x.x && inst6.y.x == inst6_2.y.x && inst6.z.x == inst6_2.z.x

type TestDup2T
    x::Vector{TestInt8T}
end

inst7 = TestDup2T([inst1, TestInt8T(2), TestInt8T(3), TestInt8T(4)])

b = FlatBuffers.Builder(TestDup2T)
FlatBuffers.build!(b, inst7)
t = FlatBuffers.Table(b)
inst7_2 = FlatBuffers.read(t)

@test all(map(x->x.x, inst7.x) .== map(x->x.x, inst7_2.x))

# self-referential type test (type has subtype of itself)
# type TestCircT
#     x::Int8
#     y::TestCircT
# end
#
# immutable TestCircI
#     x::Int8
#     y::TestCircI
# end

# simple Union (Any_)

# fbs
# table TestUnionT {
#     x::TestUnionI
# }

@union TestUnionU Union{Void,TestInt8T,TestInt8A}

type TestUnionT
    x_type::Int8
    x::TestUnionU
end

TestUnionT(x::TestUnionU) = TestUnionT(FlatBuffers.typeorder(TestUnionU, typeof(x)), x)

inst8 = TestUnionT(inst1)

b = FlatBuffers.Builder(TestUnionT)
FlatBuffers.build!(b, inst8)
t = FlatBuffers.Table(b)
inst8_2 = FlatBuffers.read(t)

@test inst8.x_type == inst8_2.x_type && inst8.x.x == inst8_2.x.x

inst9 = TestUnionT(inst3)

b = FlatBuffers.Builder(TestUnionT)
FlatBuffers.build!(b, inst9)
t = FlatBuffers.Table(b)
inst9_2 = FlatBuffers.read(t)

@test inst9.x_type == inst9_2.x_type && inst9.x.x == inst9_2.x.x

# test @struct macro
@struct immutable A
    a::Int32
end
@test sizeof(A) == 4
@test fieldnames(A) == [:a]
@test A(1) == A(1)

@struct immutable B
    a::Int8
    b::Int32
end
@test sizeof(B) == 8
@test fieldnames(B) == [:a, :_pad_a_B_0, :_pad_a_B_1, :b]
@test B(1,2) == B(1,2)

@struct immutable C
    a::Int16
    b::Int32
    c::Int16
end
@test sizeof(C) == 12
@test fieldnames(C) == [:a, :_pad_a_C_0, :b, :c, :_pad_c_C_1]
@test C(1,2,3) == C(1,2,3)

@struct immutable D
    a::Int8
    b::Int64
end
@test sizeof(D) == 16
@test fieldnames(D) == [:a, :_pad_a_D_0, :_pad_a_D_1, :_pad_a_D_2, :b]
@test D(1,2) == D(1,2)

@struct immutable E
    a::Int64
    b::Int32
end
@test sizeof(E) == 16
@test fieldnames(E) == [:a, :b, :_pad_b_E_0]
@test E(1,2) == E(1,2)

@struct immutable F
    a::Int32
    b::Int16
    c::Int32
    d::Int32
    e::Int64
end
@test sizeof(F) == 24
@test fieldnames(F) == [:a, :b, :_pad_b_F_0, :c, :d, :e]
@test F(1,2,3,4,5) == F(1,2,3,4,5)

@struct immutable G
    a::Float64
    b::Int8
    c::Int16
    d::Int32
end
@test sizeof(G) == 24
@test fieldnames(G) == [:a, :b, :_pad_b_G_0, :c, :_pad_c_G_1, :d]
@test G(1,2,3,4) == G(1,2,3,4)

@struct immutable H
    a::Float32
    b::Int8
    c::Int16
end
@test sizeof(H) == 8
@test fieldnames(H) == [:a, :b, :_pad_b_H_0, :c]
@test H(1,2,3) == H(1,2,3)

@struct immutable I
    a::Float64
    b::Int8
    c::Int32
end
@test sizeof(I) == 16
@test fieldnames(I) == [:a, :b, :_pad_b_I_0, :_pad_b_I_1, :c]
@test I(1,2,3) == I(1,2,3)

@struct immutable J
    a::Int8
    b::A
end
@test sizeof(J) == 8
@test fieldnames(J) == [:a, :_pad_a_J_0, :_pad_a_J_1, :b_A_a]
@test J(1,A(2)) == J(1,A(2))

@struct immutable K
    a::J
    b::I
    c::J
end
@test sizeof(K) == 48
@test fieldnames(K) == [:a_J_a, :a_J__pad_a_J_0, :a_J__pad__pad_a_J_0_J_0, :a_J__pad_a_J_1, :a_J__pad__pad_a_J_1_J_1, :a_J_b_A_a, :b_I_a, :b_I_b, :b_I__pad_b_I_0, :b_I__pad__pad_b_I_0_I_0, :b_I__pad_b_I_1, :b_I__pad__pad_b_I_1_I_1, :b_I_c, :c_J_a, :c_J__pad_a_J_0, :c_J__pad__pad_a_J_0_J_0, :c_J__pad_a_J_1, :c_J__pad__pad_a_J_1_J_1, :c_J_b_A_a]
@test K(J(1,A(2)), I(3.0, 4, 5), J(6, A(7))) == K(J(1,A(2)), I(3.0, 4, 5), J(6, A(7)))


text:
function logspec(S::SeisData, k::Union{Int64,String};
  nx::Int64=1024,
  ov::Float64=0.5,
  fmin::Real=0.0,
  fmax::Real=Inf,
  fmt::String="auto")

  N_ticks = 5

  if isa(k, String)
    j = findid(k, S)
    j == 0 && error("No matching ID or no channel at the specified #")
  else
    j = k
  end

  seis = deepcopy(S[j])
  T = eltype(seis.x)
  unscale!(seis)
  demean!(seis)
  taper!(seis)
  tos = seis.t[1,2]
  fs = T(seis.fs)

  # Set some image properties
  sz = Float64[8.0, 5.0]                              # Size of image in inches
  cf = 0.05                                           # Colorbar fraction
  cp = 0.02                                           # Colorbar spacing
  x0 = 0.14
  y0 = 0.12
  tsh = 0.13
  wid = 0.90-x0
  sp = (x0, y0+tsh+cp, wid, 0.98-(y0+tsh+cp))         # Spectrogram plot bounds
  tp = (x0, y0, wid*(1.0-(cf+cp)), tsh)               # Timeseries plot bounds
  fmax = T(min(fmax, 0.5*fs))                         # Maximum frequency to plot

  # Generate spectrogram
  if fmin != 0.0
    nx = nextpow(2, 4.0/fmin)
  end
  no = floor(Int64, ov*nx)
  sg = spectrogram(seis.x, nx, no, fs=fs)
  t = Int64.(round.(time(sg).*1.0e6)) .+ tos          # Time from start (s)
  f = T.(collect(freq(sg)))                           # Frequency (Hz)
  P = T(10.0)*log10.(power(sg))                       # Power (dB)

  # x ticks
  xti = first(t) : div(last(t) - first(t), N_ticks-1) : last(t)
  xtl = fill!(Array{String,1}(undef, length(xti)), "")

  # truncate P, y
  i0 = max(2, findfirst(f.≥max(0.5*fs/nx, fmin)))     # First sane frequency index
  i1 = findlast(f.≤fmax)                              # Last sane frequency index
  P = P[i0:i1, :]
  y = log10.(f)[i0:i1]

  # Remove all -∞
  mP = maximum(P)
  pτ = mP - T(60.0) - eps(T)
  P[P.<pτ] .= pτ

  # y ticks
  f0 = log10(f[i0])                                   # Minimum frequency to show
  f1 = log10(f[i1])                                   # Maximum frequency to show
  yt = collect(ceil(f0) : 0.5 : floor(f1))
  yti = T(10.0) .^ vcat(f0, yt, f1)
  if fmin < 0.1
    ytl = vcat("", [@sprintf("%0.1e", T(10.0)^i) for i in yt], "")
  else
    ytl = vcat("", [@sprintf("%0.1f", T(10.0)^i) for i in yt], "")
  end

  # Figure
  h = PyPlot.figure(figsize = sz, dpi=150)

  # Seismogram plot
  ax_ts = PyPlot.axes(tp)
  dt    = one(T)/fs
  tx    = t_expand(seis.t, seis.fs)
  xmax  = 1.0*maximum(abs.(seis.x))
  ts_xl = @sprintf("%0.2e", xmax)
  PyPlot.plot(tx, seis.x, linewidth=1.0, color="k")
  PyPlot.ylim(xmax .*(T(-1.05), T(1.05)))
  xfmt(round(Int64, first(xti)), round(Int64, last(xti)), fmt, N_ticks)
  PyPlot.xlabel(ax_ts.get_xlabel(), fontsize=12.0, color="black", fontweight="bold", family="serif")
  PyPlot.yticks([-xmax, zero(T), xmax], ["-"*ts_xl, "0.00e00", "+"*ts_xl])
  PyPlot.setp(ax_ts.get_yticklabels(), fontsize=9.0, color="black", fontweight="bold", family="serif")
  PyPlot.setp(ax_ts.get_xticklabels(), fontsize=9.0, color="black", fontweight="bold", family="serif")

  # Note: imshow won't work with unevenly-spaced y-ticks
  ax_sp = PyPlot.axes(sp)
  img = pcolormesh(t, y, P, cmap="Spectral_r", shading="goraud")
  PyPlot.ylabel("Frequency [Hz]", fontweight="bold", fontsize=12.0, family="serif", color="black")
  PyPlot.ylim(f0, f1)
  PyPlot.yticks(log10.(yti), ytl)
  PyPlot.setp(ax_sp.get_yticklabels(), fontsize=9.0, color="black", fontweight="bold", family="serif")
  PyPlot.xlim(minimum(t), maximum(t))
  PyPlot.xticks(xti, xtl)

  # Set color limits and draw colorbar
  cl = [mP-60.0, mP]
  PyPlot.clim(cl)
  cb = PyPlot.colorbar(label="Power [dB]", use_gridspec=false, fraction=cf, pad=cp)
  PyPlot.setp(cb.ax.get_yticklabels(), fontsize=9.0, color="black", fontweight="bold", family="serif")
  cb.set_label("Power [dB]", fontsize=10.0, color="black", fontweight="bold", family="serif")
  return h
end


text: Dict{String, Any}("name" => Dict{String, Any}("first" => Dict{String, Any}("value" => "Arthur", "type" => "string"), "last" => Dict{String, Any}("value" => "Dent", "type" => "string")), "arr" => Any[Dict{String, Any}("a" => Dict{String, Any}("b" => Dict{String, Any}("c" => Dict{String, Any}("value" => "1", "type" => "integer"), "d" => Dict{String, Any}("value" => "2", "type" => "integer")))), Dict{String, Any}("a" => Dict{String, Any}("b" => Dict{String, Any}("c" => Dict{String, Any}("value" => "3", "type" => "integer"), "d" => Dict{String, Any}("value" => "4", "type" => "integer"))))], "many" => Dict{String, Any}("dots" => Dict{String, Any}("here" => Dict{String, Any}("dot" => Dict{String, Any}("dot" => Dict{String, Any}("dot" => Dict{String, Any}("value" => "42", "type" => "integer")))))), "count" => Dict{String, Any}("f" => Dict{String, Any}("value" => "6", "type" => "integer"), "c" => Dict{String, Any}("value" => "3", "type" => "integer"), "e" => Dict{String, Any}("value" => "5", "type" => "integer"), "b" => Dict{String, Any}("value" => "2", "type" => "integer"), "a" => Dict{String, Any}("value" => "1", "type" => "integer"), "h" => Dict{String, Any}("value" => "8", "type" => "integer"), "i" => Dict{String, Any}("value" => "9", "type" => "integer"), "d" => Dict{String, Any}("value" => "4", "type" => "integer"), "g" => Dict{String, Any}("value" => "7", "type" => "integer"), "j" => Dict{String, Any}("value" => "10", "type" => "integer"), "k" => Dict{String, Any}("value" => "11", "type" => "integer"), "l" => Dict{String, Any}("value" => "12", "type" => "integer")), "tbl" => Dict{String, Any}("a" => Dict{String, Any}("b" => Dict{String, Any}("c" => Dict{String, Any}("value" => "42.666", "type" => "float")))), "a" => Dict{String, Any}("few" => Dict{String, Any}("dots" => Dict{String, Any}("polka" => Dict{String, Any}("dot" => Dict{String, Any}("value" => "again?", "type" => "string"), "dance-with" => Dict{String, Any}("value" => "Dot", "type" => "string"))))))


text:
mutable struct DummyState end
mutable struct DummyAct end
mutable struct DummyObs end
mutable struct DummyPOMDP <: POMDP{DummyState, DummyAct, DummyObs} end

function test_requirements()
    sim = MCVISimulator()
    solver = MCVISolver(sim, nothing, 1, 100, 8, 500, 1000, 5000, 50,
                        nothing, nothing)
    println("There should be nothing implemented but the default from POMDPs.jl: \n")
    @test_throws MethodError @requirements_info solver DummyPOMDP()

    solver = MCVISolver(sim, nothing, 1, 100, 8, 500, 1000, 5000, 50,
                        LightDark1DLowerBound(sim.rng), LightDark1DUpperBound(sim.rng))
    println("Everything should be implemented: \n")
    @requirements_info solver LightDark1D()
end


text:
# ###################################################################################################
# Simulation studies of a bimodal distribution
#
# author: weiya
# date: Dec 25, 2018
#
# refer to Section 11.5.1 of
# Liu, J. S. (2008). Monte Carlo strategies in scientific computing (2. ed). New York, NY: Springer.
#
# ###################################################################################################
using LinearAlgebra
using Distributions
using StatsBase
using Plots
using LaTeXStrings

Σ1 = Matrix(1.0I, 2, 2)
Σ2 = [1 0.9; 0.9 1]
Σ3 = [1 -0.9; -0.9 1]
invΣ2 = [1 -0.9; -0.9 1] ./ 0.19
invΣ3 = [1 0.9; 0.9 1] ./ 0.19
μ1 = [0, 0]
μ2 = [-6, -6]
μ3 = [4, 4]
w = [0.34, 0.33, 0.33]

function targetpdf(x)
    model1 = MvNormal(μ1, Σ1)
    model2 = MvNormal(μ2, Σ2)
    model3 = MvNormal(μ3, Σ3)
    return pdf(model1, x) * w[1] + pdf(model2, x) * w[2] + pdf(model3, x) * w[3]
end

function mycontour()
    # refer to https://docs.juliaplots.org/latest/examples/gr/#contours
    x = -10:0.5:10
    y = -10:0.5:10
    X = repeat(reshape(x, 1, :), length(y), 1)
    Y = repeat(y, 1, length(x))
    Z = copy(X)
    for i in 1:size(X,1)
        for j in 1:size(X,2)
            Z[i,j] = targetpdf([X[i,j], Y[i,j]])
        end
    end
    contour(x, y, Z)
    savefig("contour.png")
end

function propose(x; a = 4)
    θ = rand() * 2π
    r = rand() * 4
    return [x[1] + r*cos(θ), x[2] + r*sin(θ)]
end

function iMH(M = 200000)
    # start points
    x = [rand() - 0.5, rand() - 0.5]
    X = ones(M, 2)
    num = 0
    for m = 1:M
        xs = propose(x)
        ρ = targetpdf(xs) / targetpdf(x)
        if rand() < ρ
            x = xs
            num += 1
        end
        X[m,:] = x
    end
    return X, num/M
end

function propose2(x; a = 4)
    r = rand() * 2a - a
    return x + r
end

function iMH2(M = 200000)
    # start points
    x1 = rand() - 0.5
    x2 = rand() - 0.5
    X = ones(M, 2)
    num = 0
    for m = 1:M
        x1s = propose2(x1)
        ρ = targetpdf([x1s, x2]) / targetpdf([x1, x2])
        if rand() < ρ
            x1 = x1s
            num += 1
        end
        x2s = propose2(x2)
        ρ = targetpdf([x1, x2s]) / targetpdf([x1, x2])
        if rand() < ρ
            x2 = x2s
        end
        X[m,:] = [x1, x2]
    end
    return X, num/M # why acceptance ratio is two times than iMH. Got it, because I miss the direction random selection.
end

# derivative of targetpdf
function dm(x; step = 0.05, nmaxstep = 100)
    term1 = w[1] / 2π * exp(-0.5*sum(x.^2))
    term2 = w[2] / (2π * 0.19) * exp( -0.5 * transpose(x.+6) * invΣ2 * (x.+6) )
    term3 = w[3] / (2π * 0.19) * exp( -0.5 * transpose(x.-4) * invΣ3 * (x.-4) )
    term = term1 .+ term2 * invΣ2 + term3 * invΣ3
    # negative of gradient direction (remove the minus symbol already)
    dx1 = term[1,1]*x[1] + term[1,2]*x[2]
    dx2 = term[2,1]*x[1] + term[2,2]*x[2]
    dx = [dx1, dx2]
    # find local maximum
    bst = targetpdf(x)
    i = 0
    while i < nmaxstep
        i += 1
        cur = targetpdf(x .- step * dx)
        if cur > bst
            bst = cur
            x = x .- step * dx
        else
            break
        end
    end
    return x
end

function cgmc(M = 200000)
    # start points (each col is a stream)
    x = rand(2, 2) .- 0.5
    X = ones(M, 4)
    for m = 1:M
        # randomly choose xa
        xaIdx = sample([1,2])
        xa = x[:,xaIdx]
        # find gradient of pdf at xa (or conjugate gradient of log pdf at xa)
        y = dm(xa)
        xcIdx = 3 - xaIdx
        xc = x[:,xcIdx]
        e = y - xc
        e = e / norm(e)
        # sample r
        r, = MTM(y, e)
        # update xc
        x[:,xcIdx] = y .+ r*e
        X[m, :] .= vec(x)
    end
    return X
end

f(r, y, e) = abs(r) * targetpdf(y.+r*e)

function MTM(anchor, direction;M = 100, k = 5, σ = 10)
    # use the simplified version: OBMC algorithm
    # initial x
    x = 2*rand()-1
    num = 0
    for m = 1:M
        y = ones(k)
        wy = ones(k)
        for i = 1:k
            y[i] = rand(Normal(x, σ))
            wy[i] = f(y[i], anchor, direction)
        end
        #wy = wy ./ sum(wy)
        yl = sample(y, pweights(wy))
        # draw reference points
        xprime = ones(k)
        xprime[k] = x
        wxprime = ones(k)
        wxprime[k] = f(x, anchor, direction)
        for i = 1:k-1
            xprime[i] = rand(Normal(yl, σ))
            wxprime[i] = f(xprime[i], anchor, direction)
        end
        # accept or not
        ρ = sum(wy) / sum(wxprime)
        if rand() < ρ
            x = yl
            num += 1
        end
    end
    return x, num/M
end

# run
res, = iMH()
X = cgmc()
p1 = histogram(res[100000:end,1], xlab = L"x_1")
p2 = histogram(res[100000:end,2], xlab = L"x_2")
p3 = histogram(X[100000:end,1], xlab = L"x_1")
p4 = histogram(X[100000:end,2], xlab = L"x_2")

p5 = plot(autocor(res[100000:end,1]))
p6 = plot(autocor(X[100000:end,1]))
plot(p1, p3, p5, p6, layout = (2,2), legend = false)

max_stars_repo_name: szcf-weiya/MCMC

max_issues_repo_name: szcf-weiya/MCMC

max_forks_repo_name: szcf-weiya/MCMC


text:
using Lehmann
using Documenter

DocMeta.setdocmeta!(Lehmann, :DocTestSetup, :(using Lehmann); recursive = true)

makedocs(;
    modules = [Lehmann],
    authors = "Kun Chen, Tao Wang, Xiansheng Cai",
    repo = "https://github.com/numericalEFT/Lehmann.jl/blob/{commit}{path}#{line}",
    sitename = "Lehmann.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://numericaleft.github.io/Lehmann.jl",
        assets = String[]
    ),
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "manual/kernel.md"
        ],
        "API reference" => Any[
            "lib/dlr.md",
            "lib/spectral.md",
            "lib/discrete.md",
            "lib/functional.md",
            "lib/sample.md",
            "lib/utility.md",
            # map(s -> "lib/$(s)", sort(readdir(joinpath(@__DIR__, "src/lib"))))
            # "Internals" => map(s -> "lib/$(s)", sort(readdir(joinpath(@__DIR__, "src/lib"))))
        ]
    ]
)

deploydocs(;
    repo = "github.com/numericalEFT/Lehmann.jl",
    devbranch = "main"
)


text:
# Test struct for `rand_tangent` and `difference`.
struct Bar
    a::Float64
    b::Int
    c::Any
 end
@testset "rand_tangent" begin
    rng = MersenneTwister(123456)

    @testset "Primal: $(typeof(x)), Tangent: $T_tangent" for (x, T_tangent) in [

        # Things without sensible tangents.
        ("hi", NoTangent),
        ('a', NoTangent),
        (:a, NoTangent),
        (true, NoTangent),
        (4, NoTangent),
        (FiniteDifferences, NoTangent),  # Module object
        # Types (not instances of type)
        (Bar, NoTangent),
        (Union{Int, Bar}, NoTangent),
        (Union{Int, Bar}, NoTangent),
        (Vector, NoTangent),
        (Vector{Float64}, NoTangent),
        (Integer, NoTangent),
        (Type{<:Real}, NoTangent),

        # Numbers.
        (5.0, Float64),
        (5.0 + 0.4im, Complex{Float64}),
        (big(5.0), BigFloat),

        # StridedArrays.
        (fill(randn(Float32)), Array{Float32, 0}),
        (fill(randn(Float64)), Array{Float64, 0}),
        (randn(Float32, 3), Vector{Float32}),
        (randn(Complex{Float64}, 2), Vector{Complex{Float64}}),
        (randn(5, 4), Matrix{Float64}),
        (randn(Complex{Float32}, 5, 4), Matrix{Complex{Float32}}),
        ([randn(5, 4), 4.0], Vector{Any}),

        # Co-Arrays
        (randn(5)', Adjoint{Float64, Vector{Float64}}),  # row-vector: special
        (randn(5, 4)', Matrix{Float64}),                 # matrix: generic dense

        (transpose(randn(5)), Transpose{Float64, Vector{Float64}}),  # row-vector: special
        (transpose(randn(5, 4)), Matrix{Float64}),                   # matrix: generic dense

        # AbstactArrays of non-perturbable types
        (1:10, NoTangent),
        (1:2:10, NoTangent),
        ([false, true], NoTangent),

        # Tuples.
        ((4.0, ), Tangent{Tuple{Float64}}),
        ((5.0, randn(3)), Tangent{Tuple{Float64, Vector{Float64}}}),
        ((false, true), NoTangent),
        (Tuple{}(), NoTangent),

        # NamedTuples.
        ((a=4.0, ), Tangent{NamedTuple{(:a,), Tuple{Float64}}}),
        ((a=5.0, b=1), Tangent{NamedTuple{(:a, :b), Tuple{Float64, Int}}}),
        ((a=false, b=true), NoTangent),
        ((;), NoTangent),

        # structs.
        (Bar(5.0, 4, rand(rng, 3)), Tangent{Bar}),
        (Bar(4.0, 3, Bar(5.0, 2, 4)), Tangent{Bar}),
        (sin, NoTangent),
        # all fields NoTangent implies NoTangent
        (Pair(:a, "b"), NoTangent),
        (CartesianIndex(2, 3), NoTangent),

        # LinearAlgebra types
        (
            UpperTriangular(randn(3, 3)),
            UpperTriangular{Float64, Matrix{Float64}},
        ),
        (
            Diagonal(randn(2)),
            Diagonal{Float64, Vector{Float64}},
        ),
        (
            Symmetric(randn(2, 2)),
            Symmetric{Float64, Matrix{Float64}},
        ),
        (
            Hermitian(randn(ComplexF64, 1, 1)),
            Hermitian{ComplexF64, Matrix{ComplexF64}},
        ),
    ]
        @test rand_tangent(rng, x) isa T_tangent
        @test rand_tangent(x) isa T_tangent
    end

    @testset "erroring cases" begin
        # Ensure struct fallback errors for non-struct types.
        @test_throws ArgumentError invoke(rand_tangent, Tuple{AbstractRNG, Any}, rng, 5.0)
    end

    @testset "compsition of addition" begin
        x = Bar(1.5, 2, Bar(1.1, 3, [1.7, 1.4, 0.9]))
        @test x + rand_tangent(x) isa typeof(x)
        @test x + (rand_tangent(x) + rand_tangent(x)) isa typeof(x)
    end

    # Julia 1.6 changed to using Ryu printing algorithm and seems better at printing short
    VERSION >= v"1.6" && @testset "niceness of printing" begin
        rng = MersenneTwister(1)
        for i in 1:50
            @test length(string(rand_tangent(rng, 1.0))) <= 6
            @test length(string(rand_tangent(rng, 1.0 + 1.0im))) <= 12
            @test length(string(rand_tangent(rng, 1f0))) <= 9
            @test length(string(rand_tangent(rng, big"1.0"))) <= 9
        end
    end
end


text:
@testset "Peak" begin

    include("gaussian.jl")

    include("motif.jl")

end


text:
"""
    size(A) -> Tuple
    size(A, dim) -> Union{Int,StaticInt}

Returns the size of each dimension of `A` or along dimension `dim` of `A`. If the size of
any axes are known at compile time, these should be returned as `Static` numbers. Otherwise,
`ArrayInterface.size(A)` is identical to `Base.size(A)`

```julia
julia> using StaticArrays, ArrayInterface

julia> A = @SMatrix rand(3,4);

julia> ArrayInterface.size(A)
(static(3), static(4))
```
"""
@inline size(A) = map(static_length, axes(A))
size(x::SubArray) = eachop(_sub_size, to_parent_dims(x), x.indices)
_sub_size(x::Tuple, ::StaticInt{dim}) where {dim} = static_length(getfield(x, dim))
@inline size(B::VecAdjTrans) = (One(), length(parent(B)))
@inline size(B::MatAdjTrans) = permute(size(parent(B)), to_parent_dims(B))
@inline function size(B::PermutedDimsArray{T,N,I1}) where {T,N,I1}
    permute(size(parent(B)), static(I1))
end
function size(a::ReinterpretArray{T,N,S,A}) where {T,N,S,A}
    psize = size(parent(a))
    if _is_reshaped(typeof(a))
        if sizeof(S) === sizeof(T)
            return psize
        elseif sizeof(S) > sizeof(T)
            return (static(div(sizeof(S), sizeof(T))), psize...)
        else
            return tail(psize)
        end
    else
        return (div(first(psize) * static(sizeof(S)), static(sizeof(T))), tail(psize)...,)
    end
end
size(A::ReshapedArray) = Base.size(A)
size(A::AbstractRange) = (static_length(A),)

size(a, dim) = size(a, to_dims(a, dim))
size(a::Array, dim::Integer) = Base.arraysize(a, convert(Int, dim))
function size(a::A, dim::Integer) where {A}
    if parent_type(A) <: A
        len = known_size(A, dim)
        if len === nothing
            return Int(length(axes(a, dim)))
        else
            return StaticInt(len)
        end
    else
        return size(a)[dim]
    end
end
function size(A::SubArray, dim::Integer)
    pdim = to_parent_dims(A, dim)
    if pdim > ndims(parent_type(A))
        return size(parent(A), pdim)
    else
        return static_length(A.indices[pdim])
    end
end

"""
    known_size(::Type{T}) -> Tuple
    known_size(::Type{T}, dim) -> Union{Int,Nothing}

Returns the size of each dimension of `A` or along dimension `dim` of `A` that is known at
compile time. If a dimension does not have a known size along a dimension then `nothing` is
returned in its position.
"""
known_size(x) = known_size(typeof(x))
known_size(::Type{T}) where {T} = eachop(_known_size, nstatic(Val(ndims(T))), axes_types(T))
_known_size(::Type{T}, dim::StaticInt) where {T} = known_length(_get_tuple(T, dim))
@inline known_size(x, dim) = known_size(typeof(x), dim)
@inline known_size(::Type{T}, dim) where {T} = known_size(T, to_dims(T, dim))
@inline function known_size(::Type{T}, dim::CanonicalInt) where {T}
    if ndims(T) < dim
        return 1
    else
        return known_size(T)[dim]
    end
end


text:
"""
`SPParams` holds the algorithm parameters for a spatial pooler with nomenclature
similar to [source](https://www.frontiersin.org/articles/10.3389/fncom.2017.00111/full)

The dimension parameters are problem-specific and should be the first to be specified.

The tuning parameters have sensible defaults, which should be .
All gated features are enabled by default.

# Parameters

## Dimensions
- `szᵢₙ = (32,32)`: input dimensions
- `szₛₚ = (32,32)`: output dimensions
- `γ = 6`: receptive field radius (how large an input area an output minicolumn maps to). Must be <= min(szᵢₙ)

## Algorithm tuning
- `s = .02`: average output sparsity
- `prob_synapse = .5`: probability for each element of the `szᵢₙ × szₛₚ` space to be a synapse.
  Elements that roll below this value don't form a synapse and don't get a permanence value.
  If this is very low, the proximal synapses matrix can become sparse.
- `θ_permanence01 = .5`: synapse permanence connection threshold
- `p⁺_01 = .1 , p⁻_01 = .02`: synapse permanence adaptation rate (see [`ProximalSynapses`](@ref))
- `θ_stimulus_activate = 1`: minicolumn absolute activation threshold
- `Tboost = 200.0`: boosting mechanism's moving average filter period
- `β = 1.0`: boosting strength

## Feature gates
- `enable_local_inhibit = true`
- `enable_learning = true`
- `enable_boosting = true`
"""
@with_kw struct SPParams
  # dimensions
  szᵢₙ::Union{NTuple{N,Int} where N,Int}     = 1000; @assert all(szᵢₙ.>0)
  szₛₚ::Union{NTuple{N,Int} where N,Int}     = 5000; @assert all(szₛₚ.>0)
  γ::Int                    = 6;       @assert 0<γ<=minimum(szᵢₙ)

  # tuning
  s::Float32                = .02;     @assert s>0
  prob_synapse::Float32     = .5;      @assert 0<=prob_synapse<=1
  θ_permanence01::Float32   = .5;      @assert 0<=θ_permanence01<=1
  p⁺_01::Float32            = .1;      @assert 0<=p⁺_01<=1
  p⁻_01::Float32            = .02;     @assert 0<=p⁻_01<=1
  θ_permanence::𝕊𝕢        = round(𝕊𝕢, θ_permanence01*typemax(𝕊𝕢))
  p⁺::𝕊𝕢                  = round(𝕊𝕢, p⁺_01*typemax(𝕊𝕢))
  p⁻::𝕊𝕢                  = round(𝕊𝕢, p⁻_01*typemax(𝕊𝕢))
  θ_stimulus_activate::Int  = 1;       @assert θ_stimulus_activate>=0
  Tboost::Float32           = 200;     @assert Tboost>0
  β::Float32                = 1;       @assert β>0

  # feature gates
  enable_local_inhibit::Bool= true
  enable_learning::Bool     = true
  enable_boosting::Bool     = true
  @assert zero(𝕊𝕢)<=θ_permanence<=typemax(𝕊𝕢)
  @assert zero(𝕊𝕢)<=p⁺<=typemax(𝕊𝕢)
  @assert zero(𝕊𝕢)<=p⁻<=typemax(𝕊𝕢)
end

"""
`TMParams` holds the algorithm parameters for a Temporal Memory with nomenclature
similar to [source]()

## Dimensions
- `Nc = 2500`: number of columns
- `k = 10`: cells per column
- `Nₙ = k Nc` neurons in layer. The `Nₛ` number of dendritic segments
  is variable

## Tuning
- `p⁺_01 = .12 , p⁻_01 = .04 ∈ [0,1]`: synapse permanence adaptation rate (see [`ProximalSynapses`](@ref))
- `LTD_p⁻_01 = .002 ∈ [0,1]`: synapse long term depression rate
- `θ_permanence = .5*typemax(𝕊𝕢) ∈ 𝕊𝕢`: synapse permanence connection threshold
- `init_permanence = .4*typemax(𝕊𝕢) ∈ 𝕊𝕢`: permanence of a newly-grown synapse
- `synapseSampleSize = 25 ∈ ℕ`: target number of matching synapses per dendrite.
  Represents how many bits the dendrite targets to recognize the input.
  Dendrites with fewer synapses matching the input might grow new synapses.
- `θ_stimulus_activate = 14 ∈ ℕ`: number of matching synapses needed to depolarize the dendrite
- `θ_stimulus_learn = 12 ∈ ℕ`: number of matching synapses that are insufficient to depolarize the
  dendrite, but sufficient to trigger learning. `θ_stimulus_learn <= θ_stimulus_activate`

## Feature gates
- `enable_learning = true`
"""
@with_kw struct TMParams
  # dimensions
  Nc::Int                  = 5000;    @assert Nc>0
  k::Int                   = 10;      @assert k>0
  Nₙ::Int                  = k*Nc;    @assert Nₙ>0

  # tuning
  p⁺_01::Float32           = .12;     @assert 0<=p⁺_01<=1
  p⁻_01::Float32           = .04;     @assert 0<=p⁻_01<=1
  LTD_p⁻_01::Float32       = .002;    @assert 0<=LTD_p⁻_01<=1
  p⁺::𝕊𝕢                 = round(𝕊𝕢,p⁺_01*typemax(𝕊𝕢))
  p⁻::𝕊𝕢                 = round(𝕊𝕢,p⁻_01*typemax(𝕊𝕢))
  LTD_p⁻::𝕊𝕢             = round(𝕊𝕢,LTD_p⁻_01*typemax(𝕊𝕢))
  θ_permanence::𝕊𝕢       = round(𝕊𝕢,.5typemax(𝕊𝕢))
  init_permanence::𝕊𝕢    = round(𝕊𝕢,.4typemax(𝕊𝕢))
  synapseSampleSize::Int   = 25;      @assert synapseSampleSize>0
  θ_stimulus_activate::Int = 14;      @assert θ_stimulus_activate>0
  θ_stimulus_learn::Int    = 12;      @assert 0 < θ_stimulus_learn <= θ_stimulus_activate

  # feature gates
  enable_learning::Bool    = true

  @assert Nc*k == Nₙ
end

# Created from TMParams
@with_kw struct DistalSynapseParams
  p⁺::𝕊𝕢                 = round(𝕊𝕢,p⁺_01*typemax(𝕊𝕢))
  p⁻::𝕊𝕢                 = round(𝕊𝕢,p⁻_01*typemax(𝕊𝕢))
  LTD_p⁻::𝕊𝕢             = round(𝕊𝕢,LTD_p⁻_01*typemax(𝕊𝕢))
  θ_permanence::𝕊𝕢       = round(𝕊𝕢,.5typemax(𝕊𝕢))
  init_permanence::𝕊𝕢    = round(𝕊𝕢,.4typemax(𝕊𝕢))
  synapseSampleSize::Int   = 25;      @assert synapseSampleSize>0
  θ_stimulus_learn::Int    = 12;      @assert 0 < θ_stimulus_learn
end
DistalSynapseParams(tmParams::TMParams)= begin
  @unpack p⁺, p⁻, LTD_p⁻, θ_permanence,
      init_permanence, synapseSampleSize, θ_stimulus_learn = tmParams
  DistalSynapseParams(p⁺, p⁻, LTD_p⁻, θ_permanence,
      init_permanence, synapseSampleSize, θ_stimulus_learn)
end


text: f() = 1


text:
using Dates, Gumbo
#putting together progress reports

"""
g10_students::Dict{Int64,String}

dictionary containing g10 student names/ids
"""
g10_students = Dict(XXXXX1=>"student1",XXXXX2=>"student2",XXXXX3=>"student3",XXXXX4=>"student4 ",XXXXX5=>"student5",XXXXX6=>"student6",XXXXX7=>"student7",XXXXX8=>"student8",XXXXX9=>"student9")
"""
g11_students::Dict{Int64,String}

dictionary containing g11 student names/ids
"""
g11_students = Dict(XXXX1X=>"name1",XXXX2X=>"name2",XXXX3X=>"name3",XXXX4X=>"name4",XXXX5X=>"name5",XXXX6X=>"name6")

#running on wsl - can change
"""
header(date::String,f::IOStream)

Adds the header of the Keystone html files to IOStream.
Includes report run date.
See also: createMaster
"""
function header(date::String,f::IOStream)
write(f,"""

920 Central Road, Bloomsburg, PA 17815
Phone (800) 255-4937
Fax (570) 300-2346
info@keystonehighschool.com - www.keystonehighschool.com

	<input type="button" id="bt" value="Click to download notes about the progress reports" onclick="saveFile()"/>

""")
end
function seriesHeader(f::IOStream)
	write(f,"""

920 Central Road, Bloomsburg, PA 17815
Phone (800) 255-4937
Fax (570) 300-2346
info@keystonehighschool.com - www.keystonehighschool.com

""")
end
function imagelessHeader(f::IOStream)
	write(f,"""

""")
end

function createReport(students::Dict{Int64,String},grouping::String)

	f = open("$(grouping).html","w")

	header("12/1/2020",f)

	write(f,"
$(grouping)'s progress
")
	for value in keys(students)
		ID = "$value"
		docname = "keystone_progress_report_$ID"*".html"
		try
			file = parsehtml(read(docname,String))
			write(f,repr(Ref(file.root[2][1][2][2])))
			write(f,repr(Ref(file.root[2][1][2][4])))
			write(f,"

")
		catch e
			write(f,"data for $(students[value]) not found

")
			continue
		end
	end

	close(f)
	println("$(grouping).html succesfully generated. ")

end

function createMaster()
	#make sure the working directory is correct

	f = open("master_progress_current.html","w")
	report_log = open("master_progress_current.log","w")

	#adds the logo, info, and style information at the top
	header("12/8/2020",f)

	print(f,"
G10's Progress
")
	print(report_log,"...Writing G10's progres...")
	for value in keys(g10_students)
		ID = "$value"
		docname = "keystone_progress_report_$ID"*".html"
		try
			file = parsehtml(read(docname,String))
			println("testing if $(g10_students[value]) has data available")

			str1 = repr(file.root[2][1][2][2])[21:end]
			str2 = repr(file.root[2][1][2][4])[21:end]
			println("data found!")
			println(f,"""
$(g10_students[value]) data:
\n""")
			println(f,str1)
			println(f,str2)
			println(f,"

")
		catch e
			println(f,"data for $(g10_students[value]) not found

")
			println(report_log,"juliadata for $(g10_students[value]) not found")
			continue
		end
	end

	write(f,"
G11's Progress
")

		for value in keys(g11_students)
		ID = "$value"
		docname = "keystone_progress_report_$ID"*".html"
		try
			file = parsehtml(read(docname,String))
			println("testing if $(g10_students[value]) has data available")

			str1 = repr(file.root[2][1][2][2])[21:end]
			str2 = repr(file.root[2][1][2][4])[21:end]
			println("data found!")
			println(f,"""
$(g10_students[value]) data:
\n""")
			println(f,str1)
			println(f,str2)
			println(f,"

")
		catch e
			write(f,"data for $(g11_students[value]) not found

")
			println("data for $(g11_students[value]) not found")
			println(report_log,"juliadata for $(g11_students[value]) not found")
			continue
		end
	end
	println(f,"
See source code at github.com/thomastjdavis/K12ReportGenerator
")
	close(f)
	println("\n Master progress report assembled. ")

	end
function HTMLtables(f::HTMLDocument)
	tbls = []
	for elem in PostOrderDFS(f.root)
		#finds tables in the HTMLDocument
		if (typeof(elem) == HTMLElement{:table})
			tbl = PreOrderDFS(elem)
			for item in tbl
				row = PostOrderDFS(item)
				push!(tbls,row)
				break
			end
		end
	end
	tbls
end
function skipCols(j::Int64,f::HTMLDocument,I::Array{Int64,1},output::IOStream)
	tbls = HTMLtables(f)
	rowTag = ["td","th"]
	n=0
	vect = collect(tbls[j])

	#perhaps bad code, but finds the number of columns of the table
	for elem in vect
		if !(typeof(elem)==HTMLElement{:tr})
			if !(typeof(elem)==HTMLText)
				if(tag(elem)==:td)||(tag(elem)==:th)
					n+=1
				end
			else
				continue
			end
		else
			break
		end
	end

	for colN in I
		if (n<=colN)||(colN<0)
			@warn "Column number $(colN) cannot be taken out of the table"
		end
	end

	counter = 1
	rowCounter = 0
	colCounter = 0
	len = n - length(I)

	#prints only the correct columns into IOStream f
	#eg. I = [1,4,2]
	println(output,"""
\n
""")
	for elem in vect
		#println(typeof(elem))
		if !(typeof(elem)==HTMLElement{:tr})
			if !(typeof(elem)==HTMLText)
				if(tag(elem)==:td)||(tag(elem)==:th)
					#println(colCounter)
					if(inA(colCounter,I))
						colCounter+=1
						continue
					end
					if(colCounter == (n))
						println(output,"""
\n""")
						continue
					end
					colCounter+=1
					println(output,elem)
				end
			else
				continue
			end
		else
			colCounter =0
			println(output,"""
\n
""")
			rowCounter +=1
		end
	end

end
function generateSeries(ID::Int64)
	if Sys.iswindows()
		println("Run in Linux plz")
		return
	end

	function getReport()
		try
			doc = parsehtml(read("keystone_progress_report_$(ID).html",String))
		catch
			println("report caching failed for ID =$(ID)")
		end
	end

	filename = "test"*"$ID"*".html"
	io = open(filename,"w")
	seriesHeader(io)
	#print off appropriate header for file

	for date in dates
		dir = "/mnt/c/users/cheri/desktop/workCode/k12Progress/progress_reports/"*"$date"
		stdData = cd(getReport,dir)
		println(io,"""

$(date)'s report

Percent complete tolerance =  $(tolDates[date])

		""")
		if typeof(stdData)==Nothing
			continue
		end
		println(io,repr(HTMLtables(stdData)[2])[80:end-1])
		println(io,"""
""")
	end
	println(io,"""""")
	close(io)
	println("done: generated")
end
createMaster()


text:
using MICE
using Test
import CSV
using Glob
using FilePathsBase
using FilePathsBase: /
using DataFrames

Glob.glob(pattern, path::AbstractPath) = glob(pattern, string(path))

@testset "MICE.jl" begin
    input_file_paths = Path.(glob("*.csv", cwd() / "comparison" / "data" / "input"))
    @test length(input_file_paths) > 0 # Find input files
    input_datasets = Dict()
    for input_file_path in input_file_paths
        file_name = filename(input_file_path)
        input_datasets[file_name] = CSV.File(input_file_path, missingstring = "NA") |> DataFrame
    end

    output_datasets = Dict()
    for (dataset_name, datset) in input_datasets
        output_datasets[dataset_name] = mice(datset; max_iterations = 1)
    end
end


text:
export CoupledSimulation, step!, run!, name

"""
    AbstractSimulation

An abstract type representing a model simulation.
"""
abstract type AbstractSimulation end

abstract type AbstractAtmosSimulation <: AbstractSimulation end
name(::AbstractAtmosSimulation) = :atmos

abstract type AbstractOceanSimulation <: AbstractSimulation end
name(::AbstractOceanSimulation) = :ocean

abstract type AbstractLandSimulation <: AbstractSimulation end
name(::AbstractLandSimulation) = :land

abstract type AbstractCoupledSimulation <: AbstractSimulation end
name(::AbstractCoupledSimulation) = :coupled

struct CoupledSimulation{CS, S, CPL, L, C} <: AbstractCoupledSimulation
    "The coupled time-stepping scheme"
    coupler_solver::CS
    "The component simulations"
    simulations::S
    "The coupler"
    coupler::CPL
    "Diagnostic logger"
    logger::L
    "Clock"
    clock::C
end

"""
    run!(::CoupledSimulation)

A simple outer timestepping loop for coupled system runs.

This will be formalized when the run! functionality for component
models is implemented so to have a consistent interface.
"""
function run!(sim::CoupledSimulation)
    clock = sim.clock
    while !stop_time_exceeded(clock)
        step!(sim, clock.dt)
        tick!(clock)
    end
end

"""
    step!(sim, dt)

Advances a simulation by `dt`.

Note that `dt` is not necessarily the simulation's timestep length;
a simuation could take several shorter steps that total to `dt`.
"""
function step!(sim::AbstractSimulation, dt) end

max_stars_repo_name: CliMA/ClimaCoupler

max_issues_repo_name: CliMA/ClimaCoupler

max_forks_repo_name: CliMA/ClimaCoupler


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using ContactImplicitMPC
using LinearAlgebra

# ## Simulation
s_sim = get_simulation("hopper_2D", "piecewise1_2D_lc", "piecewise", approx = true);
s = get_simulation("hopper_2D", "flat_2D_lc", "flat");
model = s.model
env = s.env

# ## Reference Trajectory
ref_traj = deepcopy(ContactImplicitMPC.get_trajectory(s.model, s.env,
    joinpath(module_dir(), "src/dynamics/hopper_2D/gaits/gait_forward.jld2"),
    load_type = :joint_traj));
H = ref_traj.H
h = ref_traj.h

# ## MPC setup
N_sample = 5
H_mpc = 10
h_sim = h / N_sample
H_sim = 5000 #100*H*N_sample #500*H*N_sample
κ_mpc = 2.0e-4

obj = TrackingObjective(model, env, H_mpc,
    q = [Diagonal(1.0e-1 * [0.1,3,1,3])   for t = 1:H_mpc],
    u = [Diagonal(1.0e-0 * [1e-3, 1e0]) for t = 1:H_mpc],
    γ = [Diagonal(1.0e-100 * ones(model.nc)) for t = 1:H_mpc],
    b = [Diagonal(1.0e-100 * ones(model.nc * friction_dim(env))) for t = 1:H_mpc])

p = ci_mpc_policy(ref_traj, s, obj,
    H_mpc = H_mpc,
    N_sample = N_sample,
    κ_mpc = κ_mpc,
	mode = :configuration,
	ip_opts = InteriorPointOptions(
					undercut = 5.0,
					κ_tol = κ_mpc,
					r_tol = 1.0e-8,
					diff_sol = true,
					solver = :empty_solver,
					max_time = 1e5,),
    n_opts = NewtonOptions(
        r_tol = 3e-4,
        max_iter = 5),
    mpc_opts = CIMPCOptions(
        altitude_update = true,
        altitude_impact_threshold = 0.05,
        )
    );

# ## Initial conditions
q1_sim, v1_sim = initial_conditions(ref_traj);

# ## Simulator
sim = simulator(s, H_sim, h=h_sim, policy=p);

# ## Simulate
simulate!(sim, q1_sim, v1_sim);
@benchmark simulate!($sim, $q1_sim, $v1_sim)

# ## Visualizer
vis = ContactImplicitMPC.Visualizer()
ContactImplicitMPC.render(vis)

# ## Visualize
plot_surface!(vis, s_sim.env, n=200, xlims = [-1, 40]);
anim = visualize_robot!(vis, model, sim.traj, sample=5);
anim = visualize_force!(vis, model, s_sim.env, sim.traj, anim=anim, h=h_sim, sample = 5);

# ## Timing result
# Julia is [JIT-ed](https://en.wikipedia.org/wiki/Just-in-time_compilation) so re-run the MPC setup through Simulate for correct timing results.
process!(sim.stats, N_sample) # Time budget
H_sim * h_sim / sum(sim.stats.policy_time) # Speed ratio


text:
type Dlog
    fd::Cint
    started_at::UInt64
    start_ofs::Csize_t
    end_ofs::Csize_t
    cur_ofs::Csize_t
end

const RanksPerLog = 16384
const SpacePerRank = 200*1024

function Dlog(destpath::String;
              ranks_per::Int = RanksPerLog,
              space_per::Int = SpacePerRank)
    @assert space_per > 100

    mkpath(destpath)

    # which log and this rank's (0 based) index for that log
    gidx, fidx = divrem(grank(), ranks_per)
    if fidx > 0
        gidx += 1
    end
    fidx -= 1
    if fidx == -1
        fidx = ranks_per - 1
    end

    # how many logs and how many ranks in the only or last log
    ngroups, nlast = divrem(ngranks(), ranks_per)
    if nlast > 0
        ngroups += 1
    end

    # how many ranks in this rank's group
    ningroup = min(ngranks(), ranks_per)
    if gidx == ngroups && nlast > 0
        ningroup = nlast
    end

    # the first rank in the group `truncate`s the log file
    fullname = joinpath(destpath, "messages-$gidx.dlog")
    if fidx == 0
        cfd = ccall(:creat, Cint, (Ptr{UInt8}, Cint), fullname, 0o644)
        systemerror("creat()", cfd == -1)
        ccall(:close, Cint, (Cint,), cfd)
        r = ccall(:truncate, Cint, (Ptr{UInt8}, Csize_t),
                  fullname, ningroup * space_per)
        systemerror("truncate()", r == -1)
    end

    # wait till the log file is created
    sync()

    # open the log file and set up this rank's offsets into it
    fd = ccall(:open, Cint, (Ptr{UInt8}, Cint, Cint), fullname,
               Base.Filesystem.JL_O_WRONLY, 0)
    systemerror("Dlog.open()", fd == -1)

    start_ofs = fidx * space_per
    end_ofs = start_ofs + space_per

    dl = Dlog(fd, time_ns(), start_ofs, end_ofs, start_ofs)
    message(dl, "[$(grank())]")

    finalizer(dl, (x -> ccall(:close, Cint, (Cint,), dl.fd)))
    return dl
end

function message(dl::Dlog, msg...)
    msg_sec = (time_ns() - dl.started_at) / 1e9
    s = @sprintf("%.3f %s\n", msg_sec, string(msg...))
    sl = length(s)
    if dl.cur_ofs + sl >= dl.end_ofs
        dl.cur_ofs = dl.start_ofs
        message(dl, "wrap")
    end
    gc_state = ccall(:jl_gc_safe_enter, Int8, ())
    wrb = ccall(:pwrite, Cint, (Cint, Ptr{UInt8}, Csize_t, Csize_t),
                dl.fd, s, sl, dl.cur_ofs)
    ccall(:jl_gc_safe_leave, Void, (Int8,), gc_state)
    dl.cur_ofs += sl
    return
end


text:
# pretty showing
Base.show(io::IO, x::CXCursor) = print(io, "$(kind(x)): $(spelling(x))")
Base.show(io::IO, x::CXType) = print(io, "$(kind(x)): $(spelling(x))")
Base.show(io::IO, x::CXToken) = print(io, "$(kind(x))")

Base.show(io::IO, x::CLToken)  = print(io, typeof(x), "(\"", x.text, "\")")
Base.show(io::IO, x::CLType)   = print(io, "CLType (", typeof(x), ") ")
Base.show(io::IO, x::CLCursor) = print(io, "CLCursor (", typeof(x), ") ", name(x))


text:
"""
Module for statistics for common factor models (or approximate factor models), by Colin T Bowers \n
\n
This module implements the following statistics for use with common factor models: \n
    - numfactor: Estimating the number of factors using the methods from Bai, Ng (2002) Determining the Number
    of Factors in Approximate Factor Models
    - testobsfactor: Testing whether an observed factor spans the latent common factor space using the methods
    from Bai, Ng (2006) Evaluating Latent and Observed Factors in Macroeconomics and Finance
\n
This package has an MIT license. Please see associated LICENSE.md file for more detail.
"""
module CommonFactorModelStats

using Distributions

include("est_num_factor.jl")
include("test_obs_fac.jl")

# package code goes here

end # module


text:
using Pkg
mktempdir() do envdir
    Pkg.activate(envdir)
    Pkg.add(PackageSpec(path=joinpath(@__DIR__, "..")))
    Pkg.status()
    include("docchecks.jl")
end


text:
"""
    locper
tipe data untuk node dari graf expanded. Terdiri dari khazanah dan periode.
"""
struct locper
    loc::String
    per::Int
end

show(io::IO, lp::locper) = print(io, "⟦i=$(lp.loc),t=$(lp.per)⟧")

# default dataframes - DON'T FORGET TO COPY
const stock_df_def = DataFrame(id=String[], pecahan=String[], value=Float64[])
const demand_df_def = insertcols!(stock_df_def |> copy, :id, :periode => [], after=true)

"""
    Libraries(khazanah, trayek, init_stock, demand_forecast, demand_realization)
defines the basic libraries needed to load the simulation. calling it with \
empty arguments give the default dataframes.
"""
@with_kw_noshow mutable struct Libraries
    khazanah::DataFrame = DataFrame(
        id=String[], name=String[], x=Float64[], y=Float64[],
        Q=Int[], cpeti=Float64[], cjarak=Float64[]
    ) # khazanah cols
    trayek::DataFrame = DataFrame(
        u=String[], v=String[], moda=String[],
        Q=Int[], cpeti=Float64[], cjarak=Float64[]
    ) # trayek cols
    init_stock::DataFrame = copy(stock_df_def)
    demand_forecast::DataFrame = copy(demand_df_def)
    demand_realization::DataFrame = copy(demand_df_def)
end

"""
    Libraries(path_to_lib; complete)
reads the files specified on `path_to_lib` and loads it into library.
if `complete = true` trayek dataframe will be turned into complete graph.
"""
function Libraries(path_to_lib::String; complete::Bool)
    to_return = Libraries()

    filenames = ["khazanah", "trayek", "demand", "stock", "moda"]
    df = Dict{String,DataFrame}()
    for f in filenames
        df[f] = CSV.read(joinpath(path_to_lib, "$f.csv"), DataFrame)

        # trayek-related adjustments
        if f == "trayek"
            if complete
                append!(df[f], DataFrame(u=df[f].v, v=df[f].u, moda=df[f].moda))
            end
            unique!(df[f])
        end
    end

    append!(to_return.trayek, innerjoin(df["trayek"], df["moda"], on=:moda => :name))
    append!(to_return.demand_forecast, df["demand"])
    append!(to_return.khazanah, df["khazanah"])
    append!(to_return.init_stock, df["stock"])

    return to_return
end

function show(io::IO, libs::Libraries)
    khazanah_stat = "Khazanah\t\t: $(nrow(libs.khazanah)) vault entry."

    trayek_stat = "Trayek\t\t\t: $(nrow(libs.trayek)) trayek entry."

    init_stock_stat = "Initial stock (t=0)\t: \
    $(length(unique(libs.init_stock.id))) vault with \
    $(length(unique(libs.init_stock.pecahan))) pecahan."

    demand_forecast_stat = "Demand Forecast\t\t: \
    $(length(unique(libs.demand_forecast.id))) vault \
    for $(length(unique(libs.demand_forecast.periode))) periode \
        with $(length(unique(libs.demand_forecast.pecahan))) pecahan."

    demand_realization_stat = "Demand Realization\t: \
    $(length(unique(libs.demand_realization.id))) vault \
    for $(length(unique(libs.demand_realization.periode))) periode \
        with $(length(unique(libs.demand_realization.pecahan))) pecahan."

    print(io,
        "---System Libraries---\n$(khazanah_stat) \n$(trayek_stat) \n$(init_stock_stat) \n$(demand_forecast_stat) \n$(demand_realization_stat)"
    )
end

"""
    States(current_stock, dispatch_queue)
defines the states that makes up the system.
"""
@with_kw_noshow mutable struct States
    current_stock::DataFrame = copy(stock_df_def)
    dispatch_queue::MetaDigraph{locper} = MetaDigraph{locper}()
end

function show(io::IO, states::States)
    current_stock_stat = "Current Stock\t: \
    $(length(unique(states.current_stock.id))) vault \
    with $(length(unique(states.current_stock.pecahan))) pecahan."

    dispatch_queue_stat = "Dispatch Queue\t: \
    $(length(arcs(states.dispatch_queue))) trayek entry."

    print(io,
        "---System States---\n$(current_stock_stat)\n$(dispatch_queue_stat)"
    )
end

isempty(s::States) = reduce(*,
    [isempty(getfield(s, f)) for f in fieldnames(typeof(s))]
)

@with_kw_noshow mutable struct Accumulators
    executed_dispatch::MetaDigraph{locper} = MetaDigraph{locper}()
    inventory_levels::DataFrame = copy(demand_df_def)
    demand_fulfillment::DataFrame = copy(demand_df_def)
end

function show(io::IO, accumulators::Accumulators)
    executed_dispatch_stat = "Executed Dispatch \t: \
    $(accumulators.executed_dispatch.core.na) trayek entry."

    inventory_levels_stat = "Inventory Levels \t: \
    $(length(unique(accumulators.inventory_levels.id))) vault \
    for $(length(unique(accumulators.inventory_levels.periode))) periode \
        with $(length(unique(accumulators.inventory_levels.pecahan))) pecahan."

    demand_fulfillment_stat = "Demand Fulfillment \t: \
    $(length(unique(accumulators.demand_fulfillment.id))) vault \
    for $(length(unique(accumulators.demand_fulfillment.periode))) periode \
        with $(length(unique(accumulators.demand_fulfillment.pecahan))) pecahan."

    print(io,
        "---System Accumulators---\n$(executed_dispatch_stat)\n$(inventory_levels_stat)\n$(demand_fulfillment_stat)"
    )
end

isempty(a::Accumulators) = reduce(*,
    [isempty(getfield(a, f)) for f in fieldnames(typeof(a))]
)

@with_kw_noshow mutable struct Params
    T::Int # terminating timstep
    H::Int # planning horizon
    GAP::Float64 = 0.2 # MIPGap
    model::Function # model generator
end

function show(io::IO, params::Params)
    print(io,
        "---System Params---\nT\t= $(params.T) \nH\t= $(params.H) \nGAP\t= $(params.GAP)\
        \nmodel\t= $(params.model)"
    )
end

"""
    Simulation
initiating `Simulation` requires two inputs:
- libs
- params

at initiation, `current_timestep` will be set to zero.
init_stocks will be appended from lib to states and accumulators.
"""
@with_kw_noshow mutable struct Simulation
    # SCHEDULER
    t::Int = 0
    queue::Vector{Function} = Vector{Function}()

    # MODIFIERS
    fixed::Bool = false

    # GENERATED
    stt::States = States()
    acc::Accumulators = Accumulators()

    # SEED
    libs::Libraries
    params::Params
end

function show(io::IO, sim::Simulation)
    print(io,
        "\nCURRENT TIMESTEP = $(sim.t)\n\
        \n$(sim.stt)\n\
        \n$(sim.acc)\n\
        \n$(sim.params)\n"
    )
end