text:
"""
    cleanup(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)

Standardizes a molecule.

```julia
cleanup(mol)
```
"""
function cleanup(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)
    details_json::String = jsonify_details(details)
    ccall((:cleanup, librdkitcffi), Cvoid, (Ref{Cstring}, Ref{Csize_t}, Cstring), mol.mol, mol.mol_size, details_json)
end

"""
    normalize(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)

Applies a series of standard transformations to correct functional groups and recombine charges.

```julia
normalize(mol)
```
"""
function normalize(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)
    details_json::String = jsonify_details(details)
    ccall((:normalize, librdkitcffi), Cvoid, (Ref{Cstring}, Ref{Csize_t}, Cstring), mol.mol, mol.mol_size, details_json)
end

"""
    canonical_tautomer(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)

Returns the canonical tautomer for a molecule.

```julia
canonical_tautomer(mol)
```
"""
function canonical_tautomer(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)
    details_json::String = jsonify_details(details)
    ccall((:canonical_tautomer, librdkitcffi), Cvoid, (Ref{Cstring}, Ref{Csize_t}, Cstring), mol.mol, mol.mol_size, details_json)
end

"""
    charge_parent(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)

Returns the uncharged version of the largest fragment.

```julia
charge_parent(mol)
```
"""
function charge_parent(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)
    details_json::String = jsonify_details(details)
    ccall((:charge_parent, librdkitcffi), Cvoid, (Ref{Cstring}, Ref{Csize_t}, Cstring), mol.mol, mol.mol_size, details_json)
end

"""
    reionize(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)

Ensures the strongest acid groups are charged first.

```julia
reionize(mol)
```
"""
function reionize(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)
    details_json::String = jsonify_details(details)
    ccall((:reionize, librdkitcffi), Cvoid, (Ref{Cstring}, Ref{Csize_t}, Cstring), mol.mol, mol.mol_size, details_json)
end

"""
    neutralize(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)

Uncharges molecules by adding and/or removing hydrogens.

```julia
neutralize(mol)
```
"""
function neutralize(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)
    details_json::String = jsonify_details(details)
    ccall((:neutralize, librdkitcffi), Cvoid, (Ref{Cstring}, Ref{Csize_t}, Cstring), mol.mol, mol.mol_size, details_json)
end

"""
    fragment_parent(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)

Returns the largest fragment after doing a cleanup.

```julia
fragment_parent(mol)
```
"""
function fragment_parent(mol::Mol, details::Union{Dict{String,Any},Nothing}=nothing)
    details_json::String = jsonify_details(details)
    ccall((:fragment_parent, librdkitcffi), Cvoid, (Ref{Cstring}, Ref{Csize_t}, Cstring), mol.mol, mol.mol_size, details_json)
end


text:
@testset "Float64" begin
    B = IncrementBinner(blocksize=4)

    @test B.keep == 4

    for i in 1:4
        push!(B, 1.0)
    end
    for i in 1:4
        push!(B, 1.0)
        push!(B, 3.0)
    end
    for i in 1:4
        push!(B, 3.0)
        push!(B, 4.0)
        push!(B, 4.0)
        push!(B, 5.0)
    end

    xs = indices(B)
    ys = values(B)

    @test xs == Float64[1,2,3,4, 5.5,7.5,9.5,11.5, 14.5,18.5,22.5,26.5]
    @test ys == Float64[1,1,1,1, 2,2,2,2, 4,4,4,4]
end

@testset "ComplexF64" begin
    B = IncrementBinner(ComplexF64, blocksize=3)

    @test B.keep == 3

    for i in 1:3
        push!(B, 1.0+1.0im)
    end
    for i in 1:3
        push!(B, 1.0+3.0im)
        push!(B, 3.0+1.0im)
    end

    xs = indices(B)
    ys = values(B)

    @test xs == Float64[1,2,3, 4.5,6.5,8.5]
    @test ys == ComplexF64[1+1im,1+1im,1+1im, 2+2im,2+2im,2+2im]
end

@testset "Vector" begin
    B = IncrementBinner([0.0, 0.0], blocksize=2)

    @test B.keep == 2

    for i in 1:2
        push!(B, [1.0, 1.0])
    end
    for i in 1:2
        push!(B, [1.0, 3.0])
        push!(B, [3.0, 1.0])
    end

    xs = indices(B)
    ys = values(B)

    @test xs == Float64[1,2, 3.5,5.5]
    @test ys == [[1.0, 1.0], [1.0, 1.0], [2.0, 2.0], [2.0, 2.0]]
end

ext: jl

max_stars_repo_name: crstnbr/BinningAnalysis


text:
@defcomp n2ocycle begin
    e_globalN2Oemissions = Parameter(index=[time], unit="Mtonne/year")
    e_0globalN2Oemissions = Parameter(unit="Mtonne/year", default=11.22281348065454)
    c_N2Oconcentration = Variable(index=[time], unit="ppbv")
    pic_preindustconcN2O = Parameter(unit="ppbv", default=270.)
    exc_excessconcN2O = Variable(unit="ppbv")
    c0_N2Oconcbaseyr = Parameter(unit="ppbv", default=324.)
    re_remainN2O = Variable(index=[time], unit="Mtonne")
    re_remainN2Obase = Variable(unit="Mtonne")
    nte_natN2Oemissions = Variable(index=[time], unit="Mtonne/year")
    air_N2Ofractioninatm = Parameter(unit="%", default=100.)
    tea_N2Oemissionstoatm = Variable(index=[time], unit="Mtonne/year")
    teay_N2Oemissionstoatm = Variable(index=[time], unit="Mtonne/t")
    y_year = Parameter(index=[time], unit="year")
    y_year_0 = Parameter(unit="year")
    res_N2Oatmlifetime = Parameter(unit="year", default=114.)
    den_N2Odensity = Parameter(unit="Mtonne/ppbv", default=7.8)
    stim_N2Oemissionfeedback = Parameter(unit="Mtonne/degreeC", default=0.)
    rtl_g0_baselandtemp = Parameter(unit="degreeC", default=1.1683981941310047)
    rtl_g_landtemperature = Parameter(index=[time], unit="degreeC")

    function run_timestep(p, v, d, t)
        # note that Hope (2009) states that Equations 1-12 for methane also apply to N2O

        if is_first(t)
            # eq.3 from Hope (2006) - natural emissions feedback, using global temperatures calculated in ClimateTemperature component
            nte_0 = p.stim_N2Oemissionfeedback * p.rtl_g0_baselandtemp
            v.nte_natN2Oemissions[t] = p.stim_N2Oemissionfeedback * p.rtl_g0_baselandtemp
            # eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
            v.tea_N2Oemissionstoatm[t] = (p.e_globalN2Oemissions[t] + v.nte_natN2Oemissions[t]) * p.air_N2Ofractioninatm / 100
            tea_0 = (p.e_0globalN2Oemissions + nte_0) * p.air_N2Ofractioninatm / 100
            v.teay_N2Oemissionstoatm[t] = (v.tea_N2Oemissionstoatm[t] + tea_0) * (p.y_year[t] - p.y_year_0) / 2
            # adapted from eq.1 in Hope(2006) - calculate excess concentration in base year
            v.exc_excessconcN2O = p.c0_N2Oconcbaseyr - p.pic_preindustconcN2O
            # Eq. 2 from Hope (2006) - base-year remaining emissions
            v.re_remainN2Obase = v.exc_excessconcN2O * p.den_N2Odensity
            v.re_remainN2O[t] = v.re_remainN2Obase * exp(-(p.y_year[t] - p.y_year_0) / p.res_N2Oatmlifetime) +
                v.teay_N2Oemissionstoatm[t] * p.res_N2Oatmlifetime * (1 - exp(-(p.y_year[t] - p.y_year_0) / p.res_N2Oatmlifetime)) / (p.y_year[t] - p.y_year_0)
        else
            # eq.3 from Hope (2006) - natural emissions (carbon cycle) feedback, using global temperatures calculated in ClimateTemperature component
            # Here assume still using area-weighted average regional temperatures (i.e. land temperatures) for natural emissions feedback
            v.nte_natN2Oemissions[t] = p.stim_N2Oemissionfeedback * p.rtl_g_landtemperature[t - 1]
            # eq.6 from Hope (2006) - emissions to atmosphere depend on the sum of natural and anthropogenic emissions
            v.tea_N2Oemissionstoatm[t] = (p.e_globalN2Oemissions[t] + v.nte_natN2Oemissions[t]) * p.air_N2Ofractioninatm / 100
            # eq.7 from Hope (2006) - average emissions to atm over time period
            v.teay_N2Oemissionstoatm[t] = (v.tea_N2Oemissionstoatm[t] + v.tea_N2Oemissionstoatm[t - 1]) * (p.y_year[t] - p.y_year[t - 1]) / 2
            # eq.10 from Hope (2006) - remaining emissions in atmosphere
            v.re_remainN2O[t] = v.re_remainN2O[t - 1] * exp(-(p.y_year[t] - p.y_year[t - 1]) / p.res_N2Oatmlifetime) +
                v.teay_N2Oemissionstoatm[t] * p.res_N2Oatmlifetime * (1 - exp(-(p.y_year[t] - p.y_year[t - 1]) / p.res_N2Oatmlifetime)) / (p.y_year[t] - p.y_year[t - 1])
        end

        # eq.11 from Hope(2006) - N2O concentration
        v.c_N2Oconcentration[t] = p.pic_preindustconcN2O + v.exc_excessconcN2O * v.re_remainN2O[t] / v.re_remainN2Obase

    end
end


text:
struct PointTrajectory{T}
    frame::CartesianFrame3D
    trajectory::T
end

function (traj::PointTrajectory)(x, ::Val{num_derivs}) where num_derivs
    vals = traj.trajectory(x, Val(num_derivs))
    frame = traj.frame
    Point3D(frame, vals[1]), map(x -> FreeVector3D(frame, x), Base.tail(vals))...
end

(traj::PointTrajectory)(x) = Point3D(traj.frame, traj.trajectory(x))

struct FreeVectorTrajectory{T}
    frame::CartesianFrame3D
    trajectory::T
end

function (traj::FreeVectorTrajectory)(x, ::Val{num_derivs}) where num_derivs
    vals = traj.trajectory(x, Val(num_derivs))
    frame = traj.frame
    map(x -> FreeVector3D(frame, x), vals)
end

(traj::FreeVectorTrajectory)(x) = FreeVector3D(traj.frame, traj.trajectory(x))


text:
###################
# Packages ########
###################

using CSV
using DataFrames
using JuMP
using Gurobi
using ProgressMeter
using ElectronDisplay

# Loading the project module, containing essential functions and structs
include("model/src/Sesam.jl")
using .Sesam

###################
# SETS, PARAMS ####
###################

# time T
# T = collect(1)
time = collect(1)
data = "ieee_case5"

# Read nodes and load from .csv
nodes_df = CSV.read(string("data/", data, "/nodes.csv"))
load_df = CSV.read(string("data/", data, "/load.csv"))

# Create set of nodes
nodes = Nodes(nodes_df, load_df)

# Read power plants from .csv
powerplants_df = CSV.read(string("data/", data, "/powerplants.csv"))

# Create set of power plants
pp = PowerPlants(powerplants_df)

# Read renewables and availabilities from .csv
renewables_df = CSV.read(string("data/", data, "/renewables.csv"))
avail_solar_df = CSV.read(string("data/", data, "/avail_solar.csv"))
avail_wind_df = CSV.read(string("data/", data, "/avail_wind.csv"))

# Create set of renewables
res = Renewables(renewables_df, avail_solar_df, avail_wind_df)

# Read lines from .csv
lines_df = CSV.read(string("data/", data, "/lines.csv"))
lines = Lines(lines_df)

###################
# Run DCOPF #######
###################

DCOPF_mod, P_opt, Θ_opt, P_inj_opt, P_flow_opt =
	DCOPF(time, nodes, lines, pp)

# Optimisation status
JuMP.termination_status(DCOPF_mod)
JuMP.primal_status(DCOPF_mod)
JuMP.dual_status(DCOPF_mod)

# Results
JuMP.objective_value(DCOPF_mod)
P_opt
Θ_opt
P_inj_opt
P_flow_opt

###################
# Run SCOPF #######
###################

SCOPF_mod, P_opt_s, Θ_opt_s, P_inj_opt_s, P_flow_opt_s =
	SCOPF(time, nodes, lines, pp)

# Optimisation status
JuMP.termination_status(SCOPF_mod)
JuMP.primal_status(SCOPF_mod)
JuMP.dual_status(SCOPF_mod)

# Results
JuMP.objective_value(SCOPF_mod)
P_opt_s
Θ_opt_s
P_inj_opt_s
P_flow_opt_s


text:
#script that estimates the dir bin gas network model on emid data and evaluates GAS
# forecasting performances

using Utilities,AReg,StaticNets,JLD,MLBase,StatsBase,DynNets
using PyCall; pygui(:qt); using PyPlot

## Load dataj
using JLD2
 halfPeriod = false
 fold_Path =  "/home/Domenico/Dropbox/Dynamic_Networks/data/emid_data/juliaFiles/"
 loadFilePartialName = "Weekly_eMid_Data_from_"
 halfPeriod? periodEndStr =  "2012_03_12_to_2015_02_27.jld": periodEndStr =  "2009_06_22_to_2015_02_27.jld"
 @load(fold_Path*loadFilePartialName*periodEndStr, AeMidWeekly_T,banksIDs,inactiveBanks,
                YeMidWeekly_T,weekInd,datesONeMid ,degsIO_T,strIO_T)

#

#define sequence of matrices and degrees
matY_T = YeMidWeekly_T[:,:,3:end]
 matA_T = matY_T.>0
    N = length(matA_T[1,:,1])
    T = size(matA_T)[3]
    Ttrain = 100#round(Int, T/2) #70 106 #
    threshVar = 0.00#5
    degsIO_T = [sumSq(matA_T,2);sumSq(matA_T,1)]
 allFitSS =  StaticNets.estimate( StaticNets.SnapSeqNetDirBin1(degsIO_T); identPost = false,identIter= true )
 sizeLab = 30

 @load(fold_Path*"_rolling_SD_estimates_"*periodEndStr, gasFiltAndForeFitFromRollEst,allFitSS,gasEst_rolling)

initWorkers()
 using DynNets , StatsFuns
 @everywhere using DynNets, StatsFuns

#select links among largest banks in training sample
remMat =squeeze(prod(.!matA_T,3),3)#
  modAllObs = DynNets.SdErgmDirBin1(degsIO_T,"FISHER-DIAG")
  sizeLeg = 15
  legTex = []
  valsNsample = [5 10 20 30] #Vector(5:10:106)
  N_Nsample = length(valsNsample)
  maxNsteps = 25
  valsNsteps = Vector(1:maxNsteps) #[Vector(1:5); round(Int,linspace(6,maxNsteps,N_Nsteps-5))]
  N_Nsteps = length(valsNsteps)
  storeMat = false
  if storeMat
  expMat_T_means = SharedArray{Float64}(zeros(N,N,T,N_Nsteps,N_Nsample))
  expMat_T_Y = SharedArray{Float64}(zeros(N,N,T,N_Nsteps,N_Nsample))
  expMat_T_meanPar = SharedArray{Float64}(zeros(N,N,T,N_Nsteps,N_Nsample))
  else
      expMat_T_means = SharedArray{Float64}(zeros(N,N,T))
      expMat_T_Y = SharedArray{Float64}(zeros(N,N,T ))
      expMat_T_meanPar = SharedArray{Float64}(zeros(N,N,T))
  end
  aucGas = SharedArray{Float64}(zeros(N_Nsteps,N_Nsample,3))
 for n=1:N_Nsample
            @sync @parallel  for k=1:N_Nsteps
              Nsteps = valsNsteps[k]
              Nsample = valsNsample[n]
              @show((Nsteps,Nsample))

             if storeMat
             @time expMat_T_means[:,:,:,k,n],forFit ,expMat_T_Y[:,:,:,k,n],expMat_T_meanPar[:,:,:,k,n], =
                   DynNets.multiStepsForecastExpMat(modAllObs,matA_T,gasEst_rolling[1],Nsample,Nsteps,Ttrain)
               else
                   @time expMat_T_means,forFit ,expMat_T_Y,expMat_T_meanPar =
                             DynNets.multiStepsForecastExpMat(modAllObs,matA_T,gasEst_rolling[1],Nsample,Nsteps,Ttrain)
             end
             #figure(11)
              tmpTm1 =  StaticNets.nowCastEvalFitNet(   zeros(2,2) ,matA_T[:,:,Ttrain+1:end];mat2rem = remMat,expMat_T=expMat_T_means[:,:,Ttrain+1:end],shift=0,plotFlag = false)

              aucGas[k,n,1] = tmpTm1[3]
              #figure(12)
              tmpTm1 =  StaticNets.nowCastEvalFitNet(   zeros(2,2) ,matA_T[:,:,Ttrain+1:end];mat2rem = remMat,expMat_T=expMat_T_Y[:,:,Ttrain+1:end],shift=0,plotFlag = false)

              aucGas[k,n,2] = tmpTm1[3]
             # figure(13)

              tmpTm1 =  StaticNets.nowCastEvalFitNet(   zeros(2,2) ,matA_T[:,:,Ttrain+1:end];mat2rem = remMat,expMat_T=expMat_T_meanPar[:,:,Ttrain+1:end],shift=0,plotFlag = false )
              aucGas[k,n,3] = tmpTm1[3]

        end
    end
 @sync @parallel for k=1:5
 end

 #close()
 #plot(valsNsteps,aucGas[:,1,1])

aucAR =  zeros(maxNsteps)
 ## multistep ahead forecast from AR1
 for Nsteps =1:maxNsteps
 ARForeFitFromRollEst = zeros(2N,T)
 foreFitAR_roll = zeros(2N,T)
 for t=Ttrain+1:T
    tEst = t-Ttrain
    tObs = t-Nsteps
    trainFit = allFitSS[:,tEst:tEst + Ttrain-1]
    Nconst,isConIn,isConOut = StaticNets.defineConstDegs(degsIO_T[:,tEst:tEst + Ttrain],thVarConst =0.005 )
    estAR1= StaticNets.estManyAR(trainFit;isCon = [isConIn;isConOut])
    ## Nsteps  ahead forecasts for each fitness for each time in Test sample
    for n=1:2N
        #compute forecasting for all links, diagonal will be disregarded after
        foreFitAR_roll[n,t] = StaticNets.oneStepForAR1(allFitSS[n,tObs],estAR1[n,:],steps = Nsteps)
    end
 end
  inds2forecast = (Ttrain+2:T)
  @show(Nsteps)
 tmpTm1 =  StaticNets.nowCastEvalFitNet(  foreFitAR_roll[:,inds2forecast] ,matA_T[:,:,inds2forecast];mat2rem = remMat,shift = 0,plotFlag=false)
 aucAR[Nsteps] = tmpTm1[3]
 end

close()#
  remMat =squeeze(prod(.!matA_T,3),3)#
  sizeLeg = 15
  maxNsteps = 50
  legTex = []
  aucStat = zeros(maxNsteps)
    for n=1:maxNsteps
          shiftVal = n
          tmpTm1 = StaticNets.nowCastEvalFitNet(   allFitSS[:,Ttrain+1:end],matA_T[:,:,Ttrain+1:end];mat2rem = remMat,shift = shiftVal,plotFlag = false)
          aucStat[n] = tmpTm1[3]
          legTex = [legTex; "t-$(shiftVal) Par  AUC = $(round(tmpTm1[3],3))"]
    end

save_fold = "./data/multiStepForecast/"
@save(save_fold* "aucForVariousNsampleAndNsteps_gas_estimates.jld",
              aucGas, valsNsteps,valsNsample,aucAR,aucStat)

  #
  # legend(legTex,fontsize = sizeLeg)
  # title("T train = $(Ttrain) of $(T)  ",size = sizeLab)
  #
  # xlabel("False Positive Rate",size = sizeLab)
  # ylabel("False Negative Rate",size = sizeLab)

#
if false
 close()
    maxNsteps = 5
    plot(1:maxNsteps,aucStat,"r")
    legTex = "ERGM"
    plot(1:maxNsteps,aucAR)
     legTex = [legTex; "AR1" ]
    for indSam=1#:N_Nsample
    plot(valsNsteps,aucGas[:,indSam,1])
     legTex = [legTex; "SDERGM avg mean mat N sample = $(valsNsample[indSam])" ]
    end
    legend(legTex,fontsize = sizeLeg)
    title("Multi Steps Ahead Forecast for ",size = sizeLab)
    xlabel("Number Of Steps Ahead",size = sizeLab)
    ylabel("Area Under The Curve",size = sizeLab)
    grid()
end

#
# function forecastEvalGasNetDirBin1(obsNet_T::BitArray{3}, gasParEstOnTrain::Array{Array{Float64,1},1},Ttrain::Int,Nsteps::Int;Nsample::Int=100, thVarConst = 0.0005,mat2rem=falses(obsNet_T[:,:,1]) )
#     # this function evaluates the forecasting performances of the GAS model over the train sample.
#     # to do so uses gas parameters previously estimated and observations at time t of the test sample
#     # to forecast tv par at t+1, then using these estimates and t+1 observations goes to t+2 and so on
#
#     @show N2 = length(gasParEstOnTrain[1]);N = round(Int,N2/2)
#     T = length(obsNet_T[1,1,:])
#     @show Ttest = T-Ttrain
#     testObsNet_T = obsNet_T[:,:,Ttrain+1:end]
#     degsIO_T = [sumSq(obsNet_T,2); sumSq(obsNet_T,1)]
#
#
#
#     #disregard predictions for degrees that are constant in the training sample
#     # gas shoud manage but AR fitness no and I dont want to advantage gas
#     # to do so i need to count the number of links that are non constant (nnc)
#     Nc,isConIn,isConOut = StaticNets.defineConstDegs(degsIO_T[:,1:Ttrain];thVarConst = thVarConst )
#     noDiagIndnnc = putZeroDiag(((.!isConIn).*(.!isConOut')).* (.!mat2rem)    )
#
#     @show Nlinksnnc =sum(noDiagIndnnc)
#     # forecast fitnesses using Gas parameters and observations
#     foreFit,~ = score_driven_filter_or_dgp( DynNets.SdErgmDirBin1(degsIO_T),[gasParEstOnTrain[1];gasParEstOnTrain[2];gasParEstOnTrain[3]])
#
#     TRoc = Ttest-1
#     #storage variables
#     foreVals = 100 * ones(Nsteps*Nlinksnnc*(TRoc))
#     realVals = trues( Nsteps*Nlinksnnc*(TRoc))
#
#     lastInd=1
#
#     for t=(1+Nsteps):TRoc
#         adjMat = testObsNet_T[:,:,t]   #[testObsNet_T[50:end,:,t];testObsNet_T[1:49,:,t]] #a shuffle to test that the code does not not work
#         adjMat_tm1 = testObsNet_T[:,:,t-1]
#         indZeroR = sum(adjMat_tm1,2).==0
#         indZeroC = sum(adjMat_tm1,1).==0
#         indZeroMat = indZeroR.*indZeroC
#     #    println(sum(indZeroMat)/(length(indZeroC)^2))
#
#         expMat = StaticNets.expMatrix(StaticNets.fooErgmDirBin1,foreFit[:,Ttrain+1:end][:,t])
#
#          tmpAllMat,~ = multiSteps_predict_score_driven_par(modGasDirBin1_eMidTrain,N,foreFit[:,Ttrain+1:end][:,t-Nsteps],
#             gasParEstOnTrain[1],gasParEstOnTrain[2],gasParEstOnTrain[3],Nsample,Nsteps)
#
#          expMat=tmpAllMat[:,:,Nsteps+1]
#
#
#         #expMat[indZeroMat] = 0
#
#         foreVals[lastInd:lastInd+Nlinksnnc-1] = expMat[noDiagIndnnc]
#         realVals[lastInd:lastInd+Nlinksnnc-1] = adjMat[noDiagIndnnc]
#         lastInd += Nlinksnnc
#     end
#     return realVals,foreVals,foreFit
#  end

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
module AverageShiftedHistograms

import UnicodePlots, RecipesBase
using LinearAlgebra, Statistics
import StatsBase: nobs, fweights
export ash, ash!, extendrange, xy, xyz, nout, nobs, Kernels

#-----------------------------------------------------------------------# common

"""
`extendrange(x, s = .5, n = 200)`

Create a range of length `n` starting at `s` standard deviations below
`minimum(x)` and ending at `s` standard deviations above `maximum(x)`
"""
function extendrange(y::AbstractVector, s = 0.5, n = 500)
    σ = std(y)
    range(minimum(y) - s * σ, stop = maximum(y) + s * σ, length = n)
end

include("kernels.jl")
include("univariate.jl")
include("bivariate.jl")
end


text:
export bytes_to_samples

"""
Return bytes to samples dict
"""
function bytes_to_samples()

    Dict{String, Int32}(
        "TraceNumWithinLine"              =>  1  -1,
        "TraceNumWithinFile"              =>  5  -1,
        "FieldRecord"                     =>  9  -1,
        "TraceNumber"                     =>  13 -1,
        "EnergySourcePoint"               =>  17 -1,
        "CDP"                             =>  21 -1,
        "CDPTrace"                        =>  25 -1,
        "TraceIDCode"                     =>  29 -1,
        "NSummedTraces"                   =>  31 -1,
        "NStackedTraces"                  =>  33 -1,
        "DataUse"                         =>  35 -1,
        "Offset"                          =>  37 -1,
        "RecGroupElevation"               =>  41 -1,
        "SourceSurfaceElevation"          =>  45 -1,
        "SourceDepth"                     =>  49 -1,
        "RecDatumElevation"               =>  53 -1,
        "SourceDatumElevation"            =>  57 -1,
        "SourceWaterDepth"                =>  61 -1,
        "GroupWaterDepth"                 =>  65 -1,
        "ElevationScalar"                 =>  69 -1,
        "RecSourceScalar"                 =>  71 -1,
        "SourceX"                         =>  73 -1,
        "SourceY"                         =>  77 -1,
        "GroupX"                          =>  81 -1,
        "GroupY"                          =>  85 -1,
        "CoordUnits"                      =>  89 -1,
        "WeatheringVelocity"              =>  91 -1,
        "SubWeatheringVelocity"           =>  93 -1,
        "UpholeTimeSource"                =>  95 -1,
        "UpholeTimeGroup"                 =>  97 -1,
        "StaticCorrectionSource"          =>  99 -1,
        "StaticCorrectionGroup"           =>  101-1,
        "TotalStaticApplied"              =>  103-1,
        "LagTimeA"                        =>  105-1,
        "LagTimeB"                        =>  107-1,
        "DelayRecordingTime"              =>  109-1,
        "MuteTimeStart"                   =>  111-1,
        "MuteTimeEnd"                     =>  113-1,
        "ns"                              =>  115-1,
        "dt"                              =>  117-1,
        "GainType"                        =>  119-1,
        "InstrumentGainConstant"          =>  121-1,
        "InstrumntInitialGain"            =>  123-1,
        "Correlated"                      =>  125-1,
        "SweepFrequencyStart"             =>  127-1,
        "SweepFrequencyEnd"               =>  129-1,
        "SweepLength"                     =>  131-1,
        "SweepType"                       =>  133-1,
        "SweepTraceTaperLengthStart"      =>  135-1,
        "SweepTraceTaperLengthEnd"        =>  137-1,
        "TaperType"                       =>  139-1,
        "AliasFilterFrequency"            =>  141-1,
        "AliasFilterSlope"                =>  143-1,
        "NotchFilterFrequency"            =>  145-1,
        "NotchFilterSlope"                =>  147-1,
        "LowCutFrequency"                 =>  149-1,
        "HighCutFrequency"                =>  151-1,
        "LowCutSlope"                     =>  153-1,
        "HighCutSlope"                    =>  155-1,
        "Year"                            =>  157-1,
        "DayOfYear"                       =>  159-1,
        "HourOfDay"                       =>  161-1,
        "MinuteOfHour"                    =>  163-1,
        "SecondOfMinute"                  =>  165-1,
        "TimeCode"                        =>  167-1,
        "TraceWeightingFactor"            =>  169-1,
        "GeophoneGroupNumberRoll"         =>  171-1,
        "GeophoneGroupNumberTraceStart"   =>  173-1,
        "GeophoneGroupNumberTraceEnd"     =>  175-1,
        "GapSize"                         =>  177-1,
        "OverTravel"                      =>  179-1,
        "CDPX"                            =>  181-1,
        "CDPY"                            =>  185-1,
        "Inline3D"                        =>  189-1,
        "Crossline3D"                     =>  193-1,
        "ShotPoint"                       =>  197-1,
        "ShotPointScalar"                 =>  201-1,
        "TraceValueMeasurmentUnit"        =>  203-1,
        "TransductionConstnatMantissa"    =>  205-1,
        "TransductionConstantPower"       =>  209-1,
        "TransductionUnit"                =>  211-1,
        "TraceIdentifier"                 =>  213-1,
        "ScalarTraceHeader"               =>  215-1,
        "SourceType"                      =>  217-1,
        "SourceEnergyDirectionMantissa"   =>  219-1,
        "SourceEnergyDirectionExponent"   =>  223-1,
        "SourceMeasurmentMantissa"        =>  225-1,
        "SourceMeasurementExponent"       =>  229-1,
        "SourceMeasurmentUnit"            =>  231-1,
        "Unassigned1"                     =>  233-1,
        "Unassigned2"                     =>  237-1)

end


text:
module TwoDNavierStokes

export
  Problem,
  set_zeta!,
  updatevars!,

  energy,
  energy_dissipation,
  energy_work,
  energy_drag,
  enstrophy,
  enstrophy_dissipation,
  enstrophy_work,
  enstrophy_drag

  using
    CUDA,
    Reexport

@reexport using FourierFlows

using LinearAlgebra: mul!, ldiv!
using FourierFlows: parsevalsum

nothingfunction(args...) = nothing

"""
    Problem(; parameters...)

Construct a 2D turbulence problem.
"""
function Problem(dev::Device=CPU();
  # Numerical parameters
          nx = 256,
          Lx = 2π,
          ny = nx,
          Ly = Lx,
          dt = 0.01,
  # Drag and/or hyper-/hypo-viscosity
           ν = 0,
          nν = 1,
           μ = 0,
          nμ = 0,
  # Timestepper and equation options
     stepper = "RK4",
       calcF = nothingfunction,
  stochastic = false,
           T = Float64)

  grid = TwoDGrid(dev, nx, Lx, ny, Ly; T=T)

  params = Params{T}(ν, nν, μ, nμ, calcF)

  vars = calcF == nothingfunction ? Vars(dev, grid) : (stochastic ? StochasticForcedVars(dev, grid) : ForcedVars(dev, grid))

  equation = Equation(params, grid)

  return FourierFlows.Problem(equation, stepper, dt, grid, vars, params, dev)
end

# ----------
# Parameters
# ----------

"""
    Params(ν, nν, μ, nμ, calcF!)

Returns the params for two-dimensional turbulence.
"""
struct Params{T} <: AbstractParams
       ν :: T         # Vorticity viscosity
      nν :: Int       # Vorticity hyperviscous order
       μ :: T         # Bottom drag or hypoviscosity
      nμ :: Int       # Order of hypodrag
  calcF! :: Function  # Function that calculates the forcing F
end

Params(ν, nν) = Params(ν, nν, typeof(ν)(0), 0, nothingfunction)

# ---------
# Equations
# ---------

"""
    Equation(params, grid)

Returns the equation for two-dimensional turbulence with `params` and `grid`.
"""
function Equation(params::Params, grid::AbstractGrid)
  L = @. - params.ν * grid.Krsq^params.nν - params.μ * grid.Krsq^params.nμ
  CUDA.@allowscalar L[1, 1] = 0
  return FourierFlows.Equation(L, calcN!, grid)
end

# ----
# Vars
# ----

abstract type TwoDNavierStokesVars <: AbstractVars end

struct Vars{Aphys, Atrans, F, P} <: TwoDNavierStokesVars
     zeta :: Aphys
        u :: Aphys
        v :: Aphys
    zetah :: Atrans
       uh :: Atrans
       vh :: Atrans
       Fh :: F
  prevsol :: P
end

const ForcedVars = Vars{<:AbstractArray, <:AbstractArray, <:AbstractArray, Nothing}
const StochasticForcedVars = Vars{<:AbstractArray, <:AbstractArray, <:AbstractArray, <:AbstractArray}

"""
    Vars(dev, grid)

Returns the `vars` for unforced two-dimensional turbulence on device `dev` and with `grid`.
"""
function Vars(::Dev, grid::AbstractGrid) where Dev
  T = eltype(grid)
  @devzeros Dev T (grid.nx, grid.ny) zeta u v
  @devzeros Dev Complex{T} (grid.nkr, grid.nl) zetah uh vh
  return Vars(zeta, u, v, zetah, uh, vh, nothing, nothing)
end

"""
    ForcedVars(dev, grid)

Returns the vars for forced two-dimensional turbulence on device `dev` and with `grid`.
"""
function ForcedVars(dev::Dev, grid::AbstractGrid) where Dev
  T = eltype(grid)
  @devzeros Dev T (grid.nx, grid.ny) zeta u v
  @devzeros Dev Complex{T} (grid.nkr, grid.nl) zetah uh vh Fh
  return Vars(zeta, u, v, zetah, uh, vh, Fh, nothing)
end

"""
    StochasticForcedVars(dev, grid)

Returns the vars for stochastically forced two-dimensional turbulence on device `dev` and with `grid`.
"""
function StochasticForcedVars(dev::Dev, grid::AbstractGrid) where Dev
  T = eltype(grid)
  @devzeros Dev T (grid.nx, grid.ny) zeta u v
  @devzeros Dev Complex{T} (grid.nkr, grid.nl) zetah uh vh Fh prevsol
  return Vars(zeta, u, v, zetah, uh, vh, Fh, prevsol)
end

# -------
# Solvers
# -------

"""
    calcN_advection(N, sol, t, clock, vars, params, grid)

Calculates the advection term.
"""
function calcN_advection!(N, sol, t, clock, vars, params, grid)
  @. vars.uh =   im * grid.l  * grid.invKrsq * sol
  @. vars.vh = - im * grid.kr * grid.invKrsq * sol
  @. vars.zetah = sol

  ldiv!(vars.u, grid.rfftplan, vars.uh)
  ldiv!(vars.v, grid.rfftplan, vars.vh)
  ldiv!(vars.zeta, grid.rfftplan, vars.zetah)

  uζ = vars.u        # use vars.u as scratch variable
  @. uζ *= vars.zeta # u*zeta
  vζ = vars.v        # use vars.v as scratch variable
  @. vζ *= vars.zeta # v*zeta

  uζh = vars.uh
  mul!(uζh, grid.rfftplan, uζ) # \hat{u*zeta}
  vζh = vars.vh
  mul!(vζh, grid.rfftplan, vζ) # \hat{v*zeta}

  @. N = - im * grid.kr * uζh - im * grid.l * vζh
  return nothing
end

function calcN!(N, sol, t, clock, vars, params, grid)
  calcN_advection!(N, sol, t, clock, vars, params, grid)
  addforcing!(N, sol, t, clock, vars, params, grid)
  return nothing
end

addforcing!(N, sol, t, clock, vars::Vars, params, grid) = nothing

function addforcing!(N, sol, t, clock, vars::ForcedVars, params, grid)
  params.calcF!(vars.Fh, sol, t, clock, vars, params, grid)
  @. N += vars.Fh
  return nothing
end

function addforcing!(N, sol, t, clock, vars::StochasticForcedVars, params, grid)
  if t == clock.t # not a substep
    @. vars.prevsol = sol # sol at previous time-step is needed to compute budgets for stochastic forcing
    params.calcF!(vars.Fh, sol, t, clock, vars, params, grid)
  end
  @. N += vars.Fh
  return nothing
end

# ----------------
# Helper functions
# ----------------

"""
    updatevars!(prob)

Update variables in `vars` with solution in `sol`.
"""
function updatevars!(prob)
  vars, grid, sol = prob.vars, prob.grid, prob.sol
  @. vars.zetah = sol
  @. vars.uh =   im * grid.l  * grid.invKrsq * sol
  @. vars.vh = - im * grid.kr * grid.invKrsq * sol
  ldiv!(vars.zeta, grid.rfftplan, deepcopy(vars.zetah))
  ldiv!(vars.u, grid.rfftplan, deepcopy(vars.uh))
  ldiv!(vars.v, grid.rfftplan, deepcopy(vars.vh))
  return nothing
end

"""
    set_zeta!(prob, zeta)

Set the solution `sol` as the transform of `zeta` and update variables.
"""
function set_zeta!(prob, zeta)
  mul!(prob.sol, prob.grid.rfftplan, zeta)
  CUDA.@allowscalar prob.sol[1, 1] = 0 # zero domain average

  updatevars!(prob)

  return nothing
end

"""
    energy(prob)

Returns the domain-averaged kinetic energy, ∫ ½(u²+v²)dxdy / (Lx Ly), for the solution in `sol`.
"""
@inline function energy(prob)
  sol, vars, grid = prob.sol, prob.vars, prob.grid
  energyh = vars.uh # use vars.uh as scratch variable

  @. energyh = 1 / 2 * grid.invKrsq * abs2(sol)
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(energyh, grid)
end

"""
    enstrophy(prob)

Returns the domain-averaged enstrophy, ∫ ½ ζ² dxdy / (Lx Ly), for the solution in `sol`.
"""
@inline function enstrophy(prob)
  sol, grid = prob.sol, prob.grid
  return 1 / (2 * grid.Lx * grid.Ly) * parsevalsum(abs2.(sol), grid)
end

"""
    energy_dissipation(prob)

Returns the domain-averaged energy dissipation rate. `nν` must be >= 1.
"""
@inline function energy_dissipation(prob)
  sol, vars, params, grid = prob.sol, prob.vars, prob.params, prob.grid
  energy_dissipationh = vars.uh # use vars.uh as scratch variable

  @. energy_dissipationh = params.ν * grid.Krsq^(params.nν - 1) * abs2(sol)
  CUDA.@allowscalar energy_dissipationh[1, 1] = 0
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(energy_dissipationh, grid)
end

"""
    enstrophy_dissipation(prob)

Returns the domain-averaged enstrophy dissipation rate. `nν` must be >= 1.
"""
@inline function enstrophy_dissipation(prob)
  sol, vars, params, grid = prob.sol, prob.vars, prob.params, prob.grid
  enstrophy_dissipationh = vars.uh # use vars.uh as scratch variable

  @. enstrophy_dissipationh = params.ν * grid.Krsq^params.nν * abs2(sol)
  CUDA.@allowscalar enstrophy_dissipationh[1, 1] = 0
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(enstrophy_dissipationh, grid)
end

"""
    energy_work(prob)
    energy_work(sol, vars, grid)

Returns the domain-averaged rate of work of energy by the forcing `Fh`.
"""
@inline function energy_work(sol, vars::ForcedVars, grid)
  energy_workh = vars.uh # use vars.uh as scratch variable

  @. energy_workh = grid.invKrsq * sol * conj(vars.Fh)
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(energy_workh, grid)
end

@inline function energy_work(sol, vars::StochasticForcedVars, grid)
  energy_workh = vars.uh # use vars.uh as scratch variable

  @. energy_workh = grid.invKrsq * (vars.prevsol + sol) / 2 * conj(vars.Fh) # Stratonovich
  # @. energy_workh = grid.invKrsq * vars.prevsol * conj(vars.Fh)           # Ito
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(energy_workh, grid)
end

@inline energy_work(prob) = energy_work(prob.sol, prob.vars, prob.grid)

"""
    enstrophy_work(prob)
    enstrophy_work(sol, vars, grid)

Returns the domain-averaged rate of work of enstrophy by the forcing `Fh`.
"""
@inline function enstrophy_work(sol, vars::ForcedVars, grid)
  enstrophy_workh = vars.uh # use vars.uh as scratch variable

  @. enstrophy_workh = sol * conj(vars.Fh)
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(enstrophy_workh, grid)
end

@inline function enstrophy_work(sol, vars::StochasticForcedVars, grid)
  enstrophy_workh = vars.uh # use vars.uh as scratch variable

  @. enstrophy_workh = (vars.prevsol + sol) / 2 * conj(vars.Fh) # Stratonovich
  # @. enstrophy_workh = grid.invKrsq * vars.prevsol * conj(vars.Fh)           # Ito
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(enstrophy_workh, grid)
end

@inline enstrophy_work(prob) = enstrophy_work(prob.sol, prob.vars, prob.grid)

"""
    energy_drag(prob)

Returns the extraction of domain-averaged energy by drag/hypodrag μ.
"""
@inline function energy_drag(prob)
  sol, vars, params, grid = prob.sol, prob.vars, prob.params, prob.grid

  energy_dragh = vars.uh # use vars.uh as scratch variable

  @. energy_dragh = params.μ * grid.Krsq^(params.nμ - 1) * abs2(sol)
  CUDA.@allowscalar energy_dragh[1, 1] = 0

  return 1 / (grid.Lx * grid.Ly) * parsevalsum(energy_dragh, grid)
end

"""
    enstrophy_drag(prob)

Returns the extraction of domain-averaged enstrophy by drag/hypodrag μ.
"""
@inline function enstrophy_drag(prob)
  sol, vars, params, grid = prob.sol, prob.vars, prob.params, prob.grid

  enstrophy_dragh = vars.uh # use vars.uh as scratch variable

  @. enstrophy_dragh = params.μ * grid.Krsq^params.nμ * abs2(sol)
  CUDA.@allowscalar enstrophy_dragh[1, 1] = 0
  return 1 / (grid.Lx * grid.Ly) * parsevalsum(enstrophy_dragh, grid)
end

end # module


text:
using Makie
using Colors

using IntervalConstraintProgramming, IntervalArithmetic

## Constraint programming
solid_torus = @constraint (3 - sqrt(x^2 + y^2))^2 + z^2 <= 1

half_space = @constraint (x + y) + z <= 1

x = interval(-5, 5)
Y = IntervalBox(x, 3)

@time paving = pave(solid_torus ∩ half_space, Y, 0.1);

inner = paving.inner
boundary = paving.boundary;

## Makie plotting set-up
positions = Point{3, Float32}[Point3(mid(x)...) for x in vcat(inner, boundary)]
scales = Vec3f0[Vec3f0([diam(x) for x in xx]) for xx in vcat(inner, boundary)]

zs = Float32[x[3] for x in positions]
minz = minimum(zs)
maxz = maximum(zs)

xs = Float32[x[1] for x in positions]
minx = minimum(xs)
maxx = maximum(xs)

colors1 = RGBA{Float32}[RGBA( (zs[i]-minz)/(maxz-minz), (xs[i]-minx)/(maxx-minx), 0f0, 0.1f0)
        for i in 1:length(inner)];
colors2 = RGBA{Float32}[RGBA( 0.5f0, 0.5f0, 0.5f0, 0.02f0) for x in boundary];
colors = vcat(colors1, colors2);

## Makie plotting:
cube = Rect{3, Float32}(Vec3f0(-0.5, -0.5, -0.5), Vec3f0(1, 1, 1))
# centre, widths

meshscatter(positions, marker=cube, scale=scales, color=colors, transparency=true)


text:
# include("../src/LarSurf.jl")

using Revise
using LarSurf
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm, SparseArrays

using Pio3d

## Read data from file
# threshold = 4000
# pth = Pio3d.datasets_join_path("medical/orig/sample-data/nrn4.pklz")
# datap = Pio3d.read3d(pth)
#
# data3d = datap["data3d"]
# data_size = LarSurf.size_as_array(size(data3d))
# segmentation = data3d .> threshold

## Generate data
segmentation = LarSurf.generate_slope([9,10,11])
filteredFV, Flin, (V, model) = LarSurf.get_surface_grid_per_block(segmentation, [2,3,4])
# LarSurf.get_surface
# filtered_bigFV = [
#     bigFV[LarSurf.sub_grid_face_id_to_orig_grid_face_id(data_size, block_size, offset1, fid)]
#     for fid=1:length(Flin) if (Flin[fid] == 1)
# ]

Plasm.View((bigV,[bigVV, bigEV, filtered_bigFV]))
# for fid=1:size(Flin)[2]
#     if 0 < Flin[1,fid]
#         big_fid = LarSurf.sub_grid_face_id_to_orig_grid_face_id(data_size, block_size, offset1, fid)
#         bigFchar[big_fid] = 1
#     end
# end

# LarSurf.grid_face_id_to_cartesian(block_size, )
# A = Flin
# rows = rowvals(A)
# vals = nonzeros(A)
# m, n = size(A)
# for i = 1:n
#    for j in nzrange(A, i)
#       row = rows[j]
#       val = vals[j]
#       print("row: ", row)
#       print("col: ", j)
#       print("col: ", i)
#       println("val: ", val)
# #       LarSurf.sub_grid_face_id_to_orig_face_id(data_size, block_size, [0,0,0], )
#       # perform sparse wizardry...
#    end
# end


text:
function ecFPA(grid::Array{Tuple{String,String},1})

    if Fermi.CurrentOptions["precision"] == "single"
        prec = Float32
    elseif Fermi.CurrentOptions["precision"] == "double"
        prec = Float64
    else
        throw(Fermi.InvalidFermiOptions("Invalid precision: $(A)"))
    end

    ecFPA{prec}(grid)
end

function ecFPA{T}(grid::Array{Tuple{String,String},1}) where T <: AbstractFloat
    methods = ["RHF", "CAS", "ecCCSD", "ecCCSD(T)"]
    basis = String[]
    data = zeros(length(grid), 4)
    molecule = Molecule()
    x = 1
    for (b,m) in grid
        row = get_ecfpa_row(molecule, m, b, T)
        push!(basis, b)
        for i in eachindex(row)
            data[x,i] = row[i]
        end
        x += 1
    end

    extrapolated, ext_ind = extrapolate(data, basis, methods)
    push!(basis, "CBS")

    @output "\n ⇒ Final ecFPA (absolute):\n"
    fpa = ecFPA{T}(basis, methods, extrapolated, ext_ind)
    print_fpa(fpa)
    @output "⇒ Final ecFPA Energy:   {:15.10f}\n" fpa.data[end,end]
    return fpa
end

function get_ecfpa_row(molecule::Molecule, upmethod::String, basis::String, T::DataType)

    Fermi.CurrentOptions["basis"] = basis

    @output "   → Starting ecFPA row for {}\n\n" basis
    aoint = Fermi.Integrals.ConventionalAOIntegrals(molecule)
    RHFwfn = Fermi.HartreeFock.RHF(molecule, aoint)

    if lowercase(upmethod) == "rhf"
        return RHFwfn.energy
    end

    drop_occ = Fermi.CurrentOptions["drop_occ"]
    drop_vir = Fermi.CurrentOptions["drop_vir"]

    @output "Transforming Integrals..."
    moint = Fermi.Integrals.PhysRestrictedMOIntegrals{T}(RHFwfn.ndocc, RHFwfn.nvir, drop_occ, drop_vir, RHFwfn.C, aoint)
    aoint = nothing

    CCSDwfn = RCCSD{T}(RHFwfn, moint, Fermi.CoupledCluster.CTF())

    if lowercase(upmethod) == "ccsd"
        return RHFwfn.energy, CCSDwfn.GuessEnergy, CCSDwfn.CorrelationEnergy
    end

    CCSDpTwfn = RCCSDpT{T}(CCSDwfn, moint)

    if lowercase(upmethod) == "ccsd(t)"
        return RHFwfn.energy, CCSDwfn.GuessEnergy, CCSDwfn.CorrelationEnergy, CCSDpTwfn.correction+CCSDwfn.CorrelationEnergy
    end
end


text:
# ------------------------------------------------------------------------------------------
# # Linear models, loss functions, gradients, SGD
# (c) Deniz Yuret, 2019
# * Objectives: Define, train and visualize a simple model; understand gradients and SGD;
# learn to use the GPU.
# * Prerequisites: [Callable
# objects](https://docs.julialang.org/en/v1/manual/methods/#Function-like-objects-1),
# [Generator expressions](https://docs.julialang.org/en/v1/manual/arrays/#Generator-
# Expressions-1), [MNIST](20.mnist.ipynb), [Iterators](25.iterators.ipynb)
# * New functions:
# [mnistdata](https://github.com/denizyuret/Knet.jl/blob/master/data/mnist.jl),
# [accuracy](http://denizyuret.github.io/Knet.jl/latest/reference/#Knet.accuracy),
# [zeroone](http://denizyuret.github.io/Knet.jl/latest/reference/#Knet.zeroone),
# [nll](http://denizyuret.github.io/Knet.jl/latest/reference/#Knet.nll),
# [Param, @diff, value, params,
# grad](http://denizyuret.github.io/Knet.jl/latest/reference/#AutoGrad),
# [sgd](http://denizyuret.github.io/Knet.jl/latest/reference/#Knet.sgd),
# [progress,
# progress!](http://denizyuret.github.io/Knet.jl/latest/reference/#Knet.progress),
# [gpu](http://denizyuret.github.io/Knet.jl/latest/reference/#Knet.gpu),
# [KnetArray](http://denizyuret.github.io/Knet.jl/latest/reference/#Knet.KnetArray),
# [load](http://denizyuret.github.io/Knet.jl/latest/reference/#Knet.load),
# [save](http://denizyuret.github.io/Knet.jl/latest/reference/#Knet.save)
#
#
# <img src="https://www.oreilly.com/library/view/tensorflow-for-
# deep/9781491980446/assets/tfdl_0401.png" alt="A linear model" width=300/> ([image
# source](https://www.oreilly.com/library/view/tensorflow-for-deep/9781491980446/ch04.html))
#
# In Knet, a machine learning model is defined using plain Julia code. A typical model
# consists of a **prediction** and a **loss** function. The prediction function takes some
# input, returns the prediction of the model for that input. The loss function measures how
# bad the prediction is with respect to some desired output. We train a model by adjusting
# its parameters to reduce the loss.
#
# In this section we will implement a simple linear model to classify MNIST digits. The
# prediction function will return 10 scores for each of the possible labels 0..9 as a linear
# combination of the pixel values. The loss function will convert these scores to normalized
# probabilities and return the average -log probability of the correct answers. Minimizing
# this loss should maximize the scores assigned to correct answers by the model. We will
# make use of the loss gradient with respect to each parameter, which tells us the direction
# of the greatest loss increase. We will improve the model by moving the parameters in the
# opposite direction (using a GPU if available). We will visualize the model weights and
# performance over time. The final accuracy of about 92% is close to the limit of what we
# can achieve with this type of model. To improve further we must look beyond linear models.
# ------------------------------------------------------------------------------------------

# Set display width, load packages, import symbols
ENV["COLUMNS"]=72
using Pkg; for p in ("Knet","AutoGrad","Plots","Images","ImageMagick"); haskey(Pkg.installed(),p) || Pkg.add(p); end
using Statistics: mean
using Base.Iterators: flatten
import Random # seed!
using Knet: Knet, AutoGrad, dir, Data, Param, @diff, value, params, grad, progress, progress!, gpu, KnetArray, load, save
# The following are defined for instruction even though they are provided in Knet
# using Knet: accuracy, zeroone, nll, sgd

# Load data (mnistdata basically replicates mnist.ipynb)
include(Knet.dir("data","mnist.jl"))
dtrn,dtst = mnistdata(xsize=(784,:),xtype=Array)
println.(summary.((dtrn,dtst)));

# ------------------------------------------------------------------------------------------
# ## Model definition
# ------------------------------------------------------------------------------------------

# In Julia we define a new datatype using `struct`:
struct Linear; w; b; end

# The new struct comes with a default constructor:
model = Linear(0.01 * randn(10,784), zeros(10))

# We can define other constructors with different inputs:
Linear(i::Int,o::Int,scale=0.01) = Linear(scale * randn(o,i), zeros(o))

# This one allows instances to be defined using input and output sizes:
model = Linear(784,10)

# ------------------------------------------------------------------------------------------
# ## Prediction
# ------------------------------------------------------------------------------------------

# We turn Linear instances into callable objects for prediction:
(m::Linear)(x) = m.w * x .+ m.b

x,y = first(dtst) # The first minibatch from the test set
summary.((x,y))

Int.(y)' # correct answers are given as an array of integers (remember we use 10 for 0)

ypred = model(x)  # Predictions on the first minibatch: a 10x100 score matrix

# We can calculate the accuracy of our model for the first minibatch
accuracy(model,x,y) = mean(y' .== map(i->i[1], findmax(Array(model(x)),dims=1)[2]))
accuracy(model,x,y)

# We can calculate the accuracy of our model for the whole test set
accuracy(model,data) = mean(accuracy(model,x,y) for (x,y) in data)
accuracy(model,dtst)

# ZeroOne loss (or error) is defined as 1 - accuracy
zeroone(x...) = 1 - accuracy(x...)
zeroone(model,dtst)

# ------------------------------------------------------------------------------------------
# ## Loss function
# ------------------------------------------------------------------------------------------

# For classification we use negative log likelihood loss (aka cross entropy, softmax loss, NLL)
# This is the average -log probability assigned to correct answers by the model
function nll(scores, y)
    expscores = exp.(scores)
    probabilities = expscores ./ sum(expscores, dims=1)
    answerprobs = (probabilities[y[i],i] for i in 1:length(y))
    mean(-log.(answerprobs))
end

# model(x) gives predictions, let model(x,y) give the loss
(m::Linear)(x, y) = nll(m(x), y)
model(x,y)

# We can also use the Knet nll implementation for efficiency
(m::Linear)(x, y) = Knet.nll(m(x), y)
model(x,y)

# If the input is a dataset compute average loss:
(m::Linear)(data::Data) = mean(m(x,y) for (x,y) in data)

# Here is per-instance average negative log likelihood for the whole test set
model(dtst)

# ------------------------------------------------------------------------------------------
# **Bonus question:** What is special about the loss value 2.3?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Calculating the gradient using AutoGrad
# ------------------------------------------------------------------------------------------

@doc AutoGrad

# Redefine the constructor to use Param's so we can compute gradients
Linear(i::Int,o::Int,scale=0.01) =
    Linear(Param(scale * randn(o,i)), Param(zeros(o)))

# Set random seed for replicability
Random.seed!(9);

# Use a larger scale to get a large initial loss
model = Linear(784,10,1.0)

# We can still do predictions and calculate loss:
model(x,y)

# And we can do the same loss calculation also computing gradients:
J = @diff model(x,y)

# To get the actual loss value from J:
value(J)

# params(J) returns an iterator of Params J depends on (i.e. model.b, model.w):
params(J) |> collect

# To get the gradient of a parameter from J:
∇w = grad(J,model.w)

# Note that each gradient has the same size and shape as the corresponding parameter:
@show ∇b = grad(J,model.b);

# ------------------------------------------------------------------------------------------
# ## Checking the gradient using numerical approximation
#
# What does ∇b represent?
#
# ∇b[10] = 0.79 means if I increase b[10] by ϵ, loss will increase by 0.79ϵ
# ------------------------------------------------------------------------------------------

# Loss for the first minibatch with the original parameters
@show value(model.b)
model(x,y)

# To numerically check the gradient let's increase the last entry of b by +0.1.
model.b[10] = 0.1

# We see that the loss moves by ≈ +0.79*0.1 as expected.
@show value(model.b)
model(x,y)

# Reset the change.
model.b[10] = 0

# ------------------------------------------------------------------------------------------
# ## Checking the gradient using manual implementation
# ------------------------------------------------------------------------------------------

# Without AutoGrad we would have to define the gradients manually:
function nllgrad(model,x,y)
    scores = model(x)
    expscores = exp.(scores)
    probabilities = expscores ./ sum(expscores, dims=1)
    for i in 1:length(y); probabilities[y[i],i] -= 1; end
    dJds = probabilities / length(y)
    dJdw = dJds * x'
    dJdb = vec(sum(dJds,dims=2))
    dJdw,dJdb
end;

∇w2,∇b2 = nllgrad(model,x,y)

∇w2 ≈ ∇w

∇b2 ≈ ∇b

# ------------------------------------------------------------------------------------------
# ## Training with Stochastic Gradient Descent (SGD)
# ------------------------------------------------------------------------------------------

# Here is a single SGD update:
function sgdupdate!(func, args; lr=0.1)
    fval = @diff func(args...)
    for param in params(fval)
        ∇param = grad(fval, param)
        param .-= lr * ∇param
    end
    return value(fval)
end

# We define SGD for a dataset as an iterator so that:
# 1. We can monitor and report the training loss
# 2. We can take snapshots of the model during training
# 3. We can pause/terminate training when necessary
sgd(func, data; lr=0.1) =
    (sgdupdate!(func, args; lr=lr) for args in data)

# Let's train a model for 10 epochs to compare training speed on cpu vs gpu.
# progress!(itr) displays a progress bar when wrapped around an iterator like this:
# 2.94e-01  100.00%┣████████████████████┫ 6000/6000 [00:10/00:10, 592.96/s] 2.31->0.28
model = Linear(784,10)
@show model(dtst)
progress!(sgd(model, repeat(dtrn,10)))
@show model(dtst);

# ------------------------------------------------------------------------------------------
# ## Using the GPU
# ------------------------------------------------------------------------------------------

# The training would go a lot faster on a GPU:
# 2.94e-01  100.00%┣███████████████████┫ 6000/6000 [00:02/00:02, 2653.45/s]  2.31->0.28
# To work on a GPU, all we have to do is convert Arrays to KnetArrays:
if gpu() >= 0  # gpu() returns a device id >= 0 if there is a GPU, -1 otherwise
    atype = KnetArray{Float32}  # KnetArrays are stored and operated in the GPU
    dtrn,dtst = mnistdata(xsize=(784,:),xtype=atype)
    Linear(i::Int,o::Int,scale=0.01) =
        Linear(Param(atype(scale * randn(o,i))),
               Param(atype(zeros(o))))

    model = Linear(784,10)
    @show model(dtst)
    progress!(sgd(model,repeat(dtrn,10)))
    @show model(dtst)
end;

# ------------------------------------------------------------------------------------------
#
# ## Recording progress
# ------------------------------------------------------------------------------------------

function trainresults(file, model)
    if (print("Train from scratch? (~77s) "); readline()[1]=='y')
        # This will take every nth element of an iterator:
        takeevery(n,itr) = (x for (i,x) in enumerate(itr) if i % n == 1)
        # We will use it to snapshot model and results every epoch (i.e. 600 iterations)
        lin = ((deepcopy(model),model(dtrn),model(dtst),zeroone(model,dtrn),zeroone(model,dtst))
        # (progress displays a bar like progress! but returns an iterator, progress! returns nothing)
               for x in takeevery(length(dtrn), progress(sgd(model,repeat(dtrn,100)))))
        # Save it as a 5x100 array
        lin = reshape(collect(flatten(lin)),(5,:))
        # Knet.save and Knet.load can be used to store models in files
        Knet.save(file,"results",lin)
    else
        isfile(file) || download("http://people.csail.mit.edu/deniz/models/tutorial/$file", file)
        lin = Knet.load(file,"results")
    end
    return lin
end

# 2.43e-01  100.00%┣████████████████▉┫ 60000/60000 [00:44/00:44, 1349.13/s]
lin = trainresults("lin113.jld2",Linear(784,10));

# ------------------------------------------------------------------------------------------
# ## Linear model shows underfitting
# ------------------------------------------------------------------------------------------

using Plots; default(fmt = :png)

# Demonstrates underfitting: training loss not close to 0
# Also slight overfitting: test loss higher than train
plot([lin[2,:], lin[3,:]],ylim=(.0,.4),labels=[:trnloss :tstloss],xlabel="Epochs",ylabel="Loss")

# this is the error plot, we get to about 7.5% test error, i.e. 92.5% accuracy
plot([lin[4,:], lin[5,:]],ylim=(.0,.12),labels=[:trnerr :tsterr],xlabel="Epochs",ylabel="Error")

# ------------------------------------------------------------------------------------------
# ## Visualizing the learned weights
# ------------------------------------------------------------------------------------------

# Let us visualize the evolution of the weight matrix as images below
# Each row is turned into a 28x28 image with positive weights light and negative weights dark gray
using Images, ImageMagick
for t in 10 .^ range(0,stop=log10(size(lin,2)),length=20) #logspace(0,2,20)
    i = ceil(Int,t)
    f = lin[1,i]
    w1 = reshape(Array(value(f.w))', (28,28,1,10))
    w2 = clamp.(w1.+0.5,0,1)
    IJulia.clear_output(true)
    display(hcat([mnistview(w2,i) for i=1:10]...))
    display("Epoch $(i-1)")
    sleep(1) # (0.96^i)
end

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
cd(@__DIR__); import Pkg; Pkg.activate(".")
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
Pkg.develop(path="$(homedir())/github/SimpleSDMLayers.jl")
using Revise
using SimpleSDMLayers
using Test

## Step 5: investigate error with hcat/vcat
# Checkout previous commit
testpath = "../../../SimpleSDMLayers.jl/test/"
include(joinpath(testpath, "runtests.jl"))
include(joinpath(testpath, "overloads.jl"))

# Failing test
l1 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=0.0, right=10.0, bottom=0.0, top=10.0)
l2 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=0.0, right=10.0, bottom=10.0, top=20.0)
l3 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=10.0, right=20.0, bottom=0.0, top=10.0)
l4 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=10.0, right=20.0, bottom=10.0, top=20.0)

ml1 = hcat(l1, l3)
vl1 = vcat(l1, l2)
ml2 = hcat(l2, l4)
vl2 = vcat(l3, l4)

@test all(vcat(ml1, ml2).grid == hcat(vl1, vl2).grid)

# Keep values
pre = (
    l1 = l1,
    l2 = l2,
    l3 = l3,
    l4 = l4,
    ml1 = ml1,
    vl1 = vl1,
    ml2 = ml2,
    vl2 = vl2,
)
# Checkout latest commit
new = (
    l1 = l1,
    l2 = l2,
    l3 = l3,
    l4 = l4,
    ml1 = ml1,
    vl1 = vl1,
    ml2 = ml2,
    vl2 = vl2,
)

# Compare
pre
new
pre.l1[left=0.0, right=10.0, bottom=0.0, top=10.0]
#
coords_l1 = (left=0.0, right=10.0, bottom=0.0, top=10.0)
coords_l2 = (left=0.0, right=10.0, bottom=10.0, top=20.0)
coords_l3 = (left=10.0, right=20.0, bottom=0.0, top=10.0)
coords_l4 = (left=10.0, right=20.0, bottom=10.0, top=20.0)
boundingbox(l::SimpleSDMLayer) = (left = l.left, right = l.right, bottom = l.bottom, top = l.top)
function fixup(pre, new, coords)
    println(pre)
    println(new)
    fix = pre[coords]
    @test boundingbox(new) == boundingbox(fix)
    @test size(new) == size(fix)
    @test new.grid == fix.grid
    return fix
end
fixup(pre.l1, new.l1, coords_l1)
fixup(pre.l2, new.l2, coords_l2)
fixup(pre.l3, new.l3, coords_l3)
fixup(pre.l4, new.l4, coords_l4)

ml1 = hcat(pre.l1, pre.l3)
ml1 = hcat(new.l1, new.l3)

vl1 = vcat(pre.l1, pre.l2) # whut, no latitude span?
vl1 = vcat(new.l1, new.l2)

ml2 = hcat(pre.l2, pre.l4)
ml2 = hcat(new.l2, new.l4)

vl2 = vcat(pre.l3, pre.l4) # ??
vl2 = vcat(new.l3, new.l4) # ??

vcat(pre.ml1, pre.ml2)
vcat(new.ml1, new.ml2)

hcat(pre.vl1, pre.vl2)
hcat(new.vl1, new.vl2)

# First hcat
l1 = new.l1
l2 = new.l3
t1 = hcat(l1, l2)
boundingbox(t1)
latitudes(t1)
longitudes(t1)

# First vcat
l1 = new.l1
l2 = new.l2
t1 = vcat(l1, l2)
boundingbox(t1)
latitudes(t1)
longitudes(t1)
extrema(latitudes(t1))
Tuple(latitudes(t1)[[1, end]])

coords_ltot = (left = 0.0, right = 20.0, bottom = 0.0, top = 20.0)
ltot = SimpleSDMPredictor(WorldClim, BioClim, 1; coords_ltot...)
ltot.grid
ltot[0.0, 0.0]
ltot[0.0, 20.0]
ltot[20.0, 0.0]
ltot[20.0, 20.0]

# new.l1[0.0, 0.0]
# new.l2[0.0, 20.0]
# new.l3[20.0, 0.0]
# new.l4[20.0, 20.0]


text:
# Checking the sampler tables

module TablesTest  # separate namespace

using Base.Test

include("../src/randtn_zigtables.jl")

# basic properties
@test all(tnzig_yl .<= tnzig_yu)
@test all(diff(tnzig_x) .== tnzig_d)

# ncell index into x
n = length(tnzig_x)
xsimn = 10^6
for i = 1:xsimn
    x = tnzig_xmin + (tnzig_xmax - tnzig_xmin) * rand()
    i = tnzig_i0 + ifloor(x * tnzig_invh)
    # array boundaries
    @test i >= 1
    @test i <= length(tnzig_ncell)
    j = tnzig_ncell[i]
    @test j >= 1
    @test j <= n - 1
    # x[j] <= x < x[j+2] (or x[j+1] for j = n-1)
    @test tnzig_x[j] <= x
    @test tnzig_x[min(n, j+2)] > x
end

end

max_stars_repo_name: jdrugo/fasttn

max_issues_repo_name: jdrugo/fasttn

max_forks_repo_name: jdrugo/fasttn


text:
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = page_model()
    include("../src/components/RCPSSPScenario.jl")
    include("../src/components/GDP.jl")

    rcpsspscenario = addrcpsspscenario(m, "NDCs")
    gdp = add_comp!(m, GDP)

    gdp[:grw_gdpgrowthrate] = rcpsspscenario[:grw_gdpgrowthrate]
    gdp[:pop0_initpopulation] = readpagedata(m, "data/pop0_initpopulation.csv")
    gdp[:pop_population] = readpagedata(m, "test/validationdata/$valdir/pop_population.csv")
    gdp[:y_year] = Mimi.dim_keys(m.md, :time)
    gdp[:y_year_0] = 2015.

    p=load_parameters(m)
    set_leftover_params!(m,p)

    # run model
    run(m)

    # Generated data
    gdp = m[:GDP, :gdp]

    # Recorded data
    gdp_compare = readpagedata(m, "test/validationdata/$valdir/gdp.csv")

    @test gdp ≈ gdp_compare rtol=100

    cons_percap_consumption_0_compare = readpagedata(m, "test/validationdata/cons_percap_consumption_0.csv")
    @test m[:GDP, :cons_percap_consumption_0] ≈ cons_percap_consumption_0_compare rtol=1e-2
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
using CBCLFace
using Base.Test
using PyCall
@pyimport matplotlib.pyplot as plt

# write your own tests here
println("Running Tests...")

FD, NFD = traindata(hist_eq=true)

for i=1:20
  plt.imshow(reshape(FD[:,i],(19,19))';interpolation="nearest")
  plt.colorbar()
  plt.gray()
  plt.show(block=true)
end

for i=1:20
  plt.imshow(reshape(FD[:,i],(19,19))';interpolation="nearest")
  plt.colorbar()
  plt.gray()
  plt.show(block=true)
end


text:
module TestDoctest

## doctest
using MultipleTesting
using Base.Test

@unix ? begin
    using Lexicon
    dt = doctest(MultipleTesting)
    @test length(failed(dt)) == 0
end : println("! Skipping doctest on windows...")

end


text:
"""
    optimalportfolio_enumerate(mkt::SameCostsMarket) -> X, v

Produce the optimal portfolio for [`SameCostsMarket`](@ref) defined by `mkt`. Very slow;
you probably should use [`applicationorder_list`](@ref) or [`applicationorder_heap`](@ref) instead.

```julia-repl
julia> mkt = SameCostsMarket([0.2, 0.5, 0.1, 0.6, 0.1], [1, 4, 9, 1, 8], 2);

julia> optimalportfolio_enumerate(mkt)
([2, 3], 2.7)
```
"""
function optimalportfolio_enumerate(mkt::SameCostsMarket)::Tuple{Vector{Int},Float64}
    mkt.m ≥ 21 && @warn "Enumeration is slow for large markets"

    v, X = maximum(multiset_combinations(1:mkt.m, mkt.h)) do X
        return valuation_nopermute_sorted(X, mkt), X
    end

    return mkt.perm[X], v
end

"""
    optimalportfolio_enumerate(mkt::VariedCostsMarket) -> X, v

Produce the optimal portfolio for [`VariedCostsMarket`](@ref) defined by `mkt`. Very slow;
you probably should use [`optimalportfolio_dynamicprogram`](@ref) or [`optimalportfolio_branchbound`](@ref)
or [`optimalportfolio_simulatedannealing`](@ref) instead.

```julia-repl
julia> mkt = VariedCostsMarket([0.2, 0.5, 0.1, 0.6, 0.1], [1, 4, 9, 1, 8], [2, 4, 2, 5, 1], 8);

julia> optimalportfolio_enumerate(mkt)
([2, 5, 3], 3.24)
```
"""
function optimalportfolio_enumerate(mkt::VariedCostsMarket)::Tuple{Vector{Int},Float64}
    mkt.m ≥ 21 && @warn "Enumeration is slow for large markets"

    v, X = maximum(X for X in combinations(1:mkt.m) if sum(mkt.g[X]) ≤ mkt.H) do X
        return valuation_nopermute_sorted(X, mkt), X
    end

    return mkt.perm[X], v
end

max_stars_repo_name: maxkapur/OptimalApplication

max_issues_repo_name: maxkapur/OptimalApplication

max_forks_repo_name: maxkapur/OptimalApplication


text:
# # Modal analysis of simply supported/clamped beam

# Source code: [`beam_modal_tut.jl`](beam_modal_tut.jl)

# ## Description

# Vibration analysis of a beam simply supported in one plane, and clamped
# in another. The results are compared with analytical expressions.

# ## Goals

# - Introduce definition of cross-section.
# - Show generation of finite element mesh of beams.
# - Describe geometry, displacement, and rotation fields.
# - Describe application of support conditions.
# - Calculate the discrete model quantities and solve the free vibration problem.
# - Demonstrate visualization of the free vibrations.

# ## Definition of the basic inputs

# We will probably need some linear algebra functions.
using LinearAlgebra

# The finite element code relies on the basic functionality implemented in this
# package.
using FinEtools

# The material parameters may be defined with the specification of the units.
# The elastic properties are:
E = 30002.0 * phun("ksi")
nu = 0.0;

# The mass density is expressed in customary units as
rho = 0.28 * phun("lbm/in^3")
# Here are the cross-sectional dimensions and the length of the beam between supports.
b = 1.8 * phun("in"); h = 1.8 * phun("in"); L = 100 * phun("in");

# ## Analytical frequencies

# The analytical frequencies were taken from table 8-1 of Formulas for natural
# frequency and mode shape, Robert D. Blevins, Krieger publishing company,
# Malabar Florida, reprint edition 2001. Available also in FORMULAS FOR
# DYNAMICS, ACOUSTICS AND VIBRATION, by the same author, Wiley, 2016, Table 4.2.
#
# The beam is aligned with the $Y$ global Cartesian coordinate. The beam
# behaves as a simply supported beam in the vertical plane (global Cartesian
# $YZ$), while in the horizontal plane (global Cartesian $XY$) it behaves
# as a clamped beam.
#
# The cross-sectional properties are:
A = b * h;
I2 = b * h^3 / 12;
I3 = b^3 * h / 12;
# Then the analytical vibration frequencies for the first two modes are:
@show analyt_freq = [(1 * pi)^2, (4.73004074)^2] .* (sqrt(E * I2 / rho / A) / (2 * pi * L^2));

# The purpose of the numerical model is to calculate approximation to these two
# analytical natural frequencies.
neigvs = length(analyt_freq);

# ## Cross-section

# Cross-sectional properties are incorporated in the cross-section object. The
# three arguments supplied are functions. All are returning "constants", as
# appropriate for a uniform cross section beam. In particular the first two
# functions each return the dimension of the cross-section as a constant; the
# third function defines the orientation of the cross-section in the global
# Cartesian coordinates. `[1.0, 0.0, 0.0]` is the vector that together with the
# tangent to the midline curve of the beam spans the $x_1x_2$ plane of the
# local coordinates for the beam.
using FinEtoolsFlexStructures.CrossSectionModule: CrossSectionRectangle
cs = CrossSectionRectangle(s -> b, s -> h, s -> [1.0, 0.0, 0.0])

# We can compare the analytical values of the cross-section properties with
# those stored in the `cs` structure:
@show A, I2, I3
@show cs.parameters(0.0)

# ## Mesh generation

# Now we generate the mesh of the beam. The locations of its two endpoints are:
xyz = [[0 -L / 2 0]; [0 L / 2 0]]
# We will generate
n = 4
# beam elements along the member.
using FinEtoolsFlexStructures.MeshFrameMemberModule: frame_member
fens, fes = frame_member(xyz, n, cs);
# The mesh definition consists of the nodes
@show fens
# and the finite elements
@show fes
# Note  that the cross-sectional properties are incorporated through `cs`.

# ## Material

# Material properties can be now used to create a material: isotropic elasticity
# model of the `FinEtoolsDeforLinear` package is instantiated.
using FinEtoolsDeforLinear
material = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)

# ## Fields

# Now we start constructing the discrete finite element model.
# We begin by constructing the requisite fields, geometry and displacement.
# These are the so-called "configuration variables", all initialized to zero.

# This is the geometry field. It describes the locations of the nodes.
geom0 = NodalField(fens.xyz)
# This is the displacement field, three unknown displacements per node.
u0 = NodalField(zeros(size(fens.xyz, 1), 3))
# This is the rotation (orientation) field, where the three unknown rotations
# per node are represented with a rotation matrix, in total nine numbers. The
# utility function `initial_Rfield` is used to construct this initial
# orientation field where each orientation matrix is the identity.
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield
Rfield0 = initial_Rfield(fens)
# Here we verify the number of nodes and the number of degrees of freedom in the
# rotation field per node.
@show nents(Rfield0)
@show ndofs(Rfield0)

# Finally, this is the displacement and rotation field for incremental changes,
# incremental displacements and incremental rotations. In total, 6 unknowns per
# node.
dchi = NodalField(zeros(size(fens.xyz, 1), 6))

# ## Support conditions

# Now we apply the essential boundary conditions (EBCs) to enforce the action of
# the supports at the ends of the beam.

# First we select the node at the location  `[0 -L/2 0]`. We do that by
# searching for all nodes that are located inside a tiny box centered at this
# location. The tolerance of a fraction of the length of an element (i. e. the
# distance between the nodes) is used to take the location and blow it up into
# a nonzero-volume box.
l1 = selectnode(fens; box=[0 0 -L / 2 -L / 2 0 0], tolerance=L / n / 1000)
# The nodes in this list should consist of a single node (which is why the
# tolerance is so small, in order to limit the selection to a single node near
# the location given).
@show length(l1) == 1
# The boundary condition at this point dictates zero displacements (degrees of
# freedom 1, 2, and 3) and zero rotations about $Y$ (5) and $Z$ (6). This
# leaves rotation about the $X$ free.
for i in [1,2,3,5,6]
    setebc!(dchi, l1, true, i)
end
# Similarly, the node next to the other end of the beam is selected.
l1 = selectnode(fens; box=[0 0 L / 2  L / 2 0 0], tolerance=L / n / 1000)
# And an identical boundary condition combination is enforced.
for i in [1,2,3,5,6]
    setebc!(dchi, l1, true, i)
end
# These boundary conditions now need to be "applied". This simply means that the
# prescribed values of the degrees of freedom are copied into the active
# degrees of freedom.
applyebc!(dchi)
# The essential boundary conditions will also reduce the number of free
# (unknown) degrees of freedom.
numberdofs!(dchi);
# Here we inspect the degrees of freedom in the incremental
# displacement/rotation field:
@show dchi.dofnums
# Note that the degrees of freedom are actually carried by the incremental
# field, not by the displacement or the rotation fields. There are therefore 6
# degrees of freedom per node in the incremental displacement field.

# ## Assemble the global discrete system

using FinEtoolsFlexStructures.FEMMCorotBeamModule: FEMMCorotBeam
femm = FEMMCorotBeam(IntegDomain(fes, GaussRule(1, 2)), material);

# For disambiguation we will refer to the stiffness and mass functions by
# qualifying them with the corotational-beam module, `FEMMCorotBeamModule`.
# We will use the abbreviation `CB`.
using FinEtoolsFlexStructures.FEMMCorotBeamModule
CB = FEMMCorotBeamModule

# Thus we can construct the stiffness and mass matrix as follows:
# Note that the finite element machine is the first argument. This provides
# access to the integration domain. The next argument is the geometry field,
# followed by the displacement, rotations, and incremental
# displacement/rotation fields.
K = CB.stiffness(femm, geom0, u0, Rfield0, dchi);
M = CB.mass(femm, geom0, u0, Rfield0, dchi);
# We can compare the size of the stiffness matrix with the number of degrees of
# freedom that are unknown (20).
@show size(K)

# ## Solve the free-vibration problem

# The Arnoldi algorithm implemented in the well-known `Arpack` package is used
# to solve the generalized eigenvalue problem with the sparse matrices. As is
# common in structural dynamics, we request the smallest eigenvalues in
# absolute value (`:SM`). In order to help the function decide which algorithm
# is suitable, we explicitly designate the matrices as symmetric.
using Arpack
evals, evecs, nconv = eigs(Symmetric(K), Symmetric(M); nev=neigvs, which=:SM, explicittransform = :none);

# We should check that the requested eigenvalues actually converged:
@show nconv == neigvs

# The eigenvalues (i. e. the squares of the angular frequencies) are returned in
# the vector `evals`. The mode shapes constitute the columns of the matrix
# `evecs`.
@show size(evecs)

# The natural frequencies are obtained from the squares of the angular
# frequencies. We note the use of `sqrt.` which broadcasts the square root over
# the array `evals`.
fs = sqrt.(evals) / (2 * pi);

# ## Comparison of computed and analytical results

# The approximate and analytical frequencies are now reported.
println("Approximate frequencies: $fs [Hz]")
println("Analytical frequencies: $analyt_freq [Hz]")

# Close agreement between the approximate and analytical frequencies can be
# observed: The error of the numerical solution is a fraction of a percent.
errs = abs.(analyt_freq .- fs) ./ analyt_freq
println("Relative errors of frequencies: $errs [ND]")

# ## Visualize vibration modes

# The animation will show one of the vibration modes overlaid on the undeformed
# geometry.

# The visualization utilities take advantage of the PlotlyJS library.

using VisualStructures: plot_space_box, plot_solid, render, react!, default_layout_3d, save_to_json

# The configuration during the animation needs to reflect rotations.
# The function `update_rotation_field!` will update the rotation field given a
# vibration mode.
using FinEtoolsFlexStructures.RotUtilModule: update_rotation_field!

# The magnitude of the vibration modes (displacements  and rotations) will be
# amplified with this scale factor:
scale = 1.2

# This is the mode that will be animated:
mode = 1

# In order to handle variables inside loops correctly, we create a local scope
# with the `let end` block.
let
    # The extents of the box will be preserved during animation in order to
    # eliminate changes in the viewing parameters.
    tbox = plot_space_box([[-0.2 * L -L / 2 -0.2 * L]; [+0.2 * L L / 2 +0.2 * L]])
    # This is the geometry of the structure without deformation (undeformed). It
    # is displayed as gray, partially transparent.
    tenv0 = plot_solid(fens, fes; x=geom0.values, u=0.0 .* dchi.values[:, 1:3], R=Rfield0.values, facecolor="rgb(125, 155, 125)", opacity=0.3);
    # Initially the plot consists of the box and the undeformed geometry.
    plots = cat(tbox, tenv0; dims=1)
    # Create the layout for the plot. Set the size of the window.
    layout = default_layout_3d(; options = Dict(:responsive=>true))
    # Set the aspect mode to get the correct proportions.
    layout[:scene][:aspectmode] = "data"
    # Render the undeformed structure
    pl = render(plots; layout=layout, title="Mode $(mode)")
    sleep(2.115)
    # This is the animation loop.
    # - Distribute a fraction of the selected eigenvector into the incremental displacement/rotation field.
    # - Create the deformed configuration by defining displacement field `u1` and rotation field `Rfield1`.
    # - Create the plot for the deformed configuration, and add it to the list of plots.
    # - Call the `react!` function to update the display. Sleep for a brief period of time to give the display a chance to become current.
    for xscale in scale .* sin.(collect(0:1:89) .* (2 * pi / 21))
        scattersysvec!(dchi, xscale .* evecs[:, mode])
        u1 = deepcopy(u0)
        u1.values .= dchi.values[:, 1:3]
        Rfield1 = deepcopy(Rfield0)
        update_rotation_field!(Rfield1, dchi)
        tenv1 = plot_solid(fens, fes; x=geom0.values, u=dchi.values[:, 1:3], R=Rfield1.values, facecolor="rgb(50, 55, 125)");
        plots = cat(tbox, tenv0, tenv1; dims=1)
        react!(pl, plots, pl.plot.layout)
        sleep(0.115)
    end
    # Save the plot to a Json file. It can be then re-displayed later.
    save_to_json(pl, "deformed_plot.json")
end

# Load the plot from a file.
using VisualStructures: plot_from_json
plot_from_json("deformed_plot.json")

# Nothing to return.
nothing


text:
@testset "report invalid builtin call" begin
    interp, frame = analyze_call((Int, Type{Int}, Any)) do a, b, c
        isa(a, b, c)
    end
    @test length(interp.reports) === 1
    report = first(interp.reports)
    @test report isa InvalidBuiltinCallErrorReport &&
        widenconst.(report.argtypes) == [Int, Type{Int}, Any]

    @testset "constant propagation" begin
        m = gen_virtual_module()
        Core.eval(m, quote
            struct T
                v
            end
            access_field(t, sym) = getfield(t, sym)
        end)

        interp, frame = Core.eval(m, quote
            $analyze_call(t->access_field(t,:v), (T,))
        end)
        @test isempty(interp.reports)

        interp, frame = Core.eval(m, quote
            $analyze_call(t->access_field(t,:w), (T,))
        end)
        @test length(interp.reports) === 1
        er = first(interp.reports)
        @test er isa NoFieldErrorReport
        @test er.typ === m.T
        @test er.name === :w

        interp, frame = Core.eval(m, quote
            $analyze_call(t->access_field(t,:v), (T,))
        end)
        @test isempty(interp.reports)
    end
end

@testset "malformed getfield" begin
    let
        # shouldn't error
        interp, frame = analyze_call((Any,)) do a
            getfield(a)
        end
        @test length(interp.reports) == 1
        @test first(interp.reports) isa InvalidBuiltinCallErrorReport
    end
end

@testset "getfield with abstract global variable" begin
    # nested module access will be resolved as a direct call of `getfield`
    let
        res = @analyze_toplevel begin
            module foo

            const bar = sum

            module baz

            using ..foo

            foo.bar("julia") # -> NoMethodErrorReports

            end # module bar

            end # module foo
        end

        @test isempty(res.toplevel_error_reports)
        test_sum_over_string(res)
    end

    # this should work even if the accessed variable is not constant
    let
        res = @analyze_toplevel begin
            module foo

            bar = sum

            module baz

            using ..foo

            foo.bar("julia") # -> NoMethodErrorReports

            end # module bar

            end # module foo
        end

        @test isempty(res.toplevel_error_reports)
        test_sum_over_string(res)
    end
end

@testset "special case `return_type`" begin
    # don't report invalid method calls simulated in `return_type_tfunc`
    let
        interp, frame = analyze_call(()->CC.return_type(sum, Tuple{String}))
        @test isempty(interp.reports)
    end

    # report invalid call of `return_type` itself
    let
        interp, frame = analyze_call(()->CC.return_type(sum))
        @test !isempty(interp.reports)
        @test any(interp.reports) do report
            return report isa NoMethodErrorReport &&
                any(Base.Fix1(occursin, "return_type"), report.sig)
        end
    end

    # end to end
    let
        # this shouldn't report "no matching method found for call signature: Base.iterate(itr::DataType)"
        # , which otherwise will be caught in `abstract_cal` in `return_type_tfunc`
        interp, frame = analyze_call(()->Dict('a' => 1,
                                              :b => 2)
                                     )
        @test isempty(interp.reports)
    end
end


text:
import Base.Docs: Binding, @var
using Markdown: MD, Code, Paragraph, plain

# flat_content(md) = md
# flat_content(xs::Vector) = reduce((xs, x) -> vcat(xs,flat_content(x)), [], xs)
# flat_content(md::MD) = flat_content(md.content)
flatten(md::MD) = MD(flat_content(md))

# Faster version

function flat_content!(xs, out = [])
  for x in xs
    if isa(x, MD)
      flat_content!(x.content, out)
    else
      push!(out, x)
    end
  end
  return out
end

flat_content(md::MD) = flat_content!(md.content)

function hasdoc(b::Binding)
    for m in Docs.modules
        meta = Docs.meta(m)
        if haskey(meta, b) || haskey(meta, b[])
            return true
        end
    end
    false
end

hasdoc(b) = hasdoc(Docs.aliasof(b, typeof(b)))

function trygetdoc(b)
  docs = try
    Docs.doc(b)
  catch
    ""
  end
end

function fullsignature(b::Binding)
  hasdoc(b) || return
  docs = trygetdoc(b)
  isa(docs, MD) || return
  md = flatten(docs).content
  first = length(md) > 0 ? md[1] : ""
  code =
    isa(first, Code) ? first.code :
    isa(first, Paragraph) && isa(md[1], Code) ?
      md[1].code :
      ""
  if startswith(code, string(b.var))
    split(code, "\n")[1]
  end
end

function signature(b::Binding)
  sig = fullsignature(b)
  sig == nothing && return
  replace(sig, r" -> .*$" => "")
end

function returns(b::Binding)
  r = r" -> (.*)"
  sig = fullsignature(b)
  sig == nothing && return
  if occursin(r, sig)
    ret = match(r, sig).captures[1]
    if length(ret) < 10
      ret
    end
  end
end

function description(b::Binding)
  hasdoc(b) || return
  docs = trygetdoc(b)
  isa(docs, MD) || return
  md = flatten(docs).content
  length(md) > 0 || return
  first = md[1]
  if isa(first, Code)
    length(md) < 2 && return
    first = md[2]
  end
  if isa(first, Paragraph)
    desc = plain(first)
    return length(desc) > 100 ?
      desc[1:nextind(desc, 0, 99)]*"..." :
      desc
  end
end


text:
module Tutorial7
using FEniCS

T = 1.0           # final time (increase to 10 for full problem)
num_steps = 50    # number of time steps (increase for stability)
dt = T / num_steps # time step size
mu = 1             # kinematic viscosity
rho = 1            # density

# Create mesh and define function spaces
mesh = UnitSquareMesh(16, 16)
V = VectorFunctionSpace(mesh, "P", 2)
Q = FunctionSpace(mesh, "P", 1)

inflow  = "near(x[0], 0)"
outflow = "near(x[0], 1)"
walls = "near(x[1], 0) || near(x[1], 1)"
bcu_noslip  = DirichletBC(V, Constant((0, 0)), walls)
bcp_inflow  = DirichletBC(Q, Constant(8), inflow)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)

bcu = [bcu_noslip]
bcp = [bcp_inflow, bcp_outflow]

u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

u_n = FeFunction(V)
u_  = FeFunction(V)
p_n = FeFunction(Q)
p_ = FeFunction(Q)

U   = 0.5*(u_n + u)
n   = FacetNormal(mesh)
f   = Constant((0, 0))
k   = Constant(dt)
mu  = Constant(mu)
rho = Constant(rho)

function epsilon(u)
    return sym(nabla_grad(u))
end

function sigma(u, p)
    return 2*mu*epsilon(u) - p*Identity(len(u))
end

F1 = rho*dot((u - u_n) / k, v)*dx + rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx  + inner(sigma(U, p_n), epsilon(v))*dx  + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds - dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx

a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx

A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

[apply(bc,A1) for bc in bcu]
[apply(bc,A2) for bc in bcp]

# Time-stepping
global t = 0
for n = 0:(num_steps-1)

    # Update current time
    global t += dt

    # Step 1: Tentative velocity step
    b1 = assemble(L1)
    [apply(bc,b1) for bc in bcu]

    solve(A1,vector(u_),b1)

    # Step 2: Pressure correction step
    b2 = assemble(L2)
    [apply(bc,b2) for bc in bcp]
    solve(A2, vector(p_), b2)

    # Step 3: Velocity correction step
    b3 = assemble(L3)
    solve(A3, vector(u_), b3)

    # Plot solution
    # FEniCS.Plot(u_)

    # Compute error
    u_e = Expression(("4*x[1]*(1.0 - x[1])", "0"), degree=2)
    u_e = interpolate(u_e, V)

    # Update previous solution
    assign(u_n,u_)
    assign(p_n,p_)
    #println(n)
# Hold plot
end

end#module


text:
mutable struct PoissonLeapMethod{F1,F2} <: UnsafeLeapAlgorithm
  rates::Vector{Float64}
  total_rate::Float64
  stoichiometry::SparseMatrixCSC{Int,Int}
  leap_formula::F1
  validate_leap!::F2
end

generate_events!(::PoissonLeapMethod, v, rates, s) = map!(rate -> pois_rand(rate * s), v, rates)

function update!(algorithm::PoissonLeapMethod, state, model)
  # unpack
  rates = algorithm.rates
  total_rate = algorithm.total_rate
  stoichiometry = algorithm.stoichiometry

  # update jump rates
  update_jump_rates!(algorithm, state, model)

  # update auxiliary variables
  update!(algorithm.leap_formula, state, model, stoichiometry, rates, total_rate)

  return nothing
end


text:
# wrappers for functionality provided by the CUDA toolkit

# extensions to the C language
include("intrinsics/memory_shared.jl")
include("intrinsics/indexing.jl")
include("intrinsics/synchronization.jl")
include("intrinsics/warp_vote.jl")
include("intrinsics/warp_shuffle.jl")
include("intrinsics/output.jl")
include("intrinsics/assertion.jl")
include("intrinsics/memory_dynamic.jl")
include("intrinsics/atomics.jl")
include("intrinsics/misc.jl")
if VERSION >= v"1.5-"
include("intrinsics/wmma.jl")
end

# functionality from libdevice
#
# > The libdevice library is a collection of NVVM bitcode functions that implement common
# > functions for NVIDIA GPU devices, including math primitives and bit-manipulation
# > functions. These functions are optimized for particular GPU architectures, and are
# > intended to be linked with an NVVM IR module during compilation to PTX.
include("intrinsics/math.jl")
# TODO: native mathematical functions, CUDA C programming guide" > "C language extensions"
#       https://docs.nvidia.com/cuda/cuda-math-api/group__CUDA__MATH__INTRINSIC__DOUBLE.html
#       see /path/to/cuda/include/sm_20_intrinsics.h

# functionality from libcudadevrt
#
# The libcudadevrt library is a collection of PTX bitcode functions that implement parts of
# the CUDA API for execution on the device, such as device synchronization primitives,
# dynamic kernel APIs, etc.
using CEnum
include("intrinsics/libcudadevrt_common.jl")
include("intrinsics/libcudadevrt.jl")
include("intrinsics/cooperative_groups.jl")
include("intrinsics/dynamic_parallelism.jl")


text:
module MemoryAccessesTests

include("example.jl")
include("test_smoke.jl")
include("test_arrow.jl")

end  # module MemoryAccessesTests


text:
module SemiLagrangian

using Interpolations

function generate_solver(f₀, f; interpolation_order = :Cubic)
    if interpolation_order == :Linear
        interpolation_function = BSpline(Linear())
    elseif interpolation_order == :Quadratic
        interpolation_function = BSpline(Quadratic(Periodic(OnCell())))
    elseif interpolation_order == :Cubic
        interpolation_function = BSpline(Cubic(Periodic(OnCell())))
    end

    function solve!(c)
        if interpolation_order == :Linear
            itp = interpolate([f₀; f₀[1]], interpolation_function)
        else
            itp = interpolate(f₀, interpolation_function)
        end
        etp = extrapolate(itp, Periodic(OnCell()))

        for i = 1:length(f₀)
            if 1.0 ≤ i - c ≤ length(itp)
                f[i] = itp(i - c)
            else
                f[i] = etp(i - c)
            end
        end
    end

    return solve!
end

function generate_solver(f₀, f, c; interpolation_order = :Cubic)
    if interpolation_order == :Linear
        interpolation_function = BSpline(Linear())
    elseif interpolation_order == :Quadratic
        interpolation_function = BSpline(Quadratic(Periodic(OnCell())))
    elseif interpolation_order == :Cubic
        interpolation_function = BSpline(Cubic(Periodic(OnCell())))
    end

    function solve!()
        if interpolation_order == :Linear
            itp = interpolate([f₀; f₀[1]], interpolation_function)
        else
            itp = interpolate(f₀, interpolation_function)
        end
        etp = extrapolate(itp, Periodic(OnCell()))

        for i = 1:length(f₀)
            if 1.0 ≤ i - c ≤ length(itp)
                f[i] = itp(i - c)
            else
                f[i] = etp(i - c)
            end
        end
    end

    return solve!
end

end  # module SemiLagrangian


text:
# Test big leaf canopy model
@info "Testing the big leaf model...";
@testset "CanopyLayers --- big leaf model" begin
    for FT in [Float32, Float64]
        for result in [ big_leaf_partition(FT(3.0), FT(30.0), FT(1000.0)),
                        big_leaf_partition(FT(2.0), FT(30.0), FT(1000.0)),
                        big_leaf_partition(FT(1.0), FT(30.0), FT(1000.0)) ]
            @test FT_test(result, FT);
            @test NaN_test(result);
        end
    end
end

# FT and NaN tests
println();
@info "Testing the SCOPE model...";
@testset "CanopyLayers --- SCOPE model" begin
    for FT in [Float32, Float64]
        collections = initialize_rt_module(FT; nLayer=20, LAI=3);
        angles, can, can_opt, can_rad, in_rad, leaves, rt_con, rt_dim, soil, wls = collections;
        for data_set in collections
            @test FT_test(data_set, FT);
            @test NaN_test(data_set);
        end

        # add more tests
        angles.sza = 30;
        angles.raa = 0;
        can.LAI    = 3;
        VZA        = collect(FT, -89.5:0.5:89.5);
        for VZA_ in VZA
            angles.vza = VZA_;
            canopy_geometry!(can, angles, can_opt, rt_con);
            canopy_matrices!(leaves, can_opt);
            short_wave!(can, can_opt, can_rad, in_rad, soil, rt_con);
            SIF_fluxes!(leaves, can_opt, can_rad, can, soil, wls, rt_con, rt_dim);
        end
        @test true;

        angles.sza = 48;
        angles.raa = 0;
        can.LAI    = FT(3.22);
        for raa in 0:5:360
            angles.raa = raa;
            for VZA in 0:5:85
                angles.vza = VZA;
                canopy_geometry!(can, angles, can_opt, rt_con);
                canopy_matrices!(leaves, can_opt);
                short_wave!(can, can_opt, can_rad, in_rad, soil, rt_con);
                SIF_fluxes!(leaves, can_opt, can_rad, can, soil, wls, rt_con, rt_dim);
            end
        end
        @test true;

        # test warnings
        @info "Expect warnings here!";
        warn_wls   = create_wave_length(FT, collect(FT,2100:100:2600));
        warn_inrad = create_incoming_radiation(warn_wls);
        @test true;

        # add more tests that has not been used
        can.clump_a = 0.6;
        can.clump_b = 0.1;
        CanopyLayers.clumping_factor!(can, angles);
        canopy_fluxes!(can, can_opt, can_rad, in_rad, soil, leaves[1:1], wls, rt_con);
        canopy_matrices!(leaves[1:1], can_opt);
        SIF_fluxes!(leaves[1:1], can_opt, can_rad, can, soil, wls, rt_con, rt_dim);
        thermal_fluxes!(leaves[1:1], can_opt, can_rad, can, soil, [FT(400.0)], wls);
        thermal_fluxes!(leaves[1:2], can_opt, can_rad, can, soil, [FT(400.0)], wls);
        @test true;

        # utility functions
        CanopyLayers.dcum(FT(1.1), FT(1.1), FT(1.1));
        CanopyLayers.e2phot(rand(FT,10), rand(FT,10));
        CanopyLayers.volscatt!(rand(FT,4), FT(40), FT(90), FT(40), FT(0));
        @test true;
    end
end

# indices test
println();
@info "Testing the Indicies...";
@testset "CanopyLayers --- Indicies" begin
    for FT in [Float32, Float64]
        collections = initialize_rt_module(FT; nLayer=20, LAI=3);
        angles, can, can_opt, can_rad, in_rad, leaves, rt_con, rt_dim, soil, wls = collections;
        _indices = [EVI(can_rad, wls),
                    EVI2(can_rad, wls),
                    LSWI(can_rad, wls),
                    NDVI(can_rad, wls),
                    NIRv(can_rad, wls),
                    SIF_740(can_rad, wls),
                    SIF_757(can_rad, wls),
                    SIF_771(can_rad, wls)];
        @test NaN_test(_indices);
        @test FT_test(_indices, FT);

        # out bounds warnings
        @info "Expect warnings here!";
        REF_WL(can_rad, wls, FT(3000));
        SIF_WL(can_rad, wls, FT(3000));
        @test true
    end
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
using LinearAlgebra
using BenchmarkTools
using InteractiveUtils

# dimensions
T = 26
nx = 2
nu = 1
nw = 0

# objective
ft = (x, u) -> dot(x, x) + 1.0e-1 * dot(u, u)
fT = (x, u) -> 1.0 * dot(x, x)
ct = Cost(ft, nx, nu, [t for t = 1:T-1])
cT = Cost(fT, nx, 0, [T])
obj = [ct, cT]

# model
function dynamics(z, u, w)
    mass = 1.0
    lc = 0.5
    gravity = 9.81
    damping = 0.1
    [z[2], (u[1] / (mass * lc * lc) - gravity * sin(z[1]) / lc - damping * z[2] / (mass * lc * lc))]
end

h = 0.1
function discrete_dynamics(y, x, u, w)
    y - (x + h * dynamics(y, u, w))
end

dt = Dynamics(discrete_dynamics, nx, nx, nu, nw=nw);
dyn = [dt for t = 1:T-1]
model = DynamicsModel(dyn)

# constraints
x1 = [0.0; 0.0]
xT = [π; 0.0]
x_init = Bound(nx, nu, [1], xl=x1, xu=x1)
x_goal = Bound(nx, 0, [T], xl=xT, xu=xT)
cons = ConstraintSet([x_init, x_goal], [StageConstraint()])

# problem
trajopt = TrajectoryOptimizationProblem(obj, model, cons)
s = Solver(trajopt)

z = rand(s.p.num_var)
g = zeros(s.p.num_var)
c = zeros(s.p.num_con)
j = zeros(s.p.num_jac)

# MOI methods
@code_warntype MOI.eval_objective(s.p, z)
@benchmark MOI.eval_objective($s.p, $z)
@benchmark MOI.eval_objective_gradient($s.p, $g, $z)

@code_warntype MOI.eval_constraint(s.p, c, z)
@benchmark MOI.eval_constraint($s.p, $c, $z)

@code_warntype MOI.eval_constraint_jacobian(s.p, j, z)
@benchmark MOI.eval_constraint_jacobian($s.p, $j, $z)

# initialize
z0 = 0.0 * rand(s.p.num_var)
x_interp = linear_interpolation(x1, xT, T)
for (t, idx) in enumerate(s.p.trajopt.model.idx.x)
    z0[idx] = x_interp[t]
end

initialize!(s, z0)
solve!(s)