text:
function test_transition(pb1::Union{MDP, POMDP}, pb2::Union{SparseTabularMDP, SparseTabularPOMDP})
    for s in states(pb1)
        for a in actions(pb1)
            td1 = transition(pb1, s, a)
            si = stateindex(pb1, s)
            ai = actionindex(pb1, a)
            td2 = transition(pb2, si, ai)
            for (sp, p) in weighted_iterator(td1)
                spi = stateindex(pb1, sp)
                @test pdf(td2, spi) == p
            end
        end
    end
end

function test_reward(pb1::Union{MDP, POMDP}, pb2::Union{SparseTabularMDP, SparseTabularPOMDP})
    for s in states(pb1)
        for a in actions(pb1)
            si = stateindex(pb1, s)
            ai = actionindex(pb1, a)
            @test reward(pb1, s, a) == reward(pb2, si, ai)
        end
    end
end

function test_observation(pb1::POMDP, pb2::SparseTabularPOMDP)
    for s in states(pb1)
        for a in actions(pb1)
            od1 = observation(pb1, a, s)
            si = stateindex(pb1, s)
            ai = actionindex(pb1, a)
            od2 = observation(pb2, ai, si)
            for (o, p) in weighted_iterator(od1)
                oi = obsindex(pb1, o)
                @test pdf(od2, oi) == p
            end
        end
    end
end

## MDP

mdp = RandomMDP(100, 4, 0.95)

@show_requirements SparseTabularMDP(mdp)

smdp = SparseTabularMDP(mdp)

smdp2 = SparseTabularMDP(smdp)

@test smdp2.T == smdp.T
@test smdp2.R == smdp.R
@test smdp2.discount == smdp.discount

smdp3 = SparseTabularMDP(smdp, reward = zeros(length(states(mdp)), length(actions(mdp))))
@test smdp3.T == smdp.T
@test smdp3.R != smdp.R
smdp3 = SparseTabularMDP(smdp, transition = [sparse(1:length(states(mdp)), 1:length(states(mdp)), 1.0) for a in 1:length(actions(mdp))])
@test smdp3.T != smdp.T
@test smdp3.R == smdp.R
smdp3 = SparseTabularMDP(smdp, discount = 0.5)
@test smdp3.discount != smdp.discount

@test size(transition_matrix(smdp, 1)) == (length(states(smdp)), length(states(smdp)))
@test length(reward_vector(smdp, 1)) == length(states(smdp))

gw = SimpleGridWorld()
sparsegw = SparseTabularMDP(gw)
@test isterminal(sparsegw, 101)
@inferred actions(sparsegw, 101)
@test actions(sparsegw, 101) == collect(actions(sparsegw))

## POMDP

pomdp = TigerPOMDP()

@show_requirements SparseTabularPOMDP(pomdp)

spomdp = SparseTabularPOMDP(pomdp)

spomdp2 = SparseTabularPOMDP(spomdp)

@test spomdp2.T == spomdp.T
@test spomdp2.R == spomdp.R
@test spomdp2.O == spomdp.O
@test spomdp2.discount == spomdp.discount

spomdp3 = SparseTabularPOMDP(spomdp, reward = zeros(length(states(mdp)), length(actions(mdp))))
@test spomdp3.T == spomdp.T
@test spomdp3.R != spomdp.R
spomdp3 = SparseTabularPOMDP(spomdp, transition = [sparse(1:length(states(mdp)), 1:length(states(mdp)), 1.0) for a in 1:length(actions(mdp))])
@test spomdp3.T != spomdp.T
@test spomdp3.R == spomdp.R
spomdp3 = SparseTabularPOMDP(spomdp, discount = 0.5)
@test spomdp3.discount != spomdp.discount
@test size(observation_matrix(spomdp, 1)) == (length(states(spomdp)), length(observations(spomdp)))
@test observation_matrices(spomdp) == spomdp2.O
@test transition_matrices(spomdp) == spomdp2.T
@test reward_matrix(spomdp) == spomdp2.R

## Tests

@test length(states(pomdp)) == length(states(spomdp))
@test length(actions(pomdp)) == length(actions(spomdp))
@test length(observations(pomdp)) == length(observations(spomdp))
@test statetype(spomdp) == Int64
@test actiontype(spomdp) == Int64
@test obstype(spomdp) == Int64
@test discount(spomdp) == discount(pomdp)
@test isempty(spomdp.terminal_states)

test_transition(mdp, smdp)
test_transition(pomdp, spomdp)
test_reward(pomdp, spomdp)
test_observation(pomdp, spomdp)


text:
#### Precipitation component in atmosphere model
abstract type PrecipitationModel end

export PrecipitationModel, NoPrecipitation, RainModel, RainSnowModel

eq_tends(pv::PV, m::PrecipitationModel, ::Flux{O}) where {PV, O} = ()

using ..Microphysics

vars_state(::PrecipitationModel, ::AbstractStateType, FT) = @vars()

function atmos_nodal_update_auxiliary_state!(
    ::PrecipitationModel,
    m::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function compute_gradient_flux!(
    ::PrecipitationModel,
    diffusive,
    ∇transform,
    state,
    aux,
    t,
) end
function compute_gradient_argument!(
    ::PrecipitationModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
) end

"""
    NoPrecipitation <: PrecipitationModel

No precipitation.
"""
struct NoPrecipitation <: PrecipitationModel end

precompute(::NoPrecipitation, atmos::AtmosModel, args, ts, ::Source) =
    NamedTuple()

"""
    RainModel <: PrecipitationModel

Precipitation model with rain.
"""
struct RainModel <: PrecipitationModel end

vars_state(::RainModel, ::Prognostic, FT) = @vars(ρq_rai::FT)
vars_state(::RainModel, ::Gradient, FT) = @vars(q_rai::FT)
vars_state(::RainModel, ::GradientFlux, FT) = @vars(∇q_rai::SVector{3, FT})

precompute(::RainModel, atmos::AtmosModel, args, ts, ::Source) =
    (cache = warm_rain_sources(atmos, args, ts),)

function atmos_nodal_update_auxiliary_state!(
    precip::RainModel,
    atmos::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
) end

function compute_gradient_argument!(
    precip::RainModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρinv = 1 / state.ρ
    transform.precipitation.q_rai = state.precipitation.ρq_rai * ρinv
end

function compute_gradient_flux!(
    precip::RainModel,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    diffusive.precipitation.∇q_rai = ∇transform.precipitation.q_rai
end

"""
    RainSnowModel <: PrecipitationModel

Precipitation model with rain and snow.
"""
struct RainSnowModel <: PrecipitationModel end

vars_state(::RainSnowModel, ::Prognostic, FT) = @vars(ρq_rai::FT, ρq_sno::FT)
vars_state(::RainSnowModel, ::Gradient, FT) = @vars(q_rai::FT, q_sno::FT)
vars_state(::RainSnowModel, ::GradientFlux, FT) =
    @vars(∇q_rai::SVector{3, FT}, ∇q_sno::SVector{3, FT})

precompute(::RainSnowModel, atmos::AtmosModel, args, ts, ::Source) =
    (cache = rain_snow_sources(atmos, args, ts),)

function atmos_nodal_update_auxiliary_state!(
    precip::RainSnowModel,
    atmos::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
) end

function compute_gradient_argument!(
    precip::RainSnowModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρinv = 1 / state.ρ
    transform.precipitation.q_rai = state.precipitation.ρq_rai * ρinv
    transform.precipitation.q_sno = state.precipitation.ρq_sno * ρinv
end

function compute_gradient_flux!(
    precip::RainSnowModel,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    diffusive.precipitation.∇q_rai = ∇transform.precipitation.q_rai
    diffusive.precipitation.∇q_sno = ∇transform.precipitation.q_sno
end

#####
##### Tendency specifications
#####

eq_tends(::Rain, ::RainModel, ::Flux{FirstOrder}) = (PrecipitationFlux(),)
eq_tends(::Rain, ::RainModel, ::Flux{SecondOrder}) = (Diffusion(),)

eq_tends(::AbstractPrecipitationVariable, ::RainSnowModel, ::Flux{FirstOrder}) =
    (PrecipitationFlux(),)
eq_tends(
    ::AbstractPrecipitationVariable,
    ::RainSnowModel,
    ::Flux{SecondOrder},
) = (Diffusion(),)

ext: jl


text:
# Plot coordinates of the C alpha with best occupancy
@recipe function plot(residues::AbstractVector{PDBResidue})
    ca = CAmatrix(residues)
    chains = [r.id.chain for r in residues if r.id.group == "ATOM"]
    group --> chains
    ca[:,1], ca[:,2], ca[:,3]
end


text:
#Written By Johnathan Bizzano
module Activation
    import ..@Fun
    import ...Functions.CommonStaticFunctions.Single

    export Sigmoid, Tanh, ZeroMax, ZeroMin, RELU, LeakyRELU

    Sigmoid() = Single.Sigmoid()
    Tanh() = Single.Tanh()
    ZeroMax() = Single.RMax(0)
    ZeroMin() = Single.RMin(0)
    LeakyRELU(slope) = Single.LeftValuePieceWise([Single.Linear(slope, 0), Single.Linear(1, 0)], [0])
    None() = MemoryWrapper(sett -> (x -> x))

    #Aliases
    RELU() = ZeroMax()
end


text:
mutable struct SKenCarpConstantCache{N,Tab} <: StochasticDiffEqConstantCache
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::SKenCarp,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  tab = SKenCarpTableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ,tab.c3
  nlsolver = OrdinaryDiffEq.build_nlsolver(alg,u,uprev,p,t,dt,f,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,γ,c,Val(false))
  SKenCarpConstantCache(nlsolver,tab)
end

@cache mutable struct SKenCarpCache{uType,rateType,uNoUnitsType,N,Tab,kType,randType,rateNoiseType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  z₁::uType
  z₂::uType
  z₃::uType
  z₄::uType
  k1::kType
  k2::kType
  k3::kType
  k4::kType
  atmp::uNoUnitsType
  nlsolver::N
  tab::Tab
  chi2::randType
  g1::rateNoiseType
  g4::rateNoiseType
end

u_cache(c::SKenCarpCache)    = (c.z₁,c.z₂,c.z₃,c.z₄,c.nlsolver.dz)
du_cache(c::SKenCarpCache)   = (c.nlsolver.k,c.fsalfirst)

function alg_cache(alg::SKenCarp,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  tab = SKenCarpTableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ,tab.c3
  nlsolver = OrdinaryDiffEq.build_nlsolver(alg,u,uprev,p,t,dt,f,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,γ,c,Val(true))
  fsalfirst = zero(rate_prototype)

  atmp = fill!(similar(u,uEltypeNoUnits),0)
  z₁ = similar(u); z₂ = similar(u)
  z₃ = similar(u); z₄ = nlsolver.z
  if typeof(f) <: SplitSDEFunction
    k1 = zero(u); k2 = zero(u)
    k3 = zero(u); k4 = zero(u)
  else
    k1 = nothing; k2 = nothing
    k3 = nothing; k4 = nothing
  end

  if typeof(ΔW) <: Union{SArray,Number}
    chi2 = copy(ΔW)
  else
    chi2 = zero(ΔW)
  end

  g1 = zero(noise_rate_prototype); g4 = zero(noise_rate_prototype)

  SKenCarpCache{typeof(u),typeof(rate_prototype),typeof(atmp),typeof(nlsolver),
              typeof(tab),typeof(k1),typeof(chi2),typeof(g1)}(
              u,uprev,fsalfirst,z₁,z₂,z₃,z₄,k1,k2,k3,k4,atmp,nlsolver,tab,chi2,g1,g4)
end


text:
module Coverage

using CoverageTools
using LibGit2

export FileCoverage
export LCOV
export analyze_malloc
export amend_coverage_from_src!
export clean_file
export clean_folder
export get_summary
export merge_coverage_counts
export process_cov
export process_file
export process_folder

const CovCount = CoverageTools.CovCount
const FileCoverage = CoverageTools.FileCoverage
const amend_coverage_from_src! = CoverageTools.amend_coverage_from_src!
const clean_file = CoverageTools.clean_file
const clean_folder = CoverageTools.clean_folder
const get_summary = CoverageTools.get_summary
const iscovfile = CoverageTools.iscovfile
const merge_coverage_counts = CoverageTools.merge_coverage_counts
const process_cov = CoverageTools.process_cov
const process_file = CoverageTools.process_file
const process_folder = CoverageTools.process_folder

include("coveralls.jl")
include("codecovio.jl")
include("lcov.jl")
include("memalloc.jl")
include("parser.jl")

end # module


text:
"""
	$(SIGNATURES)

Stores a grid by age.
Paramterized by
- V: type of value (e.g. Float64)
- T: no of periods (Integer)
"""
struct GridByAge{V, T}
    # Grid by age
    gridV :: Vector{Vector{V}}
end

"""
Constructor using no of grid points and ranges.
Grid starts with `k0`, followed by grids of length `nPointsV` in range `[kMinV, kMaxV]`.

# Arguments
- spacing
- intervalRatio: only matters for power spaced grid.
"""
function GridByAge(k0 :: V, nPointsV :: Vector{T2},
    kMinV :: AbstractVector{V}, kMaxV :: AbstractVector{V};
    spacing :: Symbol = :linear, intervalRatio = 1.1) where  {V, T2 <: Integer}

    @assert all(nPointsV .>= 2)
    @assert all(kMinV .< kMaxV)
    if spacing == :log
        @assert all(kMinV .> 0.000001)
    end

    T = length(kMinV) + 1;
    gridV = Vector{Vector{V}}(undef, T);
    gridV[1] = [k0];

    for t = 2 : T
        if spacing == :power
            gridV[t] = make_power_spaced_grid(kMinV[t-1], kMaxV[t-1], nPointsV[t-1],
                intervalRatio);
        elseif spacing == :log
            gridV[t] = exp.(LinRange{V}(log(kMinV[t-1]), log(kMaxV[t-1]), nPointsV[t-1]));
    	elseif spacing == :linear
            gridV[t] = LinRange{V}(kMinV[t-1], kMaxV[t-1], nPointsV[t-1]);
        else
            error("Invalid spacing: $spacing")
		end
    end

    return GridByAge{V,T}(gridV)
end

function make_power_spaced_grid(kMin, kMax, nPoints, intervalRatio)
    g = PowerSpacedGrid(kMin, kMax, nPoints, intervalRatio);
    return CommonLH.grid(g)
end

# function make_linspaced_grid(kMin, kMax, nPoints, intervalRatio)
#     g = LinSpaced

## Validate grid properties
function validate(g :: GridByAge, nV :: Vector{T2}, kMinV :: Vector{Double},
    kMaxV :: Vector{Double})  where  T2 <: Integer

    T = length(nV);
    isValid =  (T == n_periods(g));
    for t = 1 : T
        gridV = grid_t(g, t);
        isValid = isValid  &&  (gridV[1] ≈ kMinV[t])  &&
            (gridV[end] ≈ kMaxV[t])  &&  (length(gridV) == nV[t])
    end
    return isValid
end

function make_test_grid(T :: Integer)
    nV = collect(3 .+ (1 : T));
    kMinV = collect(range(-2.0, -3.0, length = T));
    kMaxV = kMinV .+ 4.0;

    gridV = Vector{Vector{Double}}(undef, T);
    for t = 1 : T
        gridV[t] = LinRange{Double}(kMinV[t], kMaxV[t], nV[t]);
    end
    kg = GridByAge(gridV);

    @assert CollegeStrat.validate(kg, nV, kMinV, kMaxV)
    return kg
end

# ------------------  Retrieve

"""
	$(SIGNATURES)

Length of grid in periods.
"""
function n_periods(g :: GridByAge)
    return length(g.gridV)
end

"""
	$(SIGNATURES)

Return the date `t` grid.
"""
function grid_t(g :: GridByAge, t :: Integer)
    return g.gridV[t]
end

"""
	$(SIGNATURES)

Return points on the date `t` grid, selected by indices `idxV`.
"""
function grid_points(g :: GridByAge, t :: T1, idxV) where
    T1 <: Integer

    return @inbounds g.gridV[t][idxV]
end

grid_points(g :: GridByAge, t :: T1) where T1 <: Integer =
    g.gridV[t];

## Number of grid points at a t
function n_points(g :: GridByAge, t :: Integer)
    return length(g.gridV[t])
end

## Round to nearest grid point
function round_to_grid(g :: GridByAge, t :: Integer, inM :: Array{Double})
    return round_to_grid(inM,  g.gridV[t])
end

"Scalar input"
function round_to_grid(g :: GridByAge, t :: Integer, x :: Double)
    # Need to call function for array and return first element
    return round_to_grid([x], g.gridV[t])[1]
end

## ---------  Display

# Simple summary of the grid for each age
function show_grids(g :: GridByAge; io = stdout)
    for t = 1 : n_periods(g)
        show_grid(g, t; io = io);
    end
end

function show_grid(g :: GridByAge, t :: Integer; io = stdout)
    np = n_points(g, t);
    print(io,  "Grid for t=$t ($np points):  ");
    gridV = round.(grid_t(g, t), digits = 2);
    if np < 7
        println(io, gridV);
    else
        print(io, gridV[1 : 3]);
        print(io, " ... ");
        println(io, gridV[(np-2) : np]);
    end
end

## --------  Copied from MatrixLH

"""
round_to_grid()

Round a matrix to the nearest points on a grid
"""
function round_to_grid(xM :: Array{T}, gridV :: Vector{T}) where T <: AbstractFloat
    idxM = similar(xM, Int);
    for (i1, x) in enumerate(xM)
        iMin = findmin(abs.(gridV .- x));
        idxM[i1] = iMin[2];
    end
    return idxM
end

# ----------------


text:
# Short

@doc "f_1" ->
f_1(x) = x

@doc """f_2""" ->
f_2(x) = x

# Full

@doc "f_3" ->
function f_3(x)
end

@doc """f_4""" ->
function f_4(x)
end

# Complex type signatures.

@doc "f_5" ->
f_5{T <: Integer}(a::T, b::Vector{T}, c::Complex{T}) = ()

@doc "f_6" ->
f_6(::Type{String}, a::Array{Int8,1}) = ()

@doc "f_7" ->
f_7(n::@compat(Union{Int64,Int32})) = ()

@doc "f_8" ->
f_8{T<:Integer}(n::T, k::T) = ()

@doc "f_9" ->
f_9{T<:Real}(::Complex{T}, x::Real) = ()

@doc "f_10" ->
f_10{N}(dest::Array, src::Array, I::NTuple{N,@compat(Union{Int,AbstractVector})}...) = ()

@doc "f_11" ->
f_11(V, indexes::NTuple, dims::Dims, linindex::UnitRange{Int}) = ()

@doc "f_12" ->
f_12{T}(::Type{typejoin()}, ::Type{T}) = ()

@doc "f_13" ->
f_13{T<:@compat(Union{Float32,Float64})}(::Type{T},i::Integer) = ()

@doc "f_14" ->
f_14{T<:AbstractFloat}(f, a::T,b::T,c::T...;
                       abstol=zero(T),
                       reltol=sqrt(eps(T)),
                       maxevals=10^7,
                       order=7,
                       norm=vecnorm
                       ) = ()

@doc "f_15" ->
f_15(b, z, m = length(b) + 1) = ()

@doc "f_16" ->
f_16{T<:Integer}(::Type{T}, x::AbstractFloat; tol::Real = eps(x)) = ()

@doc "f_17" ->
f_17(io::IO, itr::AbstractArray, op, delim, cl, delim_one, compact = false) = ()

@doc "f_18" ->
f_18(io::IO, items, sep, indent::Int, prec::Int = 0, enclose_operators::Bool = false) = ()

@doc "f_19" ->
f_19{T,N,P,IV}(V::SubArray{T,N,P,IV}, m::Real) = ()

@doc "f_20" ->
f_20{T,N,P,IV}(V::SubArray{T,N,P,IV}, I::AbstractArray{Bool,N}) = ()

@doc "f_21" ->
f_21{T,N,P,IV}(V::SubArray{T,N,P,IV}, v, I::@compat(Union{Real,AbstractVector})...) = ()

@doc "f_22" ->
f_22(::Type{Ptr{Int8}}, s::String) = ()

@doc "f_23" ->
f_23(T, matvecA::Function, matvecB::Function, solveSI::Function, n::Integer,
     sym::Bool, cmplx::Bool, bmat::String,
     nev::Integer, ncv::Integer, which::String,
     tol::Real, maxiter::Integer, mode::Integer, v0::Vector) = ()

@doc "f_24" ->
f_24{T}(::Matrix{T}, A::Bidiagonal{T}) = ()

@doc "f_25" ->
f_25(f::ANY, t::Array, i, lim::Integer, matching::Array{Any,1}) = ()

@doc "f_26" ->
f_26(f::Function, types::tup(Vararg{Type})) = ()

@doc "f_27" ->
f_27(f::Base.Callable, ::tup(), ts::Tuple...) = ()

@doc "f_28" ->
f_28(x::tup(Any, Vararg{Any})) = ()

@doc "f_29" ->
f_29(f::Base.Callable, t::tup(Any,Any), s::tup(Any,Any)) = ()

# Generic functions.

@doc "g_f_1" ->
f_1

@doc """g_f_2""" ->
f_2

# Star syntax

@doc "f_30" ->
f_30() = ()

@doc+ "g_f_30" ->
f_30(x) = ()

@doc meta("f_31", returns = Bool) ->
f_31() = true

@doc+ meta("g_f_31", returns = Int) ->
f_31(x) = true

# doc_str, doc_mstr macros.

@doc doc"f_36" ->
f_36() = ()

@doc doc"""f_37""" ->
f_37() = ()

# Qualified docs.

type FooType end

@doc "f_38" ->
Base.Pkg.add(::FooType) = 1


text:
@testset "Local Level With Explanatory Model" begin
    @test has_fit_methods(LocalLevelExplanatory)
    y = [
    0.8297606696482265
    0.5151262659173474
    0.6578708499069135
    0.15479984562134974
    0.41836893137914033
    0.46381576757801835
    0.18373491281112986
    0.09712212189897196
    0.9696168042329114
    0.8623172534114631
    0.4312662075760265
    0.17938380947382027
    0.6265684534572271
    0.6395165239895426
    0.370025441453405
    ]
    X = [
    0.221409   0.646901
    0.943735   0.163734
    0.19864    0.620721
    0.462932   0.876136
    0.88467    0.673716
    0.387496   0.415591
    0.29221    0.989594
    0.487926   0.721127
    0.0321253  0.318292
    0.742561   0.582234
    0.664695   0.749821
    0.791058   0.423212
    0.0652515  0.955714
    0.842513   0.788755
    0.198105   0.79092
    ]
    model = LocalLevelExplanatory(y, X)
    fit!(model)

    # forecasting
    # For a fixed forecasting explanatory the variance must not decrease
    forec = forecast(model, ones(5, 2))
    @test monotone_forecast_variance(forec)
    kf = kalman_filter(model)
    ks = kalman_smoother(model)
end


text:
module Priority
using PyCall,Distributions,DataStructures
@pyimport networkx as nx
@pyimport matplotlib.pyplot as plt
# package code goes here
export average_bw,reach_stable_time,order_parameter,percolation,critical_threshold,correlation
export correlation_length,betweenness_distribution,average_betweenness,pathlength_and_travelingtime_distribution
export average_pathlength_and_travelingtime
using Priority_Scalefree
end # module


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

using Test, Tensors

let dtime = 0.25,
    nu = 0.3,
    R = 100.0,
    parameters = PerfectPlasticParameterState(youngs_modulus=200.0e3,
                                              poissons_ratio=nu,
                                              yield_stress=R),
    mat = PerfectPlastic(parameters=parameters),
    tostrain(vec) = fromvoigt(Symm2, vec; offdiagscale=2.0),
    tostress(vec) = fromvoigt(Symm2, vec),
    uniaxial_stress(sigma) = tostress([sigma, 0, 0, 0, 0, 0]),
    stresses_expected = [uniaxial_stress(R / 2),
                         uniaxial_stress(R),
                         uniaxial_stress(R),
                         uniaxial_stress(R / 2),
                         uniaxial_stress(-R)],
    dstrain11 = 1e-3*dtime,
    strains_expected = [tostrain(dstrain11*[1.0, -nu, -nu, 0.0, 0.0, 0.0]),
                        tostrain(dstrain11*[2, -2*nu, -2*nu, 0.0, 0.0, 0.0]),
                        tostrain(dstrain11*[3, -2*nu - 0.5, -2*nu - 0.5, 0.0, 0.0, 0.0]),
                        tostrain(dstrain11*[2, -nu - 0.5, -nu - 0.5, 0.0, 0.0, 0.0]),
                        tostrain(dstrain11*[-2, 2*nu, 2*nu, 0.0, 0.0, 0.0])],
    dtimes = [dtime, dtime, dtime, dtime, 1.0],
    dstrains11 = dstrain11*[1.0, 1.0, 1.0, -1.0, -4.0]

    for i in 1:length(dtimes)
        dstrain11 = dstrains11[i]
        dtime = dtimes[i]
        uniaxial_increment!(mat, dstrain11, dtime)
        update_material!(mat)
        @test isapprox(mat.variables.stress, stresses_expected[i])
        @test isapprox(mat.drivers.strain, strains_expected[i])
    end
end


text:
function Base.show(io::IO, dfr::DataFrameRow;
                   allcols::Bool = !get(io, :limit, false),
                   splitcols = get(io, :limit, false),
                   rowlabel::Symbol = :Row,
                   missingstring::AbstractString = "missing")
    r, c = parentindices(dfr)
    print(io, "DataFrameRow")
    _show(io, view(parent(dfr), [r], c), allcols=allcols, splitcols=splitcols,
          rowlabel=rowlabel, summary=false, missingstring=missingstring, rowid=r)
end

Base.show(dfr::DataFrameRow;
          allcols::Bool = !get(io, :limit, true),
          splitcols = get(io, :limit, true),
          rowlabel::Symbol = :Row,
          missingstring::AbstractString = "missing") =
    show(stdout, dfr, allcols=allcols, splitcols=splitcols, rowlabel=rowlabel, missingstring=missingstring)


text:
# using Pkg
# Pkg.activate(@__DIR__)

# For linear algebra
using StaticArrays: SVector
using StaticArrays
using LinearAlgebra
using SparseArrays: SparseMatrixCSC, sparse
# Fast Cholesky
using KoLesky

# multiple dispatch
import IterativeSolvers: cg!
import LinearAlgebra: mul!, ldiv!
import Base: size

# autoDiff
using ForwardDiff

# parser
using ArgParse

# logging
using Logging

# profile
# using Profile
# using BenchmarkTools

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--kernel"
            arg_type = String
            default = "Matern7half"
        "--sigma"
            help = "lengthscale"
            arg_type = Float64
            default = 0.3
        "--h"
            help = "grid size"
            arg_type = Float64
            default = 0.0125
        "--nugget"
            arg_type = Float64
            default = 1e-8
        "--GNsteps"
            arg_type = Int
            default = 3
        "--rho_big"
            arg_type = Float64
            default = 3.0
        "--rho_small"
            arg_type = Float64
            # default = -log(5*h_grid)
            default = 3.0
        "--k_neighbors"
            arg_type = Int
            default = 3
        "--compare_exact"
            arg_type = Bool
            default = false
    end
    return parse_args(s)
end

## PDEs type
abstract type AbstractPDEs end
struct MongeAmpere2d{TΩ} <: AbstractPDEs
    # minimal surface equation in [Ω[1,1],Ω[2,1]]*[Ω[1,2],Ω[2,2]]
    Ω::TΩ
    bdy::Function
    rhs::Function
end

## sample points
function sample_points_rdm(eqn::MongeAmpere2d, N_domain, N_boundary)
    Ω = eqn.Ω
    x1l = Ω[1,1]
    x1r = Ω[2,1]
    x2l = Ω[1,2]
    x2r = Ω[2,2]

    X_domain = hcat(rand(Float64,(N_domain, 1))*(x1r-x1l).+x1l,rand(Float64,(N_domain, 1))*(x2r-x2l).+x2l)

    N_bd_each=convert(Int64, N_boundary/4)
    if N_boundary != 4 * N_bd_each
        println("[sample points] N_boundary not divided by 4, replaced by ", 4 * N_bd_each)
        N_boundary = 4 * N_bd_each
    end

    X_boundary = zeros((N_boundary, 2))
    # bottom face
    X_boundary[1:N_bd_each, :] = hcat((x1r-x1l)*rand(Float64,(N_bd_each,1)).+x1l, x2l*ones(N_bd_each))
    # right face
    X_boundary[N_bd_each+1:2*N_bd_each, :] = hcat(x1r*ones(N_bd_each),(x2r-x2l)*rand(Float64,(N_bd_each,1)).+x2l)
    # top face
    X_boundary[2*N_bd_each+1:3*N_bd_each, :] = hcat((x1r-x1l)*rand(Float64,(N_bd_each,1)).+x1l, x2r*ones(N_bd_each))

    # left face
    X_boundary[3*N_bd_each+1:N_boundary, :] = hcat(x1l*ones(N_bd_each), (x2r-x2l)*rand(Float64,(N_bd_each,1)).+x2l)
    return X_domain', X_boundary'
end
function sample_points_grid(eqn::MongeAmpere2d, h_in, h_bd)
    Ω = eqn.Ω
    x1l = Ω[1,1]
    x1r = Ω[2,1]
    x2l = Ω[1,2]
    x2r = Ω[2,2]
    x = x1l + h_in:h_in:x1r-h_in
    y = x2l + h_in:h_in:x2r-h_in
    X_domain = reduce(hcat,[[x[i], y[j]] for i in 1:length(x) for j in 1:length(x)])

    l = length(x1l:h_bd:x1r-h_bd)
    X_boundary = vcat([x1l:h_bd:x1r-h_bd x2l*ones(l)], [x1r*ones(l) x2l:h_bd:x2r-h_bd], [x1r:-h_bd:x1l+h_bd x2r*ones(l)], [x1l*ones(l) x2r:-h_bd:x1l+h_bd])
    return X_domain, X_boundary'
end

## exact algorithms
# assemby Gram matrices
function get_Gram_matrices(eqn::MongeAmpere2d, cov::KoLesky.AbstractCovarianceFunction, X_domain, X_boundary, sol_xx_now, sol_xy_now,sol_yy_now)
    d = 2
    N_domain = size(X_domain,2)
    N_boundary = size(X_boundary,2)

    # get linearized PDEs correponding measurements
    meas_δ_bd = [KoLesky.PointMeasurement{d}(SVector{d,Float64}(X_boundary[:,i])) for i = 1:N_boundary]
    meas_∂∂ = [KoLesky.∂∂PointMeasurement{Float64,d}(SVector{d,Float64}(X_domain[:,i]), sol_yy_now[i], -2*sol_xy_now[i], sol_xx_now[i]) for i = 1:N_domain]
    meas_δ_int = [KoLesky.PointMeasurement{d}(SVector{d,Float64}(X_domain[:,i])) for i = 1:N_domain]

    Theta_train = zeros(N_domain+N_boundary,N_domain+N_boundary)

    measurements = Vector{Vector{<:KoLesky.AbstractPointMeasurement}}(undef,2)
    measurements[1] = meas_δ_bd; measurements[2] = meas_∂∂
    cov(Theta_train, reduce(vcat,measurements))

    measurement_test = Vector{Vector{<:KoLesky.AbstractPointMeasurement}}(undef,4)
    measurement_test[1] = meas_δ_int;
    measurement_test[2] = [KoLesky.∂∂PointMeasurement{Float64,d}(SVector{d,Float64}(X_domain[:,i]), 1.0, 0.0, 0.0) for i = 1:N_domain]
    measurement_test[3] = [KoLesky.∂∂PointMeasurement{Float64,d}(SVector{d,Float64}(X_domain[:,i]), 0.0, 1.0, 0.0) for i = 1:N_domain]
    measurement_test[4] = [KoLesky.∂∂PointMeasurement{Float64,d}(SVector{d,Float64}(X_domain[:,i]), 0.0, 0.0, 1.0) for i = 1:N_domain]

    Theta_test = zeros(4*N_domain,N_domain+N_boundary)
    cov(Theta_test, reduce(vcat,measurement_test), reduce(vcat,measurements))

    return Theta_train, Theta_test

end
# iterative GPR
function iterGPR_exact(eqn, cov, X_domain, X_boundary, sol_init, sol_init_xx, sol_init_xy, sol_init_yy, nugget, GNsteps)
    N_domain = size(X_domain,2)
    N_boundary = size(X_boundary,2)
    # get the rhs and bdy data
    rhs = [eqn.rhs(X_domain[:,i]) for i in 1:N_domain]
    bdy = [eqn.bdy(X_boundary[:,i]) for i in 1:N_boundary]
    v = sol_init
    v_xx = sol_init_xx
    v_xy = sol_init_xy
    v_yy = sol_init_yy

    for _ in 1:GNsteps
        Theta_train, Theta_test = get_Gram_matrices(eqn, cov, X_domain, X_boundary, v_xx, v_xy, v_yy)
        rhs_now = vcat(bdy, rhs+v_xx.*v_yy-v_xy.^2)

        v_all = Theta_test*(((Theta_train+nugget*diagm(diag(Theta_train))))\rhs_now)
        v = v_all[1:N_domain]

        v_xx = v_all[N_domain+1:2*N_domain]
        v_xy = v_all[2*N_domain+1:3*N_domain]
        v_yy = v_all[3*N_domain+1:4*N_domain]

        # in case that uxx or uyy becomes negative
        # eps = 1e-3
        # @inbounds for i in 1:N_domain
        #     if v_xx[i] < eps
        #         v_xx[i] = eps
        #     end
        #     if v_yy[i] < eps
        #         v_yy[i] = eps
        #     end
        #     if v_xy[i]^2 > v_xx[i]*v_yy[i]
        #         v_xy[i] = 0.5*sqrt(v_xx[i]*v_yy[i])*sign(v_xy[i])
        #     end
        # end

    end
    return v, v_xx, v_xy, v_yy
end

## algorithm using Kolesky and pcg
# struct that stores the factor of Theta_train
abstract type implicit_mtx end
struct approx_Theta_train{Tv,Ti,Tmtx<:SparseMatrixCSC{Tv,Ti}} <: implicit_mtx
    P::Vector{Ti}
    U::Tmtx
    L::Tmtx
    ∂11_coefs::Vector{Tv}
    ∂12_coefs::Vector{Tv}
    ∂22_coefs::Vector{Tv}
    N_boundary::Ti
    N_domain::Ti
end
struct precond_Theta_train{Tv,Ti,Tmtx<:SparseMatrixCSC{Tv,Ti}} <: implicit_mtx
    P::Vector{Ti}
    U::Tmtx
    L::Tmtx
end

function size(A::approx_Theta_train, num)
    return size(A.U,num)
end

function mul!(x, Θtrain::approx_Theta_train, b)
    N_domain = Θtrain.N_domain
    N_boundary = Θtrain.N_boundary
    @views temp = vcat(b[1:N_boundary], zeros(N_domain), Θtrain.∂11_coefs .* b[N_boundary+1:end], Θtrain.∂22_coefs .*b[N_boundary+1:end], -2*Θtrain.∂12_coefs.*b[N_boundary+1:end])
    temp[Θtrain.P] = Θtrain.L\(Θtrain.U\temp[Θtrain.P])

    @views x[1:Θtrain.N_boundary] = temp[1:Θtrain.N_boundary]
    @views x[Θtrain.N_boundary+1:end] .= Θtrain.∂11_coefs.* temp[N_boundary+N_domain+1:N_boundary+2*N_domain] .+ Θtrain.∂22_coefs.*temp[N_boundary+2*N_domain+1:N_boundary+3*N_domain] .- 2*Θtrain.∂12_coefs.*temp[N_boundary+3*N_domain+1:N_boundary+4*N_domain]
end

function ldiv!(x, precond::precond_Theta_train, b)
    x[precond.P] = precond.U*(precond.L*b[precond.P])
end

function iterGPR_fast_pcg(eqn, cov, X_domain, X_boundary, sol_init, sol_init_xx, sol_init_xy, sol_init_yy, nugget, GNsteps; ρ_big =4.0, ρ_small=6.0, k_neighbors = 4, lambda = 1.5, alpha = 1.0)

    @info "[Algorithm]: iterative GPR + fast cholesky factorization + pcg"
    N_domain = size(X_domain,2)
    N_boundary = size(X_boundary,2)

    # get the rhs and bdy data
    rhs = [eqn.rhs(X_domain[:,i]) for i in 1:N_domain]
    bdy = [eqn.bdy(X_boundary[:,i]) for i in 1:N_boundary]

    v = sol_init
    v_xx = sol_init_xx
    v_xy = sol_init_xy
    v_yy = sol_init_yy

    # form the fast Cholesky part that can be used to compute mtx-vct mul for Theta_test
    d = 2
    meas_δ_bd = [KoLesky.PointMeasurement{d}(SVector{d,Float64}(X_boundary[:,i])) for i = 1:N_boundary]
    meas_δ_int = [KoLesky.PointMeasurement{d}(SVector{d,Float64}(X_domain[:,i])) for i = 1:N_domain]
    meas_∂11_int = [KoLesky.∂∂PointMeasurement{Float64,d}(SVector{d,Float64}(X_domain[:,i]), 1.0, 0.0, 0.0) for i = 1:N_domain]
    meas_∂22_int = [KoLesky.∂∂PointMeasurement{Float64,d}(SVector{d,Float64}(X_domain[:,i]), 0.0, 0.0, 1.0) for i = 1:N_domain]
    meas_∂12_int = [KoLesky.∂∂PointMeasurement{Float64,d}(SVector{d,Float64}(X_domain[:,i]), 0.0, 1.0, 0.0) for i = 1:N_domain]

    measurements = Vector{Vector{<:KoLesky.AbstractPointMeasurement}}(undef,5)
    measurements[1] = meas_δ_bd
    measurements[2] = meas_δ_int
    measurements[3] = meas_∂11_int
    measurements[4] = meas_∂22_int
    measurements[5] = meas_∂12_int

    @info "[Big Theta: implicit factorization] time"
    @time implicit_bigΘ = KoLesky.ImplicitKLFactorization_DiracsFirstThenUnifScale(cov, measurements, ρ_big, k_neighbors; lambda = lambda, alpha = alpha)

    @info "[Big Theta: explicit factorization] time"
    @time explicit_bigΘ = KoLesky.ExplicitKLFactorization(implicit_bigΘ; nugget = nugget)
    U_bigΘ = explicit_bigΘ.U
    L_bigΘ = sparse(U_bigΘ')
    P_bigΘ = explicit_bigΘ.P

    Θtrain = approx_Theta_train(P_bigΘ, U_bigΘ, L_bigΘ,zeros(N_domain),zeros(N_domain),zeros(N_domain),N_boundary,N_domain)

    implicit_factor = nothing

    for step in 1:GNsteps

        @info "[Current GN step] $step"
        # get Cholesky of Theta_train

        meas_∂∂ = [KoLesky.∂∂PointMeasurement{Float64,d}(SVector{d,Float64}(X_domain[:,i]), v_yy[i], -2*v_xy[i], v_xx[i]) for i = 1:N_domain]

        measurements = Vector{Vector{<:KoLesky.AbstractPointMeasurement}}(undef,2)
        measurements[1] = meas_δ_bd; measurements[2] = meas_∂∂

        @info "[Theta Train: implicit factorization] time"
        @time if implicit_factor === nothing
            implicit_factor = KoLesky.ImplicitKLFactorization(cov, measurements, ρ_small, k_neighbors; lambda = lambda, alpha = alpha)
        else
            implicit_factor.supernodes.measurements .= reduce(vcat, collect.(measurements))[implicit_factor.P]
        end
        @info "[Theta Train: explicit factorization] time"
        @time explicit_factor = KoLesky.ExplicitKLFactorization(implicit_factor; nugget = nugget)

        U = explicit_factor.U
        L = sparse(U')
        P = explicit_factor.P
        rhs_now = vcat(bdy, rhs+v_xx.*v_yy-v_xy.^2)

        # use the approximate solution as the initial point for the pCG iteration
        Θinv_rhs = U*(L*rhs_now[P])
        Θinv_rhs[P] = Θinv_rhs

        # pcg step for Theta_train\rhs
        Θtrain.∂11_coefs .= v_yy
        Θtrain.∂12_coefs .= v_xy
        Θtrain.∂22_coefs .= v_xx
        precond = precond_Theta_train(P,U,L)

        @info "[pcg started]"
        @time Θinv_rhs, ch = cg!(Θinv_rhs, Θtrain, rhs_now; Pl = precond, log=true)
        @info "[pcg finished], $ch"

        # get approximated sol_now = Theta_test * Θinv_rhs
        tmp = vcat(Θinv_rhs[1:N_boundary], zeros(N_domain), v_yy .* Θinv_rhs[N_boundary+1:end], v_xx .*Θinv_rhs[N_boundary+1:end], -2*v_xy.*Θinv_rhs[N_boundary+1:end])
        tmp[P_bigΘ] = L_bigΘ\(U_bigΘ\tmp[P_bigΘ])
        @views v = tmp[N_boundary+1:N_domain+N_boundary]
        @views v_xx = tmp[N_domain+N_boundary+1:N_boundary+2*N_domain]
        @views v_yy = tmp[2*N_domain+N_boundary+1:N_boundary+3*N_domain]
        @views v_xy = tmp[3*N_domain+N_boundary+1:N_boundary+4*N_domain]
    end
    @info "[solver finished]"
    return v
end

function main(args)
    Ω = [[0,1] [0,1]]
    # ground truth solution
    function fun_u(x)
        return exp(0.5*sum((x.-0.5).^2))
    end

    function H_u(x)
        hessian = ForwardDiff.hessian(x -> fun_u(x),x)
        return @SVector [hessian[1],hessian[2],hessian[4]]
    end

    # right hand side
    function fun_rhs(x)
        H = H_u(x)
        return H[1]*H[3]-H[2]^2
    end

    # function fun_rhs(x)
    #     return (1+sum((x.-0.5).^2))*exp(sum((x.-0.5).^2))
    # end

    # boundary value
    function fun_bdy(x)
        return fun_u(x)
    end

    @info "[solver started] Monge-Ampere equation"
    @info "[equation] Monge-Ampere u_xxu_yy - (u_xy)² = f"
    eqn = MongeAmpere2d(Ω,fun_bdy,fun_rhs)

    h_in = args.h::Float64; h_bd = h_in
    X_domain, X_boundary = sample_points_grid(eqn, h_in, h_bd)
    # X_domain, X_boundary = sample_points_rdm(eqn, 900, 124)
    N_domain = size(X_domain,2)
    N_boundary = size(X_boundary,2)
    @info "[sample points] grid size $h_in"
    @info "[sample points] N_domain is $N_domain, N_boundary is $N_boundary"

    lengthscale = args.sigma
    if args.kernel == "Matern5half"
        cov = KoLesky.MaternCovariance5_2(lengthscale)
    elseif args.kernel == "Matern7half"
        cov = KoLesky.MaternCovariance7_2(lengthscale)
    elseif args.kernel == "Matern9half"
        cov = KoLesky.MaternCovariance9_2(lengthscale)
    elseif args.kernel == "Matern11half"
        cov = KoLesky.MaternCovariance11_2(lengthscale)
    elseif args.kernel == "Gaussian"
        cov = KoLesky.GaussianCovariance(lengthscale)
    end
    @info "[kernel] choose $(args.kernel), lengthscale $lengthscale\n"

    nugget = args.nugget::Float64
    @info "[nugget] $nugget"

    GNsteps_approximate = args.GNsteps::Int
    @info "[total GN steps] $GNsteps_approximate"

    ρ_big = args.rho_big::Float64
    ρ_small = args.rho_small::Float64
    k_neighbors = args.k_neighbors::Int
    @info "[Fast Cholesky] ρ_big = $ρ_big, ρ_small = $ρ_small, k_neighbors = $k_neighbors"

    # sol_init = [fun_u(X_domain[:,i]) for i in 1:N_domain] # initial solution
    # sol_init_derivatives = [H_u(X_domain[:,i]) for i in 1:N_domain]
    # sol_init_xx = [sol_init_derivatives[i][1] for i in 1:N_domain]
    # sol_init_xy = [sol_init_derivatives[i][2] for i in 1:N_domain]
    # sol_init_yy = [sol_init_derivatives[i][3] for i in 1:N_domain]

    sol_init = [0.0 for i in 1:N_domain] # initial solution
    sol_init_xx = [1.0 for i in 1:N_domain]
    sol_init_xy = [0.0 for i in 1:N_domain]
    sol_init_yy = [1.0 for i in 1:N_domain]

    truth = [fun_u(X_domain[:,i]) for i in 1:N_domain]

    fast_solve() = @time iterGPR_fast_pcg(eqn, cov, X_domain, X_boundary, sol_init, sol_init_xx, sol_init_xy, sol_init_yy, nugget, GNsteps_approximate; ρ_big = ρ_big, ρ_small = ρ_small, k_neighbors=k_neighbors);
    sol = fast_solve()

    pts_accuracy = sqrt(sum((truth-sol).^2)/N_domain)
    @info "[L2 accuracy: pCG method] $pts_accuracy"
    pts_max_accuracy = maximum(abs.(truth-sol))
    @info "[Linf accuracy: pCG method] $pts_max_accuracy"

    # @info  sqrt(sum((truth-sol_init).^2)/N_domain)
    # @info maximum(abs.(truth-sol_init))
    if args.compare_exact
        GNsteps_exact = 2
        @info "[comparison: exact method]"
        @time sol_exact, sol_xx, sol_xy, sol_yy = iterGPR_exact(eqn, cov, X_domain, X_boundary, sol_init, sol_init_xx, sol_init_xy, sol_init_yy, nugget, GNsteps_exact)
        pts_accuracy_exact = sqrt(sum((truth-sol_exact).^2)/N_domain)
        @info "[L2 accuracy: exact method] $pts_accuracy_exact"
        pts_max_accuracy_exact = maximum(abs.(truth-sol_exact))
        @info "[Linf accuracy: exact method] $pts_max_accuracy_exact"

        # @show sol_init_xx.*sol_yy .+ sol_init_yy.*sol_xx .- 2*sol_xy.*sol_init_xy .- 2*[eqn.rhs(X_domain[:,i]) for i in 1:N_domain]
        # @show sol_xx
        # @show sol_yy
    end

end

args = parse_commandline()
args = (; (Symbol(k) => v for (k,v) in args)...) # Named tuple from dict
main(args)


text:
using Ekztazy
using Test

include("fulltest.jl")


text:
using Logging
using SpinGlassEngine
using MetaGraphs
using LightGraphs
disable_logging(LogLevel(1))

function solve_instance(instance_path, nrows, ncols; β, sweeps, bond_dim, max_states)
    ig = ising_graph(instance_path)
    fg = factor_graph(
        ig,
        cluster_assignment_rule=super_square_lattice((nrows, ncols, 8))
    )
    network = PEPSNetwork(
        nrows,
        ncols,
        fg,
        rotation(0),
        β=β,
        bond_dim=bond_dim,
        sweeps=sweeps,
    )
    low_energy_spectrum(network, max_states, merge_branches(network, 1.0))
end

function solve_instance(instance_path, nrows, ncols, t; β, sweeps, bond_dim, max_states)
    ig = ising_graph(instance_path)
    fg = factor_graph(
        ig,
        cluster_assignment_rule=super_square_lattice((nrows, ncols, t))
    )
    network = FusedNetwork(
        nrows,
        ncols,
        fg,
        rotation(0),
        β=β,
        bond_dim=bond_dim,
        sweeps=sweeps,
    )
    low_energy_spectrum(network, max_states, merge_branches(network, 1.0))
end