__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
#
