text:
include("MarketOperator.jl")
include("System.jl")

mutable struct TechnologyPortfolio <: Container
    technologies::Array{Technology,1}
    Decomission::Function
    CheckAge::Function
    CheckProfitability::Function
    FindTechnologyDecomission::Function

    function TechnologyPortfolio(technologies::Array{Technology,1})
        self = new()
        self.technologies = technologies

        #Reach life expectancy?
        self.CheckAge = function(currentYear::Real, technology::Technology)
            retire = false
            if (technology.designProperties.constructionTimeStamp + technology.designProperties.lifeTime) <= currentYear
                retire = true
            end
            return retire
        end

        #Select the corresponding technologies that should be decomissioned
        self.FindTechnologyDecomission = function(currentYear::Real)
            technologyToDecomission = Array{Technology,1}(undef,0)
            for technology in self.technologies
                retire = self.CheckAge(currentYear,technology)
                #Retirement creteria is different in different processes
                if retire == true
                    push!(technologyToDecomission, technology)
                end
            end
            return technologyToDecomission
        end
        return self
    end
end

mutable struct GenCoHistory <: Buffer
    #The generation portfolio evolution history of a generation company
    technologyMixHistory::Array{TechnologyPortfolio,1}
    #The new investments of this generation companies in each year
    yearlyInvestment::Array{Technology,1}
    investmentHistory::Array{TechnologyPortfolio,1}
    investmentYearByType::Dict{Symbol,Array{Real,1}}
    allTechnologiesPossessed::Array{Technology,1}
    decomissionedTechnologies::Array{Technology,1}

    function GenCoHistory()
        self = new()
        self.technologyMixHistory = Array{TechnologyPortfolio,1}(undef,0)
        self.yearlyInvestment = Array{Technology,1}(undef,0)
        self.investmentHistory = Array{TechnologyPortfolio,1}(undef,0)
        self.investmentYearByType = Dict{Symbol,Array{Real,1}}()
        self.allTechnologiesPossessed = Array{Technology,1}(undef,0)
        self.decomissionedTechnologies = Array{Technology,1}(undef,0)
        return self
    end
end

mutable struct GenCo <: ActiveAgent
    name::String
    investmentHorizon::Real
    technologyPortfolio::TechnologyPortfolio
    availableTechnologyPool::Array{Technology,1}
    competitorExpansion::Array{CapacityExpansionEdge,1}
    competitorDivision::Array{CapacityDivisionEdge,1}
    history::GenCoHistory

    #NOTE:Common functions
    Decomission::Function
    CalculateTotalCapacity::Function
    EvaluateTechnology::Function
    InvestmentDecision::Function
    ExecuteInvestmentDecision::Function
    ################################
    #NOTE:for GEP investment algorithm
    BuildExistingCapacity::Function
    FindFutureDecomission::Function
    FindFutureExisting::Function
    PriceProjection::Function
    #################################
    #NOTE:for :EMLab and EMCAS investment algorithm
    CreatVirtualEnvironment::Function
    SummarizeTechnologyList::Function
    TimeSeries::Function
    VirtualAuction::Function
    CompetitorInvestment::Function
    CalculateCompetitorExpansion::Function
    SelectExpansionPath::Function
    SelectDivisionPath::Function
    #################################

    function GenCo(name::String,investmentHorizon::Real,technologyPortfolio::TechnologyPortfolio,availableTechnologyPool::Array{Technology,1};kwargs...)
        self = new()
        self.name = name
        self.investmentHorizon = investmentHorizon
        self.technologyPortfolio = technologyPortfolio
        self.availableTechnologyPool = availableTechnologyPool
        V = Dict(kwargs)
        if haskey(V,:competitorExpansion) && haskey(V,:competitorDivision)
            self.competitorExpansion = V[:competitorExpansion]
            self.competitorDivision = V[:competitorDivision]
        end
        self.history = GenCoHistory()
        push!(self.history.technologyMixHistory, deepcopy(self.technologyPortfolio))
        for technology in self.technologyPortfolio.technologies
            push!(self.history.allTechnologiesPossessed,technology)
        end

        #NOTE: Dismantle technologies of this agent
        self.Decomission = function(currentYear::Real, allPowerFlows::Array{PowerFlow,1})
            technologyToDecomission = self.technologyPortfolio.FindTechnologyDecomission(currentYear)
            powerFlowToDelete = Array{PowerFlow,1}(undef,0)
            for powerFlow in allPowerFlows
                if powerFlow.from in technologyToDecomission
                    push!(powerFlowToDelete,powerFlow)
                end
            end
            append!(self.history.decomissionedTechnologies,technologyToDecomission)
            deleteat!(self.technologyPortfolio.technologies, findall(x->x in technologyToDecomission, self.technologyPortfolio.technologies))
            deleteat!(allPowerFlows, findall(x->x in powerFlowToDelete, allPowerFlows))
            return technologyToDecomission
        end

        #NOTE: calculate total system installed capacity
        self.CalculateTotalCapacity = function(allGenCos::Array{GenCo,1})
            totalCapacity = 0.0
            for genCo in allGenCos
                for technology in genCo.technologyPortfolio.technologies
                    totalCapacity += technology.offerQuantity
                end
            end
            return totalCapacity
        end

        #################################################################################################################
        ##############################:NOTE:These functions are for GEP only #############################################
        #NOTE:initial capacity for input of the GEP
        self.BuildExistingCapacity = function(allGenCos::Array{GenCo,1},technologyPool::Array{Technology,1},year::Real,systemSettings::SystemSettings)
            longestTechnologyLifeTime = 20 #TODO:Note this is temporary for this version;
            historicalCapacitySpan = Int(longestTechnologyLifeTime / systemSettings.milestoneYear)
            initialCapacityMatrix = zeros(length(technologyPool) , historicalCapacitySpan) #-15,-10,-5,0 (where 0 refers to the real investment in the current milestone year)
            technologyType = [:Base,:Mid,:Peak]
            for genCo in allGenCos
                for technology in genCo.technologyPortfolio.technologies
                    for techCount in 1:length(technologyPool)
                        if technology.type == technologyType[techCount]
                            initialCapacityMatrix[techCount, historicalCapacitySpan - Int((year - technology.designProperties.constructionTimeStamp)/systemSettings.milestoneYear)] += technology.designProperties.installedCapacity
                        end
                    end
                end
            end
            return initialCapacityMatrix
        end

        #NOTE:return a vector storing the expected decomissioning in future years
        self.FindFutureDecomission = function(allGenCos::Array{GenCo,1},year::Real,systemSettings::SystemSettings)
            decomissionHorizon = Int(self.investmentHorizon / systemSettings.milestoneYear)
            futureDecomission = zeros(1,decomissionHorizon)
            newAllGenCos = deepcopy(allGenCos)
            firstDecomissionTime = year
            finalDecomissionTime = year + self.investmentHorizon - systemSettings.milestoneYear
            for thisYear in firstDecomissionTime:systemSettings.milestoneYear:finalDecomissionTime
                for genCo in newAllGenCos
                    technologyToDecomission = genCo.technologyPortfolio.FindTechnologyDecomission(thisYear, :virtual)
                    for technology in technologyToDecomission
                        futureDecomission[Int((thisYear-year)/systemSettings.milestoneYear)+1] += technology.designProperties.installedCapacity
                    end
                    deleteat!(genCo.technologyPortfolio.technologies, findall(x->x in technologyToDecomission, genCo.technologyPortfolio.technologies))
                end
            end
            return futureDecomission
        end

        #NOTE:return a vector storing the existing generators in future years (including those new investments)
        self.FindFutureExisting = function(allGenCos::Array{GenCo,1},year::Real,systemSettings::SystemSettings)
            futureHorizon = Int(self.investmentHorizon / systemSettings.milestoneYear)
            futureExisting = zeros(1,futureHorizon)
            newAllGenCos = deepcopy(allGenCos)
            firstExistingTime = year
            finalExistingTime = year + self.investmentHorizon -  systemSettings.milestoneYear
            for thisYear in firstExistingTime:systemSettings.milestoneYear:finalExistingTime
                for genCo in newAllGenCos
                    for technology in genCo.technologyPortfolio.technologies
                        if (technology.designProperties.constructionTimeStamp + technology.designProperties.lifeTime) > thisYear
                            futureExisting[Int((thisYear-year)/systemSettings.milestoneYear)+1] += technology.designProperties.installedCapacity
                        end
                    end
                end
            end
            return futureExisting
        end

        #NOTE: GEP investment algorithm. Run GEP to get price projection
        self.PriceProjection = function(technologyPool::Array{Technology,1},allGenCos::Array{GenCo,1},systemLoad::Load,
                                        year::Real,systemSettings::SystemSettings)
            initialCapacityMatrix = self.BuildExistingCapacity(allGenCos,technologyPool,year,systemSettings)

            variable_cost = [technology.economicProperties.VOMCost + technology.economicProperties.fuelCost for technology in technologyPool]
            fixed_cost = [(technology.economicProperties.FOMCost + technology.economicProperties.yearlyInstallment)/technology.designProperties.unitCapacity for technology in technologyPool]
            representative_day = systemSettings.simulationSettings.year2day
            representative_hour = systemSettings.simulationSettings.day2hour
            year2hour = year2day * day2hour
            VoLL = priceCap
            weight_year = systemSettings.milestoneYear
            allTechnologyLifeTime = [technology.designProperties.lifeTime for technology in technologyPool]
            scaledAllTechnologyLifeTime = Int.(allTechnologyLifeTime/systemSettings.milestoneYear)
            longestTechnologyLifeTime = 20 #TODO:Note this is temporary for this version;
                                            #NOTE: the pre-opt year information has to cover what can still exist until the end of the investment horizon
            preOpt_Year = Int(longestTechnologyLifeTime/systemSettings.milestoneYear - 1) #-15,-10,5
            optimizeYear = Int(self.investmentHorizon/systemSettings.milestoneYear)
            #gep= Model(with_optimizer(Gurobi.Optimizer,GUROBI_ENV))
            gep=JuMP.direct_model(Gurobi.Optimizer(Method=1, OutputFlag = 0,Presolve=0))
            #gep=JuMP.direct_model(Gurobi.Optimizer(Presolve=0))

            start_hour = preOpt_Year * representative_day * representative_hour + 1 #Hour 1
            initial_hour = preOpt_Year * representative_day * representative_hour #Hour 0
            stop_hour = initial_hour + optimizeYear * representative_day * representative_hour #Final hour

            #NOTE:Claim variables from past(-preOpt_Year) to future (+optimizeYear)
            @variable(gep,g_j_t[1:length(technologyPool),1:stop_hour]) #generation
            @variable(gep,x_j_t[1:length(technologyPool),1:(optimizeYear+preOpt_Year)]) #Investment
            @variable(gep,y_j_t[1:length(technologyPool),1:(optimizeYear+preOpt_Year)]) #Online capacity
            @variable(gep,z_j_t[1:length(technologyPool),1:(optimizeYear+preOpt_Year)]) #decomission
            @variable(gep,LL_t[1:stop_hour]) #Load loss

            @variable(gep,g_j_t_y[1:length(technologyPool),1:stop_hour,1:4])
            @variable(gep,y_j_t_y[1:length(technologyPool),1:optimizeYear,1:4]) #Online capacity

            #The weights of representative days, and repeat for the whole optimization horizon
            representativeDayWeightFactors = systemSettings.simulationSettings.representativeData[:Weights][:,2]
            dailyWeight = repeat(representativeDayWeightFactors,optimizeYear + preOpt_Year)
            hourlyWeight_vector = repeat(dailyWeight, inner = (day2hour,1))

            #Define objective function
            #obj = 0
            obj = AffExpr(0.0)

            for year in (preOpt_Year+1):(optimizeYear+preOpt_Year)
                for j in 1:length(technologyPool)
                    @constraint(gep, x_j_t[j,year] >= 0)
                    for yy in 1:4
                        @constraint(gep,y_j_t_y[j,year-preOpt_Year,yy] == x_j_t[j, year-yy+1])
                        add_to_expression!(obj, fixed_cost[j]* weight_year,  y_j_t_y[j,year-preOpt_Year,yy])
                    end
                end
            end

            for j = 1:length(technologyPool)
                @constraint(gep, x_j_t[j,4] >= initialCapacityMatrix[j,4])
            end

            #Assign value to previous investments, so that future decommissioning information can be derived
            for year in 1 : preOpt_Year
                for j = 1:length(technologyPool)
                    @constraint(gep, x_j_t[j,year] == initialCapacityMatrix[j,year])
                    @constraint(gep, z_j_t[j,year] == 0)
                end
            end
            for j = 1:length(technologyPool)
                @constraint(gep, z_j_t[j,preOpt_Year+1] == 0)
            end

            #Objective function: variable costs + fixed costs
            for hour in start_hour:stop_hour
                thisYear = Int(ceil(hour / year2hour))
                for j in 1:length(technologyPool)
                    for yy in 1:4
                        add_to_expression!(obj, variable_cost[j] * hourlyWeight_vector[hour] * weight_year,  g_j_t_y[j,hour,yy] )
                    end
                end
                add_to_expression!(obj, VoLL * hourlyWeight_vector[hour] * weight_year,  LL_t[hour])

            end
            @objective(gep,Min,obj)
            constraintDict = Dict{String, Array}()

            for internalYear = 1:(preOpt_Year+optimizeYear)
                for day = 1:year2day
                    for hour = 1:day2hour
                        thisHour = (internalYear-1) * year2hour + (day-1) * day2hour + hour
                        #Load loss in the past, set to zero (optional)
                        if internalYear <= preOpt_Year
                            @constraint(gep, LL_t[thisHour] == 0)
                        end
                    end
                end
            end

            demand = repeat(vec(systemLoad.representativeLoad'),optimizeYear + preOpt_Year) #load for all hours
            demandSupplyConstraints = Array{Any,1}(undef,0)
            #println(size(demand))
            for hour in start_hour:stop_hour
                thisYear = Int(ceil(hour / year2hour))
                thisDemandSupply = @constraint(gep, sum(g_j_t_y[:,hour,:]) + LL_t[hour] == demand[hour])
                push!(demandSupplyConstraints,thisDemandSupply)

                @constraint(gep, LL_t[hour] >= 0)

                for j in 1:length(technologyPool)
                    for yy in 1:4
                        @constraint(gep, g_j_t_y[j,hour,yy] <= y_j_t_y[j,thisYear-preOpt_Year,yy])
                        @constraint(gep, g_j_t_y[j,hour,yy] >= 0)
                    end
                end
            end
            status = optimize!(gep)

            futureCapacityMix = zeros(length(technologyPool),optimizeYear)

            for year in (preOpt_Year+1) : (optimizeYear+preOpt_Year)
                for j in 1:length(technologyPool)
                    for yy in 1:4
                        futureCapacityMix[j,year-preOpt_Year] += value(gep[:y_j_t_y][j,year-preOpt_Year,yy])
                    end
                end
            end
            price = (JuMP.dual.(demandSupplyConstraints))./hourlyWeight_vector[start_hour:stop_hour]/weight_year

            # priceSpikes = price[findall(x -> x>=50, price)]
            # println("Price spikes: ",priceSpikes)
            return price,futureCapacityMix
        end

        #################################################################################################################
        ##############################:NOTE:These functions are for EMLab and EMCAS #############################################
        #Create a virtual play ground w.r.t. the current system capacity mix
        #Preparation for virtual market clearing simulation
        self.CreatVirtualEnvironment = function(allGenCos::Array{GenCo,1},MO::MarketOperator)
            newAllGenCos = deepcopy(allGenCos)
            newAllPowerFlows = Array{PowerFlow,1}(undef,0)
            for genCo in newAllGenCos
                for technology in genCo.technologyPortfolio.technologies
                    newPowerFlow  = PowerFlow(technology, MO, 0,0)
                    push!(newAllPowerFlows, newPowerFlow)
                end
            end
            return newAllGenCos, newAllPowerFlows
        end

        #build a list of all generators
        self.SummarizeTechnologyList = function(allGenCos::Array{GenCo,1},allPowerFlows::Array{PowerFlow,1},MO::MarketOperator,newTechnologyUnderEvaluation::Technology,competitorInvestment::Array{Technology,1})

            allTechnologiesVirtual = Array{Technology,1}(undef,0)
            for genCo in allGenCos
                for technology in genCo.technologyPortfolio.technologies
                    push!(allTechnologiesVirtual, technology)
                end
            end
            push!(allTechnologiesVirtual,newTechnologyUnderEvaluation)
            for technology in competitorInvestment
                push!(allTechnologiesVirtual,technology)
                push!(allPowerFlows,PowerFlow(technology,MO,0,0))
            end
            return allTechnologiesVirtual
        end

        #for renewables only
        self.TimeSeries = function(day::Int64,hour::Int64, load::Load, technologyList::Array{Technology,1})
            load.TimeSeries(day,hour)
            #for renewables only, not considered in this study
            for technology in technologyList
                if typeof(technology) == RenewablePlants
                    technology.TimeSeries(day,hour)
                end
            end
        end

        #Calculate maximum supply gap
        self.CalculateCompetitorExpansion = function(allGenCos::Array{GenCo,1},systemLoad::Load,capacityExpansionEdge::CapacityExpansionEdge)
            totalSupply = self.CalculateTotalCapacity(allGenCos)

            alreadyInvested_self = 0
            for technology in self.history.yearlyInvestment
                alreadyInvested_self += technology.designProperties.installedCapacity
            end

            maxLoad = findmax(systemLoad.representativeLoad)[1]
            totalExpansion = max(capacityExpansionEdge.gapCoverPercentage * maxLoad - totalSupply + alreadyInvested_self,0)
            #println(totalExpansion)
            return totalExpansion
        end

        #how much of the gap is filled by competitors
        self.SelectExpansionPath = function()
            return self.competitorExpansion[1] #This can be changed by the user if desired
        end

        #types of technologies invested by competitors (as expected by the agent)
        self.SelectDivisionPath = function()
            return self.competitorDivision[1] #This can be changed by the user if desired
        end

        # competitors' investment (if any)
        self.CompetitorInvestment = function(allGenCos::Array{GenCo,1}, thisYear::Real,technologyPool::Array{Technology,1},systemLoad::Load,systemSettings::SystemSettings)
            competitorInvestment =Array{Technology,1}(undef,0)
            if systemSettings.simulationSettings.investmentModelSelection == :EMCAS
                capacityExpansionEdge = self.SelectExpansionPath()
                capacityDivisionEdge = self.SelectDivisionPath()
                totalExpansion = self.CalculateCompetitorExpansion(allGenCos,systemLoad,capacityExpansionEdge)

                for typeCount in 1:length(technologyPool)
                    thisType = technologyPool[typeCount].type
                    thisProportion = get(capacityDivisionEdge.typeShare, thisType,0)
                    newTechnology = deepcopy(technologyPool[typeCount])
                    thisUnitNumber = totalExpansion * thisProportion / newTechnology.designProperties.unitCapacity
                    newTechnology.MultipleUnits(thisUnitNumber,systemSettings)
                    newTechnology.designProperties.constructionTimeStamp = deepcopy(thisYear)
                    push!(competitorInvestment,newTechnology)
                end
            end
            return competitorInvestment
        end

        #virtual market clearing
        self.VirtualAuction = function(allGenCos::Array{GenCo,1},technologyPool::Array{Technology,1},newTechnologyUnderEvaluation::Technology,MO::MarketOperator,systemLoad::Load,year::Real,systemSettings::SystemSettings)
            newAllGenCos, newAllPowerFlows = self.CreatVirtualEnvironment(allGenCos,MO)
            push!(newAllPowerFlows,PowerFlow(newTechnologyUnderEvaluation,MO,0,0))
            newSystemLoad = deepcopy(systemLoad)
            year2day = systemSettings.simulationSettings.year2day
            day2hour = systemSettings.simulationSettings.day2hour
            milestoneYear = systemSettings.milestoneYear
            futureCapacityMix = zeros(length(technologyPool),Int(self.investmentHorizon/milestoneYear))
            for virtualYear in year:(year+self.investmentHorizon-1)
                if virtualYear % milestoneYear == 0
                    for genCo in newAllGenCos
                        #Decomission
                        genCo.Decomission(virtualYear,newAllPowerFlows)
                    end
                    competitorInvestment = self.CompetitorInvestment(newAllGenCos, virtualYear,technologyPool,newSystemLoad,systemSettings)

                    allTechnologiesVirtual = self.SummarizeTechnologyList(newAllGenCos,newAllPowerFlows,MO,newTechnologyUnderEvaluation,competitorInvestment)

                    for technology in allTechnologiesVirtual
                        for techCount in 1:length(technologyPool)
                            if technology.type == technologyPool[techCount].type
                                futureCapacityMix[techCount,Int((virtualYear-year)/milestoneYear)+1] += technology.designProperties.installedCapacity
                            end
                        end
                    end

                    for day in 1:year2day, hour in 1:day2hour
                        self.TimeSeries(day,hour,newSystemLoad,allTechnologiesVirtual)
                        marketPrice = MO.ClearMarket(allTechnologiesVirtual,newAllPowerFlows, newSystemLoad)
                        MO.Payment(newAllPowerFlows, day,hour, :virtual,systemSettings)
                    end
                end
            end
            return futureCapacityMix
        end
        #################################################################################################################
        #################################################################################################################

        self.EvaluateTechnology = function(allGenCos::Array{GenCo,1},technologyPool::Array{Technology,1},MO::MarketOperator,newTechnologyUnderEvaluation::Technology,systemLoad::Load,year::Real,systemSettings::SystemSettings)
            revenueComponents = Dict(:PriceSpikes => 0.0, :NormalMarket => 0.0)
            costComponents = Dict(:VariableCost => 0.0, :FixedCost => 0.0)
            futureCapacityMix = zeros(length(technologyPool),Int(self.investmentHorizon/milestoneYear))

            if systemSettings.simulationSettings.investmentModelSelection == :GEP
                priceProjection,futureCapacityMix = self.PriceProjection(technologyPool,allGenCos,systemLoad,year,systemSettings)
                Cost = zeros(newTechnologyUnderEvaluation.designProperties.lifeTime)
                Revenue = zeros(newTechnologyUnderEvaluation.designProperties.lifeTime)
                representativeDayWeightFactors = systemSettings.simulationSettings.representativeData[:Weights][:,2]
                representative_day = systemSettings.simulationSettings.year2day
                representative_hour = systemSettings.simulationSettings.day2hour
                dailyWeight = repeat(representativeDayWeightFactors,Int(newTechnologyUnderEvaluation.designProperties.lifeTime/systemSettings.milestoneYear))

                variable_cost = newTechnologyUnderEvaluation.economicProperties.VOMCost + newTechnologyUnderEvaluation.economicProperties.fuelCost
                fixed_cost = (newTechnologyUnderEvaluation.economicProperties.FOMCost + newTechnologyUnderEvaluation.economicProperties.yearlyInstallment)/newTechnologyUnderEvaluation.designProperties.unitCapacity
                hourlyFixedCost = fixed_cost/(sum(representativeDayWeightFactors)*representative_hour)

                for hour in 1:length(priceProjection)
                    thisWeightDay = dailyWeight[Int(ceil(hour/representative_hour))]
                    thisYearStart = (Int(ceil(hour/representative_day/representative_hour))-1) * systemSettings.milestoneYear + 1
                    thisYearEnd = Int(ceil(hour/representative_day/representative_hour))*systemSettings.milestoneYear

                    Cost[thisYearStart:thisYearEnd] .+= repeat([hourlyFixedCost * thisWeightDay],systemSettings.milestoneYear)
                    costComponents[:FixedCost] += hourlyFixedCost * thisWeightDay * systemSettings.milestoneYear
                    thisPrice = priceProjection[hour]
                    if variable_cost < thisPrice
                        Cost[thisYearStart:thisYearEnd] .+= repeat([variable_cost *  thisWeightDay],systemSettings.milestoneYear)
                        costComponents[:VariableCost] += variable_cost *  thisWeightDay * systemSettings.milestoneYear
                        #Revenue[thisYearStart:thisYearEnd] .+= repeat([thisPrice * markupPercent[Int(thisDayHour)] *  thisWeightDay],systemSettings.milestoneYear)
                        Revenue[thisYearStart:thisYearEnd] .+= repeat([thisPrice *  thisWeightDay],systemSettings.milestoneYear)
                        if thisPrice >= 500
                            revenueComponents[:PriceSpikes] += thisPrice *  thisWeightDay * systemSettings.milestoneYear
                        else
                            revenueComponents[:NormalMarket] += thisPrice *  thisWeightDay * systemSettings.milestoneYear
                        end
                    end
                end

                yearVector = collect(0:1:(newTechnologyUnderEvaluation.designProperties.lifeTime-1))
                lifeDiscountVector = (1-systemSettings.interestRate) .^ yearVector
                NPV = (Revenue - Cost)' * lifeDiscountVector

            else
                newTechnologyUnderEvaluation.designProperties.constructionTimeStamp = deepcopy(year)
                futureCapacityMix = self.VirtualAuction(allGenCos,technologyPool,newTechnologyUnderEvaluation,MO,systemLoad,year,systemSettings)
                newTechnologyUnderEvaluation.AggregateHistory!(newTechnologyUnderEvaluation.economicProperties.virtualHistory, systemSettings)
                NPV = newTechnologyUnderEvaluation.economicProperties.virtualHistory.cashFlow[8]
                newTechnologyUnderEvaluation.economicProperties.virtualHistory.Reinitialize!()
            end
            return NPV,revenueComponents,costComponents,futureCapacityMix
        end

        #NOTE: Make investment decision based filtered technologies
        self.InvestmentDecision = function(allGenCos::Array{GenCo,1},MO::MarketOperator,technologyPool::Array{Technology,1},systemLoad::Load,year::Real,systemSettings::SystemSettings,systemHistory::SystemHistory)

            investmentIndicator = false
            bestTechnology = Nothing()

            #title_1 = string("Electricity price evolution as seen by \n", self.name," in year ", year)

            NPV_Dict = Dict(:Base => 0.0,:Mid => 0.0,:Peak => 0.0)
            revenueComponents_Dict = Dict(:Base => Dict{Symbol,Float64}(),:Mid => Dict{Symbol,Float64}(),:Peak => Dict{Symbol,Float64}())
            costComponents_Dict = Dict(:Base => Dict{Symbol,Float64}(),:Mid => Dict{Symbol,Float64}(),:Peak => Dict{Symbol,Float64}())
            bestNPV = -Inf
            #In each portfolio, every technology has a maximum number of units that can be invested, this loop is used to evaluate all these potential combination of numbers
            for technology in technologyPool
                newTechnologyUnderEvaluation = deepcopy(technology)
                thisNPV,revenueComponents,costComponents,futureCapacityMix = self.EvaluateTechnology(allGenCos,technologyPool,MO,newTechnologyUnderEvaluation,systemLoad,year,systemSettings)
                push!(systemHistory.capacityMixProjection,futureCapacityMix)
                NPV_Dict[technology.type] = thisNPV
                revenueComponents_Dict[technology.type] = revenueComponents
                costComponents_Dict[technology.type] = costComponents
                #println(newTechnologyUnderEvaluation.type, " NPV: ",thisNPV)
                if thisNPV >= -1e-5 && thisNPV > bestNPV
                    bestNPV = deepcopy(thisNPV)
                    #Store the best one and pop up the previous best one
                    bestTechnology = deepcopy(newTechnologyUnderEvaluation)
                    #println("bestNPV: ",thisNPV)
                    #println(unitNumber)
                    #println(bestTechnology.type, " ", bestTechnology.designProperties.installedCapacity)

                end
            end
            if typeof(bestTechnology) != Nothing
                investmentIndicator = true
                #println(bestTechnology.type, " ", bestTechnology.designProperties.installedCapacity)
            end

            #return NPV_Dict,revenueComponents_Dict,costComponents_Dict,investmentIndicator,bestTechnology,priceProjection,futureCapacityMix,capacityMix_Hour,cum_futureCapacityInvestment
            return investmentIndicator,bestTechnology

        end

        #NOTE: plug the new investment to the environment
        self.ExecuteInvestmentDecision = function(bestTechnology::Technology,allPowerFlows::Array{PowerFlow,1},MO::MarketOperator,year::Real)
            push!(self.technologyPortfolio.technologies, bestTechnology)
            bestTechnology.designProperties.constructionTimeStamp = deepcopy(year)
            newPowerFlow = PowerFlow(bestTechnology, MO, 0,0)
            push!(allPowerFlows,newPowerFlow)
            push!(self.history.yearlyInvestment, bestTechnology)
            push!(self.history.allTechnologiesPossessed,bestTechnology)
        end

        return self
    end
end

mutable struct GenCoPortfolioOwnership <: Ownership
    from::Agent
    to::Container
end


text:
###
# Reflection tooling
##

using Base.Meta
using .Compiler: widenconst, argextype, Const

if VERSION >= v"1.1.0-DEV.157"
    const is_return_type = Core.Compiler.is_return_type
else
    is_return_type(f) = f === Core.Compiler.return_type
end

if VERSION >= v"1.2.0-DEV.249"
    const sptypes_from_meth_instance = Core.Compiler.sptypes_from_meth_instance
else
    sptypes_from_meth_instance(mi) = Core.Compiler.spvals_from_meth_instance(mi)
end

if VERSION >= v"1.2.0-DEV.320"
    const may_invoke_generator = Base.may_invoke_generator
else
    may_invoke_generator(meth, @nospecialize(atypes), sparams) = isdispatchtuple(atypes)
end

if VERSION < v"1.2.0-DEV.573"
    code_for_method(method, metharg, methsp, world, force=false) = Core.Compiler.code_for_method(method, metharg, methsp, world, force)
else
    code_for_method(method, metharg, methsp, world, force=false) = Core.Compiler.specialize_method(method, metharg, methsp, force)
end

using Requires
transform(::Val, callsite) = callsite
@init @require CUDAnative="be33ccc6-a3ff-5ff2-a52e-74243cff1e17" begin
    using .CUDAnative
    function transform(::Val{:CuFunction}, callsite, callexpr, CI, mi, slottypes; params=nothing, kwargs...)
        sptypes = sptypes_from_meth_instance(mi)
        tt = argextype(callexpr.args[4], CI, sptypes, slottypes)
        ft = argextype(callexpr.args[3], CI, sptypes, slottypes)
        isa(tt, Const) || return callsite
        return Callsite(callsite.id, CuCallInfo(callinfo(Tuple{widenconst(ft), tt.val.parameters...}, Nothing, params=params)))
    end
end

function find_callsites(CI, mi, slottypes; params=current_params(), kwargs...)
    sptypes = sptypes_from_meth_instance(mi)
    callsites = Callsite[]

    function process_return_type(id, c, rt)
        callinfo = nothing
        is_call = isexpr(c, :call)
        arg_base = is_call ? 0 : 1
        length(c.args) == (arg_base + 3) || return nothing
        ft = argextype(c.args[arg_base + 2], CI, sptypes, slottypes)
        argTs = argextype(c.args[arg_base + 3], CI, sptypes, slottypes)
        if isa(argTs, Const)
            sig = Tuple{widenconst(ft), argTs.val.parameters...}
            mi = first_method_instance(sig)
            if mi !== nothing
                callinfo = MICallInfo(mi, rt.val)
            else
                callinfo = FailedCallInfo(sig, rt)
            end
        else
            callinfo = FailedCallInfo((ft, argTs), rt)
        end
        return Callsite(id, ReturnTypeCallInfo(callinfo))
    end

    for (id, c) in enumerate(CI.code)
        if c isa Expr
            callsite = nothing
            if c.head === :(=)
                c = c.args[2]
                (c isa Expr) || continue
            end
            if c.head === :invoke
                rt = CI.ssavaluetypes[id]
                at = argextype(c.args[2], CI, sptypes, slottypes)
                if isa(at, Const) && is_return_type(at.val)
                    callsite = process_return_type(id, c, rt)
                else
                    callsite = Callsite(id, MICallInfo(c.args[1], rt))
                end
                mi = get_mi(callsite)
                if nameof(mi.def.module) == :CUDAnative && mi.def.name == :cufunction
                    callsite = transform(Val(:CuFunction), callsite, c, CI, mi, slottypes; params=params, kwargs...)
                end
            elseif c.head === :call
                rt = CI.ssavaluetypes[id]
                types = map(arg -> widenconst(argextype(arg, CI, sptypes, slottypes)), c.args)

                # Look through _apply
                ok = true
                while types[1] === typeof(Core._apply)
                    new_types = Any[types[2]]
                    for t in types[3:end]
                        if !(t <: Tuple) || t isa Union
                            ok = false
                            break
                        end
                        append!(new_types, t.parameters)
                    end
                    ok || break
                    types = new_types
                end
                ok || continue

                # Filter out builtin functions and intrinsic function
                if types[1] <: Core.Builtin || types[1] <: Core.IntrinsicFunction
                    continue
                end

                if isdefined(types[1], :instance) && is_return_type(types[1].instance)
                    callsite = process_return_type(id, c, rt)
                else
                    if types[1] isa Union
                        # Union{typeof(sin), typeof(cos)}
                        fts = Any[]
                        function thatcher(u)
                            if u isa Union
                                thatcher(u.a)
                                thatcher(u.b)
                            else
                                push!(fts, u)
                            end
                        end
                        sigs = map(ft-> [ft, types[2:end]], fts)
                        cis = map(types -> callinfo(Tuple{types...}, rt, params=params), sigs)
                        callsite = Callsite(id, MultiCallInfo(Tuple{types...}, rt, cis))
                    else
                        callsite = Callsite(id, callinfo(Tuple{types...}, rt, params=params))
                    end
                end
            else c.head === :foreigncall
                # special handling of jl_new_task
                length(c.args) > 0 || continue
                if c.args[1] isa QuoteNode
                    cfunc = c.args[1].value
                    if cfunc === :jl_new_task
                        func = c.args[7]
                        ftype = widenconst(argextype(func, CI, sptypes, slottypes))
                        callsite = Callsite(id, TaskCallInfo(callinfo(ftype, nothing, params=params)))
                    end
                end
            end

            if callsite !== nothing
                push!(callsites, callsite)
            end
        end
    end
    return callsites
end

if VERSION >= v"1.1.0-DEV.215"
function dce!(ci, mi)
    argtypes = Core.Compiler.matching_cache_argtypes(mi, nothing)[1]
    ir = Compiler.inflate_ir(ci, sptypes_from_meth_instance(mi),
                             argtypes)
    compact = Core.Compiler.IncrementalCompact(ir, true)
    # Just run through the iterator without any processing
    Core.Compiler.foreach(x -> nothing, compact)
    ir = Core.Compiler.finish(compact)
    Core.Compiler.replace_code_newstyle!(ci, ir, length(argtypes)-1)
end
else
function dce!(ci, mi)
end
end

function do_typeinf_slottypes(mi::Core.Compiler.MethodInstance, run_optimizer::Bool, params::Core.Compiler.Params)
    ccall(:jl_typeinf_begin, Cvoid, ())
    result = Core.Compiler.InferenceResult(mi)
    frame = Core.Compiler.InferenceState(result, false, params)
    frame === nothing && return (nothing, Any)
    if Compiler.typeinf(frame) && run_optimizer
        opt = Compiler.OptimizationState(frame)
        Compiler.optimize(opt, result.result)
        opt.src.inferred = true
    end
    ccall(:jl_typeinf_end, Cvoid, ())
    frame.inferred || return (nothing, Any)
    return (frame.src, result.result, frame.slottypes)
end

function preprocess_ci!(ci, mi, optimize)
    if optimize
        # if the optimizer hasn't run, the IR hasn't been converted
        # to SSA form yet and dce is not legal
        dce!(ci, mi)
        dce!(ci, mi)
    end
end

if :trace_inference_limits in fieldnames(Core.Compiler.Params)
    current_params() = Core.Compiler.CustomParams(ccall(:jl_get_world_counter, UInt, ()); trace_inference_limits=true)
else
    current_params() = Core.Compiler.Params(ccall(:jl_get_world_counter, UInt, ()))
end

function callinfo(sig, rt; params=current_params())
    methds = Base._methods_by_ftype(sig, -1, params.world)
    (methds === false || length(methds) < 1) && return FailedCallInfo(sig, rt)
    callinfos = CallInfo[]
    for x in methds
        meth = x[3]
        atypes = x[1]
        sparams = x[2]
        if isdefined(meth, :generator) && !may_invoke_generator(meth, atypes, sparams)
            push!(callinfos, GeneratedCallInfo(sig, rt))
        else
            mi = code_for_method(meth, atypes, sparams, params.world)
            if mi !== nothing
                push!(callinfos, MICallInfo(mi, rt))
            else
                push!(callinfos, FailedCallInfo(sig, rt))
            end
        end
    end

    @assert length(callinfos) != 0
    length(callinfos) == 1 && return first(callinfos)
    return MultiCallInfo(sig, rt, callinfos)
end

function first_method_instance(@nospecialize(F), @nospecialize(TT); params=current_params())
    sig = Base.signature_type(F, TT)
    first_method_instance(sig; params=params)
end

function first_method_instance(@nospecialize(sig); params=current_params())
    methds = Base._methods_by_ftype(sig, 1, params.world)
    (methds === false || length(methds) < 1) && return nothing
    x = methds[1]
    meth = x[3]
    if isdefined(meth, :generator) && !isdispatchtuple(Tuple{sig.parameters[2:end]...})
        return nothing
    end
    mi = code_for_method(meth, sig, x[2], params.world)
end

ext: jl


text:
using SirenSeq
using Base.Test

module TestNames

testC = 0
nextName = "dummy"
testName(name::AbstractString) = ( global testC += 1 ; global nextName = name ; testC )

handler(r::Test.Success) = println("test $(testC):  success")
handler(r::Test.Failure) = println("test $(testC):  FAILURE!\n    $(nextName)")
handler(r::Test.Error) = println("test $(testC):  ERROR!\n    $(nextName)")

end

using TestNames.testName

Test.with_handler(TestNames.handler) do

## printing tests

testName("print note")
s1 = string(N(7))
s2 = "Note:   ch1,   ofs =  0 + 0//1,  dur = 1//1,   itv =  7,  ocv = 3,  vel = 1.00,  sca = cMaj"

@test s1 == s2

testName("print sequence 1")
s1 = string(S(1,2))
s2 = "Exp:    dur = 2//1
  Note:   ch1,   ofs =  0 + 0//1,  dur = 1//1,   itv =  1,  ocv = 3,  vel = 1.00,  sca = cMaj
  Note:   ch1,   ofs =  1 + 0//1,  dur = 1//1,   itv =  2,  ocv = 3,  vel = 1.00,  sca = cMaj"

@test s1 == s2

testName("print sequence 2")
s1 = string(C(Bank(9),Prog(3)))
s2 = "Exp:    dur = 0//1
  Blank-Select:  ch1,   ofs =  0 + 0//1,  bank = 9
  Program-Select:  ch1,   ofs =  0 + 0//1,  prog = 3"

@test s1 == s2

## iseqauls tests

testName("atom equality 1")
@test isequal( N(1), N(1) )

testName("atom equality 2")
@test isequal( Bank(1113), Bank(1113) )

testName("atom equality 3")
@test isequal( Bank(1), Bank(2) ) == false

testName("atom equality 4")
@test isequal(N(1),N(2)) == false

testName("blank equality 1")
@test isequal( B(1//3), B(1//3) )

testName("blank equality 2")
@test isequal( B(1//3), B(2//3) ) == false

testName("sequence equality 1")
@test isequal( S(1,3), S(1,3) )

testName("sequence equality 2")
@test isequal( S(1,N(3)), S(N(1),3) )

testName("sequence equality 3")
@test isequal( S(1,3), S(1,4) ) == false

testName("sequence equality 4")
@test isequal( S(1,2,S(3,4)), S(1,2,3,4) )

testName("sequence equality 5")
@test isequal( S(1,2,3), S(3,2,1) ) == false

testName("chord equality 1")
@test isequal( C(1,2), C(1,2) )

testName("chord equality 2")
@test isequal( C(1,2), C(1,3) ) == false

testName("chord equality 3")
@test isequal( C(C(3,4),1,2), C(1,2,3,4) )

testName("exression equality 1")
@test isequal( C(S(1,2),S(4,5)), S(C(1,4),C(2,5)) )

testName("exression equality 2")
@test isequal( S(1,2,3,B()), C(S(1,2,3),B(4)) )

testName("accelirate equality 1")
@test isequal( A(0.5,4), A(0.5,N(4)) )

testName("accelirate equality 2")
@test isequal( N(4), A(2.0,A(0.5,4)) )

testName("accelirate equality 3")
@test isequal( A(0.7,S(1,2)), S(A(0.7,1),A(0.7,2)) )

testName("dilate equality 1")
@test isequal( D(1//3,4), D(1//3,N(4)) )

testName("dilate equality 2")
@test isequal( N(4), D(2//1,D(1//2,4)) )

testName("dilate equality 3")
@test isequal( D(3//2,S(1,2)), S(D(3//2,1),D(3//2,2)) )

testName("translate equality 1")
@test isequal( T(-1,4), N(3) )

testName("translate equality 2")
@test isequal( T(3,C(1,4)), C(4,T(-1,8)) )

testName("shift equality 1")
@test isequal( F(1//3,7), F(4//3,F(-1,7)) )

testName("shift equality 2")
@test F(3//2,D(1//5,7)).dur == 1//5

testName("shift equality 3")
@test F(3//2,S(1,2,3)).dur == 3//1

## domain errors

testName("dilate domain 1")
@test_throws DomainError D(0,N(7))

testName("dilate domain 2")
@test_throws DomainError D(0,C(7,6))

testName("dilate domain 3")
@test_throws DomainError D(-1//2,N(7))

testName("accel domain 1")
@test_throws DomainError A(0.,N(7))

testName("accel domain 2")
@test_throws DomainError A(0.,C(7,6))

testName("accel domain 3")
@test_throws DomainError A(-0.1,N(7))

end


text:
struct CartesianGrid{N} <: AbstractGrid{N}
    position::Vector    # TODO position is bottom-left corner of domain (in contrast to TreeGrid where it is the center)
    nx::NTuple{N, Int}
    periodic::Vector{Bool}

    nr_cells::Int
    nr_faces::Int
end
function CartesianGrid(position::Vector, n; periodic = fill(false, length(position)))
    N = length(position)
    if isa(periodic, Bool) periodic = fill(periodic, N) end
    nx = Tuple(fill(n, N))
    return CartesianGrid{N}(position, nx, periodic, prod(nx), prod(nx .+ 1))
end

# Simply wrap the filter such that dedicated active / at_boundary etc. are called
@inline cells(grid::CartesianGrid{N}; filter::Function = _ -> true) where N = Iterators.filter(filter, Iterators.product(CartesianIndices(size(CellVar, grid)), (grid, ), (CellVar{N, Number}, )))
@inline faces(grid::CartesianGrid{N}; filter::Function = _ -> true) where N = Iterators.filter(filter, Iterators.product(CartesianIndices(size(FaceVar, grid)), (grid, ), (FaceVar{N, Number}, )))

@inline size(type::Type{CellVar}, grid::CartesianGrid) = grid.nx
@inline size(type::Type{FaceVar}, grid::CartesianGrid{N}) where N = ((Tuple(grid.nx .+ 1))..., N)

@inline index((ind, grid, type)::Tuple{CartesianIndex, CartesianGrid, DataType}) = ind
@inline volume((ind, grid, type)::Tuple{CartesianIndex, CartesianGrid, DataType}) = 1. / prod(grid.nx)
@inline cell_volume((ind, grid, type)::Tuple{CartesianIndex, CartesianGrid, DataType}) = 1. / prod(grid.nx)

function area((ind, grid, type)::Tuple{CartesianIndex, CartesianGrid, DataType})

end

@inline direction((ind, grid, type)::Tuple{CartesianIndex, CartesianGrid{N}, DataType}) where N = ind[N+1]

# NB can't dispatch on Tuples parametrized by a type, so we must unpack the Tuple
@inline active(tup::Tuple{CartesianIndex, CartesianGrid, DataType}) = active(tup...)
function active(ind::CartesianIndex, grid::CartesianGrid, type::Type{CellVar{N, T}}) where {N, T}

end
function active(ind::CartesianIndex, grid::CartesianGrid, type::Type{FaceVar{N, T}}) where {N, T}

end

@inline centroid(tup::Tuple{CartesianIndex, CartesianGrid, DataType}) = centroid(tup...)
function centroid(ind::CartesianIndex, grid::CartesianGrid{N}, type::Type{CellVar{N, T}}) where {N, T}
    cent = copy(grid.position)
    for dim ∈ eachindex(cent)
        cent[dim] += (ind[dim] - 0.5) / grid.nx[dim]
    end
    return cent
end
function centroid(ind::CartesianIndex, grid::CartesianGrid{N}, type::Type{FaceVar{N, T}}) where {N, T}
    cent = copy(grid.position)
    direction = ind[N+1]
    for dim ∈ eachindex(cent)
        if dim == direction
            cent[dim] += (ind[dim] - 1.0) / grid.nx[dim]
        else
            cent[dim] += (ind[dim] - 1.5) / grid.nx[dim]
        end
    end
    return cent
end

@inline at_boundary(tup::Tuple{CartesianIndex, CartesianGrid, DataType}) = at_boundary(tup...)
function at_boundary(ind::CartesianIndex, grid::CartesianGrid, type::Type{FaceVar})

end

@inline at_refinement(it::Tuple{CartesianIndex, CartesianGrid, DataType}) = false


text:
using Jchemo

n = 50 ; p = 7
Xtrain = rand(n, p) ; ytrain = rand(n)
m = 3
Xtest = rand(m, p) ; ytest = rand(m)

fm = svmr(Xtrain, ytrain) ;
#fm = svmr(Xtrain, ytrain; cost = .1, gamma = 10) ;
#fm = svmr(Xtrain, ytrain; kern = "kpol", degree = 2) ;
res = Jchemo.predict(fm, Xtest)
res.pred
rmsep(res.pred, ytest)


text:
"""
    LD(A::Diagonal)
    LD(A::HBlikDiag)
    LD(A::DenseMatrix)

Return `log(det(tril(A)))` evaluated in place.
"""
LD(d::Diagonal{T}) where {T<:Number} = sum(log, d.diag)

function LD(d::UniformBlockDiagonal{T}) where {T}
    dat = d.data
    return sum(log, dat[j, j, k] for j in axes(dat, 2), k in axes(dat, 3))
end

LD(d::DenseMatrix{T}) where {T} = sum(log, d[k] for k in diagind(d))

"""
    logdet(m::LinearMixedModel)

Return the value of `log(det(Λ'Z'ZΛ + I)) + m.optsum.REML * log(det(LX*LX'))`
evaluated in place.

Here LX is the diagonal term corresponding to the fixed-effects in the blocked
lower Cholesky factor.
"""
function LinearAlgebra.logdet(m::LinearMixedModel{T}) where {T}
    L = m.L
    nre = m.dims.nretrms
    s = LD(first(L))
    for j in 2:nre
        s += LD(L[kp1choose2(j)])
    end
    if m.optsum.REML
        lastL = last(L)
        s += LD(lastL)  # this adds the log of sqrtpwrss
        s -= log(last(lastL))
    end
    return s + s
end


text:
function CauchyDistribution()
    return CauchyDistribution(())
end

function CauchyDistribution(arg0::jdouble, arg1::jdouble)
    return CauchyDistribution((jdouble, jdouble), arg0, arg1)
end

function cumulative_probability(obj::CauchyDistribution, arg0::jdouble)
    return jcall(obj, "cumulativeProbability", jdouble, (jdouble,), arg0)
end

function density(obj::CauchyDistribution, arg0::jdouble)
    return jcall(obj, "density", jdouble, (jdouble,), arg0)
end

function get_median(obj::CauchyDistribution)
    return jcall(obj, "getMedian", jdouble, ())
end

function get_numerical_mean(obj::CauchyDistribution)
    return jcall(obj, "getNumericalMean", jdouble, ())
end

function get_numerical_variance(obj::CauchyDistribution)
    return jcall(obj, "getNumericalVariance", jdouble, ())
end

function get_scale(obj::CauchyDistribution)
    return jcall(obj, "getScale", jdouble, ())
end

function get_support_lower_bound(obj::CauchyDistribution)
    return jcall(obj, "getSupportLowerBound", jdouble, ())
end

function get_support_upper_bound(obj::CauchyDistribution)
    return jcall(obj, "getSupportUpperBound", jdouble, ())
end

function inverse_cumulative_probability(obj::CauchyDistribution, arg0::jdouble)
    return jcall(obj, "inverseCumulativeProbability", jdouble, (jdouble,), arg0)
end

function is_support_connected(obj::CauchyDistribution)
    return jcall(obj, "isSupportConnected", jboolean, ())
end


text:
#!/usr/bin/env julia

using Eigenbroetler

function sunflower(x::Number, y::Number)
    if x == 0 && y == 0
        return 0.0
    end
    r = sqrt(x * x + y * y)
    ϕ = atan(y, x)
    return cos(30 * log(r) + 13 * ϕ) ^ 2 + cos(20 * log(r) - 11 * ϕ) ^ 2
end

eb = Eigenbrot(sunflower, 512, 512)
save("sunflower.png", eb, ImageSetting(RealPart, LinearScale), colours = :rainbow)


text:
"""
    Plaintext

A plaintext element, storing data as a polynomial modulo the plaintext modulus. It can be used to
create a `Ciphertext` element by encrypting it with an appropriate `Encryptor` instance. Decrypting
a `Ciphertext` with a `Decryptor` instance will again return a `Plaintext` instance.

See also: [`Ciphertext`](@ref), [`Encryptor`](@ref), [`Decryptor`](@ref)
"""
mutable struct Plaintext <: SEALObject
  handle::Ptr{Cvoid}

  function Plaintext()
    handleref = Ref{Ptr{Cvoid}}(C_NULL)
    retval = ccall((:Plaintext_Create1, libsealc), Clong,
                   (Ptr{Cvoid}, Ref{Ptr{Cvoid}}),
                   C_NULL, handleref)
    @check_return_value retval
    return Plaintext(handleref[])
  end

  function Plaintext(capacity, coeff_count)
    handleref = Ref{Ptr{Cvoid}}(C_NULL)
    retval = ccall((:Plaintext_Create3, libsealc), Clong,
                   (UInt64, UInt64, Ptr{Cvoid}, Ref{Ptr{Cvoid}}),
                   capacity, coeff_count, C_NULL, handleref)
    @check_return_value retval
    return Plaintext(handleref[])
  end

  function Plaintext(hex_poly)
    handleref = Ref{Ptr{Cvoid}}(C_NULL)
    retval = ccall((:Plaintext_Create4, libsealc), Clong,
                   (Cstring, Ptr{Cvoid}, Ref{Ptr{Cvoid}}),
                   hex_poly, C_NULL, handleref)
    @check_return_value retval
    return Plaintext(handleref[])
  end

  function Plaintext(handle::Ptr{Cvoid})
    object = new(handle)
    finalizer(destroy!, object)
    return object
  end
end

function destroy!(object::Plaintext)
  if isallocated(object)
    @check_return_value ccall((:Plaintext_Destroy, libsealc), Clong, (Ptr{Cvoid},), object)
    sethandle!(object, C_NULL)
  end

  return nothing
end

function scale(plain::Plaintext)
  value = Ref{Cdouble}(0)
  retval = ccall((:Plaintext_Scale, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{Cdouble}),
                 plain, value)
  @check_return_value retval
  return Float64(value[])
end

function scale!(plain::Plaintext, value::Float64)
  retval = ccall((:Plaintext_SetScale, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{Cdouble}),
                 plain, value)
  @check_return_value retval
  return plain
end

function parms_id(plain::Plaintext)
  parms_id_ = zeros(UInt64, 4)
  retval = ccall((:Plaintext_GetParmsId, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt64}),
                 plain, parms_id_)
  @check_return_value retval
  return parms_id_
end

function to_string(plain::Plaintext)
  len = Ref{UInt64}(0)

  # First call to obtain length (message pointer is null)
  retval = ccall((:Plaintext_ToString, libsealc), Clong,
                 (Ptr{Cvoid}, Ptr{UInt8}, Ref{UInt64}),
                 plain, C_NULL, len)
  @check_return_value retval

  # Second call to obtain message
  # Note: The "+1" is needed since the terminating NULL byte is included in the *copy* operation in
  # SEAL, but *not* in the returned length.
  message = Vector{UInt8}(undef, len[] + 1)
  retval = ccall((:Plaintext_ToString, libsealc), Clong,
                 (Ptr{Cvoid}, Ptr{UInt8}, Ref{UInt64}),
                 plain, message, len)
  @check_return_value retval

  # Return as String but without terminating NULL byte
  return String(message[1:end-1])
end

function save_size(compr_mode, plain::Plaintext)
  result = Ref{Int64}(0)
  retval = ccall((:Plaintext_SaveSize, libsealc), Clong,
                 (Ptr{Cvoid}, UInt8, Ref{Int64}),
                 plain, compr_mode, result)
  @check_return_value retval
  return Int(result[])
end
save_size(plain::Plaintext) = save_size(ComprModeType.default, plain)

function save!(buffer::DenseVector{UInt8}, length::Integer,
               compr_mode::ComprModeType.ComprModeTypeEnum, plain::Plaintext)
  out_bytes = Ref{Int64}(0)
  retval = ccall((:Plaintext_Save, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt8}, UInt64, UInt8, Ref{Int64}),
                 plain, buffer, length, compr_mode, out_bytes)
  @check_return_value retval
  return Int(out_bytes[])
end
function save!(buffer::DenseVector{UInt8}, length::Integer, plain::Plaintext)
  return save!(buffer, length, ComprModeType.default, plain)
end
function save!(buffer::DenseVector{UInt8}, plain::Plaintext)
  return save!(buffer, length(buffer), plain)
end

function Base.:(==)(plain1::Plaintext, plain2::Plaintext)
  result = Ref{UInt8}(0)
  retval = ccall((:Plaintext_Equals, libsealc), Clong,
                 (Ptr{Cvoid}, Ptr{Cvoid}, Ref{UInt8}),
                 plain1, plain2, result)
  @check_return_value retval
  return Bool(result[])
end


text:
using Random; Random.seed!(1234)
using OrdinaryDiffEq
using Statistics
using ForwardDiff, Calculus
using DiffEqSensitivity
using Test
using Zygote

@testset "LSS" begin
  @testset "Lorentz single parameter" begin
    function lorenz!(du,u,p,t)
      du[1] = 10*(u[2]-u[1])
      du[2] = u[1]*(p[1]-u[3]) - u[2]
      du[3] = u[1]*u[2] - (8//3)*u[3]
    end

    p = [28.0]
    tspan_init = (0.0,30.0)
    tspan_attractor = (30.0,50.0)
    u0 = rand(3)
    prob_init = ODEProblem(lorenz!,u0,tspan_init,p)
    sol_init = solve(prob_init,Tsit5())
    prob_attractor = ODEProblem(lorenz!,sol_init[end],tspan_attractor,p)
    sol_attractor = solve(prob_attractor,Vern9(),abstol=1e-14,reltol=1e-14)

    g(u,p,t) = u[end]
    function dg(out,u,p,t,i)
      fill!(out, zero(eltype(u)))
      out[end] = -one(eltype(u))
    end
    lss_problem1 = ForwardLSSProblem(sol_attractor, ForwardLSS(), g)
    lss_problem1a = ForwardLSSProblem(sol_attractor, ForwardLSS(), nothing, dg)
    lss_problem2 = ForwardLSSProblem(sol_attractor, ForwardLSS(alpha=DiffEqSensitivity.Cos2Windowing()), g)
    lss_problem2a = ForwardLSSProblem(sol_attractor, ForwardLSS(alpha=DiffEqSensitivity.Cos2Windowing()), nothing, dg)
    lss_problem3 = ForwardLSSProblem(sol_attractor, ForwardLSS(alpha=10), g)
    lss_problem3a = ForwardLSSProblem(sol_attractor, ForwardLSS(alpha=10), g, dg) #ForwardLSS with time dilation requires knowledge of g

    adjointlss_problem = AdjointLSSProblem(sol_attractor, AdjointLSS(alpha=10.0), g)
    adjointlss_problem_a = AdjointLSSProblem(sol_attractor, AdjointLSS(alpha=10.0), g, dg)

    res1 = shadow_forward(lss_problem1)
    res1a = shadow_forward(lss_problem1a)
    res2 = shadow_forward(lss_problem2)
    res2a = shadow_forward(lss_problem2a)
    res3 = shadow_forward(lss_problem3)
    res3a = shadow_forward(lss_problem3a)

    res4 = shadow_adjoint(adjointlss_problem)
    res4a = shadow_adjoint(adjointlss_problem_a)

    @test res1[1] ≈ 1 atol=5e-2
    @test res2[1] ≈ 1 atol=5e-2
    @test res3[1] ≈ 1 atol=5e-2

    @test res1 ≈ res1a atol=1e-10
    @test res2 ≈ res2a atol=1e-10
    @test res3 ≈ res3a atol=1e-10
    @test res3 ≈ res4 atol=1e-10
    @test res3 ≈ res4a atol=1e-10

    # fixed saveat to compare with concrete solve
    sol_attractor2 = solve(prob_attractor,Vern9(),abstol=1e-14,reltol=1e-14, saveat=0.01)
    lss_problem1 = ForwardLSSProblem(sol_attractor2, ForwardLSS(), g)
    lss_problem1a = ForwardLSSProblem(sol_attractor2, ForwardLSS(), nothing, dg)
    lss_problem2 = ForwardLSSProblem(sol_attractor2, ForwardLSS(alpha=DiffEqSensitivity.Cos2Windowing()), g)
    lss_problem2a = ForwardLSSProblem(sol_attractor2, ForwardLSS(alpha=DiffEqSensitivity.Cos2Windowing()), nothing, dg)
    lss_problem3 = ForwardLSSProblem(sol_attractor2, ForwardLSS(alpha=10), g)
    lss_problem3a = ForwardLSSProblem(sol_attractor2, ForwardLSS(alpha=10), g, dg) #ForwardLSS with time dilation requires knowledge of g

    adjointlss_problem = AdjointLSSProblem(sol_attractor2, AdjointLSS(alpha=10.0), g)
    adjointlss_problem_a = AdjointLSSProblem(sol_attractor2, AdjointLSS(alpha=10.0), g, dg)

    res1 = shadow_forward(lss_problem1)
    res1a = shadow_forward(lss_problem1a)
    res2 = shadow_forward(lss_problem2)
    res2a = shadow_forward(lss_problem2a)
    res3 = shadow_forward(lss_problem3)
    res3a = shadow_forward(lss_problem3a)

    res4 = shadow_adjoint(adjointlss_problem)
    res4a = shadow_adjoint(adjointlss_problem_a)

    @test res1[1] ≈ 1 atol=5e-2
    @test res2[1] ≈ 1 atol=5e-2
    @test res3[1] ≈ 1 atol=5e-2

    @test res1 ≈ res1a atol=1e-10
    @test res2 ≈ res2a atol=1e-10
    @test res3 ≈ res3a atol=1e-10
    @test res3 ≈ res4 atol=1e-10
    @test res3 ≈ res4a atol=1e-10

    function G(p; sensealg=ForwardLSS(), dt=0.01, g=nothing)
      _prob = remake(prob_attractor,p=p)
      _sol = solve(_prob,Vern9(),abstol=1e-14,reltol=1e-14,saveat=dt,sensealg=sensealg, g=g)
      sum(getindex.(_sol.u,3))
    end

    dp1 = Zygote.gradient((p)->G(p),p)
    @test res1 ≈ dp1[1] atol=1e-10

    dp1 = Zygote.gradient((p)->G(p, sensealg=ForwardLSS(alpha=DiffEqSensitivity.Cos2Windowing())),p)
    @test res2 ≈ dp1[1] atol=1e-10

    dp1 = Zygote.gradient((p)->G(p, sensealg=ForwardLSS(alpha=10), g=g),p)
    @test res3 ≈ dp1[1] atol=1e-10

    dp1 = Zygote.gradient((p)->G(p, sensealg=AdjointLSS(alpha=10.0), g=g),p)
    @test res4 ≈ dp1[1] atol=1e-10

    @show res1[1] res2[1] res3[1]
  end

  @testset "Lorentz" begin
    function lorenz!(du,u,p,t)
      du[1] = p[1]*(u[2]-u[1])
      du[2] = u[1]*(p[2]-u[3]) - u[2]
      du[3] = u[1]*u[2] - p[3]*u[3]
    end

    p = [10.0, 28.0, 8/3]

    tspan_init = (0.0,30.0)
    tspan_attractor = (30.0,50.0)
    u0 = rand(3)
    prob_init = ODEProblem(lorenz!,u0,tspan_init,p)
    sol_init = solve(prob_init,Tsit5())
    prob_attractor = ODEProblem(lorenz!,sol_init[end],tspan_attractor,p)
    sol_attractor = solve(prob_attractor,Vern9(),abstol=1e-14,reltol=1e-14)

    g(u,p,t) = u[end] + sum(p)
    function dgu(out,u,p,t,i)
      fill!(out, zero(eltype(u)))
      out[end] = -one(eltype(u))
    end
    function dgp(out,u,p,t,i)
      fill!(out, -one(eltype(p)))
    end

    lss_problem = ForwardLSSProblem(sol_attractor, ForwardLSS(alpha=10), g)
    lss_problem_a = ForwardLSSProblem(sol_attractor, ForwardLSS(alpha=10), g, (dgu,dgp))
    adjointlss_problem = AdjointLSSProblem(sol_attractor, AdjointLSS(alpha=10.0), g)
    adjointlss_problem_a = AdjointLSSProblem(sol_attractor, AdjointLSS(alpha=10.0), g, (dgu,dgp))

    resfw = shadow_forward(lss_problem)
    resfw_a = shadow_forward(lss_problem_a)
    resadj = shadow_adjoint(adjointlss_problem)
    resadj_a = shadow_adjoint(adjointlss_problem_a)

    @test resfw ≈ resadj rtol=1e-10
    @test resfw ≈ resfw_a rtol=1e-10
    @test resfw ≈ resadj_a rtol=1e-10

    sol_attractor2 = solve(prob_attractor,Vern9(),abstol=1e-14,reltol=1e-14, saveat=0.01)
    lss_problem = ForwardLSSProblem(sol_attractor2, ForwardLSS(alpha=10), g)
    resfw = shadow_forward(lss_problem)

    function G(p; sensealg=ForwardLSS(), dt=0.01, g=nothing)
      _prob = remake(prob_attractor,p=p)
      _sol = solve(_prob,Vern9(),abstol=1e-14,reltol=1e-14,saveat=dt,sensealg=sensealg, g=g)
      sum(getindex.(_sol.u,3)) + sum(p)
    end

    dp1 = Zygote.gradient((p)->G(p, sensealg=ForwardLSS(alpha=10), g=g),p)
    @test resfw ≈ dp1[1] atol=1e-10

    dp1 = Zygote.gradient((p)->G(p, sensealg=AdjointLSS(alpha=10.0), g=g),p)
    @test resfw ≈ dp1[1] atol=1e-10

    @show resfw
  end

  @testset "T0skip and T1skip" begin
    function lorenz!(du,u,p,t)
      du[1] = p[1]*(u[2]-u[1])
      du[2] = u[1]*(p[2]-u[3]) - u[2]
      du[3] = u[1]*u[2] - p[3]*u[3]
    end

    p = [10.0, 28.0, 8/3]

    tspan_init = (0.0,30.0)
    tspan_attractor = (30.0,50.0)
    u0 = rand(3)
    prob_init = ODEProblem(lorenz!,u0,tspan_init,p)
    sol_init = solve(prob_init,Tsit5())
    prob_attractor = ODEProblem(lorenz!,sol_init[end],tspan_attractor,p)
    sol_attractor = solve(prob_attractor,Vern9(),abstol=1e-14,reltol=1e-14, saveat=0.01)

    g(u,p,t) = u[end]^2/2 + sum(p)
    function dgu(out,u,p,t,i)
      fill!(out, zero(eltype(u)))
      out[end] = -u[end]
    end
    function dgp(out,u,p,t,i)
      fill!(out, -one(eltype(p)))
    end

    function G(p; sensealg=ForwardLSS(), dt=0.01, g=nothing, t0skip=0.0, t1skip=0.0)
      _prob = remake(prob_attractor,p=p)
      _sol = solve(_prob,Vern9(),abstol=1e-14,reltol=1e-14,saveat=dt,sensealg=sensealg, g=g, t0skip=t0skip, t1skip=t1skip)
      sum(getindex.(_sol.u,3).^2)/2 + sum(p)
    end

    ## ForwardLSS

    lss_problem = ForwardLSSProblem(sol_attractor, ForwardLSS(alpha=10), g)
    resfw = shadow_forward(lss_problem)

    res = deepcopy(resfw)

    dp1 = Zygote.gradient((p)->G(p, sensealg=ForwardLSS(alpha=10), g=g),p)
    @test res ≈ dp1[1] atol=1e-10

    resfw = shadow_forward(lss_problem; t0skip=10.0, t1skip=5.0)
    resskip = deepcopy(resfw)

    dp1 = Zygote.gradient((p)->G(p, sensealg=ForwardLSS(alpha=10), g=g, t0skip=10.0, t1skip=5.0),p)
    @test resskip ≈ dp1[1] atol=1e-10

    @show res resskip

    ## ForwardLSS with dgdu and dgdp

    lss_problem = ForwardLSSProblem(sol_attractor, ForwardLSS(alpha=10), g, (dgu,dgp))
    res2 = shadow_forward(lss_problem)
    @test res ≈ res2 atol=1e-10
    res2 = shadow_forward(lss_problem; t0skip=10.0, t1skip=5.0)
    @test resskip ≈ res2 atol=1e-10

    ## AdjointLSS

    lss_problem = AdjointLSSProblem(sol_attractor, AdjointLSS(alpha=10.0), g)
    res2 = shadow_adjoint(lss_problem)
    @test res ≈ res2 atol=1e-10
    res2 = shadow_adjoint(lss_problem; t0skip=10.0, t1skip=5.0)
    @test_broken resskip ≈ res2 atol=1e-10

    dp1 = Zygote.gradient((p)->G(p, sensealg=AdjointLSS(alpha=10.0), g=g),p)
    @test res ≈ dp1[1] atol=1e-10

    dp1 = Zygote.gradient((p)->G(p, sensealg=AdjointLSS(alpha=10), g=g, t0skip=10.0, t1skip=5.0),p)
    @test res2 ≈ dp1[1] atol=1e-10

    ## AdjointLSS with dgdu and dgd

    lss_problem = AdjointLSSProblem(sol_attractor, AdjointLSS(alpha=10.0), g, (dgu,dgp))
    res2 = shadow_adjoint(lss_problem)
    @test res ≈ res2 atol=1e-10
    res2 = shadow_adjoint(lss_problem; t0skip=10.0, t1skip=5.0)
    @test_broken resskip ≈ res2 atol=1e-10
  end
end

@testset "NILSS" begin
  @testset "Lorentz single parameter" begin
    function lorenz!(du,u,p,t)
      du[1] = 10*(u[2]-u[1])
      du[2] = u[1]*(p[1]-u[3]) - u[2]
      du[3] = u[1]*u[2] - (8//3)*u[3]
    end

    p = [28.0]
    tspan_init = (0.0,100.0)
    tspan_attractor = (100.0,120.0)
    u0 = rand(3)
    prob_init = ODEProblem(lorenz!,u0,tspan_init,p)
    sol_init = solve(prob_init,Tsit5())
    prob_attractor = ODEProblem(lorenz!,sol_init[end],tspan_attractor,p)

    g(u,p,t) = u[end]
    function dg(out,u,p,t,i)
      fill!(out, zero(eltype(u)))
      out[end] = -one(eltype(u))
    end

    nseg = 50 # number of segments on time interval
    nstep = 2001 # number of steps on each segment

    nilss_prob = NILSSProblem(prob_attractor, NILSS(nseg, nstep), g)
    res1 = DiffEqSensitivity.shadow_forward(nilss_prob,Tsit5())

    @test res1[1] ≈ 1 atol=5e-2
  end
end


text:
using LinearAlgebra
using ProgressMeter

has_real_eigvals(H) = ishermitian(H)

"""Initialize a workspace for propagation.

```julia
wrk = initpropwrk(state, tlist, method=:auto, generator...; kwargs...)
```

The resulting `wrk` can be passed to [`propagate`](@ref) or
[`propstep!`](@ref).

# Arguments

* `state`: An exemplary state for the propagation (e.g., the initial state)
* `tlist`: The time grid over which [`propagate`](@ref) will be called. Must
  include at least to points in order to determine the propagation time step
  to prepare. If the propagation will be over a `tlist` with a variable `dt`,
  the full `tlist` must be passed here.
* `generator`: An exemplary (non-time-dependent) dynamical generator. For full
  generality (if `method=:cheby`), the given `generator` should have a spectral
  range sufficiently large to encompass the entire propagation. If given
  multiple times, a spectral envelope enclosing all the generators will be
  determined automatically. In this case, you should pass the generators with
  the extremal values of all the controls.
* `method`: The propagation method to use. The default value of `:auto`
  attempts to choose the best method available, based on the properties of the
  given `state`, `tlist`, and `generator`. Alternative values are `:cheby` and
  `:newton`, and `:expprop`.
* `specrad_method`: for `method=:cheby`, method to use for estimating the
   spectral radius, see [`QuantumPropagators.SpectralRange.specrange`](@ref).
   Defaults to `:auto`.
* `tolerance`: for `method=:cheby`, a tolerance factor for the estimated
  spectral radius. That is, Chebychev coefficients will be calculated for a
  spectral radius increased by the `tolerance` factor compared to the specral
  radius estimated for the `generator`.

All other `kwargs` are filtered and passed to the contructor for returned
workspace, e.g. `limit` for `method=:cheby` or `m_max` for `method=:newton`.
For `method=:cheby`, they are also passed to
[`specrange`](@ref QuantumPropagators.SpectralRange.specrange).
"""
function initpropwrk(
    state,
    tlist,
    method::Val{:auto},
    generator...;
    verbose=false,
    kwargs...
)
    is_small = false
    if length(generator) == 0
        try
            is_small = length(state) <= 8
        catch exception
            if isa(exception, MethodError)
                @warn "Cannot determine problem dimension"
            else
                rethrow()
            end
        end
        method = is_small ? Val(:expprop) : Val(:newton)
    else # we have at least one generator
        try
            is_small = length(state) <= 8
        catch
            try
                is_small = all(size(G)[2] <= 8 for G in generator)
            catch exception
                if isa(exception, MethodError)
                    @warn "Cannot determine problem dimension"
                else
                    rethrow()
                end
            end
        end
        if is_small
            method = Val(:expprop)
        else
            method = Val(:newton)
            try
                if all(has_real_eigvals(G) for G in generator)
                    method = Val(:cheby)
                end
            catch exception
                # if we can't determine G is Hermitian, we don't do anything
                if isa(exception, MethodError)
                    @warn("Cannot determine if generator is Hermitian.")
                else
                    rethrow()
                end
            end
        end
    end
    if verbose
        @info "Auto-choosing propagation method $(repr(typeof(method).parameters[1]))"
    end
    return initpropwrk(state, tlist, method, generator...; kwargs)
end

function initpropwrk(state, tlist, method, generator...; kwargs...)
    throw(TypeError("`method` must be a known symbol, not $method"))
end

function initpropwrk(state, tlist; kwargs...)
    return initpropwrk(state, tlist, Val(:auto); kwargs...)
end

function initpropwrk(state, tlist, method::Symbol, generator...; kwargs...)
    return initpropwrk(state, tlist, Val(method), generator...; kwargs...)
end

function initpropwrk(
    state,
    tlist,
    method::Val{:cheby},
    generator...;
    specrad_method=Val(:auto),
    tolerance=0.01,
    kwargs...
)
    if length(generator) == 0
        throw(ArgumentError("The :cheby method requires at least one `generator`"))
    end
    # find a spectral envelope for all generators
    E_min, E_max = SpectralRange.specrange(generator[1], specrad_method; kwargs...)
    for G in generator[2:end]
        _E_min, _E_max = SpectralRange.specrange(G, specrad_method; kwargs...)
        E_min = (_E_min < E_min) ? _E_min : E_min
        E_max = (_E_max > E_max) ? _E_max : E_max
    end
    Δ = E_max - E_min
    δ = tolerance * Δ
    E_min = E_min - δ / 2
    Δ = Δ + δ
    dt = tlist[2] - tlist[1]
    allowed_kwargs = Set((:limit,))
    filtered_kwargs = filter(p -> p.first in allowed_kwargs, kwargs)
    return Cheby.ChebyWrk(state, Δ, E_min, dt; filtered_kwargs...)
end

function initpropwrk(state, tlist, method::Val{:newton}, generator...; kwargs...)
    allowed_kwargs = Set((:m_max,))
    filtered_kwargs = filter(p -> p.first in allowed_kwargs, kwargs)
    return Newton.NewtonWrk(state; filtered_kwargs...)
end

function initpropwrk(state, tlist, method::Val{:expprop}, generator...; kwargs...)
    # ExpPropWrk has no kwargs
    return ExpProp.ExpPropWrk(state)
end

"""Perform a single propagation step in-place.

```julia
propstep!(state, generator, dt, wrk; kwargs...)
```

The propagation method is determined by `wrk`, see [`initpropwrk`](@ref).

Generally, an in-place propagation will not be suitable for in the context of
automatic differentiation.  See [`propagate`](@ref) for a method that does not
act in-place.

The `kwargs` are forwarded to the underlying method.
"""
function propstep!(state, generator, dt, wrk::Cheby.ChebyWrk; kwargs...)
    return Cheby.cheby!(state, generator, dt, wrk; kwargs...)
end

function propstep!(state, generator, dt, wrk::Newton.NewtonWrk; kwargs...)
    return Newton.newton!(state, generator, dt, wrk; kwargs...)
end

function propstep!(state, generator, dt, wrk::ExpProp.ExpPropWrk; kwargs...)
    return ExpProp.expprop!(state, generator, dt, wrk; kwargs...)
end

"""Perform a single propagation step and return the propagated `state`.

```julia
state_out = propstep(state, generator, dt, wrk; kwargs...)
```

The propagation method is determined by `wrk`, see [`initpropwrk`](@ref).

Unlike [`propstep!`](@ref), this method does not act in place, which generally
makes it more suitable for automatic differentiation. However, there may be a
performance penalty associated with the additional memory allocations.

The `kwargs` are forwarded to the underlying method.
"""
function propstep(state, generator, dt, wrk::ExpProp.ExpPropWrk; kwargs...)
    return ExpProp.expprop(state, generator, dt, wrk; kwargs...)
end

# Default "observable" for storing the propagated state. We want to keep this
# private, as the routine is not safe unless it is the *only*
# observable. If state storage needs to be combined with other observables
# `state->copy(state)` would need to be used.
_store_state(state) = copy(state)
_store_state(state::Vector) = state

# Work around https://github.com/timholy/ProgressMeter.jl/issues/214
import ProgressMeter
struct NoProgress end
ProgressMeter.next!(p::NoProgress) = nothing

"""Propagate a state over an entire time grid.

```julia
state_out = propagate(
    state, genfunc, tlist; method=:auto,
    backwards=false; storage=nothing, observables=(<store state>, ),
    hook=nothing, showprogress=false, control_parameters=nothing,
    in_place=true, kwargs...)
```

propagates `state` over the time grid in `tlist`, using piecewise-constant
dynamical generators (Hamiltonians or Liouvillians) determined by `genfunc`,
and returns the resulting propagated state. The propagation is performed by
calling [`propstep!`](@ref) for every interval in `tlist`, or
[`propstep`](@ref) if `in_place=false`.

For the i'th time interval, `genfunc(tlist, i)` must return the generator for
that time interval. Generally, when approximating a time-continuous dynamical
generator as piecewise-constant on the time grid, it should be evaluated at the
*midpoint* of the interval. A possible exception is the first and last
interval, which may be better evaluated at `tlist[1]` and `tlist[end]` to
ensure exact boundary conditions like control fields that are exactly zero.

In addition to the two positional parameters indicating the time interval,
`genfunc` will also receive the `state` (the input state for the propagation
step), `backwards`, `storage`, `observables`, `control_parameters`, and `init`
as keyword arguments.

The `control_parameters` are an optional array of floats with parameters for
`genfunc`. This is required when `propagate` is used in the context of
automatic differentiation (AD). E.g., the
[Zygote](https://fluxml.ai/Zygote.jl/) framework can automatically calculate
gradients of a function `control_parameters -> J_T`, where `J_T` might be be a
function of the overlap between a propagated state (returned by `propagate`)
and a target state. Thus, the `control_parameters` must be *explicit* in
`propagate`.  Outside of an AD context, `control_parameters` are not generally
required: they can be implicit in `genfunc`.

The remaining keyword arguments may be used for unusual
equations of motion beyond the standard Schrödinger or Liouville-von-Neumann
equation, e.g. `state` would enter the `genfunc` for a Gross–Pitaevskii
equation. For standard equations of motion that do not use the additional
parameters, it is best to capture the keyword arguments to `genfunc` with a
definition like

```julia
genfunc(tlist, i; kwargs...) = ...
```

For valid propagation `method`s, see [`initpropwrk`](@ref).

In general, there is no requirement that `tlist` has a constant time step,
although some propagation methods (most notably [`cheby!`](@ref
QuantumPropagators.Cheby.cheby!)) only support a uniform time grid.

If `storage` is given as an Array, it will be filled with data determined by
the `observables`. The default "observable" results in the propagated states at
every point in time being stored.
The `storage` array should be created with [`init_storage`](@ref). See its
documentation for details.

The `storage` parameter may also be given as `true`, and a new storage array
will be created internally with [`init_storage`](@ref) and returned instead of
the propagated state:

```julia
data = propagate(
    state, genfunc, tlist; method=:auto
    backwards=false; storage=true, observables=observables,
    hook=nothing, showprogress=false, kwargs...)
```

If `backwards` is `true`, the input state is assumed to be at time
`tlist[end]`, and the propagation progresses backwards in time (with a negative
time step `dt`). If `storage` is given, it will be filled back-to-front during
the backwards propagation.

If `hook` is given as a callable, it will be called after each propagation
step, as `hook(state, generator, tlist, i, wrk, observables)` where `i` is the
index of the time interval on `tlist` covered by the propagation step (0 for
the initial state, respectives `lastindex(tlist)` for the backward
propagation).  The `hook` is called before calculating any observables. Example
usage includes writing data to file, or modifying `state`, e.g., removing
amplitude from the lowest and highest level to mitigate "truncation error".

If `showprogress` is given as `true`, a progress bar will be shown for
long-running propagationn. In order to customize the progress bar,
`showprogress` may also be a function that receives `length(tlist)` and returns
a `ProgressMeter.Progress` instance.

If `in_place=false` is given, the propagation avoids in-place operations by
using [`propstep`](@ref) instead of [`propstep!`](@ref). This
is often required in the context of automatic differentiation (AD), e.g., with
[Zygote](https://fluxml.ai/Zygote.jl/). That is, use `in_place=false` if
`propagate` is called inside a function to be passed to `Zygote.gradient`,
`Zygote.pullback`, or a similar function. In and AD context, `storage` and
`showprogress` should not be used.

The `propagate` routine returns the propagated state at `tlist[end]`,
respectively `tlist[1]` if `backwards=true`, or a storage array with the
stored states / observable data if `storage=true`.
"""
function propagate(state, genfunc, tlist; method=Val(:auto), kwargs...)
    return propagate(state, genfunc, tlist, method; kwargs...)
end

function propagate(state, genfunc, tlist, method::Symbol; kwargs...)
    return propagate(state, genfunc, tlist, Val(method); kwargs...)
end

# default `propagate` implementation
function propagate(state, genfunc, tlist, method::Val; kwargs...)
    # TODO: document what happens with kwargs
    backwards = get(kwargs, :backwards, false)
    storage = get(kwargs, :storage, nothing)
    observables = get(kwargs, :observables, (_store_state,))
    control_parameters = get(kwargs, :control_parameters, nothing)
    G = genfunc(
        tlist,
        1;
        state=state,
        backwards=backwards,
        storage=storage,
        observables=observables,
        control_parameters=control_parameters,
        init=true
    )
    wrk = initpropwrk(state, tlist, method, G; kwargs...)
    return _propagate(state, genfunc, tlist, wrk; kwargs...)
end

# `propagate` backend (note `wrk` argument instead of `method`)
function _propagate(
    state,
    genfunc,
    tlist,
    wrk;
    backwards=false,
    storage=nothing,
    observables=(_store_state,),
    hook=nothing,
    showprogress=false,
    control_parameters=nothing,
    in_place=true,
    kwargs...
)
    return_storage = false
    if storage === true
        storage = init_storage(state, tlist, observables)
        return_storage = true
    end
    state = copy(state)
    intervals = enumerate(tlist[2:end])
    if backwards
        intervals = Iterators.reverse(intervals)
        if storage ≠ nothing
            Storage.write_to_storage!(storage, lastindex(tlist), state, observables)
        end
    else
        if storage ≠ nothing
            Storage.write_to_storage!(storage, 1, state, observables)
        end
    end

    N = length(intervals)
    if isa(showprogress, Function)
        progressmeter = showprogress(N)
    else
        progressmeter = Progress(N, enabled=showprogress)
        if !showprogress
            # XXX: https://github.com/timholy/ProgressMeter.jl/issues/214
            progressmeter = NoProgress()
        end
    end

    for (i, t_end) in intervals
        dt = t_end - tlist[i]
        generator = genfunc(
            tlist,
            i;
            state=state,
            backwards=backwards,
            storage=storage,
            observables=observables,
            control_parameters=control_parameters,
            init=false
        )
        if in_place
            propstep!(state, generator, (backwards ? -dt : dt), wrk; kwargs...)
        else
            state = propstep(state, generator, (backwards ? -dt : dt), wrk; kwargs...)
        end
        if storage ≠ nothing
            Storage.write_to_storage!(storage, i + (backwards ? 0 : 1), state, observables)
        end
        if hook ≠ nothing
            hook(state, generator, tlist, i, wrk, observables)
        end
        next!(progressmeter)
    end
    if return_storage
        return storage
    else
        return state
    end
end

# `propagate` for Chebychev propagation
function propagate(state, genfunc, tlist, method::Val{:cheby}; kwargs...)
    backwards = get(kwargs, :backwards, false)
    storage = get(kwargs, :storage, nothing)
    observables = get(kwargs, :observables, (_store_state,))
    control_parameters = get(kwargs, :control_parameters, nothing)
    # TODO: generate multiple examplary G
    G = genfunc(
        tlist,
        1;
        state=state,
        backwards=backwards,
        storage=storage,
        observables=observables,
        control_parameters=control_parameters,
        init=true
    )
    wrk = initpropwrk(state, tlist, method, G; kwargs...)
    return _propagate(state, genfunc, tlist, wrk; kwargs...)
end


text:
### ARSState

# ARSState holds the internal state ("local variables") of the acceptance-rejection sampler

mutable struct ARSState{S<:ValueSupport, F<:VariateForm} <: MCSamplerState{F}
  pstate::ParameterState{S, F} # Parameter state used internally by ARS
  tune::MCTunerState
  logproposal::Real
  weight::Real
  diagnosticindices::Dict{Symbol, Integer}
end

ARSState(
  pstate::ParameterState{S, F},
  tune::MCTunerState=BasicMCTune()
) where {S<:ValueSupport, F<:VariateForm} =
  ARSState(pstate, tune, NaN, NaN, diagnosticindices)

### Acceptance-rejection sampler (ARS)

struct ARS <: MCSampler
  logproposal::Function # Possibly unnormalized log-proposal (envelope)
  proposalscale::Real # Scale factor for logproposal (envelope) to cover target
  jumpscale::Real # Scale factor for adapting jump size

  function ARS(logproposal::Function, proposalscale::Real, jumpscale::Real)
    @assert jumpscale > 0 "Scale factor for adapting jump size is not positive"
    new(logproposal, proposalscale, jumpscale)
  end
end

ARS(logproposal::Function; proposalscale::Real=1., jumpscale::Real=1.) = ARS(logproposal, proposalscale, jumpscale)

### Initialize ARS sampler

## Initialize parameter state

function initialize!(
  pstate::ParameterState{Continuous, F},
  parameter::Parameter{Continuous, F},
  sampler::ARS,
  outopts::Dict
) where F<:VariateForm
  parameter.logtarget!(pstate)
  @assert isfinite(pstate.logtarget) "Log-target not finite: initial value out of support"

  if !isempty(outopts[:diagnostics])
    pstate.diagnostickeys = copy(outopts[:diagnostics])
    pstate.diagnosticvalues = Array(Any, length(pstate.diagnostickeys))
  end
end

## Initialize ARSState

function sampler_state(
  parameter::Parameter{S, F},
  sampler::ARS,
  tuner::MCTuner,
  pstate::ParameterState{S, F},
  vstate::VariableStateVector,
  diagnostickeys::Vector{Symbol}
) where {S<:ValueSupport, F<:VariateForm}
  sstate = ARSState(generate_empty(pstate), tuner_state(parameter, sampler, tuner))
  set_diagnosticindices!(sstate, [:accept], diagnostickeys)
  sstate
end

## Reset parameter state

function reset!(
  pstate::ParameterState{Continuous, Univariate},
  x::Real,
  parameter::Parameter{Continuous, Univariate},
  sampler::ARS
)
  pstate.value = x
  parameter.logtarget!(pstate)
end

function reset!(
  pstate::ParameterState{Continuous, Multivariate},
  x::RealVector,
  parameter::Parameter{Continuous, Multivariate},
  sampler::ARS
)
  pstate.value = copy(x)
  parameter.logtarget!(pstate)
end

reset!(
  sstate::ARSState{S, F},
  pstate::ParameterState{S, F},
  parameter::Parameter{S, F},
  sampler::MCSampler,
  tuner::MCTuner
) where {S<:ValueSupport, F<:VariateForm} =
  reset!(sstate.tune, sampler, tuner)

show(io::IO, sampler::ARS) =
  print(io, "ARS sampler: proposal scale = $(sampler.proposalscale), jump scale = $(sampler.jumpscale)")


text:
"""
    split_on_whitespace(x)

Split a string on any amount of whitespaces.
"""
split_on_whitespace(x) = split(strip(x), r"(\s+(?=\S))")

"""
    log_error(err::Exception)

Log an exception and its backtrace to stderr.
"""
log_error(err::Exception) = @error "Function Failed: $(typeof(err))" exception=(err, catch_backtrace())

"""
    @lock expr

Locks a function.
"""
macro lock(lock, expr)
    lk = esc(lock)
    quote
        lock($lk) do
            $(esc(expr))
        end
    end
end


text:
using GFitViewer
import GFitViewer: ViewerData, viewer

function ViewerData(res::QSFitResults{T}; kw...) where T
    vd = ViewerData(res.model, res.pspec.data, res.fitres; kw...)
    vd.dict[:meta][:banner] = "QSFit (v0.1)
Date: " * string(trunc(res.fitres.timestamp, Second))

    id = 1
    m = vd.dict[:models][id][:meta]
    m[:label] = res.source.name * ", z=" * string(res.source.z) * ", E(B-V)=" * string(res.source.mw_ebv)
    m[:label_x] = "Rest frame wavelength"
    m[:unit_x]  = string(QSFit.unit_λ())
    m[:log10scale_x] = 0
    m[:label_y] = "Lum. density"

    if "UNITFUL_FANCY_EXPONENTS" in keys(ENV)
        orig = ENV["UNITFUL_FANCY_EXPONENTS"]
        ENV["UNITFUL_FANCY_EXPONENTS"] = true
        m[:unit_y]  = string(QSFit.unit_lum_density())
        ENV["UNITFUL_FANCY_EXPONENTS"] = orig
    else
        ENV["UNITFUL_FANCY_EXPONENTS"] = true
        m[:unit_y]  = string(QSFit.unit_lum_density())
        delete!(ENV, "UNITFUL_FANCY_EXPONENTS")
    end
    m[:log10scale_y] = QSFit.log10_scale_lum()
    vd.dict[:data][id][:meta][:label] = res.source.specs[id].label

    m = vd.dict[:extra][id]
    m[:EW] = GFitViewer.MDict()
    m[:EW][:label] = "Em. lines EW"
    m[:EW][:fields] = GFitViewer.MDict()
    m[:EW][:fields][:Label] = GFitViewer.MDict()
    m[:EW][:fields][:Label][:meta] = GFitViewer.MDict()
    m[:EW][:fields][:Label][:meta][:desc] = "Em. line"
    m[:EW][:fields][:Label][:data] = collect(keys(res.EW))
    m[:EW][:fields][:Value] = GFitViewer.MDict()
    m[:EW][:fields][:Value][:meta] = GFitViewer.MDict()
    m[:EW][:fields][:Label][:meta][:desc] = "EW [A]"
    m[:EW][:fields][:Value][:data] = collect(values(res.EW))

    m[:extratab_2] = GFitViewer.MDict()
    m[:extratab_2][:label] = "Second extra table"
    m[:extratab_2][:fields] = GFitViewer.MDict()
    m[:extratab_2][:fields][:fname_1] = GFitViewer.MDict()
    m[:extratab_2][:fields][:fname_1][:meta] = GFitViewer.MDict()
    m[:extratab_2][:fields][:fname_1][:meta][:desc] = "Optional. In case we want to add metedata."
    m[:extratab_2][:fields][:fname_1][:data] = [102, 203, 304]
    m[:extratab_2][:fields][:fname_2] = GFitViewer.MDict()
    m[:extratab_2][:fields][:fname_2][:meta] = GFitViewer.MDict()
    m[:extratab_2][:fields][:fname_2][:data] = [10.2, 20.3, 30.4]
    m[:extratab_2][:fields][:fname_3] = GFitViewer.MDict()
    m[:extratab_2][:fields][:fname_3][:meta] = GFitViewer.MDict()
    m[:extratab_2][:fields][:fname_3][:data] = ["String_1", "String_2", "String_3"]
    return vd
end

viewer(res::QSFitResults{T}; filename=nothing, offline=false, kw...) where T =
    viewer(ViewerData(res; kw...); filename=filename, offline=offline)

function ViewerData(res::QSFitMultiResults{T}; kw...) where T
    vd = ViewerData(res.multi, [res.pspecs[id].data for id in 1:length(res.pspecs)], res.fitres; kw...)
    vd.dict[:meta][:banner] = "QSFit (v0.1)
Date: " * string(trunc(res.fitres.timestamp, Second))

    for id in 1:length(res.multi)
        m = vd.dict[:models][id][:meta]
        m[:label] = res.source.name * ", z=" * string(res.source.z) * ", E(B-V)=" * string(res.source.mw_ebv)
        m[:label_x] = "Rest frame wavelength"
        m[:unit_x]  = string(QSFit.unit_λ())
        m[:log10scale_x] = 0
        m[:label_y] = "Lum. density"

        if "UNITFUL_FANCY_EXPONENTS" in keys(ENV)
            orig = ENV["UNITFUL_FANCY_EXPONENTS"]
            ENV["UNITFUL_FANCY_EXPONENTS"] = true
            m[:unit_y]  = string(QSFit.unit_lum_density())
            ENV["UNITFUL_FANCY_EXPONENTS"] = orig
        else
            ENV["UNITFUL_FANCY_EXPONENTS"] = true
            m[:unit_y]  = string(QSFit.unit_lum_density())
            delete!(ENV, "UNITFUL_FANCY_EXPONENTS")
        end
        m[:log10scale_y] = QSFit.log10_scale_lum()
        vd.dict[:data][id][:meta][:label] = res.source.specs[id].label

        m = vd.dict[:extra][id]
        m[:EW] = GFitViewer.MDict()
        m[:EW][:label] = "Em. lines EW"
        m[:EW][:fields] = GFitViewer.MDict()
        m[:EW][:fields][:Label] = GFitViewer.MDict()
        m[:EW][:fields][:Label][:meta] = GFitViewer.MDict()
        m[:EW][:fields][:Label][:meta][:desc] = "Em. line"
        m[:EW][:fields][:Label][:data] = collect(keys(res.EW[id]))
        m[:EW][:fields][:Value] = GFitViewer.MDict()
        m[:EW][:fields][:Value][:meta] = GFitViewer.MDict()
        m[:EW][:fields][:Label][:meta][:desc] = "EW [A]"
        m[:EW][:fields][:Value][:data] = collect(values(res.EW[id]))

        m[:extratab_2] = GFitViewer.MDict()
        m[:extratab_2][:label] = "Second extra table"
        m[:extratab_2][:fields] = GFitViewer.MDict()
        m[:extratab_2][:fields][:fname_1] = GFitViewer.MDict()
        m[:extratab_2][:fields][:fname_1][:meta] = GFitViewer.MDict()
        m[:extratab_2][:fields][:fname_1][:meta][:desc] = "Optional. In case we want to add metedata."
        m[:extratab_2][:fields][:fname_1][:data] = [102, 203, 304]
        m[:extratab_2][:fields][:fname_2] = GFitViewer.MDict()
        m[:extratab_2][:fields][:fname_2][:meta] = GFitViewer.MDict()
        m[:extratab_2][:fields][:fname_2][:data] = [10.2, 20.3, 30.4]
        m[:extratab_2][:fields][:fname_3] = GFitViewer.MDict()
        m[:extratab_2][:fields][:fname_3][:meta] = GFitViewer.MDict()
        m[:extratab_2][:fields][:fname_3][:data] = ["String_1", "String_2", "String_3"]
    end
    return vd
end

viewer(res::QSFitMultiResults{T}; filename=nothing, offline=false, kw...) where T =
    viewer(ViewerData(res; kw...); filename=filename, offline=offline)


text:
export @unit
macro unit()
  str="""
  metre=1e2;meter=1e2;cm=metre/1e2;mm=metre/1e3;um=metre/1e6;nm=metre/1e9;
  second=1e3;ms=second/1e3;
  Hz=1/second;
  voltage=1e3;mV=voltage/1e3;
  ampere=1e6;mA=ampere/1e3;uA=ampere/1e6;nA=ampere/1e9;
  farad=1e6;uF=ufarad=farad/1e6;
  siemens=1e3;mS=msiemens=siemens/1e3;nS=nsiemens=siemens/1e9;"""
  print(str,"\n")
  eval(parse(replace(str,"\n","")))
end


text:
module q2

import Evolve

steps = 16000               # number of simulation steps
samplesteps = 50            # steps between measurements
L = 2                       # width of the box
xmin = -1                   # leftmost point in the box
N = 2000                    # number of subintervals
ϵ = L/N                     # width of each subinterval
ks = [100, 200]              # wavenumbers
σs = [0.05, 0.02]            # wavepacket width
x0 = 0.0                    # wavepacket starting position

# Run the trials
function run()
for k in ks
    for σ in σs
        print("Initializing for k = $k, σ = $σ...")
        V(xs) = map(x->0.0, xs)     # zero potential energy
        AV = Evolve.A(V)            # evolution matrix constructor for square well
        ϕV = Evolve.Φ(k, σ, x0)     # initial state vector constructor for Gaussian packet

        A = AV(L, N, xmin)
        ϕ0 = ϕV(L, N, xmin)
        x = (L/N .* [0:N] + xmin)   # x positions

        print("evolving...")
        ϕ = Evolve.evolve(A, ϕ0, steps)
        phis = [ϕ0 ϕ[:,samplesteps:samplesteps:end]]

        println("done. Writing to file.")
        writedlm("phis_im_$k\_$σ.dat", imag(phis))
        writedlm("phis_re_$k\_$σ.dat", real(phis))
    end
end
end

# Run the trials
for k in ks
    for σ in σs
        print("Reverse run for k = $k, σ = $σ...")
        V(xs) = map(x->0.0, xs)     # zero potential energy
        AV = Evolve.A(V)            # evolution matrix constructor for square well
        ϕV = Evolve.Φ(k, σ, x0)     # initial state vector constructor for Gaussian packet

        A = AV(L, N, xmin)
        ϕ0 = ϕV(L, N, xmin)
        x = (L/N .* [0:N] + xmin)   # x positions

        print("evolving...")
        ϕ = Evolve.evolve(A, ϕ0, 2000)
        phis = [ϕ0 ϕ[:,samplesteps:samplesteps:end]]
        ϕr = conj(ϕ[:,end])         # reverse time evolution by taking conjugate
        ϕ = Evolve.evolve(A, ϕr, 2000)
        phis = [phis ϕ[:,samplesteps:samplesteps:end]]

        println("done. Writing to file.")
        writedlm("phis_rev_im_$k\_$σ.dat", imag(phis))
        writedlm("phis_rev_re_$k\_$σ.dat", real(phis))
    end
end

end


text:
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/build_literate/tutorial.ipynb)

# # Quick Start Tutorial

# The core of generalized Bloch model is implemented in the function [`apply_hamiltonian_gbloch!(∂m∂t, m, mfun, p, t)`](@ref), which calculates the derivative `∂m/∂t` for a given magnetization vector `m` and stores it in-place in the the variable `∂m∂t`. The function interface is written in a way that we can directly feed it into a differential equation solver of the [DifferentialEquations.jl](https://diffeq.sciml.ai/stable/) package.

# For this example, we need the following packages:
using MRIgeneralizedBloch
using DifferentialEquations
using Plots
plotlyjs(bg = RGBA(31/255,36/255,36/255,1.0), ticks=:native); #hide #!nb

# We simulate the dynamics of a coupled spin system with the following parameters:
m0s = 0.15
R1f = 0.3 # 1/s
R2f = 15 # 1/s
R1s = 2 # 1/s
T2s = 10e-6 # s
Rx = 30; # 1/s

# and the thermal equilibrium of the magnetization `m = [xf; yf; zf; zs; 1]`:
m0 = [0; 0; 1-m0s; m0s; 1];

# during a rectangular RF-pulse with the flip angle and pulse duration
α = π
TRF = 100e-6; # s

# Further, we assume a perfectly calibrated, on-resonant RF-pulse:
B1 = 1
ω0 = 0; # rad/s

# as well as a [super-Lorentzian lineshape](http://dx.doi.org/10.1002/mrm.1910330404). We interpolate the corresponding Green's function to improve performance:
G = interpolate_greens_function(greens_superlorentzian, 0, TRF / T2s);

# The generalized Bloch model is a so-called integro-differential equation where the derivative ``∂m/∂t`` at the time ``t_1`` does not just depend on ``m(t_1)``, but on ``m(t)`` for ``t \in [0, t_1]``. This is solved with a [delay differential equation (DDE) solver](https://diffeq.sciml.ai/stable/tutorials/dde_example/) that stores an interpolated *history function* `mfun(p, t)`, which we use in the [`apply_hamiltonian_gbloch!`](@ref) function to evaluate the integral. This history function has to be initialized with `mfun(p, 0) = m0`. Here, we use a slightly more complicated initialization that allows us to index the history function in [`apply_hamiltonian_gbloch!`](@ref), which improves performance:
mfun(p, t; idxs=nothing) = typeof(idxs) <: Number ? m0[idxs] : m0;

# With this, we are ready to formulate and solve the differential equation:
param = (α/TRF, B1, ω0, m0s, R1f, R2f, Rx, R1s, T2s, G) # defined by apply_hamiltonian_gbloch!
prob = DDEProblem(apply_hamiltonian_gbloch!, m0, mfun, (0, TRF), param)
sol = solve(prob)

# The plot function is implemented for such solution objects and we can plot the solution simply with
p = plot(sol, labels=["xᶠ" "yᶠ" "zᶠ" "zˢ" "1"], xlabel="t [s]", ylabel="m(t)")
#md Main.HTMLPlot(p) #hide

# More details on the interface, including the linear approximation of the generalized Bloch model can found in the following scripts that replicate all simulations, data analyses, and figures of the generalized Bloch paper.