text:
@reexport module Threshold

using Adapt
using Random
using CUDA

import ..SpikingNN: excite!, evaluate!, reset!, isactive

export AbstractThreshold, isactive

abstract type AbstractThreshold end

"""
    Ideal(vth::Real)

An ideal threshold spikes when `v > vth`.
"""
struct Ideal{T<:Real} <: AbstractThreshold
    vth::T
end

isactive(threshold::Ideal, t::Integer; dt::Real = 1.0) = false

"""
    evaluate!(threshold::Ideal, t::Real, v::Real; dt::Real = 1.0)
    (threshold::Ideal)(t::Real, v::Real; dt::Real = 1.0)
    evaluate!(thresholds::AbstractArray{<:Ideal}, t::Integer, v; dt::Real = 1.0)

Return `t` when `v > threshold.vth`.
"""
evaluate!(threshold::Ideal, t::Real, v::Real; dt::Real = 1.0) = (v >= threshold.vth) ? t : zero(t)
(threshold::Ideal)(t::Real, v::Real; dt::Real = 1.0) = evaluate!(threshold, t, v; dt = dt)
evaluate!(thresholds::T, t::Integer, v; dt::Real = 1.0) where T<:AbstractArray{<:Ideal} =
    Int.((v .>= thresholds.vth) .* adapt(typeof(v), fill(t, size(v))))

"""
    Poisson(ρ₀::Real, Θ::Real, Δᵤ::Real, rng::AbstractRNG)

Choose to output a spike based on a inhomogenous Poisson process given by

``X < \\mathrm{d}t \\: \\rho_0 \\exp\\left(\\frac{v - \\Theta}{\\Delta_u}\\right)``

where ``X \\sim \\mathrm{Unif}([0, 1])``.
Accordingly, `dt` must be set correctly so that the neuron does not always spike.

Fields:
- `ρ₀::Real`: baseline firing rate at threshold
- `Θ::Real`: firing threshold
- `Δᵤ::Real`: voltage resolution
- `rng`: random number generation
"""
struct Poisson{T<:Real,  RT<:AbstractRNG} <: AbstractThreshold
    ρ₀::T
    Θ::T
    Δᵤ::T
    rng::RT
end

function Poisson{T}(ρ₀::Real, Θ::Real, Δᵤ::Real; rng::RT = Random.GLOBAL_RNG) where {T<:Real, RT}
    Poisson{T,RT}(ρ₀, Θ, Δᵤ, rng)
end
Poisson(ρ₀::Real, Θ::Real, Δᵤ::Real; kwargs...) = Poisson{Real}(ρ₀, Θ, Δᵤ; kwargs...)

isactive(threshold::Poisson, t::Integer; dt::Real = 1.0) = true

"""
    poisson(baserate, theta, deltav, v; dt::Real, rng::AbstractRNG)
    poisson(baserate::AbstractArray{<:Real}, theta::AbstractArray{<:Real}, deltav::AbstractArray{<:Real}, v::AbstractArray{<:Real}; dt::Real, rng::AbstractRNG)
    poisson(baserate::CuVector{<:Real}, theta::CuVector{<:Real}, deltav::CuVector{<:Real}, v::CuVector{<:Real}; dt::Real, rng::AbstractRNG)

Evaluate inhomogeneous Poisson process threshold functions.
Modeled as

``X < \\mathrm{d}t \\rho_0 \\exp\\left(\\frac{v - \\Theta}{\\Delta_u}\\right)``

where ``X \\sim \\mathrm{Unif}([0, 1])``.

Use `CuVector` instead of `Vector` to evaluate on GPU.

# Fields
- `baserate`: base line firing rate
- `theta`: threshold potential
- `deltav`: potential resolution
- `v`: current membrane potential
- `dt`: simulation timestep
- `rng`: random number generation
"""
function poisson(baserate, theta, deltav, v; dt::Real, rng::AbstractRNG = Random.GLOBAL_RNG)
    rho = baserate * exp((v - theta) / deltav)

    return rand(rng) < rho * dt
end
function poisson(baserate::AbstractArray{<:Real}, theta::AbstractArray{<:Real}, deltav::AbstractArray{<:Real}, v::AbstractArray{<:Real}; dt::Real, rng::AbstractRNG = Random.GLOBAL_RNG)
    rho = baserate .* exp.((v .- theta) ./ deltav)

    return rand(rng, length(rho)) .< rho .* dt
end
function poisson(baserate::CuVector{<:Real}, theta::CuVector{<:Real}, deltav::CuVector{<:Real}, v::CuVector{<:Real}; dt::Real, rng::AbstractRNG = Random.GLOBAL_RNG)
    rho = baserate .* exp.((v .- theta) ./ deltav)

    return CuArrays.rand(rng, length(rho)) .< rho .* dt
end

"""
    evaluate!(threshold::Poisson, t::Integer, v::Real; dt::Real = 1.0)
    (::Poisson)(t::Integer, v::Real; dt::Real = 1.0)
    evaluate!(thresholds::AbstractArray{<:Poisson}, t::Integer, v; dt::Real = 1.0)

Evaluate Poisson threshold function. See [`Threshold.Poisson`](@ref).
"""
evaluate!(threshold::Poisson, t::Integer, v::Real; dt::Real = 1.0) =
    poisson(threshold.ρ₀, threshold.Θ, threshold.Δᵤ, v; dt = dt, rng = threshold.rng) ? t : zero(t)
(threshold::Poisson)(t::Integer, v::Real; dt::Real = 1.0) = evaluate!(threshold, t, v; dt = dt)
evaluate!(thresholds::T, t::Integer, v; dt::Real = 1.0) where T<:AbstractArray{<:Poisson} =
    Int.(poisson(thresholds.ρ₀, thresholds.Θ, thresholds.Δᵤ, v; dt = dt, rng = first(thresholds.rng)) .* adapt(typeof(v), fill(t, size(v))))

end


text:
module Main
export fib

fib(x::Number) = n < 2 ? n : fibo(n-1) + fibo(n-2)

end

ext: jl

max_stars_repo_name: vmchale/project-init

max_issues_repo_path: julia-lib/src/{{ Project }}.jl

max_issues_repo_name: vmchale/pi-templates

max_forks_repo_name: vmchale/project-init


text:
"""Performs a single step of the EnsembleKalmanProcess."""

# Import modules to all processes
using ArgParse
using CalibrateEDMF
using CalibrateEDMF.DistributionUtils
using CalibrateEDMF.Pipeline
using CalibrateEDMF.HelperFuncs
# Import EKP modules
using EnsembleKalmanProcesses
using EnsembleKalmanProcesses.ParameterDistributions
using JLD2

# Read iteration number of ensemble to be recovered
s = ArgParseSettings()
@add_arg_table s begin
    "--iteration"
    help = "Calibration iteration number"
    arg_type = Int
    "--job_dir"
    help = "Job output directory"
    arg_type = String
    default = "output"
end
parsed_args = parse_args(ARGS, s)
# Recover inputs for update
iteration = parsed_args["iteration"]
outdir_path = parsed_args["job_dir"]

include(joinpath(outdir_path, "config.jl"))

versions = readlines(joinpath(outdir_path, "versions_$(iteration).txt"))
priors = deserialize_prior(load(joinpath(outdir_path, "prior.jld2")))
ekobj = load(ekobj_path(outdir_path, iteration))["ekp"]
config = get_config()
ek_update(ekobj, priors, iteration, config, versions, outdir_path)


text:
using Distances

function calculate_stabilities(allSignatures, exposures, processStabAvg)

    recon = allSignatures*exposures;
    sz = size(exposures);
    totalSamples = sz[2];

    avgReconstructionError = norm(recon-originalGenomes);
    norm_genome_fro = norm(originalGenomes);
    avgReconstructionErrorPercentage = avgReconstructionError/norm_genome_fro;
    avgStability = mean(processStabAvg);

    reconError = zeros(totalSamples);
    reconErrorPercentage = zeros(totalSamples);
    simError = zeros(totalSamples);

    for j = 1 : totalSamples
        error_sample = norm(originalGenomes[:,j] - recon[:,j]);
        norm_sample_fro = norm(originalGenomes[:,j]);
        reconError[j] = error_sample;
        reconErrorPercentage[j] = error_sample/norm_sample_fro;
        accr = 1-cosine_dist(recon[:, j], originalGenomes[:, j]);
        simError[j] = accr;
    end

    return avgReconstructionError, avgReconstructionErrorPercentage, avgStability, reconError, reconErrorPercentage, simError
end

max_stars_repo_name: ypradat/sigprofilerjulia

max_issues_repo_name: ypradat/sigprofilerjulia

max_forks_repo_name: ypradat/sigprofilerjulia


text:
struct StackedArea
    time_range::Array
    data_matrix::Matrix
    labels::Array
end

struct BarPlot
    time_range::Array
    bar_data::Matrix
    labels::Array
end

struct StackedGeneration
    time_range::Array
    data_matrix::Matrix
    labels::Array
end

struct BarGeneration
    time_range::Array
    bar_data::Matrix
    labels::Array
end

"""

get_stacked_plot_data(res::OperationsProblemResults, variable::String)

This function takes in results of struct OperationsProblemResult. It takes the
dataframe from whichever variable name was given and converts it to type StackedArea.
StackedArea is the type of struct that signals the plot() function to use the
StackedArea plot recipe method.

#Example

```julia
ThermalStandard = get_stacked_plot_data(res, "P_ThermalStandard")
plot(ThermalStandard)
```
"""

function get_stacked_plot_data(res::OperationsProblemResults, variable::String; kwargs...)

    sort = get(kwargs, :sort, nothing)
    time_range = res.time_stamp[!,:Range]
    variable = res.variables[Symbol(variable)]
    alphabetical = sort!(names(variable))

    if isnothing(sort)
        variable = variable[:, alphabetical]
    else
        variable = variable[:,sort]
    end

    data_matrix = convert(Matrix, variable)
    labels = collect(names(variable))
    legend = string.(labels)

    return StackedArea(time_range, data_matrix, legend)

end

function get_bar_plot_data(res::OperationsProblemResults, variable::String; kwargs...)

    sort = get(kwargs, :sort, nothing)
    time_range = res.time_stamp[!,:Range]
    variable = res.variables[Symbol(variable)]
    alphabetical = sort!(names(variable))

    if isnothing(sort)
        variable = variable[:, alphabetical]
    else
        variable = variable[:,sort]
    end

    data = convert(Matrix, variable)
    bar_data = sum(data, dims = 1)
    labels = collect(names(variable))
    legend = string.(labels)

    return BarPlot(time_range, bar_data, legend)

end

function get_stacked_generation_data(res::OperationsProblemResults; kwargs...)

    sort = get(kwargs, :sort, nothing)
    time_range = res.time_stamp[!,:Range]
    key_name = collect(keys(res.variables))
    alphabetical = sort!(key_name)

    if !isnothing(sort)
        labels = sort
    else
        labels = alphabetical
    end

    variable = res.variables[Symbol(labels[1])]
    data_matrix = sum(convert(Matrix, variable), dims = 2)
    legend = string.(labels)

    for i in 1:length(labels)
        if i !== 1
            variable = res.variables[Symbol(labels[i])]
            data_matrix = hcat(data_matrix, sum(convert(Matrix, variable), dims = 2))
        end
    end

    return StackedGeneration(time_range, data_matrix, legend)

end

function get_bar_gen_data(res::OperationsProblemResults)

   time_range = res.time_stamp[!,:Range]
   key_name = collect(keys(res.variables))

   variable = res.variables[Symbol(key_name[1])]
   data_matrix = sum(convert(Matrix, variable), dims = 2)
   legend = string.(key_name)

    for i in 1:length(key_name)
       if i !== 1
           variable = res.variables[Symbol(key_name[i])]
           data_matrix = hcat(data_matrix, sum(convert(Matrix, variable), dims = 2))
       end
   end
   bar_data = sum(data_matrix, dims = 1)
   return BarGeneration(time_range, bar_data, legend)

end

"""
    sort_data(results::OperationsProblemResults)

This function takes in struct OperationsProblemResults,
sorts the generators in each variable, and outputs the sorted
results. The generic function sorts the generators alphabetically.

# Arguments
-`results::OperationsProblemResults`: the results of the simulation

# Accepted Key Words
-`Variables` to choose which variables to be sorted.

#Examples
```julia
new_results = sort_data(results)
```
***Note:*** only the generators included in 'my_order' will be in the
results, and consequently, only these will be plotted. This can be a nice
feature for variables with more than 5 generators.

"""
function sort_data(res::OperationsProblemResults; kwargs...)

    Variables = Dict()
    Variables[:P_ThermalStandard]  = get(kwargs, :P_ThermalStandard, nothing)
    Variables[:P_RenewableDispatch]  = get(kwargs, :P_RenewableDispatch, nothing)
    Variables[:START_ThermalStandard] = get(kwargs, :START_ThermalStandard, nothing)
    Variables[:STOP_ThermalStandard] = get(kwargs, :STOP_ThermalStandard, nothing)
    Variables[:ON_ThermalStandard] = get(kwargs, :ON_ThermalStandard, nothing)
    Variable_dict = get(kwargs, :Variables, nothing)
    key_name = collect(keys(res.variables))
    alphabetical = sort!(key_name)

    if !isnothing(Variable_dict)
        labels = Variable_dict
    else
        labels = alphabetical
    end
    variable_dict = Dict()
    for label in labels
        variable_dict[labels[label]] = res.variables[labels[label]]
    end
    for (k,v) in Variables, k in keys(variable_dict)
        variable = variable_dict[k]
        alphabetical = sort!(names(variable))
        order = Variables[k]
        if isnothing(order)
            variable = variable[:, alphabetical]
        else
            variable = variable[:, order]
        end
        variable_dict[k] = variable
    end
    res = OperationsProblemResults(variable_dict, res.total_cost, res.optimizer_log, res.time_stamp)
    return res
end


text:
abstract type SupportConstrained <: ConvexSet end

"""
#   ConstrainedBySupport(support)

Special type of convex set.

For continuous case: consists of all functions that equals zero outside some fixed area called support:
`` 𝓐_S = \\{f : f(x) = 0, x ∉ S \\} .``

For discrete case: all arrays that equals zero for indexes outside some index set:
`` 𝓐_S = \\{x : x[i] = 0, i ∉ S \\} .``

Currently supports only discrete case, with the support defined as a boolean array.

- Julia version:
- Author: Oleg Soloviev
- Date: 2019-09-01

# Examples

```jldoctest

julia> S =  ConstrainedBySupport([true, false,true])
ConstrainedBySupport(Bool[true, false, true])

julia> x = [1, 2, 3]; project(x, S)
3-element Array{Int64,1}:
 1
 0
 3

julia> S = ConstrainedBySupport([x^2 + y^2 <=1  for x in -2:.2:2, y in -2:.2:2]);
julia> x = ones(size(S.support));
julia> project(x,S)
21×21 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

```
"""
struct ConstrainedBySupport <: SupportConstrained
    support::Array{Bool}
end
export ConstrainedBySupport

# function project(x, feasset::ConstrainedBySupport)
#     return feasset.support .* x
# end

function project!(xp, x, feasset::ConstrainedBySupport)
    @inbounds for i in eachindex(xp)
        xp[i] = feasset.support[i] * x[i]
    end

    return xp
end


text:
function  generatecode(app,mesh)

    # search compilers and set options
    app = Gencode.setcompilers(app);

    # generate input files and store them in datain folder
    app, mesh, master, dmd = preprocessing(app,mesh);

    # generate source codes and store them in app folder
    Gencode.gencode(app);

    return app,mesh,master,dmd

end

max_stars_repo_name: rloekvh/Exasim

max_issues_repo_name: rloekvh/Exasim

max_forks_repo_name: rloekvh/Exasim


text:
abstract type Model end

"""
Holds the user-set parameter values.

**values set in param.json file:**

* `gamma`: CRRA coefficient
* `beta`: discount factor
* `Tbar`: maximal age
* `R`: gross interest rate (i.e. 1+r )
* `na`: number of grid points for assets
* `ny`: number of grid points for income
* `a_lowT`: lower bound on assets in final period
* `a_low`: lower bound on assets
* `a_high`: upper bound on assets
* `mu`: unconditional mean of income (iid case)
* `sigma`: unconditional variance of income (iid case)
* `rho_z`: AR1 coefficient of income (AR1 case)
* `eps_z`: standard deviation of AR1 innovation (AR1 case)
* `alpha`: disutility of work
* `inc0`: income equation
* `inc1`: income equation
* `inc2`: income equation
* `cfloor`: consumption floor
* `delta`: probability of going back to work

"""
mutable struct Param

	gamma                 :: Float64  # CRRA
	neg_gamma             :: Float64
	oneminusgamma         :: Float64
	oneover_oneminusgamma :: Float64
	neg_oneover_gamma     :: Float64
	beta                  :: Float64
	sigma                 :: Float64
	lambda                :: Float64
	R                     :: Float64
	na                    :: Int 	# asset grid
	ny                    :: Int   # income grid (support points)
	nT                    :: Int   # maximal age
	a_high                :: Float64
	a_low                 :: Float64
	a_lowT                :: Float64
	nD                    :: Int # number of discrete choices
	cfloor                :: Float64
	cfloor_plot           :: Float64
	alpha                 :: Float64
	alphaT                 :: Float64
	inc0                 :: Float64
	inc1                 :: Float64
	inc2                 :: Float64
	ρ                    :: Float64
	k                    :: Int
	retage               :: Int # mandatory retirement age
	pension              :: Float64 # pension
	delta 				 :: Float64 # probability of going back to work
	maxinterest          :: Float64 # maximal interest rate

	bbar :: Float64 # bequest luxury goods level
	ν    :: Float64 # strength of bequest motive

	# simulation
	nsims                 :: Int
	initw0                :: Float64   # low/high bounds on initial wealth from this interval
	initw1                :: Float64   # low/high bounds on initial wealth from this interval
	rseed                 :: Int     # random generator seed

	# constructor
    function Param(;par::Dict=Dict())
		f=open(joinpath(dirname(@__FILE__),"param.json"))
		j = JSON.parse(f)
		close(f)
    	this = new()
    	for (k,v) in j
    		setfield!(this,Symbol(k),v["value"])
    	end

        if length(par) > 0
            # override parameters from dict par
            for (k,v) in par
                setfield!(this,k,v)
            end
        end

    	# derived
		this.neg_gamma             = (-1.0) * this.gamma
		this.oneminusgamma         = 1.0 - this.gamma
		this.oneover_oneminusgamma = 1.0 / this.oneminusgamma
		this.neg_oneover_gamma     = (-1.0) / this.gamma

		# this.beta = 1/this.R  # hard wire that beta = 1/R
		# this is not good.

		return this
	end
end

# """
# Binary Choice GModel with iid income uncertainty

# uses cash-on-hand m=y+a as state variable
# """
# mutable struct GModel

# 	# nD is number of discrete choices: nD = 2

# 	# computation grids
# 	avec::Vector{Float64}
# 	yvec::Vector{Float64}   # income support
# 	ywgt::Vector{Float64}   # income weights

# 	# intermediate objects (na,ny,nD)
# 	m1::Dict{Int,Dict}	# a dict[it] for each period
# 	c1::Matrix{Float64}
# 	ev::Matrix{Float64}

# 	# result objects
# 	v :: Matrix{Envelope}  # a vector of Envelope objects
# 	c :: Matrix{Envelope}

# 	"""
# 	Constructor for discrete choice GModel
# 	"""
# 	function GModel(p::Param)

# 		this = new()
# 		# avec          = scaleGrid(0.0,p.a_high,p.na,2)
# 		this.avec          = collect(range(p.a_low,stop = p.a_high,length = p.na))

# 		# fedors vesion:
# 		# nodes,weights = quadpoints(p.ny,0,1)
# 		# N = Normal(0,1)
# 		# nodes = quantile.(N,nodes)
# 		# this.yvec = nodes * p.sigma
# 		# this.ywgt = weights

# 		# my version:
# 		# for y ~ N(mu,sigma), integrate y with
# 		# http://en.wikipedia.org/wiki/Gauss-Hermite_quadrature
# 		nodes,weights = gausshermite(p.ny)  # from FastGaussQuadrature
# 		this.yvec = sqrt(2.0) * p.sigma .* nodes
# 		this.ywgt = weights .* pi^(-0.5)

# 		# precompute next period's cash on hand.
# 		# (na,ny,nD)
# 		# iD = 1: no work
# 		# iD = 2: work
# 		this.m1 = Dict(it => Dict(id => Float64[this.avec[ia]*p.R + income(it,p,this.yvec[iy]) * (id-1) for ia in 1:p.na, iy in 1:p.ny  ] for id=1:p.nD) for it=1:p.nT)
# 		this.c1 = zeros(p.na,p.ny)
# 		this.ev = zeros(p.na,p.ny)

# 		this.v = [Envelope(MLine(fill(NaN,(p.na)),fill(NaN,(p.na)))) for id in 1:p.nD, it in 1:p.nT]
# 		this.c = [Envelope(MLine(fill(NaN,(p.na)),fill(NaN,(p.na)))) for id in 1:p.nD, it in 1:p.nT]
# 		# dchoice = [it => ["d" => zeros(Int,p.na), "Vzero" => 0.0, "threshold" => 0.0] for it=1:p.nT]

# 		return this
# 	end
# end

"""
Fedor's Model
"""
mutable struct FModel <: Model
	avec::Vector{Float64}
	yvec::Vector{Float64}   # income support
	ywgt::Vector{Float64}   # income support

	# intermediate objects (na,nD)
	m1::Dict{Int,Dict}	# a dict[it] for each period
	c1::Matrix{Float64}
	ev::Matrix{Float64}

	# result objects
	v :: Array{Envelope}  # arrays of Envelope objects
	c :: Array{Envelope}

	vdirty :: Array{Envelope}  # arrays of Envelope objects
	cdirty :: Array{Envelope}

	function FModel(p::Param)

		this = new()

		# fedors version:
		nodes,weights = quadpoints(p.ny,0,1)
		N = Normal(0,1)
		nodes = quantile.(N,nodes)
		this.yvec = nodes * p.sigma
		this.ywgt = weights

		# simulation shocks

		if p.a_low < 0
			@warn "FModel not built for BK and neg assets"
		end
		this.avec          = collect(range(p.a_low,stop = p.a_high,length = p.na))

		# precompute next period's cash on hand.
		# (na,ny,nD)
		# iD = 1: tomorrow work, income is work income
		# iD = 2: tomorrow no work - if delta = 0, pension is absorbing, else get penion
		# notice: you decide to work today (t), but your shock realises tomorrow (t+1) in their formulation
		this.m1 = Dict(it => Dict(id => Float64[this.avec[ia]*p.R .+ income(it+1,p,this.yvec[iy]) * (id==1) +
		                                (id==2) * (p.delta * income(it+1,p,this.yvec[iy]) + (1-p.delta) * p.pension) for iy in 1:p.ny , ia in 1:p.na ] for id=1:p.nD) for it=1:(p.nT-1))

		# result arrays: matrices of type Envelope.
		this.v = [Envelope(MLine(fill(NaN,(p.na)),fill(NaN,(p.na)))) for id in 1:p.nD, it in 1:p.nT]
		this.c = [Envelope(MLine(fill(NaN,(p.na)),fill(NaN,(p.na)))) for id in 1:p.nD ,it in 1:p.nT]

		this.vdirty = [Envelope(MLine(fill(NaN,(p.na)),fill(NaN,(p.na)))) for id in 1:p.nD, it in 1:p.nT]
		this.cdirty = [Envelope(MLine(fill(NaN,(p.na)),fill(NaN,(p.na)))) for id in 1:p.nD ,it in 1:p.nT]

		return this
	end
end

"""
# General GModel

* endogenous retirement age choice: if work get income(age), if not get fixed pension
* with bequest
* `general` because there is state dependent income, i.e an AR1 structure
* no bankruptcy
"""
mutable struct GModel <: Model

	# nD is number of discrete choices: nD = 2

	# computation grids
	avec::Vector{Float64}  # each period has its own avec
	yvec::Vector{Float64}   # income support
	ywgt::Matrix{Float64}   # income weights

	# intermediate objects (na,ny,nD)
	m1::Dict{Int,Dict}	# a dict[it] for each period
	c1::Matrix{Float64}
	ev::Matrix{Float64}

	# result objects
	v :: Array{Envelope}  # arrays of Envelope objects
	c :: Array{Envelope}

	# simulation settings
	wshocks :: Matrix{Float64}  # wage shocks
	dshocks :: Matrix{Float64}  # discrete choice shocks
	w0shocks :: Vector{Float64}  # initial wealth shock
	y0shocks :: Vector{Int}  # initial income state

	function GModel(p::Param)

		this = new()

		# fedors version:
		# nodes,weights = quadpoints(p.ny,0,1)
		# N = Normal(0,1)
		# nodes = quantile.(N,nodes)
		# this.yvec = nodes * p.sigma
		# # this.ywgt = weights'
		# this.ywgt = reshape(repeat(weights,inner = p.ny),p.ny,p.ny)  # make a matrix

		# my version:
		# for y ~ N(mu,sigma), integrate y with
		# http://en.wikipedia.org/wiki/Gauss-Hermite_quadrature

		if p.ρ == 0
			nodes,weights = gausshermite(p.ny)  # from FastGaussQuadrature
			this.yvec = sqrt(2.0) * p.sigma .* nodes
			this.ywgt = reshape(repeat(weights .* pi^(-0.5),inner = p.ny),p.ny,p.ny)  # make a matrix
		else
			# version with income persistence
			this.yvec, this.ywgt = rouwenhorst(p.ρ,0,p.sigma,p.ny)
		end

		this.avec          = scaleGrid(p.a_low,p.a_high,p.na,logorder = 2)
		# this.avec          = [collect(range(p.a_low,stop = p.a_high,length = p.na))]
		# this.avec          = [scaleGrid(p.a_low,p.a_high,p.na,logorder = 1) for it in 1:p.nT-1]
		# this.avec          = [scaleGrid(η[it],p.a_high,p.na,logorder = 1) for it in 1:p.nT-1]
		# this.avec          = [scaleGrid(p.a_low,p.a_high,p.na,logorder = 0) for it in 1:p.nT-1]
		# push!(this.avec, scaleGrid(0.0,p.a_high,p.na,logorder = 0))  # last period
		# this.avec          = collect(range(p.a_low,stop = p.a_high,length = p.na))

		# precompute next period's cash on hand.
		# (na,ny,nD)
		# iD = 1: tomorrow work
		# iD = 2: tomorrow no work
		# this.m1 = Dict(it => Dict(id => Float64[this.avec[ia]*p.R .+ income(it,p,this.yvec[iy]) * (id==1) *(it < p.nT) for iy in 1:p.ny , ia in 1:p.na ] for id=1:p.nD) for it=2:p.nT)
		this.m1 = Dict(it => Dict(id => Float64[(this.avec[ia] * p.R)  .+ (income(it,p,this.yvec[iy]) * (id==1))  .+ (p.pension * (id==2)) for iy in 1:p.ny , ia in 1:p.na ] for id=1:p.nD) for it=2:p.nT)
		this.c1 = zeros(p.na,p.ny)
		this.ev = zeros(p.na,p.ny)

		this.v = [Envelope(MLine(fill(NaN,(p.na)),fill(NaN,(p.na)))) for id in 1:p.nD, iy in 1:p.ny, it in 1:p.nT]
		this.c = [Envelope(MLine(fill(NaN,(p.na)),fill(NaN,(p.na)))) for id in 1:p.nD, iy in 1:p.ny,it in 1:p.nT]
		# dchoice = [it => ["d" => zeros(Int,p.na), "Vzero" => 0.0, "threshold" => 0.0] for it=1:p.nT]

		return this
	end
end

function gmodel()
	p = Param()
	(GModel(p),p)
end

"""
GModel with bankruptcy
"""
mutable struct BModel <: Model
	avec::Vector{Float64}
	aposvec::Vector{Float64}
	yvec::Vector{Float64}   # income support
	ywgt::Matrix{Float64}   # income support

	# intermediate objects (na,nD)
	m1::Dict{Int,Dict}	# a dict[it] for each period
	# result objects
	v :: Array{Envelope}  # arrays of Envelope objects
	c :: Array{Envelope}
	vbk :: Array{Envelope}  # arrays of Envelope objects
	cbk :: Array{Envelope}
	iazero :: Int  # index of first non-negative asset state
	rgrid :: Array{Vector{Float64}}

	function BModel(p::Param)

		this = new()

		if p.ρ == 0
			nodes,weights = gausshermite(p.ny)  # from FastGaussQuadrature
			this.yvec = sqrt(2.0) * p.sigma .* nodes
			this.ywgt = reshape(repeat(weights .* pi^(-0.5),inner = p.ny),p.ny,p.ny)  # make a matrix
		else
			# version with income persistence
			this.yvec, this.ywgt = rouwenhorst(p.ρ,0,p.sigma,p.ny)
		end

		if p.a_low >= 0
			error("need a_low < 0 for bankruptcy model")
		end

		# simulation shocks

		this.avec          = collect(range(p.a_low,stop = p.a_high,length = p.na))
		this.aposvec       = collect(range(0.0,stop = p.a_high,length = p.na))
		this.iazero = findfirst(this.avec .>= 0)

		# precompute next period's cash on hand.
		# (na,ny,nState)
		# only precompute for bankrutpcy state
		# in non-bk state cash on hand is endogenous to changing interest rate.
		this.m1 = Dict(it => Float64[p.R * this.aposvec[ia] .+ income(it,p,this.yvec[iy]) for iy in 1:p.ny , ia in 1:p.na] for it=2:(p.nT))

		# result arrays: matrices of type Envelope.
		# this allocation is only to reserve about the right amount of memory. those will be overwritten in the algo.
		this.v = [Envelope(MLine(fill(NaN,(p.na)),fill(NaN,(p.na)))) for id in 1:p.nD, iy in 1:p.ny, it in 1:p.nT]
		this.c = [Envelope(MLine(fill(NaN,(p.na)),fill(NaN,(p.na)))) for id in 1:p.nD, iy in 1:p.ny,it in 1:p.nT]
		# for bankruptcy flag on there is no discrete choice
		this.vbk = [Envelope(MLine(fill(NaN,(p.na)),fill(NaN,(p.na)))) for iy in 1:p.ny, it in 1:p.nT]
		this.cbk = [Envelope(MLine(fill(NaN,(p.na)),fill(NaN,(p.na)))) for iy in 1:p.ny,it in 1:p.nT]

		this.rgrid = [fill(NaN,p.na) for id in 1:p.nD, iy in 1:p.ny,it in 1:p.nT ]

		return this
	end
end


text:
struct Uniform{T}
    data::T
end

Uniform(x::T) where {T <:ScalarTypes} = Uniform{T}(x)
Uniform(x::T) where {T <:AbstractArray} = Uniform{T}(x)
Uniform(x) = x

@spmd Uniform(x) = Uniform(x)


text:
using PSCOPF

using Test
using Dates

@testset "test_target_ts" begin
    ts1 = DateTime("2015-01-01T11:00:00")
    TS = PSCOPF.create_target_timepoints(ts1)

    EXPECTED_TS = [DateTime("2015-01-01T11:00:00"),
                    DateTime("2015-01-01T11:15:00"),
                    DateTime("2015-01-01T11:30:00"),
                    DateTime("2015-01-01T11:45:00")]
    @test length(TS) == length(EXPECTED_TS)
    @test TS == EXPECTED_TS
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
"""
A node type that applies the swing equation and voltage dynmaics
to the frequency/angle dynamics.
In the following, we followed the implementation of the Ex11tended Model/ third
order model including voltage dynamics according to K. Schmietendorf et. al.
"Self-Organized Synchronization and Voltage Stability in Networks of Synchronous
Machines", 2013.
The third order model is a special case of the fourth order model with Ed = 0
and x11d = x11d_dash.

In addition to the third order model the  Type IEEEG1 governor was implemented.
It is used for tandem compound, double reheat steam turbine systems.
The type was implemented according to "Dynamic Models for Turbine-Governors in
Power System Studies", IEEE PES, 2013, Page 2 - 2

Additionally to u, it has the internal dynamic variable ω representing the
frequencyof the rotator relative to the grid frequency Ω, i.e. the real
frequency ωr of the rotator is given as ωr=Ω+ω.

Keyword Arguments

Third Order Model
- H: inertia constant (given in [s]
- P: active (real) power output, also called the mechanical torque applied to the shaft, given in [pu]
- D: damping coefficient, (given in [s], see P. Sauer, eq. (5.156) where the damping torque is equal Dω)
- Ω: rated frequency in [1/s] of the power grid, often 2π⋅50Hz
- T_d_dash: time constant of d-ax11is, given in [s], see P. Sauer, chapter 3.7, p. 54 for a general ex11planation on time constants
- X_d_dash: transient reactance of d-axis, given in [pu]
- X_q_dash: transient reactance of q-axis, given in [pu]
- X_d: reactance of d-, given in [pu]
- X_d: reactance of q-ax11is, given in [pu]
- E_f: scaled field voltage, which, if set equal to 1.0 pu, gives 1.0 pu open-circuit terminal voltage. The physical device that provides the value of `E_f` is called the ex11citer (according to P. Sauer, p. 65)

Governer Type
- P0:
- Pmax: Max power limit imposed by valve or gate control [pu]
- Pmin: Min power limit imposed by valve or gate control [pu]
- Pup: Max main control valve rate of change [pu/s]
- Pdown: Min main control valve rate of change [pu/s]
- T1: Controller lag compensation [s]
- T2: Controller lead compensation [s]
- T3: Governor lag [s]
- K: Total effective speed-governing system gain [pu]
"""

using PowerDynamics
using NetworkDynamics
import PowerDynamics: dimension, symbolsof, construct_vertex
import PowerDynamics: AbstractNode
import Base: @__doc__

@DynamicNode ThirdOrderEqGovernorIEEEGone(H, D, Ω, E_f,T_d_dash, X_d_dash, X_q_dash, X_d, P0, Pmax, Pmin, Pup, Pdown, T1, T2, T3, K) begin
    @assert H > 0 "inertia (H) should be >0"
    @assert D >= 0 "damping (D) should be >=0"
    @assert T_d_dash > 0 "time constant of d-axis (T_d_dash) should be >0"
    @assert X_d_dash >= 0 "transient reactance of d-ax11is (x11_d_dash) should be >=0"
    @assert X_d >= 0 "reactance of d-ax11is (x11_d_dash) should be >=0"
    Ω_H = Ω / (2 * H)
end [[θ,dθ],[ω, dω], [Pm, dPm], [x1, dx1], [z, dz], [P, dP]] begin
    i_c = 1im * i * exp(-1im * θ)
    e_c = 1im * u * exp(-1im * θ)
    p = real(u * conj(i))
    e_d = 0                         # simpilfication from K.Schmietendorf paper
    e_q = imag(e_c)
    i_d = real(i_c)
    i_q = imag(i_c)

    # Governor Type IEEEG1
    dx1 = K * (-1 / T1 * x1 + (1 - T2 /  T1) * ω) # Block Input

    dP = (1 / T1) * x1 + (T2 / T1) * ω

    y = (1 / T3) * (P0 - P - Pm)                  # Block Output
    y_temp = y                                    # temorary variable

    # Limiting the valve rate of change
    if y > Pup
        y_temp = Pup
    end

    if y < Pdown
        y_temp = Pdown
    end

    dz = y_temp
    dPm = y_temp

    # Limiting the power imposed by the valve or gate control
    if z > Pmax
        dPm = (1 - Pmax) * dPm
    end

    if z < Pmin
        dPm = (1 - Pmin) * dPm
    end

    # Third Order Model
    dθ = ω
    dω = (Pm - D * ω - p - (X_q_dash - X_d_dash) * i_d * i_q) * Ω_H
    de_q = (1 / T_d_dash) * (E_f - e_q + i_d * (X_d - X_d_dash))
    # -> u = e_q * exp(1im * θ)
    du = de_q * exp(1im * θ) + u * 1im * ω
end

export ThirdOrderEqGovernorIEEEGone

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC