text:
########################################################
####                sfmodel_spec()                  ####
########################################################

"""
sfmodel_spec(<keyword arguments>)

Provide specifications of the stochastic frontier model, including the type of model
and names of variables or matrix used in estimating the model. Two ways to
specify: Method 1: Use DataFrame as input, and Method 2: use matrix as input.

# Method 1 (DataFrame input)
Variables come from a DataFrame, and column names of the Dataframe are used in
the variable input. With this method,
equations are identified by macros but not functions (e.g., `@depvar()` but
not `depvar()`).

## Arguments of Method 1
- `sfdist(::Vararg)`: the distribution assumption of the one-sided stochastic
  variable (aka inefficiency term) of the model;
  possible choices include `truncated` (or `trun`, `t`), `half` (or `h`),
  `exponential` (or `expo`, `e`), and `trun_scale` (or `trun_scaling`, `ts`).
- `sftype(::Vararg)`: whether the model is a `production` (or `prod`) frontier
  or a `cost` frontier.
- `sfpanel(::Vararg)`: the type of panel model. Choices include `TFE_WH2010`
  (true fixed effect model of Wang and Ho 2010 JE), `TFE_CSW2014` (true fixed
  model of Chen, Schmidt, and Wang 2014 JE),  `TRE` (true random effect model
  of Greene 2004), `TimeDecay` (time decay model of Battese and Coelli 1992).
- `@depvar(::Vararg)`: the dependent variable from a DataFrame.
- `@frontier(::Vararg)`: a list of variables, separated by commas, in the frontier function.
- `@μ(::Vararg)` or `@mu(::Vararg)`: a list of variable, separated by comma,
  in the linear function of μ. (`sftype(trun)` only).
- `@σᵥ²(::Vararg)` or `@sigma_v_2(::Vararg)`: a list of variable, separated by comma, in the σᵥ²
  equation.
- `@σᵤ²(::Vararg)` or `@sigma_u_2(::Vararg)`: a list of variable, separated by comma, in the σᵤ²
  equation.
- `@σₐ²(::Vararg)` or `@sigma_a_2(::Vararg)`: a list of variable, separated by comma, in the σₐ²
  equation. `sfpanel(TRE)` only.
- `@gamma(::Vararg)`: a list of variables, separated by commas, in the gamma
  equation. `sfpanel(TimeDecay)` only.
- `@timevar(::Vararg)`: the variable containing the time period information.
  Panel data model only.
- `@idvar(::Vararg)`: the variable identifying each individual. Panel data
  model only.
- message::Bool: Whether printing (=true) or not (=false, the default) the
  confirmation message "A dictionary from sfmodel_spec() is generated."
  on the screen after `sfmodel_spec()` is successfully executed.

# Method 2 (matrix/vector input)
Data of the variables are
provided by individual matrices or vectors, and names of the mat/vec are used
in the equations. With this method, equations are identified by functions but
not macros (e.g., `depvar()` but not `@depvar()`). Note that if, for instance,
the name of `depvar` or `σᵤ²` has been used elsewhere in the program, using
these names to read in mat/vec will cause name conflict (`MethodError: objects
of type ... are not callable`). The workaround is to fully qualify the function
names, e.g., `SFrontiers.depvar`, `SFrontiers.σᵤ²`, etc. Or, use the alias (if
available), e.g., `sigma_u_2` instead of `σᵤ²`.

## Arguments of Method 2
- `sfdist(::Vararg)`: the distribution assumption on the inefficiency term;
  possible choices include `truncated` (or `trun`, `t`), `half` (or `h`),
  `exponential` (or `expo`, `e`), and `trun_scale` (or `trun_scaling`, `ts`).
- `sftype(::Vararg)`: whether the model is a `production` (or `prod`) frontier
  or a `cost` frontier.
- `sfpanel(::Vararg)`: the type of panel model. Choices include `TFE_WH2010`
  (true fixed effect model of Wang and Ho 2010 JE), `TFE_CSW2014` (true fixed
  model of Chen, Schmidt, and Wang 2014 JE),  `TRE` (true random effect model
  of Greene 2004), `TimeDecay` (time decay model of Battese and Coelli 1992).
- `depvar(::Matrix)`: Matrix or vector of the dependent variable.
- `frontier(::Matrix)`: matrix or vector for frontier function.
- `μ(::Matrix)` or `mu(::Matrix)`: matrix or vector for the (linear) μ equation (`trun` type only).
- `σᵤ²(::Matrix)` or `sigma_u_2(::Matrix)`: matrix or vector for the σᵤ² equation.
- `σₐ²(::Matrix)` or `sigma_a_2(::Matrix)`: matrix or vector for the σₐ² equation.
- message::Bool: Whether printing (=true) or not (=false, the default) the
  confirmation message "A dictionary from sfmodel_spec() is generated."
  on the screen after `sfmodel_spec()` is successfully executed.

# Examples
```julia-repl
sfmodel_spec(sftype(prod), sfdist(trun),
             @depvar(output),
             @frontier(land, , labor, bull, year, _cons),
             @μ(age, school, year, _cons),
             @σᵤ²(age, school, year, _cons),
             @σᵥ²(_cons));

sfmodel_spec(sfpanel(TRE), sftype(prod), sfdist(half),
             @timevar(yr), @idvar(id),
             @depvar(y),
             @frontier(x1, x2, _cons),
             @σₐ²(_cons),
             @σᵤ²(_cons),
             @σᵥ²(_cons),
             message = false);
```
"""
function sfmodel_spec(arg::Vararg; message::Bool=false)

    global _dicM
           _dicM = Dict{Symbol, Any}()  # nullify and initiate new dictionaries when a new model is specified

         #* -- creates default values ---

        for k in (:panel, :timevar, :idvar, :dist, :type, :depvar, :frontier, :μ, :hscale, :gamma, :σᵤ², :σᵥ², :σₐ², :hasDF, :transfer, :misc) # take out :data, :η, :λ, :τ, in this revision
            _dicM[k] = nothing
        end

        #* -- replace the defaults with user's values ---

        for d in :($(arg))
            _dicM[d[1]] = d[2]
        end

        #* ==== Method 2, matrix input (such as in simulations), create a DataFrame

           _dicM[:hasDF]    = true
           _dicM[:transfer] = false

        if typeof(_dicM[:depvar]) != Array{Symbol,1} # not a DataFrame

           _dicM[:hasDF] = false

           isa(_dicM[:depvar][1], Vector) || isa(_dicM[:depvar][1], Matrix) || throw("
           `depvar()` has to be a Vector or Matrix (e.g., Array{Float64, 1} or Array{Float64, 2}).
           Check with `isa(your_thing, Matrix)` or `isa(your_thing, Vector)`.
           Try `convert()`, `reshape()`, `Matrix()`, or something similar.")

           comDF = _dicM[:depvar][1]  # create the first data column of comDF
           varname = [:depvar]
           _dicM[:depvar] = [:depvar]

           for k in (:timevar, :idvar, :frontier, :μ, :hscale, :gamma, :σₐ², :σᵤ², :σᵥ²)
               if _dicM[k] !== nothing # if not nothing, must be Array
                  isa(_dicM[k], Vector) || isa(_dicM[k][1], Vector) || isa(_dicM[k][1], Matrix) || throw("
                     `k` has to be a Vector or Matrix (e.g., Array{Float64, 1} or Array{Float64, 2}).
                     Check with `isa(your_thing, Matrix)` or `isa(your_thing, Vector)`.
                     To convert, try `convert()`, `reshape()`, `Matrix()`, or something similar.")

                  (isa(_dicM[k], Vector)  && length(_dicM[k][1]) == 1) ?  _dicM[k] = [_dicM[k]] : nothing # ugly fix for pure vector input

                  @views comDF = hcat(comDF, _dicM[k][1]) # combine the data
                  aa = Symbol[]
                  for i in 1:size(_dicM[k][1], 2)
                      push!(aa, Symbol(String(k)*"_var$(i)")) # create name for the data
                  end
                  varname = vcat(varname, aa) # combine the dataname
                  _dicM[k] = aa
               end
           end # for k in (...)

           comDF = DataFrame(comDF, varname)
           _dicM[:sdf] = comDF

        end # if typeof(...)

        #* -- check the model identifier and the model type ---

        (_dicM[:dist] !== nothing) || throw("You need to specify dist().")
        (_dicM[:type] !== nothing) || throw("You need to specify type().")

        #* --- get the model identifier -------

        s = uppercase(String(_dicM[:dist][1])[1:1])

        global tagD
        if _dicM[:panel] === nothing  # not panel
            if (s == "T")
                tagD = Dict{Symbol, Type{Trun}}()
                tagD[:modelid] = Trun
            elseif (s == "H")
                tagD = Dict{Symbol, Type{Half}}()
                tagD[:modelid] = Half
            elseif (s == "E")
                tagD = Dict{Symbol, Type{Expo}}()
                tagD[:modelid] = Expo
            elseif (s == "S")
                tagD = Dict{Symbol, Type{Trun_Scale}}()
                tagD[:modelid] = Trun_Scale
            end
        elseif (_dicM[:panel] == [:TFE_WH2010]) && (s=="T")
            tagD = Dict{Symbol, Type{PFEWHT}}()
            tagD[:modelid] = PFEWHT
        elseif (_dicM[:panel] == [:TFE_WH2010]) && (s == "H")  # panel and dist=half
            tagD = Dict{Symbol, Type{PFEWHH}}()
            tagD[:modelid] = PFEWHH
        elseif (_dicM[:panel] == [:TFE_WH2010])
            throw("The panel TFE_WH2010 model can only have `sfdist(trun)` or `sfdist(half)`.")
        elseif (_dicM[:panel] == [:TFE_CSW2014]) && (s=="H")
            tagD = Dict{Symbol, Type{PFECSWH}}()
            tagD[:modelid] = PFECSWH
        elseif (_dicM[:panel] == [:TFE_CSW2014])
            throw("The panel TFE_CSW2014 model can only have `sfdist(half)`.")
        elseif (_dicM[:panel] == [:TRE]) && (s=="H")
            tagD = Dict{Symbol, Type{PTREH}}()
            tagD[:modelid] = PTREH
        elseif (_dicM[:panel] == [:TRE]) && (s=="T")
            tagD = Dict{Symbol, Type{PTRET}}()
            tagD[:modelid] = PTRET
        elseif (_dicM[:panel] == [:TRE])
            throw("The panel TRE model can only have `sfdist(half)` or `sfdist(trun)`.")
        elseif (_dicM[:panel] == [:TimeDecay]) && (s == "T")
            tagD = Dict{Symbol, Type{PanDecay}}()
            tagD[:modelid] = PanDecay
        elseif (_dicM[:panel] == [:TimeDecay])
            throw("The panel time-decay model can only have `sfdist(trun)`.")
        elseif (_dicM[:panel] == [:Kumbhakar1990]) && (s == "T")
            tagD = Dict{Symbol, Type{PanKumb90}}()
            tagD[:modelid] = PanKumb90
        elseif (_dicM[:panel] == [:Kumbhakar1990])
            throw("The panel Kumbhakar 1990 model can only have `sfdist(trun)`.")
        else
            throw("The `sfpanel()` and/or `sfdist()` are not specified correctly.")
        end

        #* ---- check if the model has the correct syntax ---

        SFrontiers.checksyn(tagD[:modelid])

        #* ----- make return -----------
        if message
          printstyled("A dictionary from sfmodel_spec() is generated.\n"; color = :green)
        end
        return _dicM # for debugging purpose

end  # end of sfmodel_spec()

########################################################
####                sfmodel_init()                  ####
########################################################
"""
    sfmodel_init(<keyword arguments>)

Provide initial values for the stochastic frontier model estimation. The
values could be a vector or scalars. It creates a global dictionary `_dicINI`. Optional.

# Arguments
- `all_init(::Union{Vector, Real})`: initial values of all the parameters in the model
- `frontier(::Union{Vector, Real})`: initial values of parameters in
  the `frontier()` function
- `μ(::Union{Vector, Real})` or `mu(::Union{Vector, Real})`: initial values of
  parameters in the `μ` function
- `hscale(::Union{Vector, Real})`: initial values of parameters in the `hscale()` function
- `gamma(::Union{Vector, Real})`: initial values of parameters in the `gamma()` function
- `σᵤ²(::Union{Vector, Real})` or `sigma_u_2(::Union{Vector, Real})`: initial values of parameters in the
   `σᵤ²` function
- `σᵥ²(::Union{Vector, Real})` or `sigma_v_2(::Union{Vector, Real})`: initial values of parameters in the
   `σᵥ²` function
- `σₐ²(::Union{Vector, Real})` or `sigma_a_2(::Union{Vector, Real})`: initial
  values of parameters in the `σₐ²` function
- message::Bool: Whether printing (=true) or not (=false, the default) the
  confirmation message "A dictionary from sfmodel_init() is generated."
  on the screen after `sfmodel_init()` is successfully executed.

# Remarks
- Equations do not have to follow specific orders.
- `sfmodel_init(...)` is optional but is highly recommended. If it is not
  specified or is specified as an empty set, default values are used.
- It is not necessary to specify a complete set of equations. A partial list
  or even empty lists are acceptable. Default values will be substituted for the
  missing equations.
- The generated `_dicINI` is inheritable in the sense that an exiting
  `_dicINI` (from the previous run of the same or a different model, for
  example) will be used if the current model does not have its own
  `sfmodel_init(...)`. This design has advantages in a simulations study where
  `sfmodel_init(...)` needs to be specified only once.

# Examples
```julia-repl
b_ini = ones(2)*0.2
sfmodel_init( # frontier(bb),             # may skip and use default
             μ(b_ini),                    # may use a vector
             σᵤ²(-0.1, -0.1),
             σᵥ²(-0.1) )

sfmodel_init(all_init(0.1, 0.2, 0.5, 0.0, -0.1, -0.1, -0.1),
             message = false)
```
"""
function sfmodel_init(arg::Vararg; message::Bool =false) # create a dictionary of inital vectors

   global _dicINI
          _dicINI = Dict{Symbol, Any}()

    for d in :($(arg))
        _dicINI[d[1]] = d[2]
    end

    #* If has the key, creates the alias key with the same value.

    !(haskey(_dicINI, :μ))      || (_dicINI[:eqz] = _dicINI[:μ])
    !(haskey(_dicINI, :hscale)) || (_dicINI[:eqq] = _dicINI[:hscale])
    !(haskey(_dicINI, :gamma))  || (_dicINI[:eqq] = _dicINI[:gamma])
    !(haskey(_dicINI, :σₐ²))    || (_dicINI[:eqw] = _dicINI[:σₐ²])
    !(haskey(_dicINI, :σᵤ²))    || (_dicINI[:eqw] = _dicINI[:σᵤ²])
    !(haskey(_dicINI, :σᵥ²))    || (_dicINI[:eqv] = _dicINI[:σᵥ²])

    if message
      printstyled("A dictionary from sfmodel_init() is generated.\n"; color = :green)
    end
    return _dicINI # for debugging purpose
end

########################################################
####                sfmodel_opt()                   ####
########################################################
"""
    sfmodel_opt(<keyword arguments>)

Provide options to the optimization algorithms for the maiximum likelihood
estimation. It creates a global dictionary `_dicOPT`. Optional. The `Optim`
package is used for the optimization, and a subset of
`Optim`'s keywords are directly accessible from this API.

# Arguments
- `warmstart_solver(algorithm)`: The algorithm used in the first-stage ("warmstart")
  optimization process, which serves the purpose of improving upon the initial
  values for the second-stage ("main") estimation. The default is
  `NelderMead()`. Others include `SimulatedAnnealing()`, `SAMIN()`, `ParticleSwarm()`,
  `ConjugateGradient()`, `GradientDescent()`, `BFGS()`, `LBFGS()`,
  `Newton()`, `NewtonTrustRegion()`, and `IPNewton()`. See
  http://julianlsolvers.github.io/Optim.jl/stable/ for details.
  Non-gradient based algorithms are recommended for the warmstart solver.
- `warmstart_maxIT(::Int64)`: The iteration limit for the warmstart. Default
  is 100.
- `main_solver(algorithm)`: The algorithm used in the main opimization process.
  The default is `Newton()`. Others include `SimulatedAnnealing()`, `SAMIN()`, `ParticleSwarm()`,
  `ConjugateGradient()`, `GradientDescent()`, `BFGS()`, `LBFGS()`,
  `NewtonTrustRegion()`, and `IPNewton()`. See
  http://julianlsolvers.github.io/Optim.jl/stable/ for details.
- `main_maxIT(::Int64)`: The iteration limit for the main estimation. Default
  is 2000.
- `tolerance(::Float64)`: The convergence criterion ("tolerance") based on the
  absolute value of gradients. Default is 1.0e-8. For non-gradient algorithms,
  it controls the main convergence tolerance, which is solver specific.
  See `Optim`'s `g_tol` option for more information.
- `verbose(::Bool)`: Print on screen (`true`, the default) the information of
  the model and the optimization results.
- `banner(::Bool)`: Print on screen (`true`, the default) a banner to serve as
  a visual indicator of the start of the estimation.
- `ineff_index(::Bool)`: Whether to compute the Jondrow et al. (1982)
  inefficiency index and the Battese and Coelli (1988) efficiency index. The
  defauis `true`.
- `marginal(::Bool)`: Whether to compute the marginal effects of the exogenous
  determinants of inefficiency (if any).
- `table_format()`: The format to print the coefficient tables on the screen:
  `text` (default), `html`, or `latex`. A wrapper of `PrettyTables.jl`'s
  `backend` option.
- message::Bool: Whether printing (=true) or not (=false, the default) the
  confirmation message "A dictionary from sfmodel_opt() is generated."
  on the screen after `sfmodel_opt()` is successfully executed.

# Remarks
- `sfmodel_opt(...)` is optional. It can be omitted entirely, or specifying
  only a partial list of the keywords.
- If any of the keywords are missing, default values are used.
- If warmstart is not needed, you need to give empty keyword values to
  warmstart related keys. E.g., either `warmstart_solver()` or
  `warmstart_maxIT()`, or both. Omitting the keyword entirely (i.e., not
  writing down `warmstart_solver` or `warmstart_maxIT`) will not skip the
  warmstart, but will reinstate the default.
- Users do not need to provide gradient or Hessian functions even if
  gradient-based optimization algorithms are used. The package uses automatic
  differentiation (https://en.wikipedia.org/wiki/Automatic_differentiation) to
  compute the derivatives. It is not numerical finite differentiation. It is
  fast and as accurate as the symbolic differentiation.
- The `_dicOPT` is inheritable in the sense that an exiting `_dicOPT` (from
  the previous run of the same or a different model, for example) will be used
  if the current model does not have its own `sfmodel_opt(...)`. This design
  has advantages in simulation studies where `sfmodel_opt(...)` needs to be
  specified only once.

# Examples
```julia-repl
sfmodel_opt(warmstart_solver(NelderMead()),
            warmstart_maxIT(200),
            main_solver(Newton()),
            main_maxIT(2000),
            tolerance(1e-8),
            message = false)
```
"""
function sfmodel_opt(arg::Vararg; message::Bool=false) # create a dictionary of maximization options

    global _dicOPT
           _dicOPT = Dict{Symbol, Any}()

    #* -- creates the default ---

    _dicOPT[:warmstart_solver] = :(NelderMead())
    _dicOPT[:warmstart_maxIT]  =  100
    _dicOPT[:main_solver]      = :(Newton())
    _dicOPT[:main_maxIT]       =  2000
    _dicOPT[:tolerance]        =  1.0e-8
    _dicOPT[:verbose]          =  true
    _dicOPT[:banner]           =  true
    _dicOPT[:ineff_index]      =  true
    _dicOPT[:marginal]         =  true
    _dicOPT[:table_format]     = :(text)

    #* -- replace the defaults with the user's value ---

    for d in :($(arg))
        _dicOPT[d[1]] = d[2]
    end

    #* ---- error checking --

    if (_dicOPT[:main_solver] === nothing) || (_dicOPT[:main_maxIT] === nothing) || (_dicOPT[:tolerance] === nothing)
         throw("You cannot give empty keyword values to `main_solver()`, `main_maxIT()`, or `tolerance()`. If you want to use the default, you may do so by dropping (not emptying keyword values) the keywords.")
    end

    if message
      printstyled("A dictionary from sfmodel_opt() is generated.\n"; color = :green)
    end
    return _dicOPT # for debugging purpose

end    # end of sfmodel_opt()

########################################################
###                 sfmodel_fit()                   ####
########################################################
"""
    sfmodel_fit(<keyword arguments>)

Maximum likelihood estimation of the stochastic frontier model specified
in `sfmodel_spec(...)`. Estimate the model parameters, calculate Jondrow et al.
(1982) inefficiency index and Battese and Coelli (1988) efficiency index,
compute marginal effects of inefficiency determinants (if any).
Return a dictionary with results.

# Arguments
- `useData(::DataFrame)`: The DataFrame used with the Method 1 of
  `sfmodel_spec(...)`. If use Method 2 of `sfmodel_spec(...)` (viz., data
  is supplied by individual matrices), do not need this keyword argument.

# Remarks
- Use `Optim.jl` to carry out the estimation.
- Users do not need to provide gradient or Hessian functions even if
  gradient-based optimization algorithms are used. The package uses automatic
  differentiation (https://en.wikipedia.org/wiki/Automatic_differentiation) to
  compute the derivatives. AD is not numerical finite differentiation. AD is
  fast and as accurate as the symbolic differentiation.

# Examples
```julia-repl
sfmodel_fit(useData(df))    # Method 1
sfmodel_fit()               # Method 2
```
"""
function sfmodel_fit()
     # For Method 2 of `sfmodel_spec()`.

    !(_dicM[:hasDF]) || throw("Need to specify DataFrame in `sfmodel_fit()`.")

    _dicM[:transfer] = true
    sfmodel_fit(_dicM[:sdf])

end

function sfmodel_fit(sfdat::DataFrame) #, D1::Dict = _dicM, D2::Dict = _dicINI, D3::Dict = _dicOPT)

    (_dicM[:hasDF] || _dicM[:transfer])  || throw("You provided matrix in `sfmodel_spec()` so you cannot specify a DataFrame in `sfmodel_fit()`. Leave it blank.")

   #* for simulation, add a flag
   redflag::Bool = 0

  #* ###### Check if the OPT dictionary exists #####

      @isdefined(_dicINI) || sfmodel_init()  # if not exist, create one with default values
      @isdefined(_dicOPT) || sfmodel_opt()

     if _dicOPT[:banner]
        printstyled("\n###------------------------------------###\n"; color=:yellow)
        printstyled("###  Estimating SF models using Julia  ###\n"; color=:yellow)
        printstyled("###------------------------------------###\n\n"; color=:yellow)
      end

  #* ##### Get variables from dataset #######

     # pos: (begx, endx, begz, endz, ...); variables' positions in the parameter vector.
     # num: (nofobs, nofx, ..., nofpara); number of variables in each equation
     # eqvec: ("frontier"=2, "μ"=6,...); named tuple of equation names and equation position in the table
     # eqvec2: (xeq=(1,3), zeq=(4,5),...); named tuple of equation and parameter positions, for sfmodel_predict
     # varlist: ("x1", "x2",...); variable names for making table

     (minfo1, minfo2, pos, num, eqvec, eqvec2, yvar, xvar, zvar, qvar, wvar,
      vvar,         rowIDT, varlist) = SFrontiers.getvar(tagD[:modelid], sfdat)

  #* ### print preliminary information ########

    if _dicOPT[:verbose]

      printstyled("*********************************\n "; color=:cyan)
      printstyled("      Model Specification:\n"; color=:cyan);
      printstyled("*********************************\n"; color=:cyan)

      print("Model type: "); printstyled(minfo1; color=:yellow); println();println()
      printstyled(minfo2; color=:yellow); println()
    end

  #* ##### Get the type parameter #######

     _porc::Int64 = 1

     if (_dicM[:type] == [:cost])
         _porc = -1
     end

  #* ########## Process initial value dictionary  #####
     #* --- Get OLS results and other auxiliary values. --- #

     noffixed = 0  # number of fixed effect parameter
     if _dicM[:panel] == [:TFE_WH2010] || _dicM[:panel] == [:TFE_CSW2014]  # the fixed effects and differenced off
        noffixed = size(rowIDT,1)
     end

     β0     = xvar \ yvar;  # OLS estiamte, uses a pivoted QR factorization;
     resid  = yvar - xvar*β0
     sse    = sum((resid).^2)
     ssd    = sqrt(sse/(size(resid,1)-(num.nofx + noffixed ))) # sample standard deviation; σ² = (1/(N-K))* Σ ϵ^2
     ll_ols = sum(normlogpdf.(0, ssd, resid)) # ols log-likelihood
     sk_ols = sum((resid).^3) / ((ssd^3)*(size(resid,1))) # skewnewss of ols residuals

     #* --- Create the dictionary -----------

     if (:all_init in keys(_dicINI))
         sf_init = _dicINI[:all_init]
     else
         #*  Create ini vectors from user's values; if none, use the default.--- #
         b_ini  = get(_dicINI, :frontier, β0)
         d1_ini = get(_dicINI, :eqz, ones(num.nofz) * 0.1)
         t_ini  = get(_dicINI, :eqq, ones(num.nofq) * 0.1)
         d2_ini = get(_dicINI, :eqw, ones(num.nofw) * 0.1)
         g_ini  = get(_dicINI, :eqv, ones(num.nofv) * 0.1)

         #*  Make it Array{Float64,1}; otherwise Array{Float64,2}. ---#
         #*       Could also use sf_init[:,1]. *#
         sf_init = vcat(b_ini, d1_ini, t_ini, d2_ini, g_ini)
         sf_init = vec(sf_init)
     end # if :all_init

  #* ############ Misc.  ################
     # --- check if the number of initial values is correct
        (length(sf_init) == num.nofpara) ||  throw("The number of initial values does not match the number of parameters to be estimated. Make sure the number of init values in sfmodel_init() matches the number of variabls in sfmodel_spec().")

     # --- Make sure there is no numerical issue arising from int vs. Float64.
        sf_init = convert(Array{Float64,1}, sf_init)

  #* ############# process optimization dictionary  #######

         if (_dicOPT[:warmstart_solver] === nothing) || (_dicOPT[:warmstart_maxIT] === nothing)
             do_warmstart_search = 0
         else
             do_warmstart_search = 1
             sf_ini_algo  = eval(_dicOPT[:warmstart_solver])  # warmstart search algorithms
             sf_ini_maxit = _dicOPT[:warmstart_maxIT]         # warmstart search iter limit
         end

     # ---- main maximization algorithm -----
         sf_algo  = eval(_dicOPT[:main_solver])    # main algorithm
         sf_maxit = _dicOPT[:main_maxIT]
         sf_tol   = _dicOPT[:tolerance]
         sf_table = _dicOPT[:table_format]

  #* ########  Start the Estimation  ##########

    #* ----- Define the problem's Hessian -----#

     _Hessian = TwiceDifferentiable(rho -> SFrontiers.LL_T(tagD[:modelid],
                           yvar, xvar, zvar, qvar, wvar, vvar,
                           _porc, num.nofobs, pos, rho,
                                   rowIDT, _dicM[:misc]),
                     sf_init;
                    autodiff = :forward);

    #* ---- Make placeholders for dictionary recording purposes *#

    sf_init_1st_dic  = 0
    sf_init_2nd_dic  = 0
    sf_ini_algo_dic  = nothing
    sf_ini_maxit_dic = 0
    sf_total_iter    = 0

    _run = 1  # a counter; use the -if- instead of -for- to avoid using global variables

    if (do_warmstart_search == 1) && (_run == 1)

        if _dicOPT[:verbose]
            printstyled("The warmstart run...\n\n"; color = :green)
        end

        sf_init_1st_dic  = copy(sf_init) # for dict recording
        sf_ini_algo_dic  = sf_ini_algo
        sf_ini_maxit_dic = copy(sf_ini_maxit)

        # @time
               mfun = optimize(_Hessian,
                               sf_init,         # initial values
                               sf_ini_algo,
                               Optim.Options(g_tol = sf_tol,
                                             iterations  = sf_ini_maxit,
                                             store_trace = true,
                                             show_trace  = false))

        sf_total_iter += Optim.iterations(mfun) # for later use

        sf_init = Optim.minimizer(mfun)  # save as initials for the next run
        _run    = 2                      # modify the flag

        if _dicOPT[:verbose]
            println()
            print("$mfun \n")
            print("The warmstart results are:\n"); printstyled(Optim.minimizer(mfun); color=:yellow); println("\n")
        end

   end  # if  (do_warmstart_search == 1) && (_run == 1)

   if (do_warmstart_search == 0 ) || (_run == 2) # either no warmstart run and go straight here, or the 2nd run

       sf_init_2nd_dic = copy(sf_init) # for dict recording

       if _dicOPT[:verbose]
           println()
           printstyled("Starting the optimization run...\n\n" ; color = :green)
       end

       # @time
              mfun = optimize(_Hessian,
                              sf_init,       # initial values
                              sf_algo,       # different from search run
                              Optim.Options(g_tol = sf_tol,
                                            iterations  = sf_maxit, # different from search run
                                            store_trace = true,
                                            show_trace  = false))
       sf_total_iter += Optim.iterations(mfun)

       if _dicOPT[:verbose]
             println()
             print("$mfun \n")
             print("The resulting coefficient vector is:\n"); printstyled(Optim.minimizer(mfun); color=:yellow); println("\n")
       end

       if isnan(Optim.g_residual(mfun)) || (Optim.g_residual(mfun) > 0.1)
            redflag = 1
            printstyled("Note that the estimation may not have converged properly. The gradients are problematic (too large, > 0.1, or others).\n\n", color = :red)
       end

       if Optim.iteration_limit_reached(mfun)
             redflag = 1
             printstyled("Caution: The number of iterations reached the limit.\n\n"; color= :red)
       end

   end     # if (do_warmstart_search == 0 )....

  #* ###### Post-estimation process ###############

      _coevec            = Optim.minimizer(mfun)  # coef. vec.
      numerical_hessian  = hessian!(_Hessian, _coevec)  # Hessain

     #* ------ Check if the matrix is invertible. ----

     var_cov_matrix = try
                         inv(numerical_hessian)
                      catch err
                         redflag = 1
                         checkCollinear(tagD[:modelid], xvar, zvar, qvar, wvar, vvar) # check if it is b/c of multi-collinearity in the data
                         throw("The Hessian matrix is not invertible, indicating the model does not converge properly. The estimation is abort.")
                      end

          #* In some cases the matrix is invertible but the resulting diagonal
          #*    elements are negative. Check.

          if !all( diag(var_cov_matrix) .> 0 ) # not all are positive
               redflag = 1
               printstyled("Some of the diagonal elements of the var-cov matrix are non-positive, indicating problems in the convergence. The estimation is abort.\n\n"; color = :red)
               checkCollinear(tagD[:modelid], xvar, zvar, qvar, wvar, vvar) # check if it is b/c of multi-collinearity in the data
          end

     #* ------- JLMS and BC index -------------------

     if _dicOPT[:ineff_index]
        @views (_jlms, _bc) = jlmsbc(tagD[:modelid], _porc, pos, _coevec,
                                     yvar, xvar, zvar, qvar, wvar, vvar,         rowIDT)
        _jlmsM = mean(_jlms) # sum(_jlms)/length(_jlms)
        _bcM   = mean(_bc)   # sum(_bc)/length(_bc)
     else
        _jlms  = nothing
        _bc    = nothing
        _jlmsM = nothing
        _bcM   = nothing
     end

     #* ---- marginal effect on E(u) --------------

     if _dicOPT[:marginal]
        margeff, margMinfo = get_marg(tagD[:modelid], pos, num, _coevec, zvar, qvar, wvar)
     else
        margeff, margMinfo = nothing, ()
     end

     #* ------- Make Table ------------------

     stddev  = sqrt.(diag(var_cov_matrix)) # standard error
     t_stats = _coevec ./ stddev          # t statistics
     p_value = zeros(num.nofpara)   # p values
     ci_low  = zeros(num.nofpara) # confidence interval
     ci_upp  = zeros(num.nofpara)
     tt      = cquantile(Normal(0,1), 0.025)

     for i = 1:num.nofpara
         @views p_value[i,1] = pvalue(TDist(num.nofobs - num.nofpara), t_stats[i,1]; tail=:both)
         @views ci_low[i,1] = _coevec[i,1] - tt*stddev[i,1]
         @views ci_upp[i,1] = _coevec[i,1] + tt*stddev[i,1]
     end

       #* Build the table columns *#

    table = zeros(num.nofpara, 7)  # 7 columns in the table
    table[:,2] = _coevec   # estiamted coefficients
    table[:,3] = stddev    # std deviation
    table[:,4] = t_stats   # t statistic
    table[:,5] = p_value   # p value
    table[:,6] = ci_low
    table[:,7] = ci_upp
    table      = [" " "Coef." "Std. Err." "z" "P>|z|" "95%CI_l" "95%CI_u"; table]  # add to top of the table

       #*  creating a column of function names

    table[:, 1] .= ""
    for i in 1:length(eqvec)
        @views j = eqvec[i]
        @views table[j,1] = keys(eqvec)[i]
    end

       #*  Add the column of variable names

    table = hcat(varlist, table)                      # combine the variable names column (hcat, horizontal concatenate; see also vcat)
    table[:,1], table[:,2] = table[:,2], table[:,1]   # swap the first name column and the function column
    table[1,2] = "Var."

     # * ------ Print Results ----------- *#

     if _dicOPT[:verbose]

         printstyled("*********************************\n "; color=:cyan)
         printstyled("      Estimation Results:\n"; color=:cyan);
         printstyled("*********************************\n"; color=:cyan)

         print("Model type: "); printstyled(minfo1; color=:yellow); println()
         print("Number of observations: "); printstyled(num.nofobs; color=:yellow); println()
         print("Number of total iterations: "); printstyled(sf_total_iter; color=:yellow); println()
         if Optim.converged(mfun)
             print("Converged successfully: "); printstyled(Optim.converged(mfun); color=:yellow); println()
         elseif Optim.converged(mfun) == false
             print("Converged successfully: "); printstyled(Optim.converged(mfun); color=:red); println()
             redflag = 1
         end
         print("Log-likelihood value: "); printstyled(round(-1*Optim.minimum(mfun); digits=5); color=:yellow); println()
         println()

         pretty_table(table[2:end,:],    # could print the whole table as is, but this prettier
                      header=["", "Var.", "Coef.", "Std.Err.", "z", "P>|z|",
                              "95%CI_l", "95%CI_u"],
                      formatters = ft_printf("%5.4f", 3:8),
                      compact_printing = true,
                      backend = Val(sf_table))
         println()

         # *----- Auxiliary Table, log parameters to original scales --------

         auxtable = Array{Any}(undef,2,3)
         rn = 0 # row index

         if size(wvar,2) == 1 # single variable
            varstd = sqrt(sum((wvar .- sum(wvar)/length(wvar)).^2)/length(wvar))
            if varstd  <= 1e-7 # constant, assuming =1 anyway
                rn += 1
                auxtable[rn, 1] = :σᵤ²
                auxtable[rn, 2] = exp(_coevec[pos.begw])
                auxtable[rn, 3] = exp(_coevec[pos.begw])*stddev[pos.begw]
            end
         end

         if size(vvar,2) == 1 # single variable
            varstd = sqrt(sum((vvar .- sum(vvar)/length(vvar)).^2)/length(vvar))
            if varstd  <= 1e-7 # constant
                rn += 1
                auxtable[rn, 1] = :σᵥ²
                auxtable[rn, 2] = exp(_coevec[pos.begv])
                auxtable[rn, 3] = exp(_coevec[pos.begv])*stddev[pos.begv]
            end
         end

         if rn >= 1  # table is non-empty
             println("Convert the constant log-parameter to its original scale, e.g., σ² = exp(log_σ²):")
             pretty_table(auxtable[1:rn,:],
                          header=["", "Coef.", "Std.Err."],
                          formatters = ft_printf("%5.4f", 2:3),
                          compact_printing = true,
                          backend = Val(sf_table))

             print("\nTable format: "); printstyled("$(sf_table)"; color=:yellow); println(". Use sfmodel_opt() to choose between text, html, and latex.")
             println()
         end

         printstyled("***** Additional Information *********\n"; color=:cyan)

         print("* OLS (frontier-only) log-likelihood: "); printstyled(round(ll_ols; digits=5); color=:yellow); println("")
         print("* Skewness of OLS residuals: "); printstyled(round(sk_ols; digits=5); color=:yellow); println("")
         if _dicOPT[:ineff_index]
            print("* The sample mean of the JLMS inefficiency index: "); printstyled(round(_jlmsM; digits=5); color=:yellow); println("")
            print("* The sample mean of the BC efficiency index: "); printstyled(round(_bcM; digits=5); color=:yellow); println("\n")
         end
         if length(margMinfo) >= 1
            print("* The sample mean of inefficiency determinants' marginal effects on E(u): " ); printstyled(margMinfo; color=:yellow); println("")
            println("* Marginal effects of the inefficiency determinants at the observational level are saved in the return. See the follows.\n")
         end

         println("* Use `name.list` to see saved results (keys and values) where `name` is the return specified in `name = sfmodel_fit(..)`. Values may be retrieved using the keys. For instance:")
         println("   ** `name.loglikelihood`: the log-likelihood value of the model;")
         println("   ** `name.jlms`: Jondrow et al. (1982) inefficiency index;")
         println("   ** `name.bc`: Battese and Coelli (1988) efficiency index;")
         println("   ** `name.marginal`: a DataFrame with variables' (if any) marginal effects on E(u).")
         println("* Use `keys(name)` to see available keys.")

         printstyled("**************************************\n\n\n"; color=:cyan)

     end  # if_verbose

  #* ########### create a dictionary and make a tuple for return ########### *#

      _dicRES = OrderedDict{Symbol, Any}()
      _dicRES[:converged]          = Optim.converged(mfun)
      _dicRES[:iter_limit_reached] = Optim.iteration_limit_reached(mfun)
      _dicRES[:_______________] = "___________________"  #33
      _dicRES[:n_observations]  = num.nofobs
      _dicRES[:loglikelihood]   = -Optim.minimum(mfun)
      _dicRES[:table]           = [table][1]
      _dicRES[:coeff]           = _coevec
      _dicRES[:std_err]         = stddev
      _dicRES[:var_cov_mat]     = [var_cov_matrix][1]
      _dicRES[:jlms]            = _jlms
      _dicRES[:bc]              = _bc
      _dicRES[:OLS_loglikelihood] = ll_ols
      _dicRES[:OLS_resid_skew]    = sk_ols
      _dicRES[:marginal]      = margeff
      _dicRES[:marginal_mean] = margMinfo
      _dicRES[:_____________] = "___________________"  #31
      _dicRES[:model]         = minfo1
#     _dicRES[:data]          = "$sfdat"
      _dicRES[:depvar]        = _dicM[:depvar]
      _dicRES[:frontier]      = _dicM[:frontier]
      _dicRES[:μ]             = _dicM[:μ]
      _dicRES[:hscale]        = _dicM[:hscale]
      _dicRES[:gamma]         = _dicM[:gamma]
      _dicRES[:σₐ²]           = _dicM[:σₐ²]
      _dicRES[:σᵤ²]           = _dicM[:σᵤ²]
      _dicRES[:σᵥ²]           = _dicM[:σᵥ²]
      _dicRES[:log_σₐ²]       = _dicM[:σₐ²]
      _dicRES[:log_σᵤ²]       = _dicM[:σᵤ²]
      _dicRES[:log_σᵥ²]       = _dicM[:σᵥ²]
      _dicRES[:type]          = _dicM[:type]
      _dicRES[:dist]          = _dicM[:dist]
      _dicRES[:PorC]          = _porc
      _dicRES[:timevar]       = _dicM[:timevar]
      _dicRES[:idvar]         = _dicM[:idvar] # for bootstrap marginal effect
      _dicRES[:table_format]  = _dicOPT[:table_format]
      _dicRES[:modelid]       = tagD[:modelid]
      _dicRES[:verbose]       = _dicOPT[:verbose]
      _dicRES[:hasDF]         = _dicM[:hasDF]
      _dicRES[:transfer]      = _dicM[:transfer]

    for i in 1:length(eqvec2)
        _dicRES[keys(eqvec2)[i]] = _coevec[eqvec2[i]]
    end

      _dicRES[:________________]  = "___________________" #34
      _dicRES[:Hessian]           = [numerical_hessian][1]
      _dicRES[:gradient_norm]     = Optim.g_residual(mfun)
    # _dicRES[:trace]             = Optim.trace(mfun)     # comment out because not very informative and size could be large
      _dicRES[:actual_iterations] = Optim.iterations(mfun)
      _dicRES[:______________] = "______________________" #32
      _dicRES[:warmstart_solver] = sf_ini_algo_dic
      _dicRES[:warmstart_ini]    = sf_init_1st_dic
      _dicRES[:warmstart_maxIT]  = sf_ini_maxit_dic
      _dicRES[:main_solver]      = sf_algo
      _dicRES[:main_ini]         = sf_init_2nd_dic
      _dicRES[:main_maxIT]       = sf_maxit
      _dicRES[:tolerance]        = sf_tol
      _dicRES[:eqpo]             = eqvec2

      _dicRES[:redflag]          = redflag

     #* ----- Delete optional keys that have value nothing,

         for k in (:μ, :hscale,  :gamma)
             if _dicRES[k] === nothing
                delete!(_dicRES, k)
             end
         end

     #* ----- Create a NamedTuple from the dic as the final output;
     #* -----     put the dic in the tuple.

         _ntRES = NamedTuple{Tuple(keys(_dicRES))}(values(_dicRES))
         _ntRES = (; _ntRES..., list    = _dicRES)

     #* ---- Create a gloal dictionary for sf_predict ----

        global _eqncoe
        _eqncoe = Dict{Symbol, Vector}()  # nullify and initiate new dictionaries when a new model is specified

        for i in 1:length(eqvec2)
            _eqncoe[keys(eqvec)[i]]  = _coevec[eqvec2[i]] # for sf_predict
        end

  #* ############  make returns  ############ *#

      return _ntRES

end # sfmodel_fit

#############################################################
###  bootstrapping std.err. of mean marginal effect   #######
#############################################################
"""
    sfmodel_boot_marginal(<keyword arguments>)

Bootstrap standard errors and obtain bias-corrected (BC) confidence intervals
for the mean marginal effects of inefficiency determinants. In default, return
a ``K x 2``` matrix of standard errors (1st column) and confidence intervals
(tuples, 2nd column), where ``K`` is the number of exogenous inefficiency
determinants. With `getBootData=true`, return two matrices: the first is the
same as in the default return, and the second is the ``R x K`` bootstrapped
data.

See also the help file on `sfmodel_CI()`.

# Arguments
- result=<returned result>: The returned result from `sfmodel_fit()`.
- data=<dataset>: The DataFrame dataset containing the model's data.
  Same as the one used in `sfmodel_fit()`. If the data was supplied by matrix
  (instead of DataFrame; i.e., the Method 2 of `sfmodel_spec()`), this option
  should be skipped.
- R::Integer=<number>: The number of bootstrapped samples. The default is 500.
- level::Real=<number>: The significance level (default=0.05) of the bias-corrected
  confidence intervals. If `level`>0.5, it is automatically transformed to
  `1-level`, such that `level=0.05` and `level=0.95` both return 95%
  confidene intervals at the 5% significance level.
- seed::Integer=<number>: A postive integer used to seed the random
  number generator (rng) for resampling, which ensures reproducibility.
  This rng is not global and is only effective in this function. If not
  specified, the global random number generator is used, and the bootstrap
  result may change (slightly) between different runs.
- iter::Integer=<number>: The maximum number of iterations for each bootstrapped
  sample. If the number is larger than 0, it overwrites the `main_maxIT`
  specified in `sfmodel_opt()` which is the default.
- getBootData::Bool=false: Whether to return the bootstrapped data which is
  ``R x K`` where K is the number of exogenous determinants of inefficiency.
- every::Integer=10: Print bootstrapping progress for every `every` samples.

# Remarks
- Bootstrap samples are with replacement. For panel data, it samples
  cross-sectional units with replacement.
- In the MLE estimation, estimated
  coefficients from the main result is used as initial values. There is no
  `warmstart`. The `main_solver`, `main_maxIT`, and `tolerance` specified in
  `sfmodel_opt()` are used as default, but the value of `main_maxIT` may be
  replaced by the `iter` option.

# Examples
```julia
julia> std_ci = sfmodel_boot_marginal(result=res, data=df, R=250, seed=123)
bootstrap in progress..10..20..30..40..50..60..70..80..90..100..110..120..130..140..150..160..170..180..190..200..210..220..230..240..250..Done!

┌────────┬──────────────────────┬─────────────────┬──────────────────────┐
│        │ mean of the marginal │ std.err. of the │       bias-corrected │
│        │       effect on E(u) │     mean effect │    95.0%  conf. int. │
├────────┼──────────────────────┼─────────────────┼──────────────────────┤
│    age │             -0.00264 │         0.00225 │   (-0.00734, 0.0016) │
│ school │             -0.01197 │         0.01765 │    (-0.048, 0.01224) │
│     yr │             -0.02650 │         0.01221 │ (-0.05257, -0.00447) │
└────────┴──────────────────────┴─────────────────┴──────────────────────┘

3×2 Matrix{Any}:
 0.00224786  (-0.00734, 0.0016)
 0.0176472   (-0.048, 0.01224)
 0.012213    (-0.05257, -0.00447)

julia> std_ci, bsdata = sfmodel_boot_marginal(result=res, data=df, R=250, seed=123, getBootData=true);

 (output omitted)

julia> bsdata
250×3 adjoint(::Matrix{Real}) with eltype Real:
-0.000493033  -0.00840051   -0.0339088
-0.00510175   -0.0128184    -0.0147759
-0.0024987    -0.00971428   -0.00989262
-0.00157346   -0.0265515    -0.0140001
-0.00352179   -0.00670365   -0.0246122
-0.00375162   -0.0070496    -0.0306034
-0.00153094   -0.0154201    -0.0367731
0.000149329   0.00672036   -0.0461389
⋮
-0.00373306   -0.0108364    -0.00871898
-0.00170254   -0.0393002    -0.0500638
0.000686169   0.00241594   -0.018542
0.000258745   0.000183392  -0.039621
-0.00408104   -0.014574     -0.024126
-0.00417206   -0.0192443    -0.0406959
0.00266017   -0.0396552    -0.0359759
```
"""
function sfmodel_boot_marginal(; result::Any=nothing,  data::Any=nothing,
                                 R::Integer=500,
                                 level::Real=0.05,
                                 mymisc=nothing,
                                 iter::Integer=-1,
                                 getBootData::Bool=false,
                                 seed::Integer=-1,
                                 every::Integer=10)

  if data === nothing
      if result.transfer == true
          data = _dicM[:sdf]
      else
         throw("Need the `data=<DataFrame name>` option.")
      end
  elseif typeof(data) != DataFrame
      throw("`data=` needs to be a DataFrame.")
  end

   # getvar still use info from _dicM

  if result.marginal_mean==NamedTuple() || result.marginal_mean === nothing
      printstyled("The model does not have exogenous determinants of inefficiency, ", #=
               =# "or you didn't compute the marginal effect in the main estimation. Abort the bootstrap.\n"; color = :red)
      throw("No statistics to bootstrap.")
  end

  ((level > 0.0) && (level < 1.0)) || throw("The significance level (`level`) should be between 0 and 1.")

  if level > 0.5
     level = 1-level  # 0.95 -> 0.05
  end

     # In the following lines, the integer part had been taken care of in Type.
  (seed == -1) || ( seed > 0) || throw("`seed` needs to be a positive integer.")
  (iter == -1) || ( iter > 0) || throw("`iter` needs to be a positive integer.")
  (R > 0) || throw("`R` needs to be a positive integer.")

#* ##### Get parameters #######

          _porc = result.PorC
        sf_init = result.coeff
        sf_algo = eval(result.main_solver)
       sf_maxit = result.main_maxIT
         sf_tol = result.tolerance
       sf_table = result.table_format

       if iter > 0
           sf_maxit = iter
       end

#* ########  parepare data  ##########

  (minfo1, minfo2, pos, num, eqvec, eqvec2, yvar, xvar, zvar, qvar, wvar,
   vvar,         rowIDT, varlist) = SFrontiers.getvar(result.modelid, data)

  #----- make an orginal copy; also make later code easier --------

   yvar0 = deepcopy(yvar); xvar0 = deepcopy(xvar); zvar0 = deepcopy(zvar);
   qvar0 = deepcopy(qvar); wvar0 = deepcopy(wvar); vvar0 = deepcopy(vvar);
   rowIDT0 = deepcopy(rowIDT)

#* ####### begin looping #############

 sim_res = Array{Real}(undef, ncol(result.marginal), R)

 if seed > 0
    rng = MersenneTwister(seed)
 end

 for i in 1:R

    @label start1

   # ----- sampling data for cross-sectional data --------------

   if result.idvar === nothing  # i.e., cross-sectional

     if seed == -1
        select_row = sample( 1:nrow(data), nrow(data); replace=true)  # require StatsBase.jl
     else
        select_row = sample(rng, 1:nrow(data), nrow(data); replace=true)  # require StatsBase.jl
     end

     yvar =  yvar0[select_row, :]
     xvar =  xvar0[select_row, :]
     (zvar0 == ()) ||  (zvar =  zvar0[select_row, :])
     (qvar0 == ()) ||  (qvar =  qvar0[select_row, :])
     (wvar0 == ()) ||  (wvar =  wvar0[select_row, :])
     (vvar0 == ()) ||  (vvar =  vvar0[select_row, :])

  end

   # ------ sampling panel data ------

   if result.idvar !== nothing  # i.e., panel data

      if seed == -1
        select_row = sample(     1:size(rowIDT0,1), size(rowIDT0,1); replace=true)  # require StatsBase.jl
      else
        select_row = sample(rng, 1:size(rowIDT0,1), size(rowIDT0,1); replace=true)  # require StatsBase.jl
      end

      yvar = mySample(rowIDT0, yvar0, select_row)

      xvar = mySample(rowIDT0, xvar0, select_row)
      (zvar0 == ()) ||  (zvar = mySample(rowIDT0, zvar0, select_row))
      (qvar0 == ()) ||  (qvar = mySample(rowIDT0, qvar0, select_row))
      (wvar0 == ()) ||  (wvar = mySample(rowIDT0, wvar0, select_row))
      (vvar0 == ()) ||  (vvar = mySample(rowIDT0, vvar0, select_row))

       newidvar = mySample2!(rowIDT0, ones(size(yvar,1)), select_row) # a trick to get new id; note not yvar0
       rowIDT = get_rowIDT(vec(newidvar))
       # yvar = vec(yvar)

    end

 #*  --- other small things -----

   nofobs1 = size(yvar,1)

  #* ----- Define the problem's Hessian -----#

   _Hessian = TwiceDifferentiable(rho -> SFrontiers.LL_T(result.modelid,
                         yvar, xvar, zvar, qvar, wvar, vvar,
                         _porc, nofobs1, pos, rho,
                                 rowIDT, mymisc),
                          sf_init;
                          autodiff = :forward);

  #* ---- estimate ----------------- *#

             mfun = try
                       optimize(_Hessian,
                                sf_init,       # initial values
                                sf_algo,       # different from search run
                                Optim.Options(g_tol = sf_tol,
                                              iterations  = sf_maxit, # different from search run
                                              store_trace = false,
                                              show_trace  = false))
                    catch err
                       @goto start1
                    end

#* ###### check if valid ###############

  if (Optim.iteration_limit_reached(mfun) ) ||
     (isnan(Optim.g_residual(mfun)) ) ||
     (Optim.g_residual(mfun) > 1e-1)
         @goto start1
  end

    _coevec            = Optim.minimizer(mfun)  # coef. vec.
    numerical_hessian  = hessian!(_Hessian, _coevec)  # Hessain

   #* ------ Check if the matrix is invertible. ----

   var_cov_matrix = try
                       inv(numerical_hessian)
                    catch err
                       @goto start1
                    end

        if !all( diag(var_cov_matrix) .> 0 ) # not all are positive
              @goto start1
        end

   #* ---- marginal effect on E(u) --------------

   margeff, margMinfo = get_marg(result.modelid, pos, num, _coevec, zvar, qvar, wvar)

   theM = mean.(eachcol(margeff))  # mean marginal effect

   if sum(isnan.(theM)) != 0  # in this run some of the element is NaN
      @goto start1
   end

   sim_res[:, i] = theM

   if result.verbose
      if i == 1
         printstyled("bootstrap in progress.."; color = :cyan)
      end
      if i%every == 0   # print for every 10 replications
          print("$(i)..")
      end
      if i == R
          printstyled("Done!\n\n"; color = :cyan)
      end
   end

end   # for i=1:R

 sim_res = sim_res'  # nofobs x K; ith row = ith replication, column = statistics

 #*###### compute statistics for the mean marginal effect #####*#

   theMean = mean.(eachcol(result.marginal))  # Kx1 vector
   theSTD  = sqrt.( sum((sim_res .- theMean').^2, dims=1) ./(R-1)) # 1 x K

   ci_mat = sfmodel_CI(bootdata=sim_res, observed=theMean, level=level, verbose=false);

   if result.verbose

      table = Array{Any}(undef, length(theMean), 3+1)
      vname = names(result.marginal)

      for j in 1:length(theMean)
          table[j, 1] = string(vname[j])[6:end]
      end

      table[:,2], table[:,3] = theMean, theSTD
      table[:,4] = ci_mat

       mylevel = 100*(1-level)

      table = [" " "mean of the marginal" "std.err. of the"  "bias-corrected";
               " " "effect on E(u)"       "mean effect"      "$(mylevel)%  conf. int.";
               table]

            pretty_table(table,
                         noheader= true,
                         body_hlines = [2],
                         formatters = ft_printf("%0.5f", 2:4),
                         compact_printing = true,
                         backend = Val(sf_table))
            println()

   end

   if getBootData == false   # most likely
      return hcat(theSTD', ci_mat) # K x 1
   else
      return hcat(theSTD', ci_mat), sim_res
   end

end # sfmodel_boot_marginal

# -----------------------------------------------#
# --- utilities for sfmodel_boot_marginal -------#
# -----------------------------------------------#

function mySample(rowIDT, data, select)
  temp1 = Array{Float64}(undef, 0, size(data,2))
  @inbounds for i=1:size(rowIDT,1)
       @views   cc = data[rowIDT[select[i],1], :]
                temp1 = vcat(temp1, cc)
            end
  return temp1
end

function mySample2!(rowIDT, data, select)
  temp1 = Array{Float64}(undef, 0, size(data,2))
  @inbounds for i=1:size(rowIDT,1)
       @views   cc = data[rowIDT[select[i],1], :]
                cc[:,1] .= i     # caution, it changes data
                temp1 = vcat(temp1, cc)
            end
  return temp1
end

"""
function sfmodel_CI(; bootdata::Any=nothing, observed::Union{Vector, Real, Tuple, NamedTuple}=nothing, level::Real=0.05, verbose::Bool=true)

    sfmodel_CI(<keyword arguments>)

A general purpose (not specific to stochastic frontier models) function for
obtaining bias-correctred (BC) confidence intervals from bootstrapped data.
Return a ``K x 1`` matrix of confidence intervals in the form of tuples, where
``K`` is the number of bootstrap statistics.

See also the help file on `sfmodel_boot_marginal()`.

# Arguments
- bootdata::Array=<data>: The bootstrapped data of size ``R x K``,
  where ``R`` is the number of bootstrap samples (replications) and ``K`` is
  the number of statistics. An example is the bootstrapped data from
  `sfmodel_boot_marginal(, ... getBootData=true)`.
- observed::Union{Vector, Real, Tuple, NamedTuple}=<a vector of numbers>: The
  observed values of the statistics to which the confidence intervals are to
  be calculated. The length of `observed` should be equal to ``K``. It could take
  the form of a single value (if ``K=1``), a vector, a tuple, or a NamedTuple.
- level::Real=<number>: The significance level (default=0.05) of the bias-corrected
  confidence intervals. If `level`>0.5, it is automatically transformed to
  `1-level`, such that `level=0.05` and `level=0.95` both return 95%
  confidene intervals at the 5% significance level.
- verbose::Bool=true: Print the result.

# Examples
```julia
julia> myans = sfmodel_fit(useData(df));

 (output omitted)

julia> std_ci, bsdata = sfmodel_boot_marginal(result=myans, data=df, R=250, seed=123, getBootData=true);

 (output omitted)

julia> sfmodel_CI(bootdata=bsdata, observed=myans.marginal_mean, level=0.10)

Bias-Corrected 90.0% Confidence Interval:

3×1 Matrix{Any}:
 (-0.00655, 0.0009)
 (-0.04064, 0.0086)
 (-0.04663, -0.00872)

julia> # manually input observed values

julia> sfmodel_CI(bootdata=bsdata, observed=(-0.00264, -0.01197, -0.0265), level=0.10)

Bias-Corrected 90.0% Confidence Interval:

3×1 Matrix{Any}:
 (-0.00655, 0.0009)
 (-0.04064, 0.0086)
 (-0.04663, -0.00872)
```
"""
function sfmodel_CI(; bootdata::Any=nothing, observed::Union{Vector, Real, Tuple, NamedTuple}=nothing, level::Real=0.05, verbose::Bool=true)

    # bias-corrected (but not accelerated) confidence interval
    # For the "accelerated" factor, need to estimate the SF model
    #    for every jack-knifed sample, which is expensive.

    if (observed isa NamedTuple)
        observed = values(observed)
    end

   ((level > 0.0) && (level < 1.0)) || throw("The significance level (`level`) should be between 0 and 1.")

   if level > 0.50
      level = 1-level  # 0.95 -> 0.05
   end

    nofobs, nofK = size(bootdata)  # number of statistics

    (nofK == length(observed)) || throw("The number of statistics (`observed`) does not fit the number of columns of bootstrapped data.")

    ci = Array{Any}(undef, nofK, 1)

    z1 = quantile(Normal(), level/2)
    z2 = quantile(Normal(), 1 - level/2)  #! why z1 != z2?

    for i in 1:nofK
        @views data = bootdata[:,i]

        count = sum(data .< observed[i])
        z0 = quantile(Normal(), count/nofobs) # bias corrected factor

        alpha1 = cdf(Normal(), z0 + ((z0 + z1) ))
        alpha2 = cdf(Normal(), z0 + ((z0 + z2) ))

        order_data = sort(data)

        ci[i,1] = (  round(order_data[Int(ceil(nofobs*alpha1))], digits=5),    round(order_data[Int(ceil(nofobs*alpha2))], digits=5)  )
    end

    if verbose == true
       mylevel = 100*(1-level)
       println("\nBias-Corrected $(mylevel)% Confidence Interval:\n")
    end

    return ci
end

########################################################
####            catching type error                 ####
########################################################

function sfmodel_init(arg::Vector)
    throw(ArgumentError("The initial values in sfmodel_init() are specified incorrectly. They must be supplied in macros such as all_init(0.1, 0.1)."))
end

function sfmodel_init(args::Number)
    throw(ArgumentError("The initial values in sfmodel_init() are specified incorrectly. They must be supplied in macros such as all_init(0.1, 0.1)."))
end

function sfmodel_fit(sfdat::Any)
    throw(ArgumentError("The dataset specified in sfmodel_fit() must be a DataFrame."))
end


text:
# Regressors1D is an iterable to access dynamically the context,
# i.e. x[t-1], ..., x[t-K], of a time step.

struct Regressors1D{K}
    x::PaddedView
end

function Regressors1D(K::Int64, x::Array{T}) where T <: AbstractFloat
    px = PaddedView(T(0.), x, (1-K:length(x),))
    Regressors1D{K}(px)
end

Base.length(r::Regressors1D) = length(parent(r.x))

function Base.getindex(r::Regressors1D{K}, t::Int) where K
    retval = r.x[t-K:t]
    retval[end] = 1 # "bias" regressor
    retval
end

function Base.iterate(r::Regressors1D{K}, state = nothing) where K
    if isnothing(state)
        T = eltype(r.x)
        state = (Vector{T}(undef, K+1), 1)
    end

    buffer, t = state
    if t > length(r) return nothing end

    buffer[:] = r.x[t-K:t]
    buffer[end] = 1 # bias regressor

    # we return a shallow copy of the buffer to avoid issue if the user
    # does not copy the data himself.
    return (copy(buffer), (buffer, t+1))
end

ext: jl

max_stars_repo_name: BUTSpeechFIT/ConvHMM

max_issues_repo_name: BUTSpeechFIT/ConvHMM

max_forks_repo_name: BUTSpeechFIT/ConvHMM


text:
basic_params, maneuvering_params = get_KVLCC2_L7_params()

@testset "mmg.jl KVLCC2_L7 turning" begin

    duration = 200  # [s]
    max_δ_rad = 35 * pi / 180.0  # [rad]
    n_const = 17.95  # [rps]

    sampling = duration * 10
    time_list = range(0.00, stop = duration, length = sampling)
    δ_rad_list = max_δ_rad .* ones(Float64, sampling)
    n_p_list = n_const .* ones(Float64, sampling)
    mmg_results = mmg_3dof_simulate(
        basic_params,
        maneuvering_params,
        time_list,
        δ_rad_list,
        n_p_list,
        u0 = 2.29 * 0.512,
        v0 = 0.0,
        r0 = 0.0,
    )
    u, v, r, δ, n_p = mmg_results
    x, y, ψ = calc_position(time_list, u, v, r)
end

@testset "mmg_zigzag_test" begin
    target_δ_rad = 20.0 * π / 180.0
    target_ψ_rad_deviation = 20.0 * π / 180.0
    start_time_second = 0.00
    time_second_interval = 0.01
    end_time_second = 80.00
    time_list = start_time_second:time_second_interval:end_time_second
    n_const = 17.95  # [rps]
    n_p_list = n_const * ones(Float64, length(time_list))
    u_list, v_list, r_list, ψ_list, δ_list = mmg_3dof_zigzag_test(
        basic_params,
        maneuvering_params,
        time_list,
        n_p_list,
        target_δ_rad,
        target_ψ_rad_deviation,
    )
end

@testset "mmg 3dof approx least square method func" begin
    duration = 200  # [s]
    sampling = duration * 10
    time_list = range(0.00, stop = duration, length = sampling)
    max_δ_rad = 35 * pi / 180.0  # [rad]
    n_const = 17.95  # [rps]
    δ_rad_list = max_δ_rad .* ones(Float64, sampling)
    n_p_list = n_const .* ones(Float64, sampling)
    mmg_results = mmg_3dof_simulate(
        basic_params,
        maneuvering_params,
        time_list,
        δ_rad_list,
        n_p_list,
        u0 = 2.29 * 0.512,
        v0 = 0.0,
        r0 = 0.0,
    )
    u, v, r, δ, n_p = mmg_results
    x, y, ψ = calc_position(time_list, u, v, r)
    data = ShipData(time_list, u, v, r, x, y, ψ, δ, n_p)
    estimate_mmg_approx_lsm(
        data,
        basic_params,
        maneuvering_params.k_0,
        maneuvering_params.k_1,
        maneuvering_params.k_2,
    )
end

@testset "mmg 3dof approx time window sampling least square method func" begin
    duration = 200  # [s]
    max_δ_rad = 35 * pi / 180.0  # [rad]
    n_const = 17.95  # [rps]

    sampling = duration * 10
    time_list = range(0.00, stop = duration, length = sampling)
    δ_rad_list = max_δ_rad .* ones(Float64, sampling)
    n_p_list = n_const .* ones(Float64, sampling)
    mmg_results = mmg_3dof_simulate(
        basic_params,
        maneuvering_params,
        time_list,
        δ_rad_list,
        n_p_list,
        u0 = 2.29 * 0.512,
        v0 = 0.0,
        r0 = 0.0,
    )
    u, v, r, δ, n_p = mmg_results
    noize_dist = Normal(0.0, 0.0005)
    u_obs = u + rand(noize_dist, size(u))
    v_obs = v + rand(noize_dist, size(v))
    r_obs = r + rand(noize_dist, size(r))
    x, y, ψ = calc_position(time_list, u_obs, v_obs, r_obs)
    data = ShipData(time_list, u_obs, v_obs, r_obs, x, y, ψ, δ, n_p)
    one_sample_size = 5000
    estimate_mmg_approx_lsm_time_window_sampling(
        data,
        500,
        basic_params,
        maneuvering_params.k_0,
        maneuvering_params.k_1,
        maneuvering_params.k_2,
    )
end

@testset "mmg mcmc sampling func" begin
    duration = 200  # [s]
    max_δ_rad = 35 * pi / 180.0  # [rad]
    n_const = 17.95  # [rps]

    sampling = duration * 10 + 1
    time_list = range(0.00, stop = duration, length = sampling)
    δ_rad_list = max_δ_rad .* ones(Float64, sampling)
    n_p_list = n_const .* ones(Float64, sampling)
    mmg_results = mmg_3dof_simulate(
        basic_params,
        maneuvering_params,
        time_list,
        δ_rad_list,
        n_p_list,
        u0 = 2.29 * 0.512,
        v0 = 0.0,
        r0 = 0.0,
    )
    u, v, r, δ, n_p = mmg_results
    sampling_rate = 10
    time_obs = time_list[1:sampling_rate:end]
    noize_dist = Normal(0.0, 0.01)
    u_obs = u + rand(noize_dist, size(u))
    v_obs = v + rand(noize_dist, size(v))
    r_obs = r + rand(noize_dist, size(r))
    x, y, ψ = calc_position(time_obs, u_obs, v_obs, r_obs)
    data = ShipData(
        time_obs,
        u_obs,
        v_obs,
        r_obs,
        x,
        y,
        ψ,
        δ[1:sampling_rate:end],
        n_p[1:sampling_rate:end],
    )
    n_samples = 10
    n_chains = 1
    model = create_model_for_mcmc_sample_mmg(
        data,
        basic_params,
        maneuvering_params.k_0,
        maneuvering_params.k_1,
        maneuvering_params.k_2,
    )
    chain = nuts_sampling_single_thread(model, n_samples, n_chains)
end

max_stars_repo_name: mitsuyukiLab/JuliaShipMMG

max_issues_repo_name: mitsuyukiLab/JuliaShipMMG


text:
using Test
using LightGraphs
using Serialization
include("../src/pmg_graphs.jl")
include("../src/atomgraph.jl")

@testset "AtomGraph" begin
    # build a silly little triangle graph
    g = SimpleWeightedGraph{Int32}(Float32.([0 1 1; 1 0 1; 1 1 0]))

    # add an element list that doesn't make sense
    @test_throws AssertionError AtomGraph(g, ["C"])

    # okay, now do it right, start with no features
    ag = AtomGraph(g, ["C", "C", "C"])

    # check LightGraphs fcns
    @test eltype(ag)==Int32
    @test edgetype(ag)==SimpleWeightedEdge{Int32,Float32}
    @test ne(ag)==3
    @test nv(ag)==3
    @test !is_directed(ag)
    # not sure the best way to test the ones that return iterators, e.g. edges
    @test outneighbors(ag,1)==inneighbors(ag,1)==[2,3]
    @test has_vertex(ag,1)
    @test !has_vertex(ag,4)
    @test has_edge(ag,1,2)

    # add some features
    bad_fmat = Float32.([1 2; 3 4])
    good_fmat = Float32.([1 2 3; 4 5 6])
    featurization = [AtomFeat(:feat, true, 2, false, ['a','b'])]
    @test_throws AssertionError add_features!(ag, bad_fmat, featurization)
    add_features!(ag, good_fmat, featurization)
    @test ag.features==good_fmat

    # tests for other signatures of add_features! where feature vectors are built automatically
    ag = AtomGraph(g, ["C", "C", "C"])
    ag2 = deepcopy(ag)
    ag3 = deepcopy(ag)
    feature_names = [:Block, :X]
    vecs, featurization = make_feature_vectors(feature_names, nbins=[4,3])

    add_features!(ag, vecs, featurization)
    add_features!(ag2, featurization)
    add_features!(ag3, feature_names, Int32.([4,3]))

    @test ag.features==ag2.features==ag3.features==Float32.(hcat([vecs["C"] for i in 1:3]...))
end

@testset "graph-building" begin
    ag = build_graph(joinpath(@__DIR__, "test_data", "mp-195.cif"), use_voronoi=true)
    wm_true = [0.0 1.0 1.0 1.0; 1.0 0.0 1.0 1.0; 1.0 1.0 0.0 1.0; 1.0 1.0 1.0 0.0]
    @test weights(ag) == wm_true
    @test ag.elements == ["Ho", "Pt", "Pt", "Pt"]
    ag = build_graph(joinpath(@__DIR__, "test_data", "mp-195.cif"); use_voronoi=false)
    @test weights(ag) == wm_true
    @test ag.elements == ["Ho", "Pt", "Pt", "Pt"]

    # tests for some other file formats
    graphs = AtomGraph[]
    for fp in ["mp-195.poscar", "mp-195.traj", "mp-195.xyz"]
        push!(graphs, build_graph(joinpath(@__DIR__, "test_data", fp)))
    end
    for g in graphs
        @test ag.elements == g.elements
        @test weights(ag) == weights(g)
    end

    # test for nonperiodic system
    @test_logs (:warn, "Voronoi edge weights are not supported if any direction in the structure is nonperiodic. Using cutoff weights method...") build_graph(joinpath(@__DIR__,"test_data", "methane.xyz"), use_voronoi=true)
    methane = build_graph(joinpath(@__DIR__,"test_data", "methane.xyz"))
    @test all(isapprox.(weights(methane)[2:5,1], 1.0, atol=1e-4))
    @test all(isapprox.(weights(methane)[3:2,2], 0.375, atol=1e-5))
end

@testset "save/load" begin
    g = SimpleWeightedGraph{Int32}(Float32.([0 1 1; 1 0 1; 1 1 0]))
    @test_throws AssertionError AtomGraph(g, ["C"])
    fmat = Float32.([1 2 3; 4 5 6; 0 1 0; 9 8 7])
    featurization = [AtomFeat(:feat, true, 2, false, ['a','b']), AtomFeat(:feat2, false, 2, false, [-1,0,1])]
    ag = AtomGraph(g, ["C", "C", "C"], fmat, featurization)
    serialize(joinpath(@__DIR__, "test_data", "testgraph.jls"), ag)
    ag2 = deserialize(joinpath(@__DIR__, "test_data", "testgraph.jls"))
    @test ag2.elements==ag.elements
    @test ag.lapl==ag2.lapl
    @test ag.features==ag2.features
    for i in 1:2
        for field in [:name, :categorical, :num_bins, :logspaced, :vals]
            @test getfield(ag.featurization[i], field)==getfield(ag2.featurization[i], field)
        end
    end
end

@testset "batch processing" begin
    featurization = build_atom_feats([Symbol("Atomic mass"), :Block])
    build_graphs_batch(joinpath(@__DIR__, "test_data"), joinpath(@__DIR__, "test_data", "graphs"), featurization)
    g1 = deserialize(joinpath(@__DIR__, "test_data","graphs","mp-195.jls"))
    @test size(g1)==(4,4)
    @test size(g1.features)==(14,4)
    g2 = deserialize(joinpath(@__DIR__, "test_data","graphs","mp-224.jls"))
    @test size(g2)==(6,6)
    @test size(g2.features)==(14,6)
    w = weights(g2)
    @test w[1,2]==w[1,3]==w[1,5]==w[2,4]==w[2,6]==w[3,4]==w[5,6]==0.0
    @test w[3,3]==w[4,4]==w[5,5]==w[6,6]==1.0
    @test g2.elements==["W","W","S","S","S","S"]

    # test read_graphs_batch
    gs = read_graphs_batch(joinpath(@__DIR__, "test_data", "graphs"))
    @test length(gs)>=2
end


text:
using Test
using YaoHIR
using YaoLocations
using CompilerPluginTools
using ZXCalculus
using ZXCalculus: Phase, BlockIR
using YaoHIR.IntrinsicOperation
using YaoHIR: Chain, Gate, Ctrl

ir = @make_ircode begin
    Expr(:call, :+, 1, 1)::Int
    Expr(:call, :+, 3, 3)::Int
end

circ = Chain(
    Gate(shift(1.0), Locations(1)),
    Gate(Rz(Core.SSAValue(2)), Locations(1)),
)

bir = BlockIR(ir, 4, circ)

zxd = convert_to_zxd(bir)
c = clifford_simplification(zxd)
ZXCalculus.generate_layout!(zxd)
qc_tl = convert_to_chain(phase_teleportation(zxd))
@test length(qc_tl) == 1
p = Phase(1//1)
@test p + 1 == 1 + p == p + p
@test p - 1 == 1 - p == p - p
@test p / 1 == 1 / p == p / p
@test zero(Phase) == zero(p) && one(Phase) == one(p)


text:
"""
Module: SymbolicDiff (Symbolic Operation for Arithmetic)
"""

const SymbolicEnv = Dict{Symbol,Any}
const SymbolicCache = Dict{Union{AbstractSymbolic,Tuple{AbstractSymbolic,Symbol},Tuple{AbstractSymbolic,Tuple{Symbol,Symbol}}},Any}
const globalenv = SymbolicEnv()

function Base.setindex!(h::SymbolicEnv, x::Any, k::SymbolicVariable{Tv}) where Tv
    h[k.var] = x
end


text:
immutable StatStruct
    device  :: Uint
    inode   :: Uint
    mode    :: Uint
    nlink   :: Int
    uid     :: Uint
    gid     :: Uint
    rdev    :: Uint
    size    :: Int64
    blksize :: Int64
    blocks  :: Int64
    mtime   :: Float64
    ctime   :: Float64
end

StatStruct(buf::Union(Vector{Uint8},Ptr{Uint8})) = StatStruct(
     uint(ccall(:jl_stat_dev,     Uint32,  (Ptr{Uint8},), buf)),
     uint(ccall(:jl_stat_ino,     Uint32,  (Ptr{Uint8},), buf)),
     uint(ccall(:jl_stat_mode,    Uint32,  (Ptr{Uint8},), buf)),
      int(ccall(:jl_stat_nlink,   Uint32,  (Ptr{Uint8},), buf)),
     uint(ccall(:jl_stat_uid,     Uint32,  (Ptr{Uint8},), buf)),
     uint(ccall(:jl_stat_gid,     Uint32,  (Ptr{Uint8},), buf)),
     uint(ccall(:jl_stat_rdev,    Uint32,  (Ptr{Uint8},), buf)),
    int64(ccall(:jl_stat_size,    Uint64,  (Ptr{Uint8},), buf)),
    int64(ccall(:jl_stat_blksize, Uint64,  (Ptr{Uint8},), buf)),
    int64(ccall(:jl_stat_blocks,  Uint64,  (Ptr{Uint8},), buf)),
          ccall(:jl_stat_mtime,   Float64, (Ptr{Uint8},), buf),
          ccall(:jl_stat_ctime,   Float64, (Ptr{Uint8},), buf),
)

show(io::IO, st::StatStruct) = print("StatStruct(mode=$(oct(st.mode,6)), size=$(st.size))")

# stat & lstat functions

const stat_buf = Array(Uint8, ccall(:jl_sizeof_stat, Int32, ()))
macro stat_call(sym,arg1type,arg)
    quote
        fill!(stat_buf,0)
        r = ccall($(Expr(:quote,sym)), Int32, ($arg1type,Ptr{Uint8}), $(esc(arg)), stat_buf)
        r==0 || r==UV_ENOENT || r==UV_ENOTDIR || throw(UVError("stat",r))
        st = StatStruct(stat_buf)
        if ispath(st) != (r==0)
            error("stat returned zero type for a valid path")
        end
        st
    end
end

stat(fd::RawFD)     = @stat_call jl_fstat Int32 fd.fd
stat(fd::Integer)   = @stat_call jl_fstat Int32 fd
stat(path::ByteString)  = @stat_call jl_stat  Ptr{Uint8} path
stat(path::String) = stat(bytestring(path))
lstat(path::ByteString) = @stat_call jl_lstat Ptr{Uint8} path
lstat(path::String) = lstat(bytestring(path))

stat(path...) = stat(joinpath(path...))
lstat(path...) = lstat(joinpath(path...))

# mode type predicates

    ispath(st::StatStruct) = st.mode & 0xf000 != 0x0000
    isfifo(st::StatStruct) = st.mode & 0xf000 == 0x1000
 ischardev(st::StatStruct) = st.mode & 0xf000 == 0x2000
     isdir(st::StatStruct) = st.mode & 0xf000 == 0x4000
isblockdev(st::StatStruct) = st.mode & 0xf000 == 0x6000
    isfile(st::StatStruct) = st.mode & 0xf000 == 0x8000
    islink(st::StatStruct) = st.mode & 0xf000 == 0xa000
  issocket(st::StatStruct) = st.mode & 0xf000 == 0xc000

# mode permission predicates

issetuid(st::StatStruct) = (st.mode & 0o4000) > 0
issetgid(st::StatStruct) = (st.mode & 0o2000) > 0
issticky(st::StatStruct) = (st.mode & 0o1000) > 0

  isreadable(st::StatStruct) = (st.mode & 0o444) > 0
  iswritable(st::StatStruct) = (st.mode & 0o222) > 0
isexecutable(st::StatStruct) = (st.mode & 0o111) > 0

uperm(st::StatStruct) = uint8(st.mode >> 6) & 0x7
gperm(st::StatStruct) = uint8(st.mode >> 3) & 0x7
operm(st::StatStruct) = uint8(st.mode     ) & 0x7

# mode predicate methods for file names

for f in {
    :ispath
    :isfifo
    :ischardev
    :isdir
    :isblockdev
    :isfile
    :islink
    :issocket
    :issetuid
    :issetgid
    :issticky
    :isreadable
    :iswritable
    :isexecutable
    :uperm
    :gperm
    :operm
}
    @eval ($f)(path...)  = ($f)(stat(path...))
end

islink(path...) = islink(lstat(path...))

# some convenience functions

filemode(path...) = stat(path...).mode
filesize(path...) = stat(path...).size
   mtime(path...) = stat(path...).mtime
   ctime(path...) = stat(path...).ctime

samefile(a::StatStruct, b::StatStruct) = a.device==b.device && a.inode==b.inode
samefile(a::String, b::String) = samefile(stat(a),stat(b))

max_stars_repo_name: porterjamesj/julia

max_issues_repo_name: porterjamesj/julia

max_forks_repo_name: porterjamesj/julia


text:
# ---
# title: 297. Serialize and Deserialize Binary Tree
# id: problem297
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Tree, Design
# link: <https://leetcode.com/problems/serialize-and-deserialize-binary-tree/description/>
# hidden: true
# ---
#
# Serialization is the process of converting a data structure or object into a
# sequence of bits so that it can be stored in a file or memory buffer, or
# transmitted across a network connection link to be reconstructed later in the
# same or another computer environment.
#
# Design an algorithm to serialize and deserialize a binary tree. There is no
# restriction on how your serialization/deserialization algorithm should work.
# You just need to ensure that a binary tree can be serialized to a string and
# this string can be deserialized to the original tree structure.
#
# **Clarification:** The input/output format is the same as [how LeetCode
# serializes a binary tree](https://leetcode.com/faq/#binary-tree). You do not necessarily need to
# follow this format, so please be creative and come up with different
# approaches yourself.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/15/serdeser.jpg)
#
#
#
#     Input: root = [1,2,3,null,null,4,5]
#     Output: [1,2,3,null,null,4,5]
#
#
# **Example 2:**
#
#
#
#     Input: root = []
#     Output: []
#
#
# **Example 3:**
#
#
#
#     Input: root = [1]
#     Output: [1]
#
#
# **Example 4:**
#
#
#
#     Input: root = [1,2]
#     Output: [1,2]
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is in the range `[0, 104]`.
#   * `-1000 <= Node.val <= 1000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module How2Mix

# Unassuming code in an unassuming module...
module SubFoo

function h(x)
    x = rand()
    y = rand()
    return x + y
end

function f(x)
    z = rand()
    return h(x)
end

end

using Mixtape
import Mixtape: CompilationContext, transform, allow_transform, show_after_inference,
                show_after_optimization, debug

using IRTools

# 101: How2Mix
struct MyMix <: CompilationContext end

function transform(::MyMix, ir)
    for (v, st) in ir
        st.expr isa Expr || continue
        st.expr.head == :call || continue
        st.expr.args[1] == Base.rand || continue
        ir[v] = 5
    end
    display(ir)
    return ir
end

# MyMix will only transform functions which you explicitly allow.
#allow_transform(ctx::MyMix, fn::typeof(SubFoo.h), a...) = true
allow_transform(ctx::MyMix, m::Module) = m == SubFoo
show_after_inference(ctx::MyMix) = false
show_after_optimization(ctx::MyMix) = false
debug(ctx::MyMix) = true

fn = Mixtape.jit(MyMix(), SubFoo.f, Tuple{Float64})
@time fn = Mixtape.jit(MyMix(), SubFoo.f, Tuple{Float64})

fn(5.0)
@time fn(5.0)
display(fn(5.0))
display(SubFoo.f(5.0))

end # module


text:
@testset "Default" begin
    @testset "Library" begin
        @test default_library(Polyhedra.FullDim{2}(), Float64) isa SimplePolyhedraLibrary
        @test default_library(Polyhedra.FullDim{1}(), Float64) isa IntervalLibrary
        @test default_library(Polyhedra.FullDim{1}(), Float64) isa IntervalLibrary
        @test default_library(Polyhedra.FullDim{1}(), Float64) isa IntervalLibrary
        @test default_library(Polyhedra.FullDim{2}(), AbstractFloat) isa SimplePolyhedraLibrary{Float64}
        @test default_library(Polyhedra.FullDim{1}(), AbstractFloat) isa IntervalLibrary{Float64}
    end
    @testset "LP Solver" begin
        @test Polyhedra.default_solver(hrep([HalfSpace((@SVector [1, 2]), 3)])) === nothing
        @test Polyhedra.solver(hrep([HalfSpace((@SVector [1, 2]), 3)]), lp_solver) === lp_solver
    end
end


text:
"Your optional docstring here"
function distance(a, b)
  hamming = 0

  if length(a) != length(b)
    throw(ArgumentError("Invalid length"))
  end

  for i in 1:length(a)
    if a[i] != b[i]
      hamming += 1
    end
  end

  return hamming
end


text:
# This line should have no code
f2(x) = 2x
if true
  f3(x) = 3x
end
f4(x) = 4x; f5(x) = 5x
# This line should have no code
@doc """
`f6(x)` multiplies `x` by 6
"""
f6(x) = 6x
# This line should have no code
@doc """
`f7(x)` multiplies `x` by 7
"""
function f7(x)
    7x
end

f8(x) = 8x

# This line should have no code
module TestModule
end

# This line should have no code
baremodule TestBareModule
end


text:
#TODO: seed

"""
- The multi-discrete action space consists of a series of discrete action spaces with different number of actions in eachs
- It is useful to represent game controllers or keyboards where each key can be represented as a discrete action space
- It is parametrized by passing an array of positive integers specifying number of actions for each discrete action space

Note: A value of 0 always need to represent the NOOP action.

e.g. Nintendo Game Controller
- Can be conceptualized as 3 discrete action spaces:

    1) Arrow Keys: Discrete 5  - NOOP[0], UP[1], RIGHT[2], DOWN[3], LEFT[4]  - params: min: 0, max: 4
    2) Button A:   Discrete 2  - NOOP[0], Pressed[1] - params: min: 0, max: 1
    3) Button B:   Discrete 2  - NOOP[0], Pressed[1] - params: min: 0, max: 1

- Can be initialized as

    MultiDiscrete([ 5, 2, 2 ])
    or MultiDiscrete((5, 2, 3))

"""
mutable struct MultiDiscrete# <: AbstractSpace
    nvec::NTuple{N, UInt32} where N
    dtype::DataType
    shape::Int
end

function MultiDiscrete(nvec::NTuple{N, Int} where N) # nvec: vector of counts of each categorical variable
    @assert all(nvec .> 0) "nvec (counts) have to be positive"
    MultiDiscrete(nvec, UInt32, length(nvec))
end

MultiDiscrete(nvec::Array{Int, 1}) = MultiDiscrete(Tuple(nvec))

sample(multidisc_obj::MultiDiscrete) = [multidisc_obj.dtype(rand(1:counts)) for counts in multidisc_obj.nvec]

contains(multidisc_obj::MultiDiscrete, x) = all(0 .< x .<= multidisc_obj.nvec)

Base.:(==)(multidisc_obj::MultiDiscrete, other::MultiDiscrete) = multidisc_obj.nvec == other.nvec


text:
struct MLDState{DM, DS}
    vehicle_states::Vector{SVector{3,Float64}} # last one is own state
    driver_models::Vector{DM}
    driver_states::Vector{DS}
    terminal::Nullable{Symbol}
end

MLDState(vss, dms, dss) = MLDState(vss, dms, dss, Nullable{Symbol}())

# const MLDPhysicalState = Vector{SVector{3, Float64}}
struct MLDPhysicalState
    vehicle_states::Vector{SVector{3, Float64}}
end

my_idx(s::Union{MLDState, MLDPhysicalState}) = length(s.vehicle_states)

struct MLDPOMDP{DM, DS} <: POMDP{MLDState{DM, DS}, AccelSpeed, MLDPhysicalState}
    reward::MLDRewardModel
    roadway::Highway
    vehicle_model::HighwayVehicle
    dt::Float64
    radius::Float64
    safety::MaxSafeAccel
    dm_dist::IDMGaussianCopula
    _normal_actions::Vector{AccelSpeed}
end

function MLDPOMDP(;reward=MLDRewardModel()
                  roadway::Highway=Highway(4, 4.0, 100.0),
                  vehicle_model::HighwayVehicle=HighwayVehicle(),
                  adjustment_accel::Float64=1.0,
                  dt::Float64=1.0,
                  radius::Float64=50.0,
                  safety=MaxSafeAccel(dt),
                  dm_dist::IDMGaussianCopula=IDMGaussianCopula(0.75)
                 )

    lr_speed = roadway.lanewidth/2.0 # half a lane per second
    accels = (-adjustment_accel, 0.0, adjustment_accel)
    na = vec(collect(AccelSpeed(a, s) for a in accels, s in (-lr_speed, 0.0, lr_speed)))
    DM = NoCrashMLDriver{IDM, Void, MaxSafeAccel}
    C = HighwayDrivingContext{typeof(vehicle_model), SVector{3,Float64}}
    DS = MLDriverState{C, nothing, C}
    return MLDPOMDP{DM, DS}(reward, roadway, vehicle_model, dt, safety, dm_dist, na)
end

isterminal(p::MLPOMDP, s::MLDState) = !isnull(s.terminal)

function generate_sr(p::MLDPOMDP, s::MLDState, a::AccelSpeed, rng::AbstractRNG)
    r = 0.0

    context = HighwayDrivingContext(p.vehicle_model, s.vehicle_states, p.roadway)

    # ego vehicle
    if a.long_accel < -p.reward.brake_thresh
        r -= p.reward.hard_brake_penalty
    end
    my_sp = propagate(p.vehicle_model, s.vehicle_states[end], a, context, p.dt)

    # others
    n_other = length(s.driver_models)
    vss = similar(s.vehicle_states, 0)
    dms = similar(s.driver_models, 0)
    old_is = Int[]
    own_xp = pos(my_sp)[1]
    for i in 1:n_other
        m = s.driver_models[i]
        ds = s.driver_states[i]
        da = sample_action(m, ds, i, AccelSpeed, rng)
        if da.long_accel < -p.reward.brake_thresh
            r -= p.reward.hard_brake_penalty
        end
        vs = s.vehicle_states[i]
        vsp = propagate(p.vehicle_model, vs, da, context, p.dt)
        if abs(pos(vsp)[1]-own_xp) > p.radius
            push!(vss, vsp)
            push!(dms, s.driver_models[i])
            push!(old_is, i)
        end
    end

    # insert ego vehicle
    push!(vss, my_sp)
    holdovers = HighwayDrivingContext(p.vehicle_model, vss, p.roadway)

    all_vss = vss[1:end-1] # all except ego
    newf = push_new_cars_front!(dms, all_vss, holdovers, p.dm_dist, own_xp + p.radius)
    newb = push_new_cars_back!(dms, all_vss, holdovers, p.dm_dist, own_xp - p.radius)

    push!(all_vss, my_sp)

    new_ctx = HighwayDrivingContext(p.vehicle_model, all_vss, p.roadway)::HighwayDrivingContext{typeof(p.vehicle_model), eltype(vss)}

    # see if terminal
    if pos(my_sp)[1] > p.roadway.len
        terminal = Nullable(:reached_end)
    elseif all(occupied_lanes(new_ctx, length(vss)) .== p.reward.target_lane)
        terminal = Nullable(:success)
        r += p.reward.success_reward
    else
        terminal = Nullable{Symbol}()
    end

    # update driver states
    dss = similar(s.driver_states, 0)
    for i in 1:length(vss)-1
        if i <= length(old_is)
            ds = observe(dms[i],
                         s.driver_states[old_is[i]],
                         new_ctx, p.dt)
            push!(dss, ds)
        else
            ds = initial_state(dms[i], new_ctx)
            push!(dss, ds)
        end
    end

    return MLDState(vss, s.driver_models, dss, terminal)::typeof(s), r
end

generate_o(p::MLPOMDP, sp::MLDState, rng::AbstractRNG) = MLDPhysicalState(s.vehicle_states)

function initial_state(p::MLPOMDP, rng::AbstractRNG)
    init = seed_state(p)
    burn_in_steps = 100
    my_model = NORMAL_IDM
    policy = ModelPolicy(p, my_model)
    hr = HistoryRecorder(max_steps=burn_in_steps)
    up = PrimedPreviousObservationUpdater(init)
    hist = simulate(hr, p, policy, up, init, init)
    return state_hist(hist)[end]
end

function seed_state(p::MLPOMDP{DM, DS})
    cars_per_100_m = 3
    xs = linspace(0.0, 2.0*p.radius, cars_per_100_m)
    car_inits = SVector{3, Float64}[]
    for l in 1:p.roadway.n_lanes
        y = (l-1)*p.roadway.lanewidth
        for x in xs
            if l > 1 || abs(x-p.radius) # stay away from ego at (radius,0)
                push!(car_inits, SVector(x, y, 20.0))
            end
        end
    end
    dms = DM[]
    dss = DS[]
    for 1:length(car_inits)
        dm = rand(rng, p.dm_dist)
        push!(dms, dm)
        push!(dss, initial_state(dm))
    end
    push!(car_inits, SVector(p.radius, 0.0, 20.0)) # ego
    return MLDState(vss, dms, dss)
end

# ACTIONS

function actions(p::MLDPOMDP, s::Union{MLDState, MLDPhysicalState})
    na = p._normal_actions
    c = HighwayDrivingContext(p.vehicle_model, s.vehicle_states, p.roadway)
    brake_acc = min(max_safe_acc(p.safety, c, my_idx(s)), -p.reward.brake_thresh/2.0)
    brake = AccelSpeed(brake_acc, 0.0)
    acceptable = [brake]
    for a in na
        # leaving road is always unsafe
        if !is_safe(RoadEdges(p.dt), c, my_idx(s), a)
            continue
        end

        # prevent running into the person in front or to the side
        if is_safe(p.safety, c, my_idx(s), a)
            push!(acceptable, a)
        end
    end
    return acceptable
end

struct ModelPolicy{M} <: Policy
    p::MLDPOMDP
    model::M
end

function action(policy::ModelPolicy, s::Union{MLDState, MLDPhysicalState})
    c = HighwayDrivingContext(policy.p.vehicle_model, s.vehicle_states, policy.p.roadway)
    accel = action(policy.model, c, my_idx(s), LongAccel)
    return AccelSpeed(accel, 0.0)
end

# REWARD
@with_kw struct MLDRewardModel
    brake_thresh::Float64    = 4.0 # always positive
    brake_penalty::Float64   = 1.0 # always positive
    target_lane::Int         =
    success_reward::Float64
end

# ASSUME IDM MODELS for access to sstar
function push_new_cars_front!(models, states, c::Context, dm_dist, front_x::Float64, rng::AbstractRNG)
    # the models MUST correspond to the same indices as the states in c
    for l in 1:c.roadway.n_lanes
        head = lane_head(c, l)
        if head !=0
            head_x = pos(state(c, head))[1]
            delta_x = front_x - head_x
            head_xdot = vel(state(c, head))
            delta_xdot = head_xdot - NORMAL_IDM.v0
        end
        if head == 0 || delta_x >= sstar(NORMAL_IDM, head_xdot, delta_xdot)
            push!(models, rand(rng, p.dm_dist))
            push!(states, SVector(front_x, (l-1)*c.roadway.lanewidth, NORMAL_IDM.v0))
        end
    end
end

function push_new_cars_back!(models, states, c::Context, dm_dist, back_x::Float64)
    for l in 1:c.roadway.n_lanes
        tail = lane_tail(c, l)
        if tail != 0
            tail_x = pos(state(c, tail))[1]
            delta_x = tail_x - back_x
            tail_xdot = vel(state(c, tail))
            delta_xdot = NORMAL_IDM.v0 - tail_xdot
        end
        if tail == 0 || delta_x >= sstar(NORMAL_IDM, NORMAL_IDM.v0, delta_xdot)
            push!(models, rand(rng, p.dm_dist))
            push!(states, SVector(front_x, (l-1)*c.roadway.lanewidth, NORMAL_IDM.v0))
        end
    end

    return models, states
end


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()

module Knet
using Compat

const libknet8 = Libdl.find_library(["libknet8.so"], [dirname(@__FILE__)])

using AutoGrad; export grad, gradloss, gradcheck, getval

include("gpu.jl");              export gpu
include("kptr.jl");             # KnetPtr
include("karray.jl");           export KnetArray
include("unary.jl");            export relu, sigm, invx, logp, htanh
include("broadcast.jl");        # elementwise broadcasting operations
include("reduction.jl");        export logsumexp
include("linalg.jl");           export mat # matmul, axpy!, transpose, (i)permutedims
include("conv.jl");             export conv4, pool, deconv4, unpool
include("update.jl"); 		export Sgd, Momentum, Adam, Adagrad, Adadelta, Rmsprop, update!
include("distributions.jl"); 	export gaussian, xavier, bilinear

"""
    Knet.dir(path...)

Construct a path relative to Knet root.

# Example
```julia
julia> Knet.dir("examples","mnist.jl")
"/home/dyuret/.julia/v0.5/Knet/examples/mnist.jl"
```
"""
dir(path...) = joinpath(dirname(dirname(@__FILE__)),path...)

# See if we have a gpu at initialization:
function __init__()
    try
        r = gpu(true)
        info(r >= 0 ? "Knet using GPU $r" : "No GPU found, Knet using the CPU")
    catch e
        warn("Knet using the CPU: $e")
        gpu(false)
    end
end

end # module


text:
## Ioannis Kosmdiis
## 4 Feb 2021

using ForwardDiff
using ReverseDiff
using Random
using BenchmarkTools

module logit

using Distributions

struct logistic_data
    y::Vector
    x::Array{Float64}
    m::Vector
end

function nobs(data::logistic_data)
    nx = size(data.x)[1]
    ny = length(data.y)
    nm = length(data.m)
    if (nx != ny)
        error("number of rows in of x is not equal to the length of y")
    elseif (nx != nm)
       error("number of rows in of x is not equal to the length of m")
    elseif (ny != nm)
       error("length of y is not equal to the length of m")
    end
    nx
end

function simulate(theta::Vector,
                  x::Matrix{Float64},
                  m::Vector{Int64})
    n = size(x)[2]
    y = rand.(Binomial.(m, cdf.(Logistic(), x * theta)))
    logistic_data(y, x, m);
end

function loglik1(theta::Vector,
                 data::logistic_data,
                 i::Int64)
    eta = sum(data.x[i, :] .* theta)
    mu = exp.(eta) ./ (1 .+ exp.(eta))
    data.y[i] .* log.(mu) + (data.m[i] - data.y[i]) .* log.(1 .- mu)
end

function observation(data::logistic_data, i::Int64)
    vcat(data.y[i], data.m[i], data.x[i,:])
end

function loglik2(theta::AbstractArray,
                 data::AbstractArray)
    y = data[1]
    m = data[2]
    x = data[3:end]
    eta = sum(x .* theta)
    mu = exp.(eta) ./ (1 .+ exp.(eta))
    y .* log.(mu) + (m - y) .* log.(1 .- mu)
end

end

Random.seed!(123);
n = 5000;
p = 50;
x = rand(n, p);
x[:, 1] = fill(1.0, n);
true_betas = rand(p) * 0.02;
my_data = logit.simulate(true_betas, x, fill(1, n));

## Gradient per observation as required by MEstimation
function gR(theta::Vector{T}, obs::Vector{S}) where {T <: Real,S <: Real}
    x = (theta, obs)
    f = (x0, x1) -> logit.loglik2(x0, x1)
    if !(haskey(CACHE, T))
        tape = ReverseDiff.compile(ReverseDiff.GradientTape(f, x))
        CACHE[T] = (tape, (zeros(T, length(theta)), zeros(T, length(obs))))
    end
    tape, y = CACHE[T]
    return ReverseDiff.gradient!(y, tape, x)[1]
end

function gR2(theta::Vector{T}, obs::Vector{S}) where {T <: Real,S <: Real}
    f = eta -> logit.loglik2(eta, obs)
    if !(haskey(CACHE, T))
        tape = ReverseDiff.compile(ReverseDiff.GradientTape(f, theta))
        CACHE[T] = (tape, zeros(T, length(theta)))
    end
    tape, y = CACHE[T]
    return ReverseDiff.gradient!(y, tape, theta)
end

## Same thing using forward diff
function gF(theta::Vector, data::logit.logistic_data, i::Int64)
    obs = logit.observation(data, i)
    ForwardDiff.gradient(th -> logit.loglik2(th, obs), theta)
end

const CACHE = Dict{DataType,Any}()

function runR(pars, j::Int64)
    ForwardDiff.jacobian(theta -> gR(theta, logit.observation(my_data, j)), pars)
end

function runR2(pars, j::Int64)
    ForwardDiff.jacobian(theta -> gR2(theta, logit.observation(my_data, j)), pars)
end

function runF(pars, j::Int64)
    ForwardDiff.jacobian(theta -> gF(theta, my_data, j), pars)
end

function runFhessian(pars, j::Int64)
    ForwardDiff.hessian(theta -> logit.loglik1(theta, my_data, j), pars)
end

CACHE
runR(true_betas, 1)
runR2(true_betas, 1)
CACHE
runF(true_betas, 1)
runFhessian(true_betas, 1)

## The log-likelihood hessian across all observations
ra = 1:logit.nobs(my_data)
a0 = @benchmark sum(map(j -> runR2(true_betas, j), ra))
a1 = @benchmark sum(map(j -> runR(true_betas, j), ra))
a2 = @benchmark sum(map(j -> runF(true_betas, j), ra))
a3 = @benchmark sum(map(j -> runFhessian(true_betas, j), ra))

## Forward is 15-25 times slower than reverse implementation above
ratio(mean(a2), mean(a1))
ratio(mean(a3), mean(a1))
# julia> ratio(mean(a2), mean(a1))
# BenchmarkTools.TrialRatio:
#   time:             16.721108277780537
#   gctime:           29.415474243662555
#   memory:           26.19886681971946
#   allocs:           1.0385193063358555

# julia> ratio(mean(a3), mean(a1))
# BenchmarkTools.TrialRatio:
#   time:             26.0195816865898
#   gctime:           29.945178720705652
#   memory:           32.704893760988654
#   allocs:           3.469006754781931

## Use in optimization
using NLsolve

function fR(theta::Vector)
    gr = zeros(length(theta))
    for j in 1:n
        obs = logit.observation(my_data, j)
        gr += gR(theta, obs)
    end
    gr
end

function fF(theta::Vector)
    gr = zeros(length(theta))
    for j in 1:n
        gr += gF(theta, my_data, j)
    end
    gr
end

a1 = @benchmark nlsolve(fR, true_betas)
a2 = @benchmark nlsolve(fF, true_betas)

## Optimization with forward is 3x slower and reuires 10 times the memory
ratio(mean(a2), mean(a1))
# julia> ratio(mean(a2), mean(a1))
# BenchmarkTools.TrialRatio:
#   time:             3.0670491962167503
#   gctime:           14.19187050750875
#   memory:           10.935421740548241
#   allocs:           1.1018641415575794


text:
"This script contains the implementation algorithms for the unitary transformation on discrete and finite fields. It acts as a single transformation, but it can be adecuate to be extended on a FOR loop to be applied on polychromatic fields

It provides three functions: `fnm(n, m, j, imgdata)` maps the field from position to modes. `fnmθ(θ, n, m, j, imgdat)` adds an angle dependent phase. `Fqθ(θ, q₁, q₂, j, imgdat)` returns the transformed field in the modes to the position representation. This last function is the upper layer of the transformation, it can (could, should) be used to the polychromatic duties."

using PixelsGT # Call the parent module

export fnm, fnmθ, Fqθ

"`fnm(n, m, j, imgdat)` encodes a two-dimensional position field in a modes space labeled by mode _n_ and counter _m_, who ranges are n ∈ [0, 2j] and |m| ≤ n. 'imgdat' is a position field, that is a N×N matrix, for N = 2j."
function fnm(n, m, j, imgdat)
    vnm = zeros(ComplexF64, (2 * j + 1, 2 * j + 1))
    for q₁ in -j:j
        for q₂ in -j:j
            vnm[q₁ + j + 1, q₂ + j + 1] = LK(q₁, q₂, n, m, j) * imgdat[q₁ + j + 1, q₂ + j + 1]
        end
    end
    return sum(vnm)
end

"`fnmθ(θ, n, m, j, imgdat)` is a function who applies an angle-dependent phase to the `fnm` function. It is the preamble to perform the transformation, given the value of θ."
fnmθ(θ, n, m, j, imgdat) = exp(-1im * m * θ) * fnm(n, m, j, imgdat)

"`Fqθ(θ, q₁, q₂, j, imgdat)` is the top function who perform the whole transformation on the position field. Given and N×N field, an angle θ, the values of the transformed pair (q₁, q₂) are given by this function. Note: It is extremely slow, so slow."
function Fqθ(θ, q₁, q₂, j, imgdat)
    rd = sum([conj(LK(q₁, q₂, n, m, j)) * fnmθ(θ, n, m, j, imgdat) for n in 0:2*j for m in -n:2:n])
    ru = sum([conj(LK(q₁, q₂, n, m, j)) * fnmθ(θ, n, m, j, imgdat) for n in (2*j + 1):4*j for m in -(4 * j - n):2:(4 * j - n)])
    return rd + ru
end