text:
__precompile__()
"""
    EasySearch

searching steady states & transition path;
concise design;
requires modules/packages in LOAD_PATH
"""
module EasySearch
    import House  # household life-cycle optimization problems
    import EasyEcon  # economic functions, e.g. CD production function
    import EasyMath: vecExpand, diagr, diagw!  # generalized mathematical and linear-algebra methods

# ==============================================================================
## SECTION 0: Assistant Functions 辅助函数
# NOTE: in this section, we define some assistant functions which are used in searching
#       本节我们定义一批在求解时很有用的辅助函数
#       they are:
#       1. SummaryYear(), prints a summary of the search of steady states
#       2. SummaryTransition(), prints a summary of the search of transition path
#       3. DatSlice4Household(), slices data packages for the life-cycle optimization problem of household department
# --------------
    """
        SummaryYear( t::Int, Dt::Dict, Dst::Dict, Pt::Dict, Ps::Dict, S::Int, Sr::Int )

    Summarizes a specific cross-sectional; usually used in the summary of steady states;
    It returns nothing but prints to console;
    """
    function SummaryYear( t::Int, Dt::Dict, Dst::Dict, Pt::Dict, Ps::Dict, S::Int, Sr::Int )
        println("\t","."^40, " Factor Markets")
        println("\t+ Aggregated Capital: ", Dt[:K][t] )
        println("\t+ Aggregated Labor  : ", Dt[:L][t] )
        println("\t+ Net Intereste Rate: ", Dt[:r][t] )
        println("\t+ Average Wage Level: ", Dt[:w̄][t] )
        println("\t","-"^40, " Good Market")
        println("\t+ GDP               : ", Dt[:Y][t] )
        println("\t+ C, I, G           : ", [ Dt[:C][t], Dt[:I][t], Dt[:G][t] ]  )
        println("\t+ C, I, G Shares    : ", [ Dt[:C][t], Dt[:I][t], Dt[:G][t] ] ./ (Dt[:Y][t] + eps())  )
        println("\t","."^40, " General Fiscal")
        println("\t+ TaxRev Consumption: ", Dt[:TRc][t] )
        println("\t+ TaxRev Wage Income: ", Dt[:TRw][t] )
        println("\t+ Debt & Debt/GDP   : ", [ Dt[:D][t], Dt[:D2Y][t] ] )
        println("\t","."^40, " PAYG Pension")
        println("\t+ Benefit Amount    : ", Dt[:Λ][t] )
        println("\t+ Substitution Rate : ", Dt[:Λ][t] / (Dt[:w̄][t] * Dt[:L][t] / sum(Ps[:N][t, 1:Sr]) + eps())  )
        println("\t","."^40, " UE-BMI")
        println("\t+ Gap of Pooling Acc: ", Dt[:LI][t] )
        println("\t+ Gap/GDP           : ", Dt[:LI][t] / (Dt[:Y][t] + eps()) )
        println("\t+ Gap/TaxRev        : ", Dt[:LI][t] / (Dt[:TRc][t] + Dt[:TRw][t] + eps()) )
        println("\t","."^40, " Household Department")
        println("\t+ Agg Wealth : ", sum(Ps[:N][t, 1:S] .* Dst[:𝒜][t, 1:S]) )
        println("\t+ Agg Personal Asset: ", sum(Ps[:N][t, 1:S] .* Dst[:a][t, 1:S]) )
        println("\t+ Agg UEBMI Indi-Acc: ", sum(Ps[:N][t, 1:S] .* Dst[:Φ][t, 1:S]) )
        println("\t+ Max Wealth & Loc  : ", findmax(Dst[:𝒜][t, 1:S]) )  # prepared for possible DP, which can help check if DP works well
        println("\t","-"^40)

        return nothing
    end
    # ------------------------------
    """
    DatSlice4Household( t::Int, Dt::Dict, Dst::Dict, Pt::Dict, Ps::Dict ; S::Int = 80, Sr::Int = 40, a1::Float64 = 0.0, Φ1::Float64 = 0.0, SearchType::String = "SteadyState" )

    Slices data, generates a Dict **d** for LifeDecision() in EasyHousehold module (just input **Pc** to LifeDecision());
    based on data collections generated by proc_InitPars.jl & proc_VarsDeclare.jl;
    the String parameter **SearchType** controls the type of data slices. It should be one of:
    1. "SteadyState":  The Search of Steady State (S > Sr > 1), i.e. starting from s = 1
    2. "AliveInYear0Working": For those alive in initial steady state but re-optimize working & retired paths in year 1 (S > Sr >= 1)

    The function now can automatically handle the case that there are several years toching possible final steady state;
    returns a Dict **OriginData**.

    用于为稳态、转轨路径的家户最优决策问题切片数据；已经能够自动集中处理各类异常（尤其是touch到最终稳态时的数据拓展问题）

    # Depends on:
    1. EasyMath.vecExpand() : expand vectors for those born in laste years of transition path. some years of these men have reached the final steady state
    2. EasyMath.diagr(): read diagonal elements from a matrix, automatically add/drop extra elements
    """
    function DatSlice4Household( t::Int, Dt::Dict, Dst::Dict, Pt::Dict, Ps::Dict, Pc::Dict, env::NamedTuple ;
        S::Int = 80, Sr::Int = 40, a1::Float64 = 0.0, Φ1::Float64 = 0.0,
        SearchType::String = "SteadyState" )
        # ----------------------
        local OriginData = Dict{Symbol,Union{Float64,Int,Vector{Float64}}}()  # use strict type assertions

        if SearchType == "SteadyState"
        # CASE 1: The Search of Steady State (S > Sr >= 1), i.e. starting from s = 1
        # CASE 1: 搜索稳态（S > Sr >= 1),即从出生开始到死
        # NOTE: :t now indicates which year to search steady state on
            # assert: S > Sr > 1
            @assert( S > Sr > 1, "requiring: S > Sr > 1; CASE SteadyState" )
            # --------- Constants
            OriginData[:Smax] = S; OriginData[:Sret] = Sr
            OriginData[:alpha] = Pc[:α]; OriginData[:gamma] = Pc[:γ]
            OriginData[:k1] = a1 + Φ1
            # --------- Constants in this paper but converted to vectors in a standard problem
            OriginData[:σ] = fill( Pc[:σ] , Sr )
            OriginData[:μ] = fill( Pc[:μ] , S  )
            OriginData[:δ] = fill( Pc[:δ] , S  )
            # --------- Vectors
                # -------- len = S
            OriginData[:Survival] = 1.0 .- Ps[:F][t,1:S] # Survival probs
            OriginData[:r] = fill(Dt[:r][t],S) # interest rates
            OriginData[:cpB] = fill(Pt[:cpB][t],S) # co-payment rate of inpatient expenditure
            OriginData[:p] = Ps[:p][1:S] # MA/MB ratio
                # OriginData[:q] = fill(Pt[:q][t],S) # m2c ratio, the ratio of health expenditure on consumption
            # NOTE: bug fixed, q_t has been upgraded to q_{s,t} which is saved in Ps::Dict rather than Pt::Dict
            OriginData[:q] = Ps[:q][t,1:S]
                # -------- len = Sr
            OriginData[:w] = Dst[:w][t,1:Sr] # wage level
            OriginData[:z] = fill(Pt[:z][t],Sr) # collection rate of pension
            OriginData[:ϕ] = fill(Pt[:ϕ][t],Sr) # contribution: agent → UEBMI
            OriginData[:ζ] = fill(Pt[:ζ][t],Sr) # contribution: firm → UEBMI
            OriginData[:η] = fill(Pt[:η][t],Sr) # contribution: firm → pension
            OriginData[:θ] = fill(Pt[:θ][t],Sr) # contribution: agent → pension
            OriginData[:𝕒] = fill(Pt[:𝕒][t],Sr) # transfer rate from firm contribution of UEBMI to working agents
                # -------- len = S - Sr
            OriginData[:Λ] = fill(Dt[:Λ][t],S-Sr) # pension benefit amounts
            OriginData[:𝕡] = fill(Dt[:𝕡][t],S-Sr) # transfer amounts from firm contribution of UEBMI to retired generations

        elseif SearchType == "AliveInYear0Working"
        # CASE 2: For those alive in initial steady state but re-optimize working & retired paths in year 1 (S > Sr >= 1)
        # CASE 2: 为那些在初始稳态时仍存活且【还未】退休的人，这些人由于第1年的冲击会重新决策剩下年岁里的资产、消费、劳动；此时有 S > Sr >= 1
        # NOTE: now, Sr is the left years of working, and S is the left years to live!
        # NOTE: using the REAL :S in :env to locate the starting years (to slice data from our data packages)
        #       使用env环境变量Tuple里真正的最大年龄env.S来确定此时出发的真正位置（用于在矩阵里取数据）
        # NOTE: the principle to compute offset= : 计算diagr的offset参数的规则
        #       Princp: the agent is s=2 years old at the BEGINNING of year 1 (Steady State)
        #               he/she at most lives for S=4 years (dies at the END of age 4);
        #               and both shock & life-cycle-re-optimization happens at the BEGINNING of year 1;
        #               Therefore, he/she has S`=3=S-s+1 years to live, i.e. the :S parameter in this function;
        #               to slice his/her data, we need an offset=s-1=1, i.e. offset=S-S`
        #       Formula: s = S - S` + 1
        # NOTE: in fact, this formula assumes that we always slice data from the top row of matrix data structures! (i.e. CASE 2,3 ONLY works for the search of first years of transition path!)
        #       事实上，使用这样一个公式意味着CASE 2,3只能用于搜索转轨路径的最开始几年！（因为它意味着总是从矩阵的第一行，最顶端开始取值）
        # NOTE: some lines are specially tendented to indicate that they are specially modified !
        # --------------------
            # assert  S > Sr >= 1
            @assert( S > Sr >= 1, "requiring: S > Sr >= 1; CASE AliveInYear0Working" )
            # --------- Constants
            OriginData[:Smax] = S; OriginData[:Sret] = Sr
            OriginData[:alpha] = Pc[:α]; OriginData[:gamma] = Pc[:γ]
            OriginData[:k1] = a1 + Φ1
            # --------- Constants in this paper but converted to vectors in a standard problem
            OriginData[:σ] = fill( Pc[:σ] , Sr )
            OriginData[:μ] = fill( Pc[:μ] , S  )
            OriginData[:δ] = fill( Pc[:δ] , S  )
            # --------- Vectors
                # -------- len = S
                OriginData[:Survival] = 1.0 .- diagr(Ps[:F], offset = env.S - S , LEN = S)  # please refer to the notes of this CASE
                    # NOTE: bug fixed, Mar 14, 2019; see below
                    # OriginData[:q] = vecExpand(Pt[:q][t:end], S) # m2c ratio
            OriginData[:r] = vecExpand(Dt[:r][t:end], S) # interest rates
            OriginData[:cpB] = vecExpand(Pt[:cpB][t:end], S) # co-payment rate of inpatient expenditure
                OriginData[:p] = Ps[:p][env.S-S+1:env.S] # NOTE: s -> real maximum age
                # NOTE: bug fixed, q_t has been upgraded to q_{s,t} which is saved in Ps::Dict rather than Pt::Dict
                OriginData[:q] = diagr(Ps[:q], offset = env.S - S, LEN = S)
                # -------- len = Sr
                OriginData[:w] = diagr(Dst[:w], offset = env.S - S, LEN = Sr) # wage, (sliced)
            OriginData[:z] = vecExpand(Pt[:z][t:end], Sr) # collection rate of pension
            OriginData[:ϕ] = vecExpand(Pt[:ϕ][t:end], Sr) # contribution: agent → UEBMI
            OriginData[:ζ] = vecExpand(Pt[:ζ][t:end], Sr) # contribution: firm → UEBMI
            OriginData[:η] = vecExpand(Pt[:η][t:end], Sr) # contribution: firm → pension
            OriginData[:θ] = vecExpand(Pt[:θ][t:end], Sr) # contribution: agent → pension
            OriginData[:𝕒] = vecExpand(Pt[:𝕒][t:end], Sr) # transfer rate from firm contribution of UEBMI to working agents
                # -------- len = S - Sr
            OriginData[:Λ] = vecExpand(Dt[:Λ][t:end], S-Sr) # pension benefit amounts
            OriginData[:𝕡] = vecExpand(Dt[:𝕡][t:end], S-Sr) # transfer amounts from firm contribution of UEBMI to retired generations

        elseif SearchType == "AliveInYear0Retired"
        # CASE 3: For those alive in initial steady state but re-optimize ONLY retired paths in year 1 (S >= 1)
        #         为那些在初始稳态时仍存活且【已经】退休的人，这些人由于第1年的冲击会重新决策剩下年岁里的资产、消费、劳动；此时有 S >= 1
        # NOTE: now, Sr is no more required, and S is the left years to live!
        # NOTE: using the REAL :S in :env to locate the starting years (to slice data from our data packages)
        #       使用env环境变量Tuple里真正的最大年龄env.S来确定此时出发的真正位置（用于在矩阵里取数据）
        # NOTE: the principle to compute offset= is the same as that of CASE 2 计算diagr的offset参数的规则和上一种CASE相同
        # -----------
            # assert  S >= 1
            @assert( S >= 1, "requiring: S > Sr >= 1; CASE AliveInYear0Retired" )
            # --------- Constants
            OriginData[:Smax] = S
            OriginData[:alpha] = Pc[:α]; OriginData[:gamma] = Pc[:γ]
            OriginData[:k1] = a1 + Φ1
            # --------- Constants in this paper but converted to vectors in a standard problem
            OriginData[:μ] = fill( Pc[:μ] , S  )
            OriginData[:δ] = fill( Pc[:δ] , S  )
            # --------- Vectors
                # -------- len = S
                OriginData[:Survival] = 1.0 .- diagr(Ps[:F], offset = env.S - S , LEN = S)  # please refer to the notes of this CASE
                    # NOTE: bug fixed, Mar 14, 2019; see below
                    # OriginData[:q] = vecExpand(Pt[:q][t:end], S)
                # NOTE: bug fixed, q_t has been upgraded to q_{s,t} which is saved in Ps::Dict rather than Pt::Dict
                OriginData[:q] = diagr(Ps[:q], offset = env.S - S, LEN = S)
            OriginData[:r] = vecExpand(Dt[:r][t:end], S)
            OriginData[:cpB] = vecExpand(Pt[:cpB][t:end], S)
                OriginData[:p] = Ps[:p][env.S-S+1:env.S] # NOTE: s -> real maximum age
                # -------- len = S (NOTE: :S is the left years to live!, and the indices of Λ & 𝕡 have been changed)
                OriginData[:Λ] = vecExpand(Dt[:Λ][t:end], S )
                OriginData[:𝕡] = vecExpand(Dt[:𝕡][t:end], S )

        elseif SearchType == "BornInTransition"
        # CASE 4: For those born in transition path (including those who have reached final steady state before natural death)
        #         为那些在转轨路径上出生的人（包括有某些年已经到了最终稳态的那些人）
        # NOTE: now, we always compute FULL life-cycle paths, i.e. S > Sr >= 1 must be satisfied
        #       现在我们总是计算完整的生命期消费劳动资产等路径，必有S>Sr>=1成立
        #       and, because we use special-designed functions EasyMath.vecExpand() & EasyMath.diagr()
        #       并且，由于我们使用了特别设计的两个函数
        #       the function can automatically construct full-length vectors for those who have reached final steady state before natural death
        #       我们可以自动构建全长度的数据向量，为那些在死之前就已经到达最终稳态的人(i.e. t_born + S > T_max)
        #       where actually in our data structures, the years after T_max (the moment reaching final steady state) are not saved
        #       这样设计是因为在我们的数据结构里，当到达最终稳态即超过最大年数T（转轨路径长度）之后的数据并不会特别储存（因为它们和最终稳态一样）
        # NOTE: different from CASE 2,3, we use a different principle to determine the offset= in diagr()
        #       不同于CASE 2,3，我们现在使用不同的规则来决定diagr()中offset=参数的值：
        #       Principle: offset = 1 - t, where t is the year when agents born, t = 1,...,T
        # -----------
            # assert  S > Sr >= 1
            @assert( S > Sr >= 1, "requiring: S > Sr >= 1; CASE AliveInYear0Working" )
            # --------- Constants
            OriginData[:Smax] = S; OriginData[:Sret] = Sr
            OriginData[:alpha] = Pc[:α]; OriginData[:gamma] = Pc[:γ]
            OriginData[:k1] = a1 + Φ1
            # --------- Constants in this paper but converted to vectors in a standard problem
            OriginData[:σ] = fill( Pc[:σ] , Sr )
            OriginData[:μ] = fill( Pc[:μ] , S  )
            OriginData[:δ] = fill( Pc[:δ] , S  )
            # --------- Vectors
                # -------- len = S
                OriginData[:Survival] = 1.0 .- diagr(Ps[:F], offset = 1 - t, LEN = S)
                    # NOTE: bug fixed, Mar 14, 2019; see below
                    # OriginData[:q] = vecExpand(Pt[:q][t:end], S) # m2c ratio
                OriginData[:q] = diagr(Ps[:q], offset = 1 - t, LEN = S)
            OriginData[:r] = vecExpand(Dt[:r][t:end], S) # interest rates
            OriginData[:cpB] = vecExpand(Pt[:cpB][t:end], S) # co-payment rate of inpatient expenditure
                OriginData[:p] = Ps[:p][env.S-S+1:env.S] # NOTE: s -> real maximum age
                # -------- len = Sr
                OriginData[:w] = diagr(Dst[:w], offset = 1 - t, LEN = Sr) # wage, (sliced)
            OriginData[:z] = vecExpand(Pt[:z][t:end], Sr) # collection rate of pension
            OriginData[:ϕ] = vecExpand(Pt[:ϕ][t:end], Sr) # contribution: agent → UEBMI
            OriginData[:ζ] = vecExpand(Pt[:ζ][t:end], Sr) # contribution: firm → UEBMI
            OriginData[:η] = vecExpand(Pt[:η][t:end], Sr) # contribution: firm → pension
            OriginData[:θ] = vecExpand(Pt[:θ][t:end], Sr) # contribution: agent → pension
            OriginData[:𝕒] = vecExpand(Pt[:𝕒][t:end], Sr) # transfer rate from firm contribution of UEBMI to working agents
                # -------- len = S - Sr
            OriginData[:Λ] = vecExpand(Dt[:Λ][t:end], S-Sr) # pension benefit amounts
            OriginData[:𝕡] = vecExpand(Dt[:𝕡][t:end], S-Sr) # transfer amounts from firm contribution of UEBMI to retired generations

        else
        # CASE : wrong SearchType= input
            throw(ErrorException("in valid keyword parameter SearchType=, please use one of : SteadyState, BornInTransition, AliveInYear0Retired, AliveInYear0Working"))
        end

        # return
        return OriginData::Dict
    end
    # ----------------

# ==============================================================================
## SECTION 1: Steady State 稳态（长期均衡）搜索
# NOTE: in this section, we define an in-place function SteadyState!(t)
#       本节我们定义一个直接作用于数据包上的函数SteadyState!(t)
#       this function searches steady states (SS) based on input datasets, where t is the year to search a steady state
#       这个函数基于输入的数据集，在特定年上搜索稳态(SS)
#       the results are made through modifying :Dt & :Dst
#       搜索结果直接反映在Dt和Dst这两个数据包里，直接修改而不是返回新的数据包
#       this kind of design allows us to simutaneously solve multiple different models which are defined by data
#       这样的设计允许我们同时求解一揽子参数/数据不同的模型
# --------------
    """
        SteadyState!( t::Int, Guess::NamedTuple, Dt::Dict, Dst::Dict, Pt::Dict, Ps::Dict, Pc::Dict, env::NamedTuple ; atol::Float64 = 1E-8, MaxIter::Int = 100, PrintMode::String = "full", MagicNum::Real = 2.0, StepLen::Real = 0.5 )

    Searches steady state on a specific **t** slice, using guesses of interest rate & labor (**Guess**);
    the results will be directly/in-place written in data collections (**Dt**, **Dst**);
    returns nothing.

    ## Inputs:
    1. t : the year to search steady state on, it should be less than or equal to *env.MAX_YEAR*
    2. Guess : a NamedTuple consisting of two elements: ( r = interest rate in digits, L = labor factor )
    3. Dt : a data collection to modify/write, containing year data; defined in proc_VarsDeclare
    4. Dst : a data collection to modify/write, containing age × year data; defined in proc_VarsDeclare
    5. env : a parameter collection containing basic parameters like MAX_AGE, RETIRE_AGE and others
    6. Pt : a parameter collection containing year data, won't be modified
    7. Ps : a parameter collection containing age specific data & demographic data (population & mortality); won't be modified
    8. Pc : a parameter collection containing constant data, won't be modified
    9. atol : tolerance of Gauss-Seidel iteration
    10. MaxIter : maximum loops
    11. PrintMode : a string to decide how to display results & loop information
        1. "full" : print all information, including round, error, summary of economy in each round; print a detailed summary when exit
        2. "concise" : print round, error, but no summary of economy for each round; print a detailed summary when ends
        3. "final": print nothing for each round; print a detailed summary when iteration ends
        4. "silent" : nothing printed for each round; nothing printed when iteration ends
    12. MagicNum : a lower bound of capital per labor (K/L) to prevent divergence out of interest rate (a function of capital per labor)
    13. StepLen : relative step length to update guesses of Gauss-Seidel algorithm, in range (0,1]

    ## Modified:
    1. Dt: update economic variables (e.g. GDP, interest rate)
    2. Dst: update consumption, labor and other age-specific data

    ## Iteration Steps:
    1. Firm department
    2. Pension
    3. Household problems
    4. Fiscal, Aggregation & Update

    ## Depends on:
    1.

    """
    function SteadyState!( t::Int, Guess::NamedTuple, Dt::Dict, Dst::Dict,
        Pt::Dict, Ps::Dict, Pc::Dict, env::NamedTuple ;
        atol::Float64 = 1E-8,  # tolerance of Gauss-Seidel iteration
        MaxIter::Int = 100,  # maximum loops
        PrintMode::String = "full",  # mode of printing, one of ["full","concise","final","silent"]
        MagicNum::Real = 2.0,  # magic number, the lower bound of K/L (capital per labor)
        StepLen::Real = 0.5  # relative step length to update guesses, in range (0,1]
    )
        ## Section: validation
        @assert( 0 < t <= env.T , string("t not in range [1, ",env.T,"], received :",t) )
        @assert( isa(Guess.r, Real) & (Guess.L > 0.0), "please check your Guess::NamedTuple" )  # simutaneously check existence of r, L in Guess and their range
        @assert( PrintMode in ["full","concise","final","silent"] , "undefined PrintMode string" )
        @assert( 0.0 < StepLen <= 1.0 , "invalid StepLen, it should be in range (0,1]" )
        @assert( MagicNum > 0.0 , "invalid MagicNum (capital per capita) which must be greater than 0" )
        Ps[:F][t,env.S] = 0.0  # forcely refresh mortality in the very last year 强制刷新最后一岁的死亡率为0
        ## Section: malloc 内存预分配
        local tmpK1::Float64, tmpK2::Float64  # temporary variables of capital factor in iterations
        local tmpL1::Float64, tmpL2::Float64  # temporary variables of labor factor
            tmpL1 = Guess.L  # save the guess of aggregated labor (initialization)
            Dt[:r][t] = Guess.r  # save the guess of interest rate (initialiation)
        ## Section: intermediate variables & abbreviations 快捷变量
        local tmpVal::Float64 = 1.0 + Pt[:z][t] * Pt[:η][t] + Pt[:ζ][t] # a temporary variable
        local πCoef::Float64 = Pt[:z][t] * (Pt[:θ][t] + Pt[:η][t]) / tmpVal # total contribution to pension (on nomial wage level)
        local πMf::Float64 = Pt[:ζ][t] / tmpVal # firm contribution rate to UE-BMI
        local πMp::Float64 = Pt[:ϕ][t] / tmpVal # employee/personal contribution rate to UE-BMI
        local πM = πMf + πMp # total contribution rate to UE-BMI on nomial wage level
        ## Section: fill labor with an even distribution to initialize labor 使用平坦的劳动力分布填充以初始化
        Dst[:Lab][t,1:env.Sr] .= tmpL1 / env.Sr

        ## Section: Gauss-Seidel Iterations
        for idx in 1:MaxIter
            # SubSection 1: firm department 厂商部门
                # 1. production function & firm optimal decisions 生产函数 & 厂商最优决策
                if idx == 1   # NOTE: use interest rate & labor to get capital in the first round  第一轮使用利率和劳动力反推出资本存量
                    tmpK1 = tmpL1 * ( (Guess.r + Pc[:κ]) / (Pt[:β][t] * Pt[:A][t]) ) ^ (1.0 / (Pt[:β][t] - 1.0))
                    @assert( tmpK1 > 0.0 , string("negative K found in round: ", idx) ) # check capital
                    Dt[:Y][t], tmpVal, Dt[:w̄][t] = EasyEcon.CDProdFunc( Pt[:A][t], tmpK1, tmpL1, Pt[:β][t], κ = Pc[:κ], GetPrice = true, TechType = "Hicks"  ) # optimal decision 最优决策，返回价格
                else
                    Dt[:Y][t], Dt[:r][t], Dt[:w̄][t] = EasyEcon.CDProdFunc( Pt[:A][t], tmpK1, tmpL1, Pt[:β][t], κ = Pc[:κ], GetPrice = true, TechType = "Hicks"  ) # optimal decision 最优决策，返回价格
                end
                # 2. wage profiling 工资曲线
                Dst[:w][t,:], Dt[:o][t] = EasyEcon.WageProfile( Dt[:w̄][t], Ps[:ε][1:env.Sr], Ps[:N][t,1:env.Sr], Dst[:Lab][t,1:env.Sr], GetScalingCoef = true  )

            # SubSection 2: PAYG pension benefits 养老金给付
                Dt[:Λ][t] = EasyEcon.PAYGPension( πCoef, Dst[:w][t,:], Ps[:N][t,:], Dst[:Lab][t,:], env.Sr )
            # SubSection 3: average transfer amount from firm medical contribution to retired generations 当期企业医保缴纳对退休人员的平均补贴量
                Dt[:𝕡][t] = EasyEcon.Get𝕡( Pt[:𝕓][t], πMf, Dst[:w][t,:], Dst[:Lab][t,1:env.Sr], Ps[:N][t,1:env.S], env.Sr )

            # SubSection 4: household lifetime optimization problem 家庭部门生命期效用最大化问题
                # 1. prepare sliced data 准备输入数据切片
                local tmpOriginData = DatSlice4Household(t, Dt,Dst,Pt,Ps,Pc,env,
                    S = env.S, Sr = env.Sr, a1 = 0.0, Φ1 = 0.0, SearchType = "SteadyState"  )
                # 2. solve the optimization, get paths (wealth/capital, asset, UEBMI-indi, consumption, leisure) 求解，得到路径（财富、资产、个人医保、消费、闲暇）
                local tmpRet = House.HHSolve( tmpOriginData, ReturnData = true )
                # 2.5 extract asset (a_{s}) & UEBMI-indi (Φ_{s}) 提取个人资产和UEBMI个人账户
                House.ExtractAPhi!( tmpRet, tmpOriginData, a1 = 0.0 )
                # 3. distribute results 存储结果
                # NOTE: k,a,Φ (len=S+1) have an extra element: bequest, so, drop it
                Dst[:𝒜][t, 1:env.S] = tmpRet[:ks][1:env.S]
                Dst[:a][t, 1:env.S] = tmpRet[:as][1:env.S]
                Dst[:Φ][t, 1:env.S] = tmpRet[:Φs][1:env.S]
                # NOTE: c,M,MA,MB (len=S)
                Dst[:c][t, 1:env.S] = tmpRet[:cs][1:env.S]
                Dst[:Lab][t, 1:env.Sr] = tmpRet[:ls][1:env.Sr]  # labor supply (not leisure!)
                Dst[:m][t, 1:env.S] .= tmpRet[:Ms][1:env.S]   # medical expenditure
                Dst[:MA][t, 1:env.S] .= tmpRet[:MAs][1:env.S] # outpatient expenditure
                Dst[:MB][t, 1:env.S] .= tmpRet[:MBs][1:env.S] # inpatient expenditure
                Dst[:ΦGaps][t, 1:env.S] .= tmpRet[:ΦGaps][1:env.S] # inpatient expenditure

            # SubSection 5: fiscal & UE-BMI 财政 &　UE-BMI
                # 1. tax revenues 税收
                Dt[:TRc][t] = Pc[:μ] * sum( Dst[:c][t, 1:env.S]   .* Ps[:N][t,1:env.S] )  # consumption tax 消费税
                Dt[:TRw][t] = Pc[:σ] * sum( Dst[:Lab][t, 1:env.Sr] .* Ps[:N][t,1:env.Sr] .* Dst[:w][t, 1:env.Sr] )  # wage tax 工资税
                # 2. the gaps of the social pooling account of UE-BMI (positive for gap, negative for surplus)
                Dt[:LI][t] = sum( ( 1 .- Pt[:cpB][t] ) .* Dst[:MB][t, 1:env.S] .* Ps[:N][t,1:env.S] )
                Dt[:LI][t] -= ( 1 .- Pt[:𝕒][t] .- Pt[:𝕓][t] ) .* Pt[:ζ][t] ./ ( 1 .+ Pt[:η][t] + Pt[:ζ][t] ) .* sum( Ps[:N][t,1:env.Sr] .* Dst[:w][t, 1:env.Sr] .* Dst[:Lab][t, 1:env.Sr] )

            # SubSection 6: update aggregated labor supply 更新劳动力供应
                tmpL2 = sum( Ps[:N][t,1:env.Sr] .* Dst[:Lab][t, 1:env.Sr] )
            # SubSection 7: update GDP with the updated labor supply (through production function) 使用更新的劳动要素更新GDP估计
                Dt[:Y][t] = EasyEcon.CDProdFunc( Pt[:A][t], tmpK1, tmpL2, Pt[:β][t], κ = Pc[:κ], GetPrice = true, TechType = "Hicks"  )[1]

            # SubSection 8: update updating aggregated capital (K)
                # NOTE: there are two ways to update K: one uses the capital market clearing condition (Method A);
                # and the other one uses the capital growth dynamics (Method B); the two methods are equivalent; and I use Method A here.
                # 1. get aggregated investment (I) through capita dynamics 通过资本动态条件得到投资
                Dt[:I][t] = tmpK1 * Pc[:κ]
                # 2. aggregate consumption 汇总社会总消费
                Dt[:C][t] = sum( Dst[:c][t, 1:env.S] .* Ps[:N][t,1:env.S] )
                # 3. get government purchase (G) through the good market clearing condition 通过商品市场出清得到政府购买
                Dt[:G][t] = Dt[:Y][t] - Dt[:I][t] - Dt[:C][t]
                # 4. get government outstanding debt through the fiscal budget 通过政府预算约束得到政府债务余额
                Dt[:D][t] = ( Dt[:G][t] + Dt[:LI][t] - Dt[:TRc][t] - Dt[:TRw][t] ) / (1 - Dt[:r][t])
                # 5. the upper bound of government outstanding debt 政府未偿债务软约束
                Dt[:D][t] = max( 0.0, min( Dt[:D][t], Dt[:Y][t] * Pt[:D2Ycap][t] ) )
                # 6. record the ratio of government outstanding debt on GDP 记录未偿债务与GDP的比
                Dt[:D2Y][t] = Dt[:D][t] / Dt[:Y][t]
                # 7. aggregate/update capital 汇总/更新（年初）资本存量
                tmpK2 = sum( Dst[:𝒜][t, 1:env.S] .* Ps[:N][t,1:env.S] )

            # Convergence check & go to the next loop 收敛检查
                # 1. compute errors, using a minor number to avoid exact zeros
                local Err = ( K = abs(tmpK2 / (tmpK1 + eps()) - 1 ), L = abs(tmpL2 / (tmpL1 + eps()) - 1 ) )
                # 2. record original K, L
                local OriGuess = [tmpK1, tmpL1]
                # 3. check
                if (Err.K < atol) & (Err.L < atol)  # converged
                    # 1. save the converged K, L; using mean values
                    Dt[:K][t] = ( tmpK1 + tmpK2 ) / 2; Dt[:L][t] = ( tmpL1 + tmpL2 ) / 2;
                    # 2. print final summary (if not silently solved)
                    if PrintMode != "silent"
                        println("\n\t+ Status: Converged")  # status
                        println("\t+ Round: ", idx)
                        println("\t+ Relative Errors: ",Err)  # print error first
                        println("\t+ Updating of Kap: ",round.([ OriGuess[1], tmpK2, Dt[:K][t] ],digits=6) ) # then print the updating process of K,L
                        println("\t+ Updating of Lab: ",round.([ OriGuess[2], tmpL2, Dt[:L][t] ],digits=6) )
                        println("\t+ NOTE: guess -> re-aggregated -> updated")
                        SummaryYear( t, Dt, Dst, Pt, Ps, env.S, env.Sr )  # print the details of the steady state
                    end
                    break  # ends iteration
                elseif idx == MaxIter  # diverged or too few maxmimum rounds set
                    Dt[:K][t] = ( tmpK1 + tmpK2 ) / 2; Dt[:L][t] = ( tmpL1 + tmpL2 ) / 2;
                    if PrintMode != "silent"
                        println("\n\t+ Status: Maximum iteration reached, not converged")
                        println("\t+ Round: ", idx)
                        println("\t+ Relative Errors: ",Err)
                        println("\t+ Updating of Kap: ",round.([ OriGuess[1], tmpK2, Dt[:K][t] ],digits=6) ) # then print the updating process of K,L
                        println("\t+ Updating of Lab: ",round.([ OriGuess[2], tmpL2, Dt[:L][t] ],digits=6) )
                        println("\t+ NOTE: guess -> re-aggregated -> updated")
                        SummaryYear( t, Dt, Dst, Pt, Ps, env.S, env.Sr )
                    end
                    break  # ends iteration
                else  # go to the next loop
                    # 1. check if the updated labor supply touches the bottom of zero
                    (tmpL2 <     0 || ~isreal(tmpL2))    &&    begin @warn("Labor lower than 0 or complex!"); tmpL2 = 0.01; end
                    # 2. update labor supply
                    tmpL1 += StepLen * (tmpL2 - tmpL1)
                    # 3. use the magic number to set the bottom of capital, which bounds the interest rate in a reasonable range
                    local tmpKfloor = MagicNum * tmpL1
                    # 4. check & update capital supply
                    (tmpK2 < tmpKfloor || ~isreal(tmpK2))    &&    begin @warn("Capital lower than the lower bound or complex!"); tmpK2 = tmpKfloor; end
                    tmpK1 += StepLen * (tmpK2 - tmpK1)
                    # 5. print information, according to :PrintMode
                    if PrintMode in ["full", "concise"]
                        println("\n\t+ Round: ", idx)
                        println("\t+ Relative Errors: ", Err)
                        println("\t+ Updating of Kap: ",round.([ OriGuess[1], tmpK2, Dt[:K][t] ],digits=6) ) # then print the updating process of K,L
                        println("\t+ Updating of Lab: ",round.([ OriGuess[2], tmpL2, Dt[:L][t] ],digits=6) )
                        println("\t+ NOTE: guess -> re-aggregated -> updated")
                    end
                    if PrintMode == "full"
                        SummaryYear( t, Dt, Dst, Pt, Ps, env.S, env.Sr )
                    end
                    # 6. explicitly write continue (as a reminder without practical jobs)
                    continue
                end  # check ends

        end # Gauss-Seidel Iteration ends
        # nomial return
        return nothing
    end  # function ends

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

# ==============================================================================
## SECTION 2: Transition Search 转轨路径搜索
# NOTE: in this section, we define an in-place function Transition!()
#       本节我们定义一个直接作用于数据包上的函数Transition!()
#       it is used to search a transition path on input datasets
#       using Gauss-Seidel iterations
#       returns nothing
# ------------
    """
        Transition!( Dt::Dict, Dst::Dict, Pt::Dict, Ps::Dict, Pc::Dict, env::NamedTuple ; atol::Float64 = 1E-8, MaxIter::Int = 100, PrintMode::String = "full", MagicNum::Real = 2.0, StepLen::Real = 0.5 )

    Searches a transition path based on the input datasets;
    requires Initial & Final Steady States searched;
    returns nothing;

    And, the definitions of input parameters are the same as those of SteadyState!()
    """
    function Transition!( Dt::Dict, Dst::Dict, Pt::Dict, Ps::Dict, Pc::Dict, env::NamedTuple ;
        atol::Float64 = 1E-8,  # tolerance of Gauss-Seidel iteration
        MaxIter::Int = 100,  # maximum loops
        PrintMode::String = "full",  # mode of printing, one of ["full","concise","final","silent"]
        MagicNum::Real = 2.0,  # magic number, the lower bound of K/L (capital per labor)
        StepLen::Real = 0.5, # relative step length to update guesses, in range (0,1]
        ReturnLog::Bool = false ) # whether to return a Dict of convergence (errors per loop)
        # ------------
        ## Section: validation
            @assert( PrintMode in ["full","concise","final","silent"] , "undefined PrintMode string" )
            @assert( 0.0 < StepLen <= 1.0 , "invalid StepLen, it should be in range (0,1]" )
            @assert( MagicNum > 0.0 , "invalid MagicNum (capital per capita) which must be greater than 0" )
            Ps[:F][:,env.S] .= 0.0  # forcely refresh mortality in the very last year 强制刷新最后一岁的死亡率为0
        ## Section: malloc 内存预分配
            # a. linearly fill capital & labor 线性填充资本和劳动供应
            local tmpK1::Vector{Float64} = LinRange( Dt[:K][1], Dt[:K][env.T], env.T )
            local tmpK2::Vector{Float64} = LinRange( Dt[:K][1], Dt[:K][env.T], env.T )
            local tmpL1::Vector{Float64} = LinRange( Dt[:L][1], Dt[:L][env.T], env.T )
            local tmpL2::Vector{Float64} = LinRange( Dt[:L][1], Dt[:L][env.T], env.T )
            # b. initialization, evenly fill every year's labor among generations 平均填充每年每岁人的劳动力供应
            for t in 2:env.T-1
                Dst[:Lab][t,1:env.Sr] .= tmpL1[t] ./ env.Sr
            end
            # c. other convenient variables 其他快捷变量
            local tmpVal = 1.0    .+  Pt[:z] .* Pt[:η]  .+ Pt[:ζ]
            local πCoef  = Pt[:z] .* (Pt[:θ] .+ Pt[:η]) ./ tmpVal
            local πMf    = Pt[:ζ] ./  tmpVal
            local πMp    = Pt[:ϕ] ./  tmpVal
            local πM     = πMf    .+  πMp
            # d. convenient index 快捷索引
            local idxS  = 1:env.S
            local idxWorking = 1:env.Sr
            local idxRetired = env.Sr+1:env.S
            local idx2toT = 2:env.T
            local idx1toT = 1:env.T
            local idx2toTminus1 = 2:env.T-1
            local idx1toTminus1 = 1:env.T-1

            # e. save a copy of k_{s}, a_{s}, Φ_{s} for both init & final steady states
            # NOTE: because we may over-write them, if we want better readiability of the code
            local copy_InitSSk = ( k = Dst[:𝒜][1,:],     a = Dst[:a][1,:],     Φ = Dst[:Φ][1,:] )
            local copy_FinaSSk = ( k = Dst[:𝒜][env.T,:], a = Dst[:a][env.T,:], Φ = Dst[:Φ][env.T,:] )
            # f. logs of performance/convergence 迭代日志
            local PerfLog = Dict( :K => Array{Float64,1}(), :L => Array{Float64,1}() )

        ## Section: Gauss-Seidel Iterations 搜索
        for idxiter in 1:MaxIter
            # check: NaN of Capital & Labor 资本和劳动的NaN值检查（必要，防止可能的数值错误后仍继续狂飙）
            @assert( !any(isnan.(tmpK1)) , "NaN found in agg capital" )
            @assert( !any(isnan.(tmpL1)) , "NaN found in agg labor" )

            # Section: Firm Department & PAYG Pension 厂商部门与养老金
            for t in idx1toTminus1
                # 1. GDP, r, avg w GDP，利率，平均工资
                Dt[:Y][t], Dt[:r][t], Dt[:w̄][t] =
                    EasyEcon.CDProdFunc( Pt[:A][t], tmpK1[t], tmpL1[t], Pt[:β][t],
                    κ = Pc[:κ], GetPrice = true, TechType = "Hicks"  )
                # 2. wage profile 工资曲线
                Dst[:w][t,idxWorking], Dt[:o][t] =
                    EasyEcon.WageProfile( Dt[:w̄][t], Ps[:ε][idxWorking],
                    Ps[:N][t,idxWorking], Dst[:Lab][t,idxWorking], GetScalingCoef = true  )
                # 3. the benefits of PAYG pension 养老金给付
                Dt[:Λ][t] = EasyEcon.PAYGPension( πCoef[t], # NOTE: diff
                    Dst[:w][t,idxWorking], Ps[:N][t,idxS], Dst[:Lab][t,idxWorking], env.Sr )
                # 4. average transfer amount from firm medical contribution to retired generations 当期企业医保缴纳对退休人员的平均补贴量
                Dt[:𝕡][t] = EasyEcon.Get𝕡( Pt[:𝕓][t], πMf[t],  # NOTE: diff
                    Dst[:w][t,idxWorking], Dst[:Lab][t,idxWorking], Ps[:N][t,idxS], env.Sr )
            end # firm & pension ends

            # Section: Household Department & Health Expenditure 家庭部门与医疗支出
            # NOTE: there are three cases to compute in this section:
            #       本节计算三种情形：
            #       1. for those alive in t=1, still working, but change their paths in t>2,...
            #       2. for those alive in t=1, have retired, but change their paths in t>2,...
            #       3. for those born in t=1,...,T (i.e. on the transition path)
            for s in 2:env.Sr
            # NOTE: CASE 1: alive in t=1, working
                # 1. prepare data
                # NOTE: t = 1 (mark year); S = S-s+1 (left years to live); Sr = Sr-s+1 (left years to work)
                #       remember, when s=Sr, we are standing at the beginning of the year of age Sr
                local YearsToLive::Int = env.S  - s + 1
                local YearsToWork::Int = env.Sr - s + 1
                local tmpOriginData = DatSlice4Household(1, Dt,Dst,Pt,Ps,Pc,env, S = YearsToLive, Sr = YearsToWork,
                    a1 = Dst[:a][1,s], Φ1 = Dst[:Φ][1,s], SearchType = "AliveInYear0Working"  )
                # 2. optimization
                local tmpRet = House.HHSolve( tmpOriginData, ReturnData = true )
                # 3. extract asset (a_{s}) & UEBMI-indi (Φ_{s})
                House.ExtractAPhi!( tmpRet, tmpOriginData, a1 = Dst[:a][1,s] )
                # 4. distribute results
                # NOTE: k,a,Φ (len=S+1) have an extra element: bequest, so, drop it
                diagw!( Dst[:𝒜], tmpRet[:ks][1:YearsToLive], offset = s - 1 )
                diagw!( Dst[:a], tmpRet[:as][1:YearsToLive], offset = s - 1 )
                diagw!( Dst[:Φ], tmpRet[:Φs][1:YearsToLive], offset = s - 1 )
                # NOTE: c,M,MA,MB (len=S)
                diagw!( Dst[:c], tmpRet[:cs][1:YearsToLive], offset = s - 1 )
                diagw!( Dst[:Lab], tmpRet[:ls][1:YearsToWork], offset = s - 1 )
                diagw!( Dst[:m], tmpRet[:Ms][1:YearsToLive], offset = s - 1 )
                diagw!( Dst[:MA], tmpRet[:MAs][1:YearsToLive], offset = s - 1 )
                diagw!( Dst[:MB], tmpRet[:MBs][1:YearsToLive], offset = s - 1 )
                diagw!( Dst[:ΦGaps], tmpRet[:ΦGaps][1:YearsToLive], offset = s - 1 )
            end
            for s in env.Sr+1:env.S
            # NOTE: CASE 2: alive in t=1, retired
                # 1. prepare data
                # NOTE: t = 1 (mark year); S = S-s+1 (left years to live)
                local YearsToLive::Int = env.S  - s + 1
                local tmpOriginData = DatSlice4Household(1, Dt,Dst,Pt,Ps,Pc,env, S = YearsToLive, Sr = 0,
                    a1 = Dst[:a][1,s], Φ1 = Dst[:Φ][1,s], SearchType = "AliveInYear0Retired"  )
                # 2. optimization
                local tmpRet = House.HHSolve_Retired( tmpOriginData, ReturnData = true )
                # 3. extract asset (a_{s}) & UEBMI-indi (Φ_{s})
                House.ExtractAPhi_Retired!( tmpRet, tmpOriginData, a1 = Dst[:a][1,s] )
                # 4. distribute results
                # NOTE: k,a,Φ (len=S+1) have an extra element: bequest, so, drop it
                diagw!( Dst[:𝒜], tmpRet[:ks][1:YearsToLive], offset = s - 1 )
                diagw!( Dst[:a], tmpRet[:as][1:YearsToLive], offset = s - 1 )
                diagw!( Dst[:Φ], tmpRet[:Φs][1:YearsToLive], offset = s - 1 )
                # NOTE: c,M,MA,MB (len=S)
                diagw!( Dst[:c], tmpRet[:cs][1:YearsToLive], offset = s - 1 )
                diagw!( Dst[:m], tmpRet[:Ms][1:YearsToLive], offset = s - 1 )
                diagw!( Dst[:MA], tmpRet[:MAs][1:YearsToLive], offset = s - 1 )
                diagw!( Dst[:MB], tmpRet[:MBs][1:YearsToLive], offset = s - 1 )
                diagw!( Dst[:ΦGaps], tmpRet[:ΦGaps][1:YearsToLive], offset = s - 1 )
            end
            # NOTE: refresh the capital/asset/UEBMI-indi in our initial steady state
                Dst[:𝒜][1,idxS] = copy_InitSSk.k
                Dst[:a][1,idxS] = copy_InitSSk.a
                Dst[:Φ][1,idxS] = copy_InitSSk.Φ
            for t in 1:env.T-1
            # NOTE: CASE 3: born on transition path
                # 1. prepare data
                # NOTE: S = S, Sr = Sr, t = t
                #       everyone has a complete life-cycle
                local tmpOriginData = DatSlice4Household(t, Dt,Dst,Pt,Ps,Pc,env, S = env.S, Sr = env.Sr,
                    a1 = 0.0, Φ1 = 0.0, SearchType = "BornInTransition"  )
                # 2. optimization
                local tmpRet = House.HHSolve( tmpOriginData, ReturnData = true )
                # 3. extract asset (a_{s}) & UEBMI-indi (Φ_{s})
                House.ExtractAPhi!( tmpRet, tmpOriginData, a1 = 0.0 )
                # 4. distribute results
                # NOTE: k,a,Φ (len=S+1) have an extra element: bequest, so, drop it
                diagw!( Dst[:𝒜], tmpRet[:ks][idxS], offset = 1 - t )
                diagw!( Dst[:a], tmpRet[:as][idxS], offset = 1 - t )
                diagw!( Dst[:Φ], tmpRet[:Φs][idxS], offset = 1 - t )
                # NOTE: c,M,MA,MB (len=S)
                diagw!( Dst[:c], tmpRet[:cs][idxS], offset = 1 - t )
                diagw!( Dst[:Lab], tmpRet[:ls][idxWorking], offset = 1 - t )
                diagw!( Dst[:m], tmpRet[:Ms][idxS], offset = 1 - t )
                diagw!( Dst[:MA], tmpRet[:MAs][idxS], offset = 1 - t )
                diagw!( Dst[:MB], tmpRet[:MBs][idxS], offset = 1 - t )
                diagw!( Dst[:ΦGaps], tmpRet[:ΦGaps][idxS], offset = 1 - t )
            end
            # NOTE: refresh the capital/asset/UEBMI-indi in our final steady state
                Dst[:𝒜][env.T,idxS] = copy_FinaSSk.k
                Dst[:a][env.T,idxS] = copy_FinaSSk.a
                Dst[:Φ][env.T,idxS] = copy_FinaSSk.Φ

            # Section: Fiscal, Consumption & Labor Aggregation 财政，消费与劳动力市场
            for t in idx1toTminus1
                # 1. aggregated consumption 总消费
                Dt[:C][t] = sum( Dst[:c][t,idxS] .* Ps[:N][t,idxS] )
                # 2. tax revenues 财政收入
                Dt[:TRc][t] = Pc[:μ] * Dt[:C][t]
                Dt[:TRw][t] = Pc[:σ] * sum( Dst[:Lab][t, idxWorking] .* Ps[:N][t,idxWorking] .* Dst[:w][t, idxWorking] )
                # 3. the gaps of the social pooling account of UE-BMI (positive for gap, negative for surplus)
                Dt[:LI][t] = sum( ( 1 .- Pt[:cpB][t] ) .* Dst[:MB][t, idxS] .* Ps[:N][t,idxS] )
                Dt[:LI][t] -= ( 1 .- Pt[:𝕒][t] .- Pt[:𝕓][t] ) .* Pt[:ζ][t] ./
                    ( 1 .+ Pt[:η][t] + Pt[:ζ][t] ) .*
                    sum( Ps[:N][t,idxWorking] .* Dst[:w][t,idxWorking] .* Dst[:Lab][t,idxWorking] )
                # 4. update aggregated labor supply 更新劳动力供应
                tmpL2[t] = sum( Ps[:N][t,idxWorking] .* Dst[:Lab][t,idxWorking] )
                # 5. update GDP using new L 更新GDP
                Dt[:Y][t] = EasyEcon.CDProdFunc( Pt[:A][t], tmpK1[t], tmpL2[t], Pt[:β][t], κ = Pc[:κ], GetPrice = true, TechType = "Hicks"  )[1]
            end # fiscal, consumption & labor-agg end

            # Section: Government debt 政府债务
            # NOTE: in this paper, D is always zero (self-supported fiscal budget)
            for t in idx1toTminus1
                # 1. government outstanding debt, through the capital market 政府负债，通过资本市场均衡
                Dt[:D][t] = tmpK1[t] - sum( Ps[:N][t,idxS] .* ( Dst[:a][t,idxS] .+ Dst[:Φ][t,idxS] ) )
                Dt[:D][t] = max( 0.0, min( Dt[:D][t], Dt[:Y][t] * Pt[:D2Ycap][t] ) )
                Dt[:D2Y][t] = Dt[:D][t] / Dt[:Y][t]
            end # government

            # Section: Government purchase, Investment & Capital 政府购买，投资与资本市场
            # NOTE: because we use D_{t+1} to compute G_{t} (well ... though doesnot matter in this paper but for generality)
            #       we seperate D_{t} from other variables, in a single loop
            for t in idx1toTminus1
                # 1. government purchase 政府购买
                Dt[:G][t] = Dt[:TRw][t] + Dt[:TRc][t] + Dt[:D][t+1] -
                            Dt[:LI][t] - Dt[:r][t] * Dt[:D][t]
                # 2. investment 投资
                Dt[:I][t] = Dt[:Y][t] - Dt[:G][t] - Dt[:C][t]

            end # gov-purchase, investment & capital-agg end

            # Section: Reaggregate Capital 汇总资本
            # NOTE: notice the index of looping! :)
            # NOTE: re-aggregate capital through the dynamics of capital growth
            for t in env.T-1:-1:2
                tmpK2[t]  = tmpK1[t+1] - Dt[:I][t]
                tmpK2[t] /= 1.0 - Pc[:κ]
            end

            # Convergence check & go to the next loop 收敛检查
            # 1. compute errors, using a minor number to avoid exact zeros
            # NOTE: we care the maximum of the errors of a complete path 我们关心整条序列误差的最大值
            local Err = ( K = findmax( abs.(tmpK2 ./ tmpK1 .- 1.0 )[idx2toTminus1] )[1],
                          L = findmax( abs.(tmpL2 ./ tmpL1 .- 1.0 )[idx2toTminus1] )[1]   )
            # performance log 记录误差
                append!(PerfLog[:K], Err.K)
                append!(PerfLog[:L], Err.L)
            # 2. check
            if ( Err.K < atol ) & (Err.L < atol)  # converged
                # 1. save the converged K,L, using mean values
                Dt[:K][idx2toTminus1] = ( tmpK1 .+ tmpK2 )[idx2toTminus1] ./ 2
                Dt[:L][idx2toTminus1] = ( tmpL1 .+ tmpL2 )[idx2toTminus1] ./ 2
                # 2. print final summary (if not silently solved)
                if PrintMode != "silent"
                    println("\n\t+ Status: Converged")  # status
                    println("\t+ Round: ", idxiter)
                    println("\t+ Max Relative Errors: ",Err)  # print error first
                    # SummaryYear( t, Dt, Dst, Pt, Ps, env.S, env.Sr )  # print the details of the steady state
                end
                break  # ends iteration
            elseif idxiter == MaxIter # diverged or too few maximum rounds of loops
                # 1. save the converged K,L, using mean values
                Dt[:K][idx2toTminus1] = ( tmpK1 .+ tmpK2 )[idx2toTminus1] ./ 2
                Dt[:L][idx2toTminus1] = ( tmpL1 .+ tmpL2 )[idx2toTminus1] ./ 2
                # 2. print final summary (if not silently solved)
                if PrintMode != "silent"
                    println("\n\t+ Status: Maximum iteration reached, not converged")
                    println("\t+ Round: ", idxiter)
                    println("\t+ Max Relative Errors: ",Err)  # print error first
                    # SummaryYear( t, Dt, Dst, Pt, Ps, env.S, env.Sr )  # print the details of the steady state
                end
                break  # ends iteration

            else  # go to the next loop
                # 1. check if the updated labor supply touches the bottom of zero
                local chkidx = (tmpL2 .< 0.0) .| (tmpL2 .!= real.(tmpL2))
                tmpL2[chkidx] .= 0.001
                # 2. update labor supply
                tmpL1[idx2toTminus1] .+= StepLen .* ( tmpL2 .- tmpL1 )[idx2toTminus1]
                # 3. use the magic number to set the bottom of capital, which bounds the interest rate in a reasonable range
                local tmpKfloor = MagicNum .* findmin(tmpL1)[1]
                # 4. check & update capital supply
                chkidx = (tmpK2 .< tmpKfloor) .| (tmpK2 .!= real.(tmpK2))
                tmpK2[chkidx] .= tmpKfloor
                tmpK1[idx2toTminus1] .+= StepLen .* ( tmpK2 .- tmpK1 )[idx2toTminus1]
                # 5. print information, according to :PrintMode
                if PrintMode in ["full", "concise"]
                    println("\n\t+ Round: ", idxiter)
                    println("\t+ Max Relative Errors: ", Err)
                    # println("\t+ Updating of Kap: ",round.([ OriGuess[1], tmpK2, Dt[:K][t] ],digits=6) ) # then print the updating process of K,L
                    # println("\t+ Updating of Lab: ",round.([ OriGuess[2], tmpL2, Dt[:L][t] ],digits=6) )
                    # println("\t+ NOTE: guess -> re-aggregated -> updated")
                end
                # if PrintMode == "full"
                #     # SummaryYear( t, Dt, Dst, Pt, Ps, env.S, env.Sr )
                # end
                # 6. explicitly write continue (as a reminder without practical jobs)
                continue

            end  # branch ends

        end  # Gauss-Seidel ends
        # nominal returns
        if ReturnLog
            return PerfLog::Dict
        else
            return nothing
        end
    end
# -----------------

# ==============================================================================
# SECTION 3: Assistant Functions After-Transition-Search 转轨路径搜索后的辅助函数
# NOTE: in this section, we define some functions used when the search of transition path finished
#       本节我们定义一些完成转轨路径搜索后使用的辅助函数
#       these functions are used to do economic analysis
#       这些函数用于分析具体的经济问题
#       they compute variables based on a converged/realized transition path
#       这些函数的计算基于已经求解完的转轨路径，属于对结果的应用
# --------------
"""
    ProcAfterTransition!( t::Int, Dt::Dict, Dst::Dict, Pt::Dict, Ps::Dict, Pc::Dict, env::NamedTuple )

Computes variables based on the datasets after the search of transition path;
returns nothing;

Computes:
1. Aggregated Utility (Dt[:U], len = T)
2. GDP Growth Rate % (Dt[:GDPgrowth], len = T-1) [in percentile]
3. The rate of UEBMI-pool gap to GDP (Dt[:LI2Y], len = T)
4. The proportion of labor population (Dt[:WorkNRat], len = T)
5. Substitution Rate of PAYG pension (Dt[:SubstiRat], len = T)
6. Population (Dt[:N], len = T)

"""
function ProcAfterTransition!( Dt::Dict, Dst::Dict, Pt::Dict, Ps::Dict, Pc::Dict, env::NamedTuple )
    # 1. Aggregated Utility (Dt[:U], len = T)
    Dt[:U] = Array{Float64,1}()
    local tmpus = zeros( env.S )  # temp utility vector for a year (populi cruxes)
    for t in 1:env.T
        local tmplongLab = cat( Dst[:Lab][t,1:env.Sr], ones(env.S - env.Sr), dims = 1 )
        for s in 1:env.S
            tmpus[s] = House.u( Dst[:c][t,s], tmplongLab[s], q = Ps[:q][t,s], alpha = Pc[:α], gamma = Pc[:γ] )
        end
        push!( Dt[:U], sum(tmpus) )
    end

    # 2. GDP Growth Rate % (Dt[:GDPgrowth], len = T-1), in percentile
    # NOTE: Dt[:GDPgrowth][t] = the GDP growth rate in year t (from time t to time t+1)
    Dt[:GDPgrowth] = Array{Float64,1}()
    for t in 1:env.T-1
        push!( Dt[:GDPgrowth], 100 * ( Dt[:Y][t+1] / Dt[:Y][t] - 1.0 ) )
    end

    # 3. The rate of UEBMI-pool gap to GDP (Dt[:LI2Y], len = T)
    Dt[:LI2Y] = Dt[:LI] ./ Dt[:Y]

    # 4. The proportion of labor population (Dt[:WorkNRat], len = T)
    # NOTE: using [:,1] to convert a Array{Float64,2} to Array{Float64,1}
    Dt[:WorkNRat] = ( sum( Ps[:N][:,1:env.Sr], dims = 2 ) ./ sum( Ps[:N][:,1:env.S], dims = 2 ) )[:,1]

    # 5. Substitution Rate of PAYG Pension (Dt[:SubstiRat], len = T)
    # NOTE: cross-sectional
    #       替代率是截面的，即当年退休人的平均养老金给付（Λ_{t}）/当年平均总工资收入（L_{t} * w̄_{t} / N_{t}）
    Dt[:SubstiRat] = Array{Float64,1}()
    for t in 1:env.T
        push!( Dt[:SubstiRat], Dt[:Λ][t] / ( Dt[:L][t] * Dt[:w̄][t] / sum(Ps[:N][t,1:env.Sr]) )   )
    end

    # 6. Population in year t (Dt[:N], len = T)
    Dt[:N] = Array{Float64,1}()
    for t in 1:env.T
        push!( Dt[:N], sum(Ps[:N][t,1:env.S]) )
    end

    # 7. (Social) aggregated outpatient expenditure & inpatient expenditure
    Dt[:AggMA] = Array{Float64,1}()
    Dt[:AggMB] = Array{Float64,1}()
    for t in 1:env.T
        push!( Dt[:AggMA], sum( Ps[:N][t,1:env.S] .* Dst[:MA][t,1:env.S] ) )
        push!( Dt[:AggMB], sum( Ps[:N][t,1:env.S] .* Dst[:MB][t,1:env.S] ) )
    end

    # 8. Aggregated Income & Expenditure of the pooling account of UE-BMI
    Dt[:AggPoolExp] = Array{Float64,1}()
    Dt[:AggPoolIn] = Array{Float64,1}()
    for t in 1:env.T
        push!( Dt[:AggPoolExp], Dt[:AggMB][t] * ( 1.0 - Pt[:cpB][t] ) )
        push!( Dt[:AggPoolIn], Dt[:AggPoolExp][t] - Dt[:LI][t] )
    end

    return nothing
end

# ==============================================================================
end  # module ends
#


text:
using GeometricIntegrators.BasisFunctions
using GeometricIntegrators.CommonFunctions
using GeometricIntegrators.Config
using GeometricIntegrators.Integrators
using GeometricIntegrators.Quadratures
using GeometricIntegrators.Solvers
using GeometricIntegrators.TestProblems.HarmonicOscillatorProblem
using GeometricIntegrators.Utils
using Test

using GeometricIntegrators.TestProblems.HarmonicOscillatorProblem: Δt, nt, refx, refq, refp

iode = harmonic_oscillator_iode()

QGau4 = GaussLegendreQuadrature(4)
BGau4 = LagrangeBasis(nodes(QGau4))

### CGVI Integrators ###

cgint = IntegratorCGVI(iode, BGau4, QGau4, Δt)
cgsol = integrate(cgint, nt)

@test rel_err(cgsol.q, refx) < 1E-7

### DGVI Integrators ###

dgint = IntegratorDGVI(iode, BGau4, QGau4, Δt)
dgsol = integrate(dgint, nt)

@test rel_err(dgsol.q, refx) < 1E-7

ext: jl


text:
#___  ____       ___ _____   _   _       _ _
#|  \/  | |     |_  /  __ \ | | | |     (_| |
#| .  . | |       | | /  \/ | | | |_ __  _| |_ ___
#| |\/| | |       | | |     | | | | '_ \| | __/ _ \
#| |  | | |___/\__/ | \__/\ | |_| | | | | | || (_) |
#_______\_____\____/ \____/  _____|_|___|_____\_____ _____ _____
#| ___ \        (_)         | | \ \ / / / __  |  _  / __  |  _  |
#| |_/ _ __ ___  _  ___  ___| |_ \ V /  `' / /| |/' `' / /| |/' |
#|  __| '__/ _ \| |/ _ \/ __| __|/   \    / / |  /| | / / |  /| |
#| |  | | | (_) | |  __| (__| |_/ /^\ \ ./ /__\ |_/ ./ /__\ |_/ /
#\_|  |_|  \___/| |\___|\___|\__\/   \/ \_____/\___/\_____/\___/
#              _/ |
#             |__/
#
# This code is part of the proposal of the team "MLJC UniTo" - University of Turin
# for "ProjectX 2020" Climate Change for AI.
# The code is licensed under MIT 3.0
# Please read readme or comments for credits and further information.

# Compiler: Julia 1.5

# Short description of this file: Burger's equation Neural vs Numerical
#Implementation

using ModelingToolkit, NeuralPDE, DiffEqFlux, Flux, Plots, GalacticOptim, Test, Optim

cb = function (p,l)
    println("Current loss is: $l")
    return false
end

@parameters t, x, θ
@variables u(..)
@derivatives Dt'~t
@derivatives Dx'~x
@derivatives Dxx''~x

#2D PDE
eq  = Dt(u(t,x,θ)) + u(t,x,θ)*Dx(u(t,x,θ)) - (0.01/pi)*Dxx(u(t,x,θ)) ~ 0

# Initial and boundary conditions
bcs = [u(x,0,θ) ~ 2*x]
#u(0,x,θ) ~ -sin(pi*x),
       #u(t,-1,θ) ~ 0.,
       #u(t,1,θ) ~ 0.]

# Space and time domains
domains = [t ∈ IntervalDomain(0.0,10.0),
           x ∈ IntervalDomain(-1.0,1.0)]
# Discretization
dx = 0.1

# Neural network
chain = FastChain(FastDense(2,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1))

strategy = GridTraining()
discretization = PhysicsInformedNN(dx,chain,strategy=strategy)

indvars = [t,x]
depvars = [u]
dim = length(domains)

expr_pde_loss_function = build_loss_function(eq,indvars,depvars)
expr_bc_loss_functions = [build_loss_function(bc,indvars,depvars) for bc in bcs]

train_sets = generate_training_sets(domains,dx,bcs,indvars,depvars)

train_domain_set, train_bound_set, train_set= train_sets

phi = discretization.phi
autodiff = discretization.autodiff
derivative = discretization.derivative

initθ = discretization.initθ

pde_loss_function = get_loss_function(eval(expr_pde_loss_function),
                                      train_domain_set,
                                      phi,
                                      derivative,
                                      strategy)

bc_loss_function = get_loss_function(eval.(expr_bc_loss_functions),
                                     train_bound_set,
                                     phi,
                                     derivative,
                                     strategy)

function loss_function(θ,p)
    return pde_loss_function(θ) + bc_loss_function(θ)
end

f = OptimizationFunction(loss_function, GalacticOptim.AutoZygote())

prob = GalacticOptim.OptimizationProblem(f, initθ)

# optimizer
opt = Optim.BFGS()

res = GalacticOptim.solve(prob, opt; cb = cb, maxiters=15)

phi = discretization.phi

ts,xs = [domain.domain.lower:dx:domain.domain.upper for domain in domains]
u_predict_contourf = reshape([first(phi([t,x],res.minimizer)) for t in ts for x in xs] ,length(xs),length(ts))
plot(ts, xs, u_predict_contourf, linetype=:surface,title = "predict")

u_predict = [[first(phi([t,x],res.minimizer)) for x in xs] for t in ts ]
p1= plot(xs, u_predict[2],title = "t = 0.1");
p2= plot(xs, u_predict[6],title = "t = 0.5");
p3= plot(xs, u_predict[end],title = "t = 1");
plot(p1,p2,p3)

analytic_sol_func(x,t) = 2*x/(1+2*t)
#analytic_sol_func(x,t) = [exp(-t)*sin(pi*x), exp(-t)*cos(pi*x), (1+pi^2)*exp(-t)]
#u_real  = [[analytic_sol_func(x,t)[i] for t in ts for x in xs] for i in 1:3]
u_real  = [analytic_sol_func(x,t) for t in ts for x in xs]

#=
for i in 1:3
    p1 = Plots.plot(xs, ts, u_real[i], st=:surface,title = "u$i, analytic");
    Plots.plot(p1)
    Plots.savefig("sol_real$i")
end
=#

plot(ts, xs, u_real, linetype=:surface,title = "real_contour")


text:
import Term.Trees: Tree
import Term: hLine, vLine

function pprint(tree::Tree)
    print(" " * hLine(tree.measure.w; style = "red"))
    print(vLine(tree.measure.h; style = "red") * tree)
    return println(tree.measure)
end

tree_dict = Dict(
    "nested\nasdasdsada\nsdasdsad\nsadasdasdasd\nsadasasdasdas\nsfsdfdsfds\nfdsfsdfsdf\ndsfsdfd" =>
        Dict("n1" => 1, "n2" => 2),
)

tree_dict_1 =
    Dict("nested" => Dict("n1" => 1, "n2" => 2), "nested2" => Dict("n1" => "a", "n2" => 2))

tree_dict_2 = Dict(
    "nested" => Dict("n1" => 1, "n2" => 2),
    "leaf2" => 2,
    "leaf" => 2,
    "leafme" => "v",
    "canopy" => "test",
)

tree_dict_3 = Dict(
    "nested" => Dict(
        "deeper" => Dict("aleaf" => "unbeliefable", "leaflet" => "level 3"),
        "n2" => Int,
        "n3" => 1 + 2,
    ),
    "nested2" => Dict("n1" => "a", "n2" => 2),
)

tree_dict_4 = Dict(
    "nested" => Dict(
        "deeper" => Dict(
            "aleaf" => "unbeliefable",
            "leaflet" => "level 3",
            "sodeep" => Dict("a" => 4),
        ),
        "n2" => Int,
        "n3" => 1 + 2,
        "adict" => Dict("x" => 2),
    ),
    "nested2" => Dict("n1" => "a", "n2" => 2),
)

print("\n"^20)

for dd in (tree_dict, tree_dict_1, tree_dict_2, tree_dict_3, tree_dict_4)
    pprint(Tree(dd))
    print("\n")
end


text:
using TDD
using Test

@testset "TDD.jl" begin
    include("graph.jl")
    include("interval.jl")
end


text:
module ImageContainers

export
    ImageContainer,
    storeimage

include("show.jl")

end # module


text:
using StaticArrays: SArray, @SMatrix
using LinearAlgebra: I, norm

# This matrix-based implementation is very cool (IMO)
# and only ~25% slower than the reference implementation for addition, multiplication, subtraction
#
# We get a 10x speedup for division and abs (although, for abs in particular, base widens to guard against overflow)
struct ComplexNumber{T<:Real} <: Number
    mat::SArray{Tuple{2,2},T,2,4}
    ComplexNumber(mat) = new{eltype(mat)}(mat)
    ComplexNumber(real::T1,imag::T2) where {T1 <: Real, T2 <: Real} = new{promote_type(T1,T2)}(@SMatrix([real  -imag; imag real]))
end

Base.:*(l::ComplexNumber,r::ComplexNumber) = ComplexNumber(l.mat * r.mat)
Base.:+(l::ComplexNumber,r::ComplexNumber) = ComplexNumber(l.mat + r.mat)
Base.:-(l::ComplexNumber,r::ComplexNumber) = ComplexNumber(l.mat - r.mat)
Base.:/(l::ComplexNumber,r::ComplexNumber) = ComplexNumber(l.mat * inv(r.mat))
Base.:*(z::ComplexNumber,r::Real) = ComplexNumber(r*z.mat)
Base.:*(r::Real, z::ComplexNumber) = z*r
Base.:+(z::ComplexNumber,r::Real) = ComplexNumber(r*I + z.mat)
Base.:+(r::Real, z::ComplexNumber) = z + r

Base.show(io::IO, z::ComplexNumber) = print(io, "ComplexNumber(", Re(z), ", ", Im(z), ")")
Base.:≈(l::ComplexNumber,r::ComplexNumber) = l.mat ≈ r.mat # I had to make one of the division tests ≈ to pass it
Base.:(==)(l::ComplexNumber,r::ComplexNumber) = l.mat == r.mat

Base.abs(z::ComplexNumber) = norm(z.mat[:,1])

# Sorry about the capital letter. That's what maths does.
@inline Re(z::ComplexNumber) = @inbounds z.mat[1,1]
@inline Im(z::ComplexNumber) = @inbounds z.mat[2,1]

Base.real(z::ComplexNumber) = Re(z)
Base.imag(z::ComplexNumber) = Im(z)

Base.conj(z::ComplexNumber) = ComplexNumber(transpose(z.mat))

Base.exp(z::ComplexNumber) = ComplexNumber(exp(z.mat))

# There's not an iota of a chance that ι will take off as the sensible symbol to use for the imaginary unit
const ι = jm = ComplexNumber(0,1)


text:
@kwargs function emul!(
    w::GBVector,
    u::GBVector,
    v::GBVector;
    op::MonoidSemiringOrBinary = BinaryOps.TIMES
)
    size(w) == size(u) == size(v) || throw(DimensionMismatch())
    op = getoperator(op, eltype(w))
    if op isa libgb.GrB_Semiring
        libgb.GrB_Vector_eWiseMult_Semiring(w, mask, accum, op, u, v, desc)
        return w
    elseif op isa libgb.GrB_Monoid
        libgb.GrB_Vector_eWiseMult_Monoid(w, mask, accum, op, u, v, desc)
        return w
    elseif op isa libgb.GrB_BinaryOp
        libgb.GrB_Vector_eWiseMult_BinaryOp(w, mask, accum, op, u, v, desc)
        return w
    else
        error("Unreachable")
    end
end

@kwargs function emul(
    u::GBVector,
    v::GBVector;
    op::MonoidSemiringOrBinary = BinaryOps.TIMES
)
    if op isa OpUnion
        t = ztype(op)
    else
        t = optype(eltype(u), eltype(v))
    end
    w = GBVector{t}(size(u))
    return emul!(w, u, v; op, mask , accum, desc)
end

@kwargs function emul!(
    C::GBMatrix,
    A::GBMatOrTrans,
    B::GBMatOrTrans;
    op::MonoidSemiringOrBinary = BinaryOps.TIMES
)
    size(C) == size(A) == size(B) || throw(DimensionMismatch())
    A, desc, B = _handletranspose(A, desc, B)
    op = getoperator(op, eltype(C))
    if op isa libgb.GrB_Semiring
        libgb.GrB_Matrix_eWiseMult_Semiring(C, mask, accum, op, A, B, desc)
        return C
    elseif op isa libgb.GrB_Monoid
        libgb.GrB_Matrix_eWiseMult_Monoid(C, mask, accum, op, A, B, desc)
        return C
    elseif op isa libgb.GrB_BinaryOp
        libgb.GrB_Matrix_eWiseMult_BinaryOp(C, mask, accum, op, A, B, desc)
        return C
    end
end

@kwargs function emul(
    A::GBMatOrTrans,
    B::GBMatOrTrans;
    op::MonoidSemiringOrBinary = BinaryOps.TIMES
)
    if op isa OpUnion
        t = ztype(op)
    else
        t = optype(eltype(A), eltype(B))
    end
    C = GBMatrix{t}(size(A))
    return emul!(C, A, B; op, mask, accum, desc)
end

@kwargs function eadd!(
    w::GBVector,
    u::GBVector,
    v::GBVector;
    op::MonoidSemiringOrBinary = BinaryOps.PLUS
)
    size(w) == size(u) == size(v) || throw(DimensionMismatch())
    op = getoperator(op, eltype(w))
    if op isa libgb.GrB_Semiring
        libgb.GrB_Vector_eWiseAdd_Semiring(w, mask, accum, op, u, v, desc)
        return w
    elseif op isa libgb.GrB_Monoid
        libgb.GrB_Vector_eWiseAdd_Monoid(w, mask, accum, op, u, v, desc)
        return w
    elseif op isa libgb.GrB_BinaryOp
        libgb.GrB_Vector_eWiseAdd_BinaryOp(w, mask, accum, op, u, v, desc)
        return w
    else
        error("Unreachable")
    end
end

@kwargs function eadd(
    u::GBVector,
    v::GBVector;
    op::MonoidSemiringOrBinary = BinaryOps.PLUS
)
    if op isa OpUnion
        t = ztype(op)
    else
        t = optype(eltype(u), eltype(v))
    end
    w = GBVector{t}(size(u))
    return eadd!(w, u, v; op, mask, accum, desc)
end

@kwargs function eadd!(
    C::GBMatrix,
    A::GBMatOrTrans,
    B::GBMatOrTrans;
    op::MonoidSemiringOrBinary = BinaryOps.PLUS
)
    size(C) == size(A) == size(B) || throw(DimensionMismatch())
    A, desc, B = _handletranspose(A, desc, B)
    op = getoperator(op, eltype(C))
    if op isa libgb.GrB_Semiring
        libgb.GrB_Matrix_eWiseAdd_Semiring(C, mask, accum, op, A, B, desc)
        return C
    elseif op isa libgb.GrB_Monoid
        libgb.GrB_Matrix_eWiseAdd_Monoid(C, mask, accum, op, A, B, desc)
        return C
    elseif op isa libgb.GrB_BinaryOp
        libgb.GrB_Matrix_eWiseAdd_BinaryOp(C, mask, accum, op, A, B, desc)
        return C
    else
        error("Unreachable")
    end
end

@kwargs function eadd(
    A::GBMatOrTrans,
    B::GBMatOrTrans;
    op::MonoidSemiringOrBinary = BinaryOps.PLUS
)
    if op isa OpUnion
        t = ztype(op)
    else
        t = optype(eltype(A), eltype(B))
    end
    C = GBMatrix{t}(size(A))
    return eadd!(C, A, B; op, mask, accum, desc)
end

@kwargs Base.broadcasted(::typeof(+), u::GBVector, v::GBVector) = eadd(u, v)
@kwargs Base.broadcasted(::typeof(*), u::GBVector, v::GBVector) = emul(u, v)

@kwargs Base.broadcasted(::typeof(+), A::GBMatOrTrans, B::GBMatOrTrans) = eadd(A, B)
@kwargs Base.broadcasted(::typeof(*), A::GBMatOrTrans, B::GBMatOrTrans) = emul(A, B)


text:
@testset "mesh" begin
    eM = eMesh{Tri,Tet}()
    push!(eM.point, SVector{3,Float64}(NaN, NaN, NaN))
    @test_throws AssertionError verify_mesh(eM)
    push!(eM.ϵ, NaN)

    append!(eM, eMesh_half_plane())
    mesh_remove_unused_points!(eM)
    @test length(eM.point) == 4

    append!(eM, eMesh_half_plane())
    mesh_inplace_rekey!(eM)
    @test length(eM.point) == 4

    @test !isempty(eM)
    @test_throws MethodError empty(eM)
    eM = as_tri_eMesh(eM)
    empty!(eM)
    @test isempty(eM)

    eM = as_tri_eMesh(eMesh_box())
    @test area(crop_mesh(eM, SMatrix{1,4,Float64,4}(0.0, 0.0, +1.0, 0.0))) ≈ 12.0
    @test area(crop_mesh(eM, SMatrix{1,4,Float64,4}(0.0, 0.0, -1.0, 0.0))) ≈ 12.0
    @test volume(eMesh_box()) ≈ 8.0
end

@testset "verify_mesh" begin
    p1 = SVector(0.0, 0.0, 0.0)
    p2 = SVector(1.0, 0.0, 0.0)
    p3 = SVector(0.0, 1.0, 0.0)
    p4 = SVector(0.0, 0.0, 1.0)

    ### TRI
    eM_tri = eMesh{Tri,Nothing}()
    @test verify_mesh(eM_tri) == nothing

    push!(eM_tri.tri, SVector(1, 2, 3))
    @test_throws BoundsError verify_mesh(eM_tri)  # no points

    push!(eM_tri.point, p1, p2, p3)
    @test verify_mesh(eM_tri) == nothing

    push!(eM_tri.tri, SVector(1, 2, 0))
    @test_throws BoundsError verify_mesh(eM_tri)  # index less than 1

    ### TET
    eM_tet = eMesh{Nothing, Tet}()
    @test verify_mesh(eM_tet) == nothing

    push!(eM_tet.tet, SVector(1, 2, 3, 4))
    @test_throws BoundsError verify_mesh(eM_tet)  # no points

    eM_tet = eMesh{Nothing, Tet}()
    push!(eM_tet.point, p1, p2, p3, p4)
    @test_throws ErrorException verify_mesh(eM_tet)  # no ϵ

    push!(eM_tet.ϵ, 0.0, 0.0, 0.0, -1.0)
    @test_throws ErrorException verify_mesh(eM_tet)  # negative ϵ

    eM_tet.ϵ[end] = 1.0
    @test verify_mesh(eM_tet) == nothing

    push!(eM_tet.tet, SVector(1, 2, 4, 3))
    @test_throws ErrorException verify_mesh(eM_tet)  # negative volume)

    eM_tet.tet[1] = SVector(1, 2, 3, 4)
    @test verify_mesh(eM_tet) == nothing
end

@testset "delete_triangles" begin
	random_points = [randn(SVector{3,Float64}) for _ = 1:3]
	i3_tri = SVector{3,Int64}(1,2,3)
	i3_tri_oppose = SVector{3,Int64}(i3_tri[2], i3_tri[3], i3_tri[1])

	eM_delete = eMesh(random_points, [i3_tri, i3_tri_oppose])
	@test 2 == delete_triangles!(eM_delete)
	@test isempty(eM_delete)

	eM_delete = eMesh(random_points, [i3_tri, SVector{3,Int64}(2,3,4)])
	@test 0 == delete_triangles!(eM_delete)
	@test !isempty(eM_delete)

	eM_delete = eMesh(random_points, [i3_tri, i3_tri])
	@test_throws ErrorException delete_triangles!(eM_delete)

	eM_delete = eMesh(random_points, [i3_tri, i3_tri_oppose, i3_tri])
	@test_throws ErrorException delete_triangles!(eM_delete)
end

two = 2.0
@testset "half_plane" begin
    eM_hp_f = eMesh_half_plane(two, false)
    eM_hp_t = eMesh_half_plane(two, true)
    z_val_f = [k[3] for k = get_point(eM_hp_f)]
    z_val_t = [k[3] for k = get_point(eM_hp_t)]
    @test z_val_f == z_val_t
    @test -two == minimum(z_val_f)
    @test 0.0 == maximum(z_val_f)
    @test 4 == n_point(eM_hp_f)
    for k = 1:3
        @test -0.5 ≈ dot(get_point(eM_hp_f)[k], get_point(eM_hp_f)[mod1(k + 1, 3)])
    end
    @test n_tri(eM_hp_f) == 1
    @test n_tri(eM_hp_t) == 4
    @test n_tet(eM_hp_f) == n_tet(eM_hp_t) == 1
	@test verify_mesh(eM_hp_f) == nothing
	@test verify_mesh(as_tri_eMesh(eM_hp_t)) == nothing
	@test verify_mesh(as_tet_eMesh(eM_hp_t)) == nothing
end

@testset "sphere" begin
    for k = 1:4
        eM_sphere = eMesh_sphere(two, k)
        n_zero = 0
        @test n_tri(eM_sphere) == n_tet(eM_sphere) == 20 * k^2
        for k = 1:n_point(eM_sphere)
            norm_point = norm(get_point(eM_sphere)[k])
            if norm_point == 0.0
                n_zero += 1
            else
                @test two ≈ norm_point
            end
        end
        @test n_zero == 1
        (k == 1) && @test n_point(eM_sphere) == 1 + 12
        (k == 2) && @test n_point(eM_sphere) == 1 + 12 + 30  # each edge bisected
        (k == 3) && @test n_point(eM_sphere) == 1 + 12 + 2 * 30 + 20  # each edge bisected + point in center
        (k == 4) && @test n_point(eM_sphere) == 1 + 12 + 3 * 30 + 3 * 20  # each edge bisected + two more points in center
		@test verify_mesh(eM_sphere) == nothing
    end
end

@testset "polygon circle" begin
	for n = 3:64
		eM_circle = PressureFieldContact.Geometry.create_2D_circle(1.0, n=n)
		b = sin(pi / n)  # distance from center to bisector
		a = sqrt(1.0^2 - b^2)  # distance from top of triangle to point on circle
		my_area = n * a * b
		@test my_area ≈ area(eM_circle)
		@test verify_mesh(eM_circle) == nothing
	end
end

@testset "prism" begin
	thick = 0.1
	p1 = SVector(0.0, 0.0, 0.0)
	p2 = SVector(1.0, 0.0, 0.0)
	p3 = SVector(0.0, 1.0, 0.0)
	eM_tri = eMesh([p1, p2, p3], [SVector(1, 2, 3)])
	eM_prism = extrude_mesh(eM_tri, thick)
	vol_ = volume(eM_prism)
	area_ = area(eM_prism)
	@test vol_ ≈ thick * 0.5
	@test area_ ≈ 1.0 + thick * (2.0 + sqrt(2))
	@test verify_mesh(eM_prism) == nothing
end

@testset "cylinder" begin
	for n = 3:64
		rad = 1.0
		height = 1000.0
		eM_cylinder = eMesh_cylinder(rad, height, n=n)
		eM_circle = PressureFieldContact.Geometry.create_2D_circle(rad, n=n)
		area_circle = area(eM_circle)
		@test volume(eM_cylinder) ≈ height * area_circle
		@test area(eM_cylinder) ≈ 2 * area_circle + height * 2 * n * sin(pi / n)
		@test verify_mesh(eM_cylinder) == nothing
	end
end


text:
export Domain,GlobalData,updateStates!,updateDomainStateBoundary!,getExternalForce,
    setNeumannBoundary!, setGeometryPoints!

@doc """
GlobalData

Store data for finite element update, assume the problem has n freedoms

- `state`: Float64[n],  displacement array at the current time, only for active freedoms,
   the ordering is based on the equation number, here ndims=2 is the dimension of the problem space
- `state`: Float64[n],  displacement array at the previous time
- `velo`: Float64[n],  velocity array at the current
- `acce`: Float64[n],  acceleration array at the current
- `time`: float, current time
- `M`: Float64[n,n] spares mass matrix
- `Mlumped`: Float64[n] lumped mass array
- `MID`: Float64[n, nd1] off-diagonal part of the mass matrix, between the active freedoms and the time-dependent Dirichlet freedoms, assume there are nd time-dependent Dirichlet freedoms
- `EBC_func`: function Float64:t-> Float64[n_d1] float array, time-dependent Dirichlet boundary condition (ordering is direction first then node number, u1, u3, ... v1, v4 ...)
- `FBC_func`: function Float64:t-> Float64[n_t1] float array, time-dependent load boundary condition (ordering is direction first then node number, u1, u3, ... v1, v4 ...)
"""->
mutable struct GlobalData
    state::Array{Float64}    #u
    Dstate::Array{Float64}   #uk
    velo::Array{Float64}     #∂u
    acce::Array{Float64}     #∂∂u
    time::Float64
    M::Union{SparseMatrixCSC{Float64,Int64},Array{Float64}}
    Mlumped::Array{Float64}
    MID::Array{Float64}

    EBC_func::Union{Function,Nothing}  #time dependent Dirichlet boundary condition
    FBC_func::Union{Function,Nothing}  #time force load boundary condition

end

function GlobalData(state::Array{Float64},Dstate::Array{Float64},velo::Array{Float64},acce::Array{Float64}, neqs::Int64,
        EBC_func::Union{Function, Nothing}=nothing, FBC_func::Union{Function, Nothing}=nothing)
    time = 0.0
    M = Float64[]
    Mlumped = Float64[]
    MID = Float64[]
    GlobalData(state, Dstate, velo, acce, time, M, Mlumped, MID, EBC_func, FBC_func)
end

@doc """
Domain

Date structure for the computatational domain

- `nnodes`: Int64, number of nodes (each quadratical quad element has 9 nodes)
- `nodes`: Float64[nnodes, ndims], coordinate array of all nodes
- `neles`: number of elements
- `elements`: Element[neles], element array, each element is a struct
- `ndims`: Int64, dimension of the problem space
- `state`: Float64[nnodes*ndims] current displacement of all nodal freedoms, Float64[1:nnodes] are for the first direction.
- `Dstate`: Float64[nnodes*ndims] previous displacement of all nodal freedoms, Float64[1:nnodes] are for the first direction.
- 'LM':  Int64[neles][ndims], LM(e,d) is the global equation number(active freedom number) of element e's d th freedom,
         -1 means fixed (time-independent) Dirichlet
         -2 means time-dependent Dirichlet
- 'DOF': Int64[neles][ndims], DOF(e,d) is the global freedom number of element e's d th freedom
- 'ID':  Int64[nnodes, ndims], ID(n,d) is the equation number(active freedom number) of node n's dth freedom,
         -1 means fixed (time-independent) Dirichlet
         -2 means time-dependent Dirichlet
- 'neqs':  Int64,  number of equations or active freedoms
- 'eq_to_dof':  Int64[neqs], map from to equation number(active freedom number) to the freedom number (Int64[1:nnodes] are for the first direction)
- 'dof_to_eq':  Bool[nnodes*ndims], map from freedom number(Int64[1:nnodes] are for the first direction) to booleans (active freedoms(equation number) are true)
- 'EBC':  Int64[nnodes, ndims], EBC[n,d] is the displacement boundary condition of node n's dth freedom,
           -1 means fixed(time-independent) Dirichlet boundary nodes
           -2 means time-dependent Dirichlet boundary nodes
- 'g':  Float64[nnodes, ndims], values for fixed(time-independent) Dirichlet boundary conditions of node n's dth freedom,
- 'FBC': Int64[nnodes, ndims], FBC[n,d] is the force load boundary condition of node n's dth freedom,
           -1 means constant(time-independent) force load boundary nodes
           -2 means time-dependent force load boundary nodes
- 'fext':  Float64[neqs], constant(time-independent) force load boundary conditions for these freedoms
- 'time': Float64, current time
- 'npoints': Int64, number of points (each quadratical quad element has 4 points, npoints==nnodes, when porder==1)
- 'node_to_point': Int64[nnodes]:map from node number to point point, -1 means the node is not a geometry point

- 'ii_stiff': Int64[], first index of the sparse matrix representation of the stiffness matrix
- 'jj_stiff': Int64[], second index of the sparse matrix representation of the stiffness matrix
- 'vv_stiff_ele_indptr': Int64[], Int64[e] is the first index entry for the e's element of the sparse matrix representation of the stiffness matrix
- 'vv_stiff': Float64[], values of the sparse matrix representation of the stiffness matrix

- 'ii_dfint_dstress': Int64[], first index of the sparse matrix representation of the dfint_dstress matrix
- 'jj_dfint_dstress': Int64[], second index of the sparse matrix representation of the dfint_dstress matrix
- 'vv_dfint_dstress_ele_indptr': Int64[], Int64[e] is the first index entry for the e's element of the sparse matrix representation of the dfint_dstress matrix
- 'vv_dfint_dstress': Float64[], values of the sparse matrix representation of the dfint_dstress matrix

- 'ii_dstrain_dstate': Int64[], first index of the sparse matrix representation of the dstrain_dstate matrix
- 'jj_dstrain_dstate': Int64[], second index of the sparse matrix representation of the dstrain_dstate matrix
- 'vv_dstrain_dstate_ele_indptr': Int64[], Int64[e] is the first index entry for the e's element of the sparse matrix representation of the stiffness matrix
- 'vv_dstrain_dstate': Float64[], values of the sparse matrix representation of the dstrain_dstate matrix

- 'history': Dict{String, Array{Array{Float64}}}, dictionary between string and its time-histories quantity Float64[ntime][]
"""->
mutable struct Domain
    nnodes::Int64
    nodes::Array{Float64}
    neles::Int64
    elements::Array
    ndims::Int64
    state::Array{Float64}
    Dstate::Array{Float64}
    LM::Array{Array{Int64}}
    DOF::Array{Array{Int64}}
    ID::Array{Int64}
    neqs::Int64
    eq_to_dof::Array{Int64}
    dof_to_eq::Array{Bool}
    EBC::Array{Int64}  # Dirichlet boundary condition
    g::Array{Float64}  # Value for Dirichlet boundary condition
    FBC::Array{Int64}  # Nodal force boundary condition
    fext::Array{Float64}  # Value for Nodal force boundary condition
    time::Float64

    npoints::Int64     # number of mesh points(the same as nodes, when porder==1)
    node_to_point::Array{Int64} # map from node to point, -1 means the node is not a geometry point

    ii_stiff::Array{Int64}
    jj_stiff::Array{Int64}
    vv_stiff_ele_indptr::Array{Int64}
    vv_stiff::Array{Float64}

    ii_dfint_dstress::Array{Int64}
    jj_dfint_dstress::Array{Int64}
    vv_dfint_dstress_ele_indptr::Array{Int64}
    vv_dfint_dstress::Array{Float64}

    ii_dstrain_dstate::Array{Int64}
    jj_dstrain_dstate::Array{Int64}
    vv_dstrain_dstate_ele_indptr::Array{Int64}
    vv_dstrain_dstate::Array{Float64}

    history::Dict{String, Array{Array{Float64}}}

end

function Base.:copy(g::Union{GlobalData, Domain})
    names = fieldnames(g)
    args = [copy(getproperty(g, n)) for n in names]
    GlobalData(args...)
end

@doc """
    Creating a finite element domain

    - `nodes`: Float64[nnodes, ndims], coordinate array of all nodes
    - `elements`: Element[neles], element array, each element is a struct
    - `ndims`: Int64, dimension of the problem space
    - 'EBC':  Int64[nnodes, ndims], EBC[n,d] is the displacement boundary condition of node n's dth freedom,
           -1 means fixed(time-independent) Dirichlet boundary nodes
           -2 means time-dependent Dirichlet boundary nodes
    - 'g':  Float64[nnodes, ndims], values for fixed(time-independent) Dirichlet boundary conditions of node n's dth freedom,
    - 'FBC': Int64[nnodes, ndims], FBC[n,d] is the force load boundary condition of node n's dth freedom,
           -1 means constant(time-independent) force load boundary nodes
           -2 means time-dependent force load boundary nodes
    - 'f':  Float64[nnodes, ndims], values for constant(time-independent) force load boundary conditions of node n's dth freedom,

    Return: Domain, the finite element domain
"""->
function Domain(nodes::Array{Float64}, elements::Array, ndims::Int64, EBC::Array{Int64}, g::Array{Float64}, FBC::Array{Int64}, f::Array{Float64})
    nnodes = size(nodes,1)
    neles = size(elements,1)
    state = zeros(nnodes * ndims)
    Dstate = zeros(nnodes * ndims)
    LM = Array{Int64}[]
    DOF = Array{Int64}[]
    ID = Int64[]
    neqs = 0
    eq_to_dof = Int64[]
    dof_to_eq = zeros(Bool, nnodes * ndims)
    fext = Float64[]

    npoints = -1
    node_to_point = Int64[]

    history = Dict("state"=>Array{Float64}[], "acc"=>Array{Float64}[], "fint"=>Array{Float64}[],
                "fext"=>Array{Float64}[], "strain"=>[], "stress"=>[], "time"=>Array{Float64}[])

    domain = Domain(nnodes, nodes, neles, elements, ndims, state, Dstate,
    LM, DOF, ID, neqs, eq_to_dof, dof_to_eq,
    EBC, g, FBC, fext, 0.0, npoints, node_to_point,
    Int64[], Int64[], Int64[], Float64[],
    Int64[], Int64[], Int64[], Float64[],
    Int64[], Int64[], Int64[], Float64[], history)

    #set fixed(time-independent) Dirichlet boundary conditions
    setDirichletBoundary!(domain, EBC, g)
    #set constant(time-independent) force load boundary conditions
    setNeumannBoundary!(domain, FBC, f)

    assembleSparseMatrixPattern!(domain)

    domain
end

@doc """
    In the constructor
    Update the node_to_point map

    - `self`: Domain, finit element domain
    - 'npoints': Int64, number of points (each quadratical quad element has 4 points, npoints==nnodes, when porder==1)
    - 'node_to_point': Int64[nnodes]:map from node number to point point, -1 means the node is not a geometry point

""" ->
function setGeometryPoints!(self::Domain, npoints::Int64, node_to_point::Array{Int64})
    self.npoints = npoints
    self.node_to_point = node_to_point
end

@doc """
    Update current step strain, stress in the history map of the Domain
    strain and stress are both Float[ngp, nstrain], ngp is the number of Gaussian quadrature point
    nstrain=1 for 1D and nstrain = 3 for 2D
    The strain is in Voigt notation
""" ->
function commitHistory(domain::Domain)
    for e in domain.elements
        commitHistory(e)
    end

    # 1D, nstrain=1; 2D, nstrain=3
    eledim = domain.elements[1].eledim
    nstrain = div((eledim + 1)*eledim, 2)
    ngp = domain.neles * length(domain.elements[1].weights)
    if nstrain==1
        strain = zeros(ngp)
        stress = zeros(ngp)
        k = 1
        for e in domain.elements
            for igp in e.mat
                strain[k] = igp.ε0
                stress[k] = igp.σ0
                k += 1
            end
        end
    else
        strain = zeros(ngp, nstrain)
        stress = zeros(ngp, nstrain)
        k = 1
        for e in domain.elements
            for igp in e.mat
                strain[k,:] = igp.ε0
                stress[k,:] = igp.σ0
                k += 1
            end
        end
    end

    push!(domain.history["strain"], strain)
    push!(domain.history["stress"], stress)
end

@doc """
    In the constructor
    Update the fixed (time-independent Dirichlet boundary) state entries
    Build LM, and DOF array
    - 'self': Domain
    - 'EBC':  Int64[nnodes, ndims], EBC[n,d] is the displacement boundary condition of node n's dth freedom,
           -1 means fixed(time-independent) Dirichlet boundary nodes
           -2 means time-dependent Dirichlet boundary nodes
    - 'g':  Float64[nnodes, ndims], values for fixed(time-independent) Dirichlet boundary conditions of node n's dth freedom,
    ''
""" ->
function setDirichletBoundary!(self::Domain, EBC::Array{Int64}, g::Array{Float64})

    # ID(n,d) is the global equation number of node n's dth freedom,
    # -1 means fixed (time-independent) Dirichlet
    # -2 means time-dependent Dirichlet

    nnodes, ndims = self.nnodes, self.ndims
    neles, elements = self.neles, self.elements
    #ID = zeros(Int64, nnodes, ndims) .- 1

    ID = copy(EBC)

    eq_to_dof, dof_to_eq = Int64[], zeros(Bool, nnodes * ndims)
    neqs = 0
    for idof = 1:ndims
      for inode = 1:nnodes
          if (EBC[inode, idof] == 0)
              neqs += 1
              ID[inode, idof] = neqs
              push!(eq_to_dof,inode + (idof-1)*nnodes)
              dof_to_eq[(idof - 1)*nnodes + inode] = true
          elseif (EBC[inode, idof] == -1)
              #update state fixed (time-independent) Dirichlet boundary conditions
              self.state[inode + (idof-1)*nnodes] = g[inode, idof]
          end
      end
    end

    self.ID, self.neqs, self.eq_to_dof, self.dof_to_eq = ID, neqs, eq_to_dof, dof_to_eq

    # LM(e,d) is the global equation number of element e's d th freedom
    LM = Array{Array{Int64}}(undef, neles)
    for iele = 1:neles
      el_nodes = getNodes(elements[iele])
      ieqns = ID[el_nodes, :][:]
      LM[iele] = ieqns
    end
    self.LM = LM

    # DOF(e,d) is the global dof number of element e's d th freedom

    DOF = Array{Array{Int64}}(undef, neles)
    for iele = 1:neles
      el_nodes = getNodes(elements[iele])
      DOF[iele] = [el_nodes;[idof + nnodes for idof in el_nodes]]
    end
    self.DOF = DOF

end

@doc """
    In the constructor
    Update the external force vector fext
    - 'self': Domain
    - 'FBC': Int64[nnodes, ndims], FBC[n,d] is the force load boundary condition of node n's dth freedom,
           -1 means constant(time-independent) force load boundary nodes
           -2 means time-dependent force load boundary nodes
    - 'f':  Float64[nnodes, ndims], values for constant(time-independent) force load boundary conditions of node n's dth freedom,

""" ->
function setNeumannBoundary!(self::Domain, FBC::Array{Int64}, f::Array{Float64})

    fext = zeros(Float64, self.neqs)
    # ID(n,d) is the global equation number of node n's dth freedom, -1 means no freedom

    nnodes, ndims, ID = self.nnodes, self.ndims, self.ID
    for idof = 1:ndims
      for inode = 1:nnodes
          if (FBC[inode, idof] == -1)
              @assert ID[inode, idof] > 0
              fext[ID[inode, idof]] += f[inode, idof]
          end
        end
    end
    self.fext = fext
end

@doc """
    At each time step
    Update the state and Dstate in Domain from GlobalData
    state and Dstate in GlobalData are only for active freedoms (equations)
    state and Dstate in Domain are only for all freedoms, they are used for constructing the internal force and/or stiffness matrix
    update the state and acc history of the Domain
    - 'self': Domain
    - 'globaldat': GlobalData
""" ->
function updateStates!(self::Domain, globaldat::GlobalData)

    self.state[self.eq_to_dof] = globaldat.state[:]

    self.time = globaldat.time
    push!(self.history["state"], copy(self.state))
    push!(self.history["acc"], copy(globaldat.acce))

    updateDomainStateBoundary!(self, globaldat)

    self.Dstate = self.state[:]
end

# @doc """
#     :param self: Domain
#     :param state : 1D array to convert
#     :param compress_or_expand  : "Compress" or "Expand"

#     "Compress", the state has all freedoms on all nodes, remove these freedoms on EBC
#     "Expand",   the state has only active freedoms on active nodes (active means not prescribed),
#                 set these freedoms on EBC to 0

#     :return:
# """ ->
# function convertState(self::Domain, state::Array{Float64}, compress_or_expand::String)

#     if compress_or_expand == "Expand"
#         new_state = zeros(Float64, self.nnodes*self.ndims)
#         new_state[self.eq_to_dof] = state[:]
#         return new_state

#     elseif compress_or_expand == "Compress"
#         return state[self.eq_to_dof]

#     else
#         error("convertStats error, compress_or_expand is ", compress_or_expand)
#     end

# end

@doc """
    Update time dependent boundary of state fext in Domain based on the time-dependent boundary functions in GlobalData.
    - 'self': Domain
    - 'globaldat': GlobalData
    """ ->
function updateDomainStateBoundary!(self::Domain, globaldat::GlobalData)
    if globaldat.EBC_func != nothing
        disp, acce = globaldat.EBC_func(globaldat.time) # user defined time-dependent boundary
        dof_id = 0

        #update state of all nodes
        for idof = 1:self.ndims
            for inode = 1:self.nnodes
                if (self.EBC[inode, idof] == -2)
                    dof_id += 1
                    self.state[inode + (idof-1)*self.nnodes] = disp[dof_id]
                end
            end
        end
    end

    if globaldat.FBC_func != nothing
        ID = self.ID
        nodal_force = globaldat.FBC_func(globaldat.time) # user defined time-dependent boundary
        # @info nodal_force
        dof_id = 0
        #update fext for active nodes (length of neqs)
        for idof = 1:self.ndims
            for inode = 1:self.nnodes
                if (self.FBC[inode, idof] == -2)
                    dof_id += 1
                    @assert ID[inode, idof] > 0
                    self.fext[ID[inode, idof]] = nodal_force[dof_id]
                end
            end
        end
    end
end

@doc """
    Compute external force vector, including external force load and time-dependent Dirichlet boundary conditions
    The function needs to be called after
    function updateDomainStateBoundary!(self::Domain, globaldat::GlobalData)
    which computes the external force vector from external force load

    - 'self': Domain
    - 'globaldat': GlobalData
    - 'fext': Float64[neqs], container for the external force vector

    return external force vector at the present step
""" ->
function getExternalForce!(self::Domain, globaldat::GlobalData, fext::Array{Float64})
    fext[:] = self.fext
    if globaldat.EBC_func != nothing
        MID = globaldat.MID
        _, acce = globaldat.EBC_func(globaldat.time)

        fext -= MID * acce
    end
    fext
end

@doc """
    Get the coordinates of several nodes (possibly in one element)
    - 'self': Domain
    - 'el_nodes': Int64[n], node array

    Return: Float64[n, ndims], the coordinates of these nodes
""" ->
function getCoords(self::Domain, el_nodes::Array{Int64})
    return self.nodes[el_nodes, :]
end

@doc """
    Get the global freedom numbers of the element
    - 'self': Domain
    - 'iele': Int64, element number

    Return: Int64[], the global freedom numbers of the element (ordering in local element ordering)

""" ->
function getDofs(self::Domain, iele::Int64)
    return self.DOF[iele]
end

@doc """
    Get the equation numbers(active freedom numbers) of the element
    - 'self': Domain
    - 'iele': Int64, element number

    Return: Int64[], the equation numbers(active freedom numbers) of the element (ordering in local element ordering)

""" ->
function getEqns(self::Domain, iele::Int64)
    return self.LM[iele]
end

@doc """
    Get the displacements of several nodes (possibly in one element)
    - 'self': Domain
    - 'el_nodes': Int64[n], node array

    Return: Float64[n, ndims], the displacements of these nodes
""" ->
function getState(self::Domain, el_dofs::Array{Int64})
    return self.state[el_dofs]
end

@doc """
    Get the displacements of several nodes (possibly in one element) at the previous time step
    - 'self': Domain
    - 'el_nodes': Int64[n], node array

    Return: Float64[n, ndims], the displacements of these nodes at the previous time step
""" ->

function getDstate(self::Domain, el_dofs::Array{Int64})
    return self.Dstate[el_dofs]
end

@doc """
    Compute constant stiff, dfint_dstress, dstrain_dstate matrix patterns
    - 'self': Domain
"""->
function assembleSparseMatrixPattern!(self::Domain)

    neles = self.neles
    eledim = self.elements[1].eledim
    nstrain = div((eledim + 1)*eledim, 2)
    ngps_per_elem = length(self.elements[1].weights)
    neqs = self.neqs

    ii_stiff = Int64[]; jj_stiff = Int64[]; vv_stiff_ele_indptr = ones(Int64, neles+1);
    ii_dfint_dstress = Int64[]; jj_dfint_dstress = Int64[]; vv_dfint_dstress_ele_indptr = ones(Int64, neles+1);
    ii_dstrain_dstate = Int64[]; jj_dstrain_dstate = Int64[]; vv_dstrain_dstate_ele_indptr = ones(Int64, neles+1);

    neles = self.neles

    # Loop over the elements in the elementGroup
    for iele  = 1:neles
      element = self.elements[iele]

      el_eqns = getEqns(self,iele)

      el_dofs = getDofs(self,iele)

      el_state  = getState(self, el_dofs)

      gp_ids = (iele-1)*ngps_per_elem+1 : iele*ngps_per_elem

      # Assemble in the global array
      el_eqns_active = el_eqns .>= 1
      el_eqns_active_idx = el_eqns[el_eqns_active]
      # K[el_eqns[el_eqns_active], el_eqns[el_eqns_active]] += stiff[el_eqns_active,el_eqns_active]

      el_eqns_active_idx = el_eqns[el_eqns_active]

      for j = 1:length(el_eqns_active_idx)
        for i = 1:length(el_eqns_active_idx)
          push!(ii_stiff, el_eqns_active_idx[i])
          push!(jj_stiff, el_eqns_active_idx[j])
          #push!(vv_stiff, stiff_active[i,j])
        end
      end
      vv_stiff_ele_indptr[iele+1] = vv_stiff_ele_indptr[iele] + length(el_eqns_active_idx)*length(el_eqns_active_idx)

      for j = 1:ngps_per_elem*nstrain
        for i = 1:length(el_eqns_active_idx)
          push!(ii_dfint_dstress, el_eqns_active_idx[i])
          push!(jj_dfint_dstress, (iele-1)*ngps_per_elem*nstrain+j)
          #push!(vv_dfint_dstress, dfint_dstress_active[i,j])
        end
      end
      vv_dfint_dstress_ele_indptr[iele+1] = vv_dfint_dstress_ele_indptr[iele] + ngps_per_elem*nstrain*length(el_eqns_active_idx)

      for j = 1:length(el_eqns_active_idx)
        for i = 1:ngps_per_elem*nstrain

          push!(ii_dstrain_dstate, (iele-1)*ngps_per_elem*nstrain+i)
          push!(jj_dstrain_dstate, el_eqns_active_idx[j])
          #push!(vv_dstrain_dstate, dstrain_dstate_active[i,j])
        end
      end
      vv_dstrain_dstate_ele_indptr[iele+1] = vv_dstrain_dstate_ele_indptr[iele] + ngps_per_elem*nstrain*length(el_eqns_active_idx)

    end

    self.ii_stiff = ii_stiff; self.jj_stiff = jj_stiff;
    self.vv_stiff_ele_indptr = vv_stiff_ele_indptr; self.vv_stiff = similar(ii_stiff)

    self.ii_dfint_dstress = ii_dfint_dstress; self.jj_dfint_dstress = jj_dfint_dstress;
    self.vv_dfint_dstress_ele_indptr = vv_dfint_dstress_ele_indptr; self.vv_dfint_dstress= similar(ii_dfint_dstress)

    self.ii_dstrain_dstate = ii_dstrain_dstate; self.jj_dstrain_dstate = jj_dstrain_dstate;
    self.vv_dstrain_dstate_ele_indptr = vv_dstrain_dstate_ele_indptr; self.vv_dstrain_dstate = similar(ii_dstrain_dstate)

  end


text:
include("../src/Dataclass.jl")
using .DataClass
using Test

@dataclass mutable struct Register
    value::Int64
    valuestr::String
    _initvars = [_init=true, _repr=true, _eq=true, _order=true, _unsafe_hash=false, _frozen=true]
end

@testset "Dataclass.jl" begin
    a = Register(1, "1")
    b = Register(10, "4")
    setfield!(a, :value, 10) # Tests if struct fields can still be edited
    __init__(a, 10, "2")
    @test __eq__(a, Register(10, "2")) # Test init
    @test __repr__(a) == "Register(10, 2)"
    @test __eq__(a,b) == false
    @test __lt__(a, b) == true
    @test __le__(a, b) == true
    @test __gt__(a, b) == false
    @test __ge__(a, b) == false
    @test __key(a) == (10, "2")
    @test __hash__(a) == 2845903754746135373
end


text:
#=
Module for solving the Leitner Queue Network static planning problem
@author Siddharth Reddy <sgr45@cornell.edu>
=#

using JuMP
using Ipopt

eps = 1e-9

function allocate_work_rates(num_decks, work_rate_budget, difficulty)
    m = Model(solver=IpoptSolver()) # Uses IP-OPT solver by default
    @variable(m, u[1:num_decks] >= 0)
    @variable(m, a >= 0)
    @variable(m, l[1:num_decks] >= 0)
    @constraint(m, a + sum(u) <= work_rate_budget)
    for i = 1:num_decks
        @constraint(m, l[i] <= u[i])
    end
    @NLconstraint(m, l[1] == a + (1 - (u[1] - l[1]) / (u[1] - l[1] + difficulty)) * l[1] +
            (1 - (u[2] - l[2]) / (u[2] - l[2] + difficulty / 2)) * l[2])
    for i = 2:num_decks-1
        @NLconstraint(m, l[i] == (u[i-1] - l[i-1]) / (u[i-1] - l[i-1] + difficulty / (i-1)) *
                l[i-1] + (1 - (u[i+1] - l[i+1]) / (u[i+1] - l[i+1] + difficulty / (i+1))) * l[i+1])
    end
    @NLconstraint(m, l[num_decks] == (u[num_decks-1] - l[num_decks-1]) /
            (u[num_decks-1] - l[num_decks-1] + difficulty / (num_decks-1)) * l[num_decks-1])
    @NLobjective(m, Max, a)
    TT = STDOUT
    redirect_stdout()
    status = solve(m)
    redirect_stdout(TT)
    return getobjectivevalue(m), getvalue(u), getvalue(l)
end

# Assuming fixed work rates, what is the corresponding equilibrium?
function throughput_for_work_rates(num_decks, u, difficulty)
    m = Model(solver=IpoptSolver())
    @variable(m, a >= 0)
    @variable(m, l[1:num_decks] >= 0)
    for i = 1:num_decks
        @constraint(m, l[i] <= u[i])
    end
    @NLconstraint(m, l[1] == a + (1 - (u[1] - l[1]) / (u[1] - l[1] + difficulty)) * l[1] +
            (1 - (u[2] - l[2]) / (u[2] - l[2] + difficulty / 2)) * l[2])
    for i = 2:num_decks-1
        @NLconstraint(m, l[i] == (u[i-1] - l[i-1]) / (u[i-1] - l[i-1] + difficulty / (i-1)) *
                l[i-1] + (1 - (u[i+1] - l[i+1]) / (u[i+1] - l[i+1] + difficulty / (i+1))) * l[i+1])
    end
    @NLconstraint(m, l[num_decks] == (u[num_decks-1] - l[num_decks-1]) /
            (u[num_decks-1] - l[num_decks-1] + difficulty / (num_decks-1)) * l[num_decks-1])
    @NLobjective(m, Max, a)
    TT = STDOUT
    redirect_stdout()
    status = solve(m)
    redirect_stdout(TT)
    return getobjectivevalue(m)
end

# Assuming fixed work rates and arrival rate, what is the corresponding equilibrium?
function eq_flow_rates_for_work_rates_and_arrival_rate(num_decks, u, difficulty, a)
    m = Model(solver=IpoptSolver())
    @variable(m, l[1:num_decks] >= 0)
    for i = 1:num_decks
        @constraint(m, l[i] <= u[i])
    end
    @NLconstraint(m, l[1] == a + (1 - (u[1] - l[1]) / (u[1] - l[1] + difficulty)) * l[1] +
            (1 - (u[2] - l[2]) / (u[2] - l[2] + difficulty / 2)) * l[2])
    for i = 2:num_decks-1
        @NLconstraint(m, l[i] == (u[i-1] - l[i-1]) / (u[i-1] - l[i-1] + difficulty / (i-1)) *
                l[i-1] + (1 - (u[i+1] - l[i+1]) / (u[i+1] - l[i+1] + difficulty / (i+1))) * l[i+1])
    end
    @NLconstraint(m, l[num_decks] == (u[num_decks-1] - l[num_decks-1]) /
            (u[num_decks-1] - l[num_decks-1] + difficulty / (num_decks-1)) * l[num_decks-1])
    @NLobjective(m, Max, l[num_decks])
    TT = STDOUT
    redirect_stdout()
    status = solve(m)
    redirect_stdout(TT)
    return getobjectivevalue(m), getvalue(l)
end

# Use a heuristic to generate a work rate allocation that satisfies a budget
function work_rates_of_heuristic(num_decks, work_rate_budget, f)
    work_rates = map(f, 1:num_decks)
    return work_rates / sum(work_rates) * work_rate_budget
end

max_stars_repo_name: exKAZUu-Research/SmartMotivator

max_issues_repo_name: exKAZUu-Research/SmartMotivator

max_forks_repo_name: exKAZUu-Research/WillingQuiz


text:
InputFolder = './Images/Yeasts/';
OutputFolder = './Results/Images/Yeasts2/';
Rescale = 2;

@iA = '*.tif';

@fxg_mBlockOtsuThr [iA] > [M];
params.OpenRad = 2;
params.TopHatRad = 50;
params.MinRatio = 5;
params.BlckSize = 64;
params.BlckShft = 8;
params.Lvl = 0.5;
params.Bck = 0.5;
params.AbsBck = 0.1;
/endf

@fxm_mModBinWat [M] > [M];
params.SmallHolesArea = 0;
params.DistStdRad = 7;
params.MinDistLocVar = 9;
/endf

/show iA > M;
/keep M > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER