text:
# generated
SummonerSpells = Dict{Int,NamedTuple}(
    21 => (id = "SummonerBarrier", name = "방어막", description = "2초 동안 방어막으로 감싸 피해를 115~455(챔피언 레벨에 따라 변동)만큼 흡수합니다.", tooltip = "2초간 방어막이 챔피언을 감싸 {{ f1 }}의 피해를 흡수합니다.", maxrank = 1, cooldown = Any[180], cooldownBurn = "180", cost = Any[0], costBurn = "0", datavalues = NamedTuple(), effect = Any[nothing, Any[95], Any[20], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0]], effectBurn = Any[nothing, "95", "20", "0", "0", "0", "0", "0", "0", "0", "0"], vars = Any[], key = "21", summonerLevel = 4, modes = Any["NEXUSBLITZ", "URF", "PRACTICETOOL", "TUTORIAL", "CLASSIC", "ARAM", "DOOMBOTSTEEMO", "ULTBOOK", "ONEFORALL", "ARSR", "ASSASSINATE", "FIRSTBLOOD", "PROJECT", "STARGUARDIAN"], costType = "소모값 없음", maxammo = "-1", range = Any[1200], rangeBurn = "1200", image = (full = "SummonerBarrier.png", sprite = "spell0.png", group = "spell", x = 0, y = 0, w = 48, h = 48), resource = "소모값 없음"),
    1 => (id = "SummonerBoost", name = "정화", description = "챔피언에 걸린 모든 이동 불가와 (제압 및 공중으로 띄우는 효과 제외) 소환사 주문에 의한 해로운 효과를 제거하고 새로 적용되는 이동 불가 효과들의 지속시간을 3초간 65% 감소시킵니다.", tooltip = "챔피언에 적용되는 모든 방해 효과와 (제압 및 공중으로 띄우는 효과 제외) 해로운 소환사 주문 효과를 제거하고 다음 {{ f1 }}초 동안 방해 효과의 지속시간을 65% 감소시킵니다.", maxrank = 1, cooldown = Any[210], cooldownBurn = "210", cost = Any[0], costBurn = "0", datavalues = NamedTuple(), effect = Any[nothing, Any[0.65], Any[3], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0]], effectBurn = Any[nothing, "0.65", "3", "0", "0", "0", "0", "0", "0", "0", "0"], vars = Any[], key = "1", summonerLevel = 9, modes = Any["URF", "NEXUSBLITZ", "ARSR", "ONEFORALL", "ARAM", "CLASSIC", "PRACTICETOOL", "DOOMBOTSTEEMO", "TUTORIAL", "ULTBOOK", "FIRSTBLOOD", "PROJECT"], costType = "소모값 없음", maxammo = "-1", range = Any[200], rangeBurn = "200", image = (full = "SummonerBoost.png", sprite = "spell0.png", group = "spell", x = 48, y = 0, w = 48, h = 48), resource = "소모값 없음"),
    14 => (id = "SummonerDot", name = "점화", description = "적 챔피언을 불태워 5초 동안 70~410의 고정 피해(챔피언 레벨에 따라 변동)를 입히고 모습을 드러내며 치료 효과를 감소시킵니다.", tooltip = "대상 적 챔피언을 불태워, 5초에 걸쳐 {{ f1 }}의 고정 피해를 입히고 그동안 적의 위치를 드러내며 {{ grievousamount*100 }}%의 고통스러운 상처를 적용합니다.

<rules>(고통스러운 상처는 치유 및 회복 효과를 가장 높은 수치만큼 감소시킵니다. 또한, 은신 중인 적은 위치가 드러나지 않습니다.)</rules>", maxrank = 1, cooldown = Any[180], cooldownBurn = "180", cost = Any[0], costBurn = "0", datavalues = NamedTuple(), effect = Any[nothing, Any[5], Any[10], Any[4], Any[100], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0]], effectBurn = Any[nothing, "5", "10", "4", "100", "0", "0", "0", "0", "0", "0"], vars = Any[], key = "14", summonerLevel = 9, modes = Any["NEXUSBLITZ", "URF", "PRACTICETOOL", "TUTORIAL", "CLASSIC", "ARAM", "DOOMBOTSTEEMO", "ULTBOOK", "ONEFORALL", "ARSR", "ASSASSINATE", "FIRSTBLOOD"], costType = "소모값 없음", maxammo = "-1", range = Any[600], rangeBurn = "600", image = (full = "SummonerDot.png", sprite = "spell0.png", group = "spell", x = 96, y = 0, w = 48, h = 48), resource = "소모값 없음"),
    3 => (id = "SummonerExhaust", name = "탈진", description = "적 챔피언을 지치게 만들어 3초 동안 이동 속도를 30% 느리게 하고 이 동안 가하는 피해량을 40% 낮춥니다.", tooltip = "적 챔피언을 지치게 만들어 3초 동안 이동 속도를 {{ f3 }}% 느리게 하고 이 동안 가하는 피해량을 {{ f2 }}% 낮춥니다.", maxrank = 1, cooldown = Any[210], cooldownBurn = "210", cost = Any[0], costBurn = "0", datavalues = NamedTuple(), effect = Any[nothing, Any[3], Any[40], Any[0], Any[0], Any[30], Any[0], Any[0], Any[0], Any[0], Any[0]], effectBurn = Any[nothing, "3", "40", "0", "0", "30", "0", "0", "0", "0", "0"], vars = Any[], key = "3", summonerLevel = 4, modes = Any["NEXUSBLITZ", "URF", "PRACTICETOOL", "TUTORIAL", "CLASSIC", "ARAM", "DOOMBOTSTEEMO", "ULTBOOK", "ONEFORALL", "ARSR"], costType = "소모값 없음", maxammo = "-1", range = Any[650], rangeBurn = "650", image = (full = "SummonerExhaust.png", sprite = "spell0.png", group = "spell", x = 144, y = 0, w = 48, h = 48), resource = "소모값 없음"),
    4 => (id = "SummonerFlash", name = "점멸", description = "커서 방향으로 챔피언이 짧은 거리를 순간이동합니다.", tooltip = "커서 방향으로 챔피언이 짧은 거리를 순간이동합니다.", maxrank = 1, cooldown = Any[300], cooldownBurn = "300", cost = Any[0], costBurn = "0", datavalues = NamedTuple(), effect = Any[nothing, Any[400], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0]], effectBurn = Any[nothing, "400", "0", "0", "0", "0", "0", "0", "0", "0", "0"], vars = Any[], key = "4", summonerLevel = 7, modes = Any["NEXUSBLITZ", "URF", "PRACTICETOOL", "SNOWURF", "TUTORIAL", "CLASSIC", "ARAM", "DOOMBOTSTEEMO", "ULTBOOK", "ONEFORALL", "ARSR", "ASSASSINATE", "FIRSTBLOOD", "PROJECT", "STARGUARDIAN"], costType = "소모값 없음", maxammo = "-1", range = Any[425], rangeBurn = "425", image = (full = "SummonerFlash.png", sprite = "spell0.png", group = "spell", x = 192, y = 0, w = 48, h = 48), resource = "소모값 없음"),
    6 => (id = "SummonerHaste", name = "유체화", description = "챔피언이 10초 동안 유닛과 충돌하지 않게 되며 챔피언 레벨에 따라 이동 속도가 24~48% 증가합니다. 처치 관여 시 유체화 지속시간이 늘어납니다.", tooltip = "챔피언이 {{ duration }}초 동안 유닛과 충돌하지 않게 되며 챔피언 레벨에 따라 이동 속도가 {{ movespeedmod }} 증가합니다.

처치 관여 시 유체화 지속시간이 {{ takedownextension }}초 늘어납니다.", maxrank = 1, cooldown = Any[210], cooldownBurn = "210", cost = Any[0], costBurn = "0", datavalues = NamedTuple(), effect = Any[nothing, Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0]], effectBurn = Any[nothing, "0", "0", "0", "0", "0", "0", "0", "0", "0", "0"], vars = Any[], key = "6", summonerLevel = 1, modes = Any["NEXUSBLITZ", "TUTORIAL_MODULE_2", "URF", "PRACTICETOOL", "TUTORIAL", "CLASSIC", "ARAM", "DOOMBOTSTEEMO", "ULTBOOK", "TUTORIAL_MODULE_1", "ONEFORALL", "ARSR", "ASSASSINATE", "FIRSTBLOOD", "PROJECT", "STARGUARDIAN"], costType = "소모값 없음", maxammo = "-1", range = Any[200], rangeBurn = "200", image = (full = "SummonerHaste.png", sprite = "spell0.png", group = "spell", x = 240, y = 0, w = 48, h = 48), resource = "소모값 없음"),
    7 => (id = "SummonerHeal", name = "회복", description = "자신과 대상 아군 챔피언의 체력을 90~345만큼 (챔피언 레벨에 따라 변동) 회복시키고 1초 동안 이동 속도가 30% 증가합니다. 최근 소환사 주문 회복의 영향을 받은 유닛의 경우 치유량이 절반만 적용됩니다.", tooltip = "자신과 대상 아군 챔피언의 체력을 {{ f1 }} 회복하고 1초 동안 이동 속도를 30% 올려줍니다. 최근 소환사 주문 회복의 영향을 받은 유닛의 경우 치유량이 절반만 적용됩니다.

이 주문은 대상을 정하지 않은 경우, 범위 내에서 가장 큰 부상을 입은 아군 챔피언에게 시전됩니다.", maxrank = 1, cooldown = Any[240], cooldownBurn = "240", cost = Any[0], costBurn = "0", datavalues = NamedTuple(), effect = Any[nothing, Any[0.3], Any[75], Any[15], Any[0.5], Any[826], Any[0.5], Any[0], Any[0], Any[0], Any[0]], effectBurn = Any[nothing, "0.3", "75", "15", "0.5", "826", "0.5", "0", "0", "0", "0"], vars = Any[], key = "7", summonerLevel = 1, modes = Any["NEXUSBLITZ", "TUTORIAL_MODULE_2", "URF", "PRACTICETOOL", "TUTORIAL", "CLASSIC", "ARAM", "DOOMBOTSTEEMO", "ULTBOOK", "TUTORIAL_MODULE_1", "ONEFORALL", "ARSR", "ASSASSINATE", "PROJECT", "STARGUARDIAN"], costType = "소모값 없음", maxammo = "-1", range = Any[850], rangeBurn = "850", image = (full = "SummonerHeal.png", sprite = "spell0.png", group = "spell", x = 288, y = 0, w = 48, h = 48), resource = "소모값 없음"),
    13 => (id = "SummonerMana", name = "총명", description = "최대 마나량의 50%를 회복합니다. 주변 아군도 최대 마나량의 25%가 회복됩니다.", tooltip = "최대 마나의 {{ f1 }}%를 회복하고 주변 아군의 마나를 최대 마나의 {{ f2 }}%만큼 회복시켜줍니다.", maxrank = 1, cooldown = Any[240], cooldownBurn = "240", cost = Any[0], costBurn = "0", datavalues = NamedTuple(), effect = Any[nothing, Any[50], Any[25], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0]], effectBurn = Any[nothing, "50", "25", "0", "0", "0", "0", "0", "0", "0", "0"], vars = Any[], key = "13", summonerLevel = 6, modes = Any["ARAM", "FIRSTBLOOD"], costType = "소모값 없음", maxammo = "-1", range = Any[600], rangeBurn = "600", image = (full = "SummonerMana.png", sprite = "spell0.png", group = "spell", x = 336, y = 0, w = 48, h = 48), resource = "소모값 없음"),
    30 => (id = "SummonerPoroRecall", name = "왕을 향해!", description = "포로 왕의 곁으로 빠르게 이동합니다.", tooltip = "기본 지속 효과: 포로로 적 챔피언을 맞히면 아군에게 포로 표식이 생깁니다. 포로 표식이 10개가 되면 아군은 포로 왕을 소환하여 함께 싸울 수 있습니다. 포로 왕이 활성화된 동안에는 어느 팀도 포로 표식을 획득할 수 없습니다.

사용 시: 포로 왕 쪽으로 빠르게 질주합니다. 아군이 포로 왕을 소환한 동안에만 사용할 수 있습니다.

''포로들이 줄을 당깁니다. 여러분은 몸을 맡기기만 하면 됩니다.''</mainText>", maxrank = 1, cooldown = Any[10], cooldownBurn = "10", cost = Any[0], costBurn = "0", datavalues = NamedTuple(), effect = Any[nothing, Any[3000], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0]], effectBurn = Any[nothing, "3000", "0", "0", "0", "0", "0", "0", "0", "0", "0"], vars = Any[], key = "30", summonerLevel = 1, modes = Any["KINGPORO"], costType = "소모값 없음", maxammo = "-1", range = Any[200], rangeBurn = "200", image = (full = "SummonerPoroRecall.png", sprite = "spell0.png", group = "spell", x = 384, y = 0, w = 48, h = 48), resource = "소모값 없음"),
    31 => (id = "SummonerPoroThrow", name = "포로 던지기", description = "포로를 적에게 던집니다. 적이 맞으면 해당 적에게 빠르게 이동할 수 있습니다.", tooltip = "포로를 멀리 던져서 첫 번째로 맞는 적에게 {{ f2 }}의 고정 피해를 입히고 대상에 대한 절대 시야를 얻습니다.

적을 맞히면 3초 안에 이 스킬을 다시 사용하여 맞힌 대상에게 질주할 수 있습니다. 대상에게 질주해 가면 {{ f2 }}의 추가 고정 피해를 입히고 포로 던지기의 재사용 대기시간이 {{ e4 }}초 감소합니다.

포로는 마법이 아니라 동물이므로, 주문 보호막이나 바람 장벽으로 막을 수 없습니다!

''룬테라 공기 역학의 결정판, 포로!''</mainText>", maxrank = 1, cooldown = Any[20], cooldownBurn = "20", cost = Any[0], costBurn = "0", datavalues = NamedTuple(), effect = Any[nothing, Any[20], Any[10], Any[3], Any[5], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0]], effectBurn = Any[nothing, "20", "10", "3", "5", "0", "0", "0", "0", "0", "0"], vars = Any[], key = "31", summonerLevel = 1, modes = Any["KINGPORO"], costType = "소모값 없음", maxammo = "-1", range = Any[2500], rangeBurn = "2500", image = (full = "SummonerPoroThrow.png", sprite = "spell0.png", group = "spell", x = 432, y = 0, w = 48, h = 48), resource = "소모값 없음"),
    11 => (id = "SummonerSmite", name = "강타", description = "대상 에픽 및 대형/중형 몬스터 또는 미니언에게 @SmiteBaseDamage@의 고정 피해를 입힙니다.", tooltip = "대상 몬스터 또는 미니언에게 {{ f1 }}의 고정 피해를 입힙니다. 몬스터에게 사용하면 체력도 {{ f6 }}(+{{ f7 }})만큼 회복됩니다.

강타는 {{ ammorechargetime }}초에 한 번씩 충전되며, 최대 2회까지 충전됩니다.

<rules>작은 몬스터에게는 사용할 수 없습니다.</rules>", maxrank = 1, cooldown = Any[15], cooldownBurn = "15", cost = Any[0], costBurn = "0", datavalues = NamedTuple(), effect = Any[nothing, Any[15], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0]], effectBurn = Any[nothing, "15", "0", "0", "0", "0", "0", "0", "0", "0", "0"], vars = Any[], key = "11", summonerLevel = 9, modes = Any["NEXUSBLITZ", "URF", "PRACTICETOOL", "TUTORIAL", "CLASSIC", "DOOMBOTSTEEMO", "ULTBOOK", "ONEFORALL", "ARSR"], costType = "소모값 없음", maxammo = "2", range = Any[500], rangeBurn = "500", image = (full = "SummonerSmite.png", sprite = "spell0.png", group = "spell", x = 0, y = 48, w = 48, h = 48), resource = "소모값 없음"),
    39 => (id = "SummonerSnowURFSnowball_Mark", name = "표식", description = "적을 향해 직선으로 눈덩이를 던집니다. 눈덩이가 적을 맞히면 표식이 남아 절대 시야가 생기고, 표식이 붙은 대상을 향해 챔피언이 빠르게 이동할 수 있습니다.", tooltip = "멀리 눈덩이를 던져 첫 번째 맞은 적에게 {{ f1 }}의 고정 피해를 입히고 대상에 대한 절대 시야를 얻습니다. 적을 맞히면 이 스킬을 {{ e3 }}초 내에 재사용하여 표식이 남은 유닛에게 돌진하여 추가로 {{ f1 }}의 고정 피해를 입힙니다. 대상에게 돌진하면 표식의 재사용 대기시간이 {{ e4 }}% 감소합니다.

표식의 투사체는 주문 보호막이나 투사체 약화 효과로 막을 수 없습니다.", maxrank = 1, cooldown = Any[80], cooldownBurn = "80", cost = Any[0], costBurn = "0", datavalues = NamedTuple(), effect = Any[nothing, Any[10], Any[5], Any[3], Any[25], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0]], effectBurn = Any[nothing, "10", "5", "3", "25", "0", "0", "0", "0", "0", "0"], vars = Any[], key = "39", summonerLevel = 6, modes = Any["SNOWURF"], costType = "소모값 없음", maxammo = "-1", range = Any[8000], rangeBurn = "8000", image = (full = "SummonerSnowURFSnowball_Mark.png", sprite = "spell0.png", group = "spell", x = 48, y = 48, w = 48, h = 48), resource = "소모값 없음"),
    32 => (id = "SummonerSnowball", name = "표식", description = "적을 향해 직선으로 눈덩이를 던집니다. 눈덩이가 적을 맞히면 표식이 남아 절대 시야가 생기고, 표식이 붙은 대상을 향해 챔피언이 빠르게 이동할 수 있습니다.", tooltip = "멀리 눈덩이를 던져 첫 번째 맞은 적에게 {{ f1 }}의 고정 피해를 입히고 대상에 대한 절대 시야를 얻습니다. 적을 맞히면 이 스킬을 {{ e3 }}초 내에 재사용하여 표식이 남은 유닛에게 돌진하여 추가로 {{ f1 }}의 고정 피해를 입힙니다. 대상에게 돌진하면 표식의 재사용 대기시간이 {{ e4 }}% 감소합니다.

표식의 투사체는 주문 보호막이나 투사체 약화 효과로 막을 수 없습니다.", maxrank = 1, cooldown = Any[80], cooldownBurn = "80", cost = Any[0], costBurn = "0", datavalues = NamedTuple(), effect = Any[nothing, Any[10], Any[5], Any[3], Any[25], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0]], effectBurn = Any[nothing, "10", "5", "3", "25", "0", "0", "0", "0", "0", "0"], vars = Any[], key = "32", summonerLevel = 6, modes = Any["ARAM", "FIRSTBLOOD"], costType = "소모값 없음", maxammo = "-1", range = Any[1600], rangeBurn = "1600", image = (full = "SummonerSnowball.png", sprite = "spell0.png", group = "spell", x = 96, y = 48, w = 48, h = 48), resource = "소모값 없음"),
    12 => (id = "SummonerTeleport", name = "순간이동", description = "4초 동안 정신을 집중한 다음, 챔피언이 지정한 아군 구조물, 미니언, 혹은 와드로 순간이동하고 이동 속도가 증가합니다. 순간이동의 재사용 대기시간은 챔피언 레벨에 따라 420~210초입니다.", tooltip = "{{ channelduration }}초 동안 정신 집중을 한 후 아군 구조물, 미니언, 혹은 와드로 순간이동하고 {{ msduration }}초 동안 이동 속도가 {{ msamount*100 }}% 증가합니다.", maxrank = 1, cooldown = Any[0], cooldownBurn = "0", cost = Any[0], costBurn = "0", datavalues = NamedTuple(), effect = Any[nothing, Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0]], effectBurn = Any[nothing, "0", "0", "0", "0", "0", "0", "0", "0", "0", "0"], vars = Any[], key = "12", summonerLevel = 7, modes = Any["ONEFORALL", "PRACTICETOOL", "CLASSIC", "TUTORIAL", "ULTBOOK", "ARSR", "ASSASSINATE", "DOOMBOTSTEEMO"], costType = "소모값 없음", maxammo = "-1", range = Any[25000], rangeBurn = "25000", image = (full = "SummonerTeleport.png", sprite = "spell0.png", group = "spell", x = 144, y = 48, w = 48, h = 48), resource = "소모값 없음"),
    54 => (id = "Summoner_UltBook_Placeholder", name = "게임 시작 후 결정", description = "", tooltip = "선택한 궁극기 소환사 주문으로 대체됩니다.", maxrank = 1, cooldown = Any[0], cooldownBurn = "0", cost = Any[0], costBurn = "0", datavalues = NamedTuple(), effect = Any[nothing, Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0], Any[0]], effectBurn = Any[nothing, "0", "0", "0", "0", "0", "0", "0", "0", "0", "0"], vars = Any[], key = "54", summonerLevel = 1, modes = Any["ULTBOOK"], costType = " ", maxammo = "-1", range = Any[400], rangeBurn = "400", image = (full = "Summoner_UltBook_Placeholder.png", sprite = "spell0.png", group = "spell", x = 192, y = 48, w = 48, h = 48), resource = " "),
)


text:
struct RouletteWheelSelection end
# selection(pop::Vector{<:GACreature}, n::Integer, rng)
function selection(::RouletteWheelSelection,
                   pop::Vector{<:GACreature}, n::Integer, rng=Base.GLOBAL_RNG)
    wmin,wmax = extrema(fitness(c) for c in pop)
    weight = wmax - wmin
    function stochasticpick()
        while true
            i = rand(rng,1:length(pop))
            if wmin + weight * rand(rng,typeof(weight)) <= fitness(pop[i])
                return i
            end
        end
    end
    parents = Vector{Tuple{Int,Int}}(undef, n)
    for k = 1:n
        i = stochasticpick()
        j = i
        while i==j
            j = stochasticpick()
        end
        parents[k] = (i,j)
    end
    parents
end

struct TournamentSelection
    k::Int # if k=2 then binary tournament selection
    TournamentSelection(k=2) = new(k)
end
function selection(sel::TournamentSelection,
                   pop::Vector{<:GACreature}, n::Integer, rng=Base.GLOBAL_RNG)
    function stochasticpick()
        si = rand(rng,1:length(pop))
        weighti = fitness(pop[si])
        for _ = 1:sel.k-1
            sj = rand(rng,1:length(pop))
            weightj = fitness(pop[sj])
            if weightj > weighti
                si = sj
                weighti = weightj
            end
        end
        si
    end
    parents = Vector{Tuple{Int,Int}}(undef, n)
    for k = 1:n
        i = stochasticpick()
        j = i
        while i==j
            j = stochasticpick()
        end
        parents[k] = (i,j)
    end
    parents
end

ext: jl


text:
function GenParam(generators, setI, setI_hat, wind_LBset, wind_OSset, pvset, pv_BTMset, hydroset, nucset, coalset, oilset, gasset, σ0_pv, σ0_wind, υ0)
    filename_windcurve = string(pwd(),"/data/hourly_wind_gen_2018.csv")
    filename_pvcurve = string(pwd(),"/data/hourly_solar_gen_2018.csv")
    # pvcurve = GenCurveLoad(filename_pvcurve)[2:31:end,:]/2883.8  # pick one day a month over the year / nameplate capacity = 2883MW
    # windcurve = GenCurveLoad(filename_windcurve)[2:31:end,:]/1300 # pick one day a month over the year / nameplate capacity = 1300MW

    WTcf = 0.29
    PVcf = 0.14
    # https://www.nyiso.com/power-trends

    temp_pv = GenCurveLoad(filename_pvcurve)
    pvcurve_year = temp_pv / mean(temp_pv) * PVcf
    temp_wind = GenCurveLoad(filename_windcurve)
    windcurve_year = temp_wind / mean(temp_wind) * WTcf

    pvcurve = sum(pvcurve_year, dims=1)/365
    windcurve = sum(windcurve_year, dims=1)/365

    # DCRF = 6.3281e-4 # h=5years / 5%
    # DCRF = 3.5481e-4 * length(setT) / 24 # h=10years / 5%
    # DCRF = 2.1984e-4 # h=20years / 5%
    # DCRF = 3.05e-4 # h=10years / 2%
    # DCRF = 1.6755e-4 # h=20years / 2%

    Cinv = Dict()
    # Investment cost source:  https://www.eia.gov/electricity/generatorcosts/
    # PV & WT ref: https://www.eia.gov/outlooks/aeo/section_issue_renewables.php

    Cinv["Wind_LB"] = DCRF * 1260 * 1e3
    Cinv["Wind_OS"] = DCRF * 5446 * 1e3
    Cinv["PV"] = DCRF * 1307 * 1e3
    Cinv["PV_BTM"] = DCRF * 2019 * 1e3
    # Cinv["PV_BTM"] = DCRF * 1539 * 1e3
    # Cinv["PV_BTM"] = DCRF * 1307 * 1e3
    Cinv["Hydro"] = DCRF * 5312 * 1e3
    Cinv["Nuclear"] = DCRF * 1000 * 1e3
    Cinv["Coal"] = DCRF * 1000 * 1e3
    Cinv["Oil"] = DCRF * 1672 * 1e3
    # Cinv["Gas"] = DCRF * 1563 * 1e3
    Cinv["Gas"] = DCRF * 895 * 1e3

    Coper = Dict()
    Coper["Wind_LB"] = [0, 1.1, 0]
    Coper["Wind_OS"] = [0, 1.1, 0]
    Coper["PV"] = [0, 0.4, 0]
    Coper["PV_BTM"] = [0, 0.8, 0]
    # Coper["PV_BTM"] = [0, 0, 0]
    Coper["Hydro"] = [0, 4, 0]
    Coper["Nuclear"] = [0, 7, 1250]
    Coper["Coal"] = [0.0, 19.98, 1535.0]
    Coper["Oil"] = [0.0, 192.06, 5000.0]
    Coper["Gas"] = [0.0, 44.92, 612.84] # 20

    Cinvg = Vector{Float64}(undef, length(setI))
    Coperg = Vector{Vector{Float64}}(undef, length(setI))
    for g in setI_hat
        if g in wind_LBset
            Cinvg[g] = Cinv["Wind_LB"]
            Coperg[g] = Coper["Wind_LB"]
        elseif g in wind_OSset
            Cinvg[g] = Cinv["Wind_OS"]
            Coperg[g] = Coper["Wind_OS"]
        elseif g in pvset
            Cinvg[g] = Cinv["PV"]
            Coperg[g] = Coper["PV"]
        elseif g in pv_BTMset
            Cinvg[g] = Cinv["PV_BTM"]
            Coperg[g] = Coper["PV_BTM"]
        elseif g in hydroset
            Cinvg[g] = Cinv["Hydro"]
            Coperg[g] = Coper["Hydro"]
        elseif g in nucset
            Cinvg[g] = Cinv["Nuclear"]
            Coperg[g] = Coper["Nuclear"]
        elseif g in coalset
            Cinvg[g] = Cinv["Coal"]
            Coperg[g] = Coper["Coal"]
        elseif g in oilset
            Cinvg[g] = Cinv["Oil"]
            Coperg[g] = Coper["Oil"]
        elseif g in gasset
            Cinvg[g] = Cinv["Gas"]
            Coperg[g] = Coper["Gas"]
        end
    end

    for g in setdiff(setI, setI_hat)
        if g in wind_LBset
            Cinvg[g] = Cinv["Wind_LB"]
            Coperg[g] = generators[g].cost
        elseif g in wind_OSset
            Cinvg[g] = Cinv["Wind_OS"]
            Coperg[g] = generators[g].cost
        elseif g in pvset
            Cinvg[g] = Cinv["PV"]
            Coperg[g] = generators[g].cost
        elseif g in pv_BTMset
            Cinvg[g] = Cinv["PV_BTM"]
            Coperg[g] = generators[g].cost
        elseif g in hydroset
            Cinvg[g] = Cinv["Hydro"]
            Coperg[g] = generators[g].cost
        elseif g in nucset
            Cinvg[g] = Cinv["Nuclear"]
            Coperg[g] = generators[g].cost
        elseif g in coalset
            Cinvg[g] = Cinv["Coal"]
            Coperg[g] = generators[g].cost
        elseif g in oilset
            Cinvg[g] = Cinv["Oil"]
            Coperg[g] = generators[g].cost
        elseif g in gasset
            Cinvg[g] = Cinv["Gas"]
            Coperg[g] = generators[g].cost
        else
            Coperg[g] = generators[g].cost
        end
    end

    βFC = 1.0*ones(length(setI),setT[end])
    for g in union(pvset, pv_BTMset)
        for t in setT
            βFC[g,t] = pvcurve[1,t]
        end
    end
    for g in union(wind_LBset, wind_OSset)
        for t in setT
            βFC[g,t] = windcurve[1,t]
        end
    end

    σ = zeros(length(setI))
    for g in union(pvset, pv_BTMset)
        σ[g] = σ0_pv
    end
    for g in union(wind_LBset, wind_OSset)
        σ[g] = σ0_wind
    end

    υ = υ0*ones(length(setI), setT[end])
    # kappa = [1; 0.252; 0.75; 0.411; 0.44; 0.385]
    # PolicyB = 1e8  * ones(length(stateset)) # Policy Budget
    return Cinvg, Coperg, βFC, σ, υ
end

max_stars_repo_name: JipKim/IA

max_issues_repo_name: JipKim/IA

max_forks_repo_name: JipKim/IA


text:
# Projected Alternating Least Squared method
#
#  Solve the following problem via alternate updating:
#
#     (1/2) * ||X - WH||^2
#   + (lambda_w/2) * ||W||^2
#   + (lambda_h/2) * ||H||^2
#
#  At each iteration, the algorithm updates H and W, holding
#  the other fixed. Particularly, it obtains H and W by solving
#  an unconstrained least square problem and then casting
#  negative entries back to zeros.
#

using LinearAlgebra: copytri!
using LinearAlgebra.LAPACK: potrs!, potrf!, potri!

mutable struct ProjectedALS{T}
    maxiter::Int            # maximum number of iterations (in main procedure)
    verbose::Bool           # whether to show procedural information
    tol::T                  # tolerance of changes on W and H upon convergence
    update_H::Bool          # whether to update H
    lambda_w::T             # L2 regularization coefficient for W
    lambda_h::T             # L2 regularization coefficient for H

    function ProjectedALS{T}(;maxiter::Integer=100,
                              verbose::Bool=false,
                              tol::Real=cbrt(eps(T)),
                              update_H::Bool=true,
                              lambda_w::Real=cbrt(eps(T)),
                              lambda_h::Real=cbrt(eps(T))) where T

        new{T}(maxiter, verbose, tol, update_H, lambda_w, lambda_h)
    end
end

solve!(alg::ProjectedALS, X, W, H) =
    nmf_skeleton!(ProjectedALSUpd(alg.update_H, alg.lambda_w, alg.lambda_h),
                  X, W, H, alg.maxiter, alg.verbose, alg.tol)

struct ProjectedALSUpd{T} <: NMFUpdater{T}
    update_H::Bool
    lambda_w::T
    lambda_h::T
end

struct ProjectedALSUpd_State{T}
    WH::Matrix{T}
    WtW::Matrix{T}
    HHt::Matrix{T}
    XHt::Matrix{T}

    function ProjectedALSUpd_State{T}(X, W::Matrix{T}, H::Matrix{T}) where T
        p, n, k = nmf_checksize(X, W, H)
        new{T}(W * H,
               Matrix{T}(undef, k, k),
               Matrix{T}(undef, k, k),
               Matrix{T}(undef, p, k))
    end
end

prepare_state(::ProjectedALSUpd{T}, X, W, H) where {T} = ProjectedALSUpd_State{T}(X, W, H)

function evaluate_objv(u::ProjectedALSUpd{T}, s::ProjectedALSUpd_State{T}, X, W, H) where T
    r = convert(T, 0.5) * sqL2dist(X, s.WH)
    if u.lambda_w > 0
        r += (convert(T, 0.5) * u.lambda_w) * abs2(norm(W))
    end
    if u.lambda_h > 0
        r += (convert(T, 0.5) * u.lambda_h) * abs2(norm(H))
    end
    return r
end

function update_wh!(upd::ProjectedALSUpd{T}, s::ProjectedALSUpd_State{T},
                    X,
                    W::Matrix{T},
                    H::Matrix{T}) where T

    # fields
    WH = s.WH
    WtW = s.WtW
    HHt = s.HHt
    XHt = s.XHt
    lambda_w = upd.lambda_w
    lambda_h = upd.lambda_h

    # update H
    if upd.update_H
        Wt = transpose(W)
        adddiag!(mul!(WtW, Wt, W), lambda_h)
        mul!(H, Wt, X)       # H <- W'X
        pdsolve!(WtW, H)     # H <- inv(WtW) * H
        projectnn!(H)
    end

    # update W
    Ht = transpose(H)
    adddiag!(mul!(HHt, H, Ht), lambda_w)
    mul!(XHt, X, Ht)    # XHt <- XH'
    pdrsolve!(XHt, HHt, W)  # W <- XHt * inv(HHt)
    projectnn!(W)

    # update WH
    mul!(WH, W, H)
end


text:
# This file contains routines to generate MDF files
export saveasMDF, loadDataset, loadMetadata, loadMetadataOnline, setparam!, compressCalibMDF
export saveasMDFHacking # temporary Hack

function setparam!(params::Dict, parameter, value)
  if value != nothing
    params[parameter] = value
  end
end

# we do not support all conversion possibilities
function loadDataset(f::MPIFile; frames=1:acqNumFrames(f), applyCalibPostprocessing=false)
  params = loadMetadata(f)

  # call API function and store result in a parameter Dict
  if experimentHasMeasurement(f)
    loadMeasParams(f, params, skipMeasData = true)
    if !applyCalibPostprocessing
      if frames!=1:acqNumFrames(f)
        setparam!(params, "measData", measData(f,frames))
        setparam!(params, "acqNumFrames", length(frames))
        setparam!(params, "measIsBGFrame", measIsBGFrame(f)[frames])
      else
        setparam!(params, "measData", measData(f))
      end
    else
        @info "load measurement data"
        data = getMeasurementsFD(f, false, frames=1:acqNumFrames(f), sortFrames=true,
              spectralLeakageCorrection=false, transposed=true, tfCorrection=false)
        setparam!(params, "measData", data)
        setparam!(params, "measIsFourierTransformed", true)
        setparam!(params, "measIsFastFrameAxis", true)
        setparam!(params, "measIsFramePermutation", true)
        setparam!(params, "measFramePermutation", fullFramePermutation(f))

        setparam!(params, "measIsBGFrame",
          cat(zeros(Bool,acqNumFGFrames(f)),ones(Bool,acqNumBGFrames(f)), dims=1))

        snr = calibSNR(f)
	if snr == nothing
          @info "calculate SNR"
          snr = calculateSystemMatrixSNR(f, data)
        end
        setparam!(params, "calibSNR", snr)
    end
  end

  loadCalibParams(f, params)
  loadRecoParams(f, params)

  return params
end

const defaultParams =[:version, :uuid, :time, :dfStrength, :acqGradient, :studyName, :studyNumber, :studyUuid, :studyTime, :studyDescription,
          :experimentName, :experimentNumber, :experimentUuid, :experimentDescription,
          :experimentSubject,
          :experimentIsSimulation, :experimentIsCalibration,
          :tracerName, :tracerBatch, :tracerVendor, :tracerVolume, :tracerConcentration,
          :tracerSolute, :tracerInjectionTime,
          :scannerFacility, :scannerOperator, :scannerManufacturer, :scannerName,
          :scannerTopology, :acqNumPeriodsPerFrame, :acqNumAverages,
          :acqStartTime, :acqOffsetField, :acqNumFrames,
          :dfNumChannels, :dfPhase, :dfBaseFrequency, :dfDivider,
          :dfCycle, :dfWaveform, :rxNumChannels, :rxBandwidth,
          :rxNumSamplingPoints, :rxTransferFunction, :rxInductionFactor,
          :rxUnit, :rxDataConversionFactor]

function loadMetadata(f, inputParams = MPIFiles.defaultParams)
  params = Dict{String,Any}()
  # call API function and store result in a parameter Dict
  for op in inputParams
    setparam!(params, string(op), eval(op)(f))
  end
  return params
end

function loadRecoParams(f, params = Dict{String,Any}())
  if experimentHasReconstruction(f)
    for op in [:recoData, :recoSize, :recoFov, :recoFovCenter, :recoOrder,
           :recoPositions, :recoParameters]
      setparam!(params, string(op), eval(op)(f))
    end
  end

  return params
end

function loadCalibParams(f, params = Dict{String,Any}())
  if experimentIsCalibration(f)
    for op in [:calibFov, :calibFovCenter,
               :calibSize, :calibOrder, :calibPositions, :calibOffsetField,
             :calibDeltaSampleSize, :calibMethod]
      setparam!(params, string(op), eval(op)(f))
    end
    if !haskey(params, "calibSNR")
      setparam!(params, "calibSNR", calibSNR(f))
    end
  end
  return params
end

function loadMeasParams(f, params = Dict{String,Any}(); skipMeasData = false)
  if experimentHasMeasurement(f)
    for op in [:measIsFourierTransformed, :measIsTFCorrected,
                 :measIsBGCorrected,
                 :measIsFastFrameAxis, :measIsFramePermutation, :measIsFrequencySelection,
                 :measIsSpectralLeakageCorrected,
                 :measFramePermutation, :measIsBGFrame]
      setparam!(params, string(op), eval(op)(f))
    end
  end

  if !skipMeasData
    setparam!(params, "measData", measData(f))
  end

  return params
end

function appendBGDataset(params::Dict, filenameBG::String; kargs...)
  fBG = MPIFile(filenameBG)
  return appendBGDataset(params, fBG; kargs...)
end

function appendBGDataset(params::Dict, fBG::MPIFile; frames=1:acqNumFrames(fBG))
  paramsBG = loadDataset(fBG, frames=frames)
  paramsBG["measIsBGFrame"][:] = true

  params["measData"] = cat(4, params["measData"], paramsBG["measData"])
  params["measIsBGFrame"] = cat(params["measIsBGFrame"], paramsBG["measIsBGFrame"], dims=1)
  params["acqNumFrames"] += paramsBG["acqNumFrames"]

  return params
end

isConvertibleToMDF(f::MPIFile) = true
function isConvertibleToMDF(f::BrukerFile)
  # check if raw data is consistent
  if !rawDataLengthConsistent(f::BrukerFile)
    return false
  end
  # check if conversion of method is supported
  #TODO use regex matching if list grows too large
  whitelist = ["User:ukeMPI337", "User:mkaul_MPI", "User:ukeMPINew", "User:ukeMPI335", "User:uke_mpi", "Bruker:MPICalibration", "Bruker:MPI"]
  if !(f["ACQ_method"] in whitelist)
    @warn "conversion of method not supported" f["ACQ_method"]
    return false
  end
  return true
end

function saveasMDF(filenameOut::String, filenameIn::String; kargs...)
  saveasMDF(filenameOut, MPIFile(filenameIn); kargs...)
end

function saveasMDF(filenameOut::String, f::MPIFile; filenameBG = nothing, kargs...)
  if isConvertibleToMDF(f)
    params = loadDataset(f;kargs...)
    if filenameBG != nothing
      appendBGDataset(params, filenameBG)
    end
    saveasMDF(filenameOut, params)
  else
    error("File not Convertible")
  end
end

function saveasMDFHacking(filenameOut::String, f::MPIFile)
    dataSet=loadDataset(f)
    dataSet["acqNumFrames"]=dataSet["acqNumPeriodsPerFrame"]*dataSet["acqNumFrames"]
    dataSet["acqNumPeriodsPerFrame"]=1
    dataSet["measData"]=reshape(dataSet["measData"],size(dataSet["measData"],1),size(dataSet["measData"],2),1,size(dataSet["measData"],3)*size(dataSet["measData"],4))
    dataSet["dfStrength"]=dataSet["dfStrength"][:,:,1:1]
    dataSet["acqOffsetField"]=dataSet["acqOffsetField"]
    #dataSet["acqOffsetFieldShift"]=dataSet["acqOffsetFieldShift"][:,1:1]
    dataSet["dfPhase"]=dataSet["dfPhase"][:,:,1:1]
    saveasMDF(filenameOut, dataSet)
    return dataSet
end

function saveasMDF(filename::String, params::Dict)
  # file has to be removed if exists. Otherwise h5create fails.
  isfile(filename) && rm(filename)
  h5open(filename, "w") do file
    saveasMDF(file, params)
  end
end

function compressCalibMDF(filenameOut::String, f::MPIFile; SNRThresh=2.0, kargs...)
  idx = Int64[]

  SNR = calibSNR(f)[:,:,1]
  for k=1:size(SNR,1)
    if maximum(SNR[k,:]) > SNRThresh
      push!(idx, k)
    end
  end

  compressCalibMDF(filenameOut, f, idx; kargs...)
end

function compressCalibMDF(filenamesOut::Vector{String}, f::MultiMPIFile; SNRThresh=2.0, kargs...)
  idx = Int64[]

  # We take the SNR from the first SF
  SNR = calibSNR(f[1])[:,:,1]
  for k=1:size(SNR,1)
    if maximum(SNR[k,:]) > SNRThresh
      push!(idx, k)
    end
  end

  for (i,f_) in enumerate(f)
    compressCalibMDF(filenamesOut[i], f_, idx; kargs...)
  end
end

function compressCalibMDF(filenameOut::String, f::MPIFile, idx::Vector{Int64};
                          sparsityTrafoRedFactor=1.0, sparsityTrafo="DCT-IV")
  params = loadMetadata(f)
  loadMeasParams(f, params, skipMeasData = true)
  loadCalibParams(f, params)
  params["calibSNR"] = calibSNR(f)
  loadRecoParams(f, params)

  data = systemMatrixWithBG(f, idx)

  params["calibSNR"] = params["calibSNR"][idx,:,:]
  if haskey(params, "rxTransferFunction")
    params["rxTransferFunction"] = params["rxTransferFunction"][idx,:]
  end
  params["measIsFrequencySelection"] = true
  params["measFrequencySelection"] = idx

  if sparsityTrafoRedFactor == 1.0
    params["measData"] = data
  else
    B = linearOperator(sparsityTrafo, calibSize(f))
    N = prod(calibSize(f))
    NBG = size(data,1) - N
    D = size(data,3)
    P = size(data,4)
    NRed = max(1, floor(Int, sparsityTrafoRedFactor*N))
    dataOut = similar(data, NBG+NRed, length(idx), D, P)
    subsamplingIndices = zeros(Int32, NRed, length(idx), D, P)

    fgdata = data[measFGFrameIdx(f),:,:,:]
    bgdata = data[measBGFrameIdx(f),:,:,:]

    dataOut[(NRed+1):end,:,:,:] = bgdata

    for k=1:length(idx), d=1:D, p=1:P
      I = B * fgdata[:,k,d,p]
      subsamplingIndices[:,k,d,p] = round.(Int32,reverse(sortperm(abs.(I)),dims=1)[1:NRed])
      dataOut[1:NRed,k,d,p] = I[vec(subsamplingIndices[:,k,d,p])]
    end

    params["measData"] = dataOut
    params["measIsSparsityTransformed"] = true
    params["measSparsityTransformation"] = sparsityTrafo
    params["measSubsamplingIndices"] = subsamplingIndices

    bgFrame = zeros(Bool, NRed+NBG)
    bgFrame[(NRed+1):end] .= true
    params["measIsBGFrame"] = bgFrame
    params["acqNumFrames"] = NRed+NBG
  end

  saveasMDF(filenameOut, params)
end

function loadAndProcessFFData(f::BrukerFile, nAverages::Int64, skipSwitchingFrames::Int64;addToEnd=0)
  dataFilename = joinpath(f.path,"rawdata.job0")
  #FileSize = stat(dataFilename).size
  #AllFrames = convert(Int,round(Int,FileSize/26928/4/3))
  #AllFrames = acqNumPeriodsPerPatch(f)*acqNumPeriodsPerFrame(f)
  AllFrames = acqNumPeriodsPerFrame(f)
  (Nx,Ny,Nz) = [length(union(acqOffsetField(f)[ll,1,:])) for ll in collect(1:3)]

  data_ = zeros(ComplexF64,Nx,Ny,Nz,rxNumFrequencies(f),rxNumChannels(f));

#  Pos_=acqOffsetField(f).*[-1,-1,-1]; #Field to position
  Pos=acqOffsetField(f).*[-1,-1,-1]; #Field to position

#  Pos = zeros(3,1,AllFrames);
#  Pos[1,1,:] = kron(vec(Pos_[1,1,:]),ones(2))
#  Pos[2,1,:] = kron(vec(Pos_[2,1,:]),ones(2))
#  Pos[3,1,:] = kron(vec(Pos_[3,1,:]),ones(2))

  dType = Int16 #acqNumAverages(f) == 1 ? Int16 : Int32  # Data are not averaged

  (PosNx,PosNy,PosNz)=[union(Pos[ll,1,:]) for ll in collect(1:3)]

  ds = open(dataFilename)
  raw = Mmap.mmap(ds, Array{dType,4},(rxNumSamplingPoints(f),1,rxNumChannels(f),AllFrames));

  for p = collect(1:Nx*Ny*Nz)
#    p = Nx*Ny*(pz-1)+pxy
    st = p*skipSwitchingFrames+(p-1)*(nAverages+addToEnd)+1

    ind = collect(st:st+nAverages-1)

    if any([length(union(Pos[ll,1,ind])) for ll in collect(1:3)].>1)
      error("Averaged data are not at the same FF-Position! Check Parameter Averages and/or skipSwitchingFrames (MPIFiles function loadAndProcessFFData)")
    end
    # sequence is not relevant
    x = collect(Nx:-1:1)[PosNx.==union(Pos[1,1,ind])]
    y = collect(Ny:-1:1)[PosNy.==union(Pos[2,1,ind])]
    z = collect(Nz:-1:1)[PosNz.==union(Pos[3,1,ind])]
    #for ch =1:3
      #dataSC = spectralLeakageCorrectedData(raw[:,1,ch,ind]);
      #data_[x,y,z,:,ch] = rfft(vec(mean(dataSC,dims=2)));
      data_[x,y,z,:,:] = rfft(mean(raw[:,1,:,ind],dims=3),1);
    #end
  end
  data = reshape(data_,:,rxNumFrequencies(f),rxNumChannels(f),1)
  return convert(Array{Complex{Float32},4},data); # allPos x freq x channel x 1
end

function convertCustomSF(filenameOut::String, f::BrukerFile, fBG::BrukerFile,nAverages::Int64,skipSwitchingFrames::Int64; nAveragesBG = nAverages,skipSwitchingFramesBG=skipSwitchingFrames)

  params = loadMetadata(f)
  loadMeasParams(f, params, skipMeasData = true)
  loadCalibParams(f, params)

  paramsBG = loadMetadata(fBG)

  params["calibSize"] = [length(union(params["acqOffsetField"][ll,1,:])) for ll in collect(1:3)]
  numFGFrames = prod(params["calibSize"])

  params["acqGradient"] = reshape(params["acqGradient"][:,:,1,1],3,3,1,1)
  params["dfPhase"] =  reshape(params["dfPhase"][:,:,1],1,3,1)
  params["calibDeltaSampleSize"] = [0.0, 0.0, 0.0] #Todo
  #params["time"] = now()
  params["calibMethod"] = "BrukerCustom"
  #params["version"] = v"2.0.0"
  params["dfStrength"] = reshape(params["dfStrength"][:,:,1,1],1,3,1)
  params["experimentIsCalibration"] = true
  #params["uuid"] = uuid4()
  params["measIsFourierTransformed"] = true
  params["measIsFastFrameAxis"] = true
  params["calibFovCenter"] = [mean(extrema(params["acqOffsetField"][ll,1,:])) for ll in collect(1:3)]
  params["acqNumPeriodsPerFrame"] = 1

  params["measIsFramePermutation"] = 1

  if rxHasTransferFunction(f) &&!measIsTFCorrected(f)
    params["transferFunction"] = rxTransferFunction(f)
  end

  params["calibFov"] = [sum(abs.(extrema(params["acqOffsetField"][ll,1,:])))./abs.(params["acqGradient"][ll,ll]) for ll in collect(1:3)]
  params["calibOrder"] = "xyz"

  params["acqOffsetField"] = reshape([mean(extrema(params["acqOffsetField"][ll,1,:])) for ll in collect(1:3)],3,1,1)
  params["acqNumAverages"] = nAverages

println("Part1")
  # Daten Laden
@time  fgFrames = loadAndProcessFFData(f,nAverages, skipSwitchingFrames) # allePos x freqs x r x 1
@time  bgFramesFull = loadAndProcessFFData(fBG,nAveragesBG, skipSwitchingFramesBG) # allePos x freqs x r x 1
  nAveragesBGHalf = nAveragesBG >1 ? div(nAveragesBG,2) : nAveragesBG
@time  bgFramesHalf = loadAndProcessFFData(fBG,div(nAveragesBG,2), skipSwitchingFramesBG+nAveragesBG-div(nAveragesBG,2)) # allePos x freqs x r x 1

  params["calibSNR"] = calculateSNRCustomSF(f,fgFrames,bgFramesFull,bgFramesHalf)

  (xBG,yBG,zBG) = [length(union(paramsBG["acqOffsetField"][ll,1,:])) for ll in collect(1:3)]
  (xFG,yFG,zFG) = params["calibSize"]

  bgFullReshaped = reshape(bgFramesFull,xBG,yBG,zBG,rxNumFrequencies(fBG),rxNumChannels(fBG))

  itp = interpolate(bgFullReshaped,(NoInterp(),BSpline(Linear()),BSpline(Linear()),NoInterp(),NoInterp())); #Interpolire bgFramesFull auf

  bgFramesFullInterp = itp(collect(1:xBG),range(1,yBG,length=yFG),range(1,zBG,length=zFG),collect(1:size(bgFullReshaped,4)),collect(1:size(bgFullReshaped,5)))

  params["acqNumFrames"] = numFGFrames + prod(size(bgFramesFullInterp)[1:3])

  tt = [round.(Int,collect(range(1,xBG+xFG,length=xBG))).+(kk-1)*(xBG+xFG) for kk =1:yFG*yFG];
  idxBGFrames = vcat(tt...);
  idxAllFrames = collect(1:params["acqNumFrames"])
  idxAllFrames[idxBGFrames] .= 0
  idxFGFrames = idxAllFrames[idxAllFrames.!=0];

  params["measFramePermutation"] =  vcat(idxFGFrames,idxBGFrames)
  params["measIsBGFrame"] = vcat(zeros(numFGFrames),ones(length(idxBGFrames)))
  dataTemp_ = vcat(fgFrames,reshape(bgFramesFullInterp,:,size(bgFullReshaped)[4:5]...,1))  #size(params["measData"])(22621, 26929, 3, 1)
  params["measData"] = convert(Array{Complex{Float32},4},dataTemp_);  #size(params["measData"])(22621, 26929, 3, 1)
  saveasMDF(filenameOut, params)
end

function convertCustomSF(filenameOut::String, f::Array{BrukerFileMeas}, fBG::BrukerFile,nAverages::Int64,skipSwitchingFrames::Int64; nAveragesBG = nAverages,skipSwitchingFramesBG=skipSwitchingFrames)

  paramsArray =[];
  paramsBG = loadMetadata(fBG)

  params = loadMetadata(f[1])
  loadMeasParams(f[1], params, skipMeasData = true)
  loadCalibParams(f[1], params)

  for fi in f[1:end]

  params_ = loadMetadata(fi)
  loadMeasParams(fi, params_, skipMeasData = true)
  loadCalibParams(fi, params_)
  push!(paramsArray,params_);
  end

  for i =2:length(f)
  params["acqOffsetField"] = cat(params["acqOffsetField"],paramsArray[i]["acqOffsetField"],dims=3)
  end

  params["calibSize"] = [length(union(params["acqOffsetField"][ll,1,:])) for ll in collect(1:3)]
  numFGFrames = prod(params["calibSize"])

  println(params["calibSize"])

  params["acqGradient"] = reshape(params["acqGradient"][:,:,1,1],3,3,1,1)
  params["dfPhase"] =  reshape(params["dfPhase"][:,:,1],1,3,1)
  params["calibDeltaSampleSize"] = [0.0, 0.0, 0.0] #Todo
  #params["time"] = now()
  params["calibMethod"] = "BrukerCustom"
  #params["version"] = v"2.0.0"
  params["dfStrength"] = reshape(params["dfStrength"][:,:,1,1],1,3,1)
  params["experimentIsCalibration"] = true
  #params["uuid"] = uuid4()
  params["measIsFourierTransformed"] = true
  params["measIsFastFrameAxis"] = true
  params["calibFovCenter"] = [mean(extrema(params["acqOffsetField"][ll,1,:])) for ll in collect(1:3)]
  params["acqNumPeriodsPerFrame"] = 1

  params["measIsFramePermutation"] = 1

  if rxHasTransferFunction(f[1]) &&!measIsTFCorrected(f[1])
    params["transferFunction"] = sampleTF(TransferFunction("/opt/data/TFs/TFHDF5/TF_3DMouse_NM0_07102019.h5"), f[1])
#    params["transferFunction"] = rxTransferFunction(f[1])
  end

  params["calibFov"] = [sum(abs.(extrema(params["acqOffsetField"][ll,1,:])))./abs.(params["acqGradient"][ll,ll]) for ll in collect(1:3)]
  params["calibOrder"] = "xyz"

  params["acqOffsetField"] = reshape([mean(extrema(params["acqOffsetField"][ll,1,:])) for ll in collect(1:3)],3,1,1)
  params["acqNumAverages"] = nAverages

println("Part1")
  # Daten Laden
#fgFrames =zeros(ComplexF64,prod(Int64,params["calibSize"]),26929,3,1)
@time global fgFrames =loadAndProcessFFData(f[1],nAverages, skipSwitchingFrames)
#global tmp =1
for (ii,ffi)  in enumerate(f[2:end])
println(ffi)
#pos = prod(Int64,[length(union(paramsArray[ii]["acqOffsetField"][ll,1,:])) for ll in collect(1:3)])
#@time fgFrames[tmp:tmp+pos-1,:,:,:] = loadAndProcessFFData(ffi,nAverages, skipSwitchingFrames) # allePos x freqs x r x 1
@time global fgFrames = cat(fgFrames,loadAndProcessFFData(ffi,nAverages, skipSwitchingFrames),dims=1) # allePos x freqs x r x 1
#global tmp = tmp+pos
end

@time  bgFramesFull = loadAndProcessFFData(fBG,nAveragesBG, skipSwitchingFramesBG) # allePos x freqs x r x 1
  nAveragesBGHalf = nAveragesBG >1 ? div(nAveragesBG,2) : nAveragesBG
@time  bgFramesHalf = loadAndProcessFFData(fBG,div(nAveragesBG,2), skipSwitchingFramesBG+nAveragesBG-div(nAveragesBG,2)) # allePos x freqs x r x 1

  params["calibSNR"] = calculateSNRCustomSF(f[1],fgFrames,bgFramesFull,bgFramesHalf)

  (xBG,yBG,zBG) = [length(union(paramsBG["acqOffsetField"][ll,1,:])) for ll in collect(1:3)]
  (xFG,yFG,zFG) = params["calibSize"]

  bgFullReshaped = reshape(bgFramesFull,xBG,yBG,zBG,rxNumFrequencies(fBG),rxNumChannels(fBG))

  itp = interpolate(bgFullReshaped,(NoInterp(),BSpline(Linear()),BSpline(Linear()),NoInterp(),NoInterp())); #Interpolire bgFramesFull auf

  bgFramesFullInterp = itp(collect(1:xBG),range(1,yBG,length=yFG),range(1,zBG,length=zFG),collect(1:size(bgFullReshaped,4)),collect(1:size(bgFullReshaped,5)))

  params["acqNumFrames"] = numFGFrames + prod(size(bgFramesFullInterp)[1:3])

  tt = [round.(Int,collect(range(1,xBG+xFG,length=xBG))).+(kk-1)*(xBG+xFG) for kk =1:yFG*yFG];
  idxBGFrames = vcat(tt...);
  idxAllFrames = collect(1:params["acqNumFrames"])
  idxAllFrames[idxBGFrames] .= 0
  idxFGFrames = idxAllFrames[idxAllFrames.!=0];

  params["measFramePermutation"] =  vcat(idxFGFrames,idxBGFrames)
  params["measIsBGFrame"] = vcat(zeros(numFGFrames),ones(length(idxBGFrames)))
  dataTemp_ = vcat(fgFrames,reshape(bgFramesFullInterp,:,size(bgFullReshaped)[4:5]...,1))  #size(params["measData"])(22621, 26929, 3, 1)
  params["measData"] = convert(Array{Complex{Float32},4},dataTemp_);  #size(params["measData"])(22621, 26929, 3, 1)
  saveasMDF(filenameOut, params)
end

hasKeyAndValue(paramDict,param) = haskey(paramDict, param) && paramDict[param] != nothing

function writeIfAvailable(file, paramOut, paramDict, paramIn )
  if hasKeyAndValue(paramDict, paramIn)
    write(file, paramOut, paramDict[paramIn])
  end
end

function saveasMDF(file::HDF5File, params::Dict)
  # general parameters
  write(file, "/version", "2.0.1")
  write(file, "/uuid", string(get(params,"uuid",uuid4() )))
  write(file, "/time", "$( get(params,"time", Dates.unix2datetime(time())) )")

  # study parameters
  write(file, "/study/name", get(params,"studyName","default") )
  write(file, "/study/number", get(params,"studyNumber",0))
  if hasKeyAndValue(params,"studyUuid")
    studyUuid = params["studyUuid"]
  else
    studyUuid = uuid4()
  end
  write(file, "/study/uuid", string(studyUuid))
  write(file, "/study/description", get(params,"studyDescription","n.a."))
  if hasKeyAndValue(params,"studyTime")
    write(file, "/study/time", string(params["studyTime"]))
  end

  # experiment parameters
  write(file, "/experiment/name", get(params,"experimentName","default") )
  write(file, "/experiment/number", get(params,"experimentNumber",0))
  if hasKeyAndValue(params,"experimentUuid")
    expUuid = params["experimentUuid"]
  else
    expUuid = uuid4()
  end
  write(file, "/experiment/uuid", string(expUuid))
  write(file, "/experiment/description", get(params,"experimentDescription","n.a."))
  write(file, "/experiment/subject", get(params,"experimentSubject","n.a."))
  write(file, "/experiment/isSimulation", Int8(get(params,"experimentIsSimulation",false)))

  # tracer parameters
  write(file, "/tracer/name", get(params,"tracerName","n.a") )
  write(file, "/tracer/batch", get(params,"tracerBatch","n.a") )
  write(file, "/tracer/vendor", get(params,"tracerVendor","n.a") )
  write(file, "/tracer/volume", get(params,"tracerVolume",0.0))
  write(file, "/tracer/concentration", get(params,"tracerConcentration",0.0) )
  write(file, "/tracer/solute", get(params,"tracerSolute","Fe") )
  tr = [string(t) for t in get(params,"tracerInjectionTime", [Dates.unix2datetime(time())]) ]
  write(file, "/tracer/injectionTime", tr)

  # scanner parameters
  write(file, "/scanner/facility", get(params,"scannerFacility","n.a") )
  write(file, "/scanner/operator", get(params,"scannerOperator","n.a") )
  write(file, "/scanner/manufacturer", get(params,"scannerManufacturer","n.a"))
  write(file, "/scanner/name", get(params,"scannerName","n.a"))
  write(file, "/scanner/topology", get(params,"scannerTopology","FFP"))

  # acquisition parameters
  write(file, "/acquisition/numAverages",  params["acqNumAverages"])
  write(file, "/acquisition/numFrames", get(params,"acqNumFrames",1))
  write(file, "/acquisition/numPeriods", get(params,"acqNumPeriodsPerFrame",1))
  write(file, "/acquisition/startTime", "$( get(params,"acqStartTime", Dates.unix2datetime(time())) )")

  writeIfAvailable(file, "/acquisition/gradient", params, "acqGradient")
  writeIfAvailable(file, "/acquisition/offsetField", params, "acqOffsetField")

  # drivefield parameters
  write(file, "/acquisition/drivefield/numChannels", size(params["dfStrength"],2) )
  write(file, "/acquisition/drivefield/strength", params["dfStrength"])
  write(file, "/acquisition/drivefield/phase", params["dfPhase"])
  write(file, "/acquisition/drivefield/baseFrequency", params["dfBaseFrequency"])
  write(file, "/acquisition/drivefield/divider", params["dfDivider"])
  write(file, "/acquisition/drivefield/cycle", params["dfCycle"])
  if !haskey(params, "dfWaveform")
    params["dfWaveform"] = "sine"
  end
  write(file, "/acquisition/drivefield/waveform", params["dfWaveform"])

  # receiver parameters
  write(file, "/acquisition/receiver/numChannels", params["rxNumChannels"])
  write(file, "/acquisition/receiver/bandwidth", params["rxBandwidth"])
  write(file, "/acquisition/receiver/numSamplingPoints", params["rxNumSamplingPoints"])
  if !haskey(params, "rxUnit")
    params["rxUnit"] = "V"
  end
  write(file, "/acquisition/receiver/unit",  params["rxUnit"])
  write(file, "/acquisition/receiver/dataConversionFactor",  params["rxDataConversionFactor"])
  if hasKeyAndValue(params,"rxTransferFunction")
    tf = params["rxTransferFunction"]
    group = file["/acquisition/receiver"]
    writeComplexArray(group, "transferFunction", tf)
  end
  writeIfAvailable(file, "/acquisition/receiver/inductionFactor",  params, "rxInductionFactor")

  # measurements
  if hasKeyAndValue(params, "measData")
    meas = params["measData"]
    if eltype(meas) <: Complex
      group = g_create(file,"/measurement")
      writeComplexArray(group, "/measurement/data", meas)
    else
      write(file, "/measurement/data", meas)
    end
    write(file, "/measurement/isFourierTransformed", Int8(params["measIsFourierTransformed"]))
    write(file, "/measurement/isSpectralLeakageCorrected", Int8(params["measIsSpectralLeakageCorrected"]))
    write(file, "/measurement/isTransferFunctionCorrected", Int8(params["measIsTFCorrected"]))
    write(file, "/measurement/isFrequencySelection", Int8(params["measIsFrequencySelection"]))
    write(file, "/measurement/isBackgroundCorrected",  Int8(params["measIsBGCorrected"]))
    write(file, "/measurement/isFastFrameAxis",  Int8(params["measIsFastFrameAxis"]))
    write(file, "/measurement/isFramePermutation",  Int8(params["measIsFramePermutation"]))
    writeIfAvailable(file, "/measurement/frequencySelection",  params, "measFrequencySelection")

    if hasKeyAndValue(params, "measFramePermutation")
      write(file, "/measurement/framePermutation", params["measFramePermutation"] )
    end
    if hasKeyAndValue(params, "measIsBGFrame")
      write(file, "/measurement/isBackgroundFrame", convert(Array{Int8}, params["measIsBGFrame"]) )
    end
    if hasKeyAndValue(params, "measIsSparsityTransformed")
      write(file, "/measurement/isSparsityTransformed", params["measIsSparsityTransformed"] )
      write(file, "/measurement/subsamplingIndices", params["measSubsamplingIndices"] )
      write(file, "/measurement/sparsityTransformation", params["measSparsityTransformation"] )
    end
  end

  # calibrations
  writeIfAvailable(file, "/calibration/snr",  params, "calibSNR")
  writeIfAvailable(file, "/calibration/fieldOfView",  params, "calibFov")
  writeIfAvailable(file, "/calibration/fieldOfViewCenter",  params, "calibFovCenter")
  writeIfAvailable(file, "/calibration/size",  params, "calibSize")
  writeIfAvailable(file, "/calibration/order",  params, "calibOrder")
  writeIfAvailable(file, "/calibration/positions",  params, "calibPositions")
  writeIfAvailable(file, "/calibration/offsetField",  params, "calibOffsetField")
  writeIfAvailable(file, "/calibration/deltaSampleSize",  params, "calibDeltaSampleSize")
  writeIfAvailable(file, "/calibration/method",  params, "calibMethod")
  if hasKeyAndValue(params, "calibIsMeanderingGrid")
    write(file, "/calibration/isMeanderingGrid", Int8(params["calibIsMeanderingGrid"]))
  end
  # reconstruction
  if hasKeyAndValue(params, "recoData")
    write(file, "/reconstruction/data", params["recoData"])
    write(file, "/reconstruction/fieldOfView", params["recoFov"])
    write(file, "/reconstruction/fieldOfViewCenter", params["recoFovCenter"])
    write(file, "/reconstruction/size", params["recoSize"])
    write(file, "/reconstruction/order", get(params,"recoOrder", "xyz"))
    if hasKeyAndValue(params,"recoPositions")
      write(file, "/reconstruction/positions", params["recoPositions"])
    end
    if hasKeyAndValue(params,"recoParameters")
      saveParams(file, "/reconstruction/_parameters", params["recoParameters"])
    end
  end

  writeIfAvailable(file, "/custom/auxiliaryData", params, "auxiliaryData")
end


text:
# Trivial BFV, scalar encoding small plaintext space, pow2 cyclotomic ring
include("bfv_triv.jl")

# Trivial BGV
include("bgv_triv.jl")

# BFV with p=65537 SIMD over cyclotomic ring
include("bfv_simd.jl")

# No encryption, just thest the PolyCRT encoding code
include("polycrt_encoding.jl")

# BFV with p=256 SIMD over non-cyclotomic ring
include("bfv_uint8.jl")

# Keyswitching for BFV
include("bfv_keyswitch.jl")

# Noise measurement for bfv
include("bfv_noise.jl")

# BFV with CRT representation
include("bfv_crt.jl")

# CKKS
include("ckks_triv.jl")

include("ckks_modswitch.jl")

include("ckks_rotate.jl")

include("ckks_matmul.jl")

include("ckks_modraise.jl")


text:
using DEBTrait
using Test

@testset "DEBTrait.jl" begin
    # Write your tests here.
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "capnproto"
version = v"0.7.0"

# Collection of sources required to build capnproto
sources = [
    ArchiveSource("https://capnproto.org/capnproto-c++-$(version).tar.gz",
                  "c9a4c0bd88123064d483ab46ecee777f14d933359e23bff6fb4f4dbd28b4cd41"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/capnproto-*/
atomic_patch -p2 ${WORKSPACE}/srcdir/patches/aligned-alloc.patch
(
    # Do native build to get a capnp that we can run
    mkdir build_native && cd build_native
    export CC=${CC_BUILD}
    export CXX=${CXX_BUILD}
    export LD=${LD_BUILD}
    ../configure --host=${MACHTYPE} --with-pic \
        lt_cv_prog_compiler_pic_works=yes \
        lt_cv_prog_compiler_pic_works_CXX=yes
    make -j${nproc}
)
export CAPNP=build_native/capnp
./configure --prefix=$prefix --host=$target --with-external-capnp
make -j${nproc}
make install
if [[ "${target}" == *-mingw* ]]; then
    # The build for Windows creates "${bindir}/capnpc" as a broken link to
    # "capnp": let's remove it and copy capnp.exe to capnpc.exe
    rm "${bindir}/capnpc"
    cp "${bindir}/capnp${exeext}" "${bindir}/capnpc${exeext}"
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct(["libcapnp", "libcapnp-$(version.major)-$(version.minor)"], :libcapnp),
    ExecutableProduct("capnp", :capnp)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies,
              preferred_gcc_version=v"5")

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
# ---
# title: 1154. Day of the Year
# id: problem1154
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/day-of-the-year/description/>
# hidden: true
# ---
#
# Given a string `date` representing a [Gregorian
# calendar](https://en.wikipedia.org/wiki/Gregorian_calendar) date formatted as
# `YYYY-MM-DD`, return the day number of the year.
#
#
#
# **Example 1:**
#
#
#
#     Input: date = "2019-01-09"
#     Output: 9
#     Explanation: Given date is the 9th day of the year in 2019.
#
#
# **Example 2:**
#
#
#
#     Input: date = "2019-02-10"
#     Output: 41
#
#
# **Example 3:**
#
#
#
#     Input: date = "2003-03-01"
#     Output: 60
#
#
# **Example 4:**
#
#
#
#     Input: date = "2004-03-01"
#     Output: 61
#
#
#
#
# **Constraints:**
#
#   * `date.length == 10`
#   * `date[4] == date[7] == '-'`, and all other `date[i]`'s are digits
#   * `date` represents a calendar date between Jan 1st, 1900 and Dec 31, 2019.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# ---
# title: 1422. Maximum Score After Splitting a String
# id: problem1422
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/maximum-score-after-splitting-a-string/description/>
# hidden: true
# ---
#
# Given a string `s` of zeros and ones, _return the maximum score after
# splitting the string into two **non-empty** substrings_ (i.e. **left**
# substring and **right** substring).
#
# The score after splitting a string is the number of **zeros** in the **left**
# substring plus the number of **ones** in the **right** substring.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "011101"
#     Output: 5
#     Explanation:
#     All possible ways of splitting s into two non-empty substrings are:
#     left = "0" and right = "11101", score = 1 + 4 = 5
#     left = "01" and right = "1101", score = 1 + 3 = 4
#     left = "011" and right = "101", score = 1 + 2 = 3
#     left = "0111" and right = "01", score = 1 + 1 = 2
#     left = "01110" and right = "1", score = 2 + 1 = 3
#
#
# **Example 2:**
#
#
#
#     Input: s = "00111"
#     Output: 5
#     Explanation: When left = "00" and right = "111", we get the maximum score = 2 + 3 = 5
#
#
# **Example 3:**
#
#
#
#     Input: s = "1111"
#     Output: 3
#
#
#
#
# **Constraints:**
#
#   * `2 <= s.length <= 500`
#   * The string `s` consists of characters '0' and '1' only.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module CoupledConditionalSMC

using RNGPool
import Statistics.mean
using Random

include("structures.jl")
include("subroutines.jl")
include("algorithms.jl")
include("unbiased.jl")
include("couplingTime.jl")
include("jls.jl")

export CCSMCIO, ccXpf!, initializeCCSMC

end # module


text:
function ode_addsteps!{calcVal,calcVal2,calcVal3}(k,t,uprev,u,dt,f,cache::Rosenbrock23ConstantCache,always_calc_begin::Type{Val{calcVal}} = Val{false},allow_calc_end::Type{Val{calcVal2}} = Val{true},force_calc_end::Type{Val{calcVal3}} = Val{false})
  if length(k)<2 || calcVal
    dT = ForwardDiff.derivative(tf,t)
    J = ForwardDiff.derivative(uf,uprev)
    W = 1-dt*d*J
    k₁ = W\(integrator.fsalfirst + dt*d*dT)
    f₁ = f(t+dt/2,uprev+dt*k₁/2)
    k₂ = W\(f₁-k₁) + k₁
    copyat_or_push!(k,1,k₁)
    copyat_or_push!(k,2,k₂)
  end
  nothing
end

function ode_addsteps!{calcVal,calcVal2,calcVal3}(k,t,uprev,u,dt,f,cache::Rosenbrock32ConstantCache,always_calc_begin::Type{Val{calcVal}} = Val{false},allow_calc_end::Type{Val{calcVal2}} = Val{true},force_calc_end::Type{Val{calcVal3}} = Val{false})
  if length(k)<2 || calcVal
    dT = ForwardDiff.derivative(tf,t)
    J = ForwardDiff.derivative(uf,uprev)
    W = 1-dt*d*J
    k₁ = W\(integrator.fsalfirst + dt*d*dT)
    f₁ = f(t+dt/2,uprev+dt*k₁/2)
    k₂ = W\(f₁-k₁) + k₁
    copyat_or_push!(k,1,k₁)
    copyat_or_push!(k,2,k₂)
  end
  nothing
end

function ode_addsteps!{calcVal,calcVal2,calcVal3}(k,t,uprev,u,dt,f,cache::Rosenbrock23Cache,always_calc_begin::Type{Val{calcVal}} = Val{false},allow_calc_end::Type{Val{calcVal2}} = Val{true},force_calc_end::Type{Val{calcVal3}} = Val{false})
  if length(k)<2 || calcVal
    @unpack k₁,k₂,k₃,du1,du2,f₁,vectmp,vectmp2,vectmp3,fsalfirst,fsallast,dT,J,W,tmp,uf,tf,linsolve_tmp = cache
    @unpack c₃₂,d = cache.tab
    uidx = eachindex(uprev)

    #=

    ### Jacobian does not need to be re-evaluated after an event
    ### Since it's unchanged

    # Setup Jacobian Calc
    sizeu  = size(u)
    tf.vf.sizeu = sizeu
    tf.uprev = uprev
    uf.vfr.sizeu = sizeu
    uf.t = t

    ForwardDiff.derivative!(dT,tf,t) # Time derivative of each component
    ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev))
    =#

    a = -dt*d
    for i in 1:length(uprev), j in 1:length(uprev)
      W[i,j] = I[i,j]+a*J[i,j]
    end

    Wfact = cache.factorization(W)

    a = -a
    for i in uidx
      linsolve_tmp[i] = @muladd fsalfirst[i] + a*dT[i]
    end

    @into! vectmp = Wfact\linsolve_tmp
    recursivecopy!(k₁,reshape(vectmp,size(u)...))
    for i in uidx
      tmp[i]=uprev[i]+dt*k₁[i]/2
    end
    f(t+dt/2,tmp,f₁)

    for i in uidx
      linsolve_tmp[i] = f₁[i]-k₁[i]
    end

    @into! vectmp2 = Wfact\linsolve_tmp
    for i in uidx
      k₂[i] = tmp[i] + k₁[i]
    end
    copyat_or_push!(k,1,k₁)
    copyat_or_push!(k,2,k₂)
  end
  nothing
end

function ode_addsteps!{calcVal,calcVal2,calcVal3}(k,t,uprev,u,dt,f,cache::Rosenbrock32Cache,always_calc_begin::Type{Val{calcVal}} = Val{false},allow_calc_end::Type{Val{calcVal2}} = Val{true},force_calc_end::Type{Val{calcVal3}} = Val{false})
  if length(k)<2 || calcVal
    @unpack k₁,k₂,k₃,du1,du2,f₁,vectmp,vectmp2,vectmp3,fsalfirst,fsallast,dT,J,W,tmp,uf,tf,linsolve_tmp = cache
    @unpack c₃₂,d = cache.tab
    uidx = eachindex(uprev)

    #=

    ### Jacobian does not need to be re-evaluated after an event
    ### Since it's unchanged

    # Setup Jacobian Calc
    sizeu  = size(u)
    tf.vf.sizeu = sizeu
    tf.uprev = uprev
    uf.vfr.sizeu = sizeu
    uf.t = t

    ForwardDiff.derivative!(dT,tf,t) # Time derivative of each component
    ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev))
    =#

    a = -dt*d
    for i in 1:length(u), j in 1:length(u)
      W[i,j] = @muladd I[i,j]+a*J[i,j]
    end

    Wfact = cache.factorization(W)

    a = -a
    for i in uidx
      linsolve_tmp[i] = @muladd fsalfirst[i] + a*dT[i]
    end

    @into! vectmp = Wfact\linsolve_tmp
    recursivecopy!(k₁,reshape(vectmp,size(u)...))
    for i in uidx
      tmp[i]=uprev[i]+dt*k₁[i]/2
    end
    f(t+dt/2,tmp,f₁)

    for i in uidx
      linsolve_tmp[i] = f₁[i]-k₁[i]
    end

    @into! vectmp2 = Wfact\linsolve_tmp
    for i in uidx
      k₂[i] = tmp[i] + k₁[i]
    end
    copyat_or_push!(k,1,k₁)
    copyat_or_push!(k,2,k₂)
  end
  nothing
end


text:
using GeoData, Test, ArchGDAL, Dates, Statistics

include(joinpath(dirname(pathof(GeoData)), "../test/test_utils.jl"))

A = [missing 7.0f0; 2.0f0 missing]
ga = GeoArray(A, (X, Y); missingval=missing)
ga99 = replace_missing(ga, -9999)
@test all(ga99 .=== [-9999.0f0 7.0f0; 2.0f0 -9999.0f0])
gaNaN = replace_missing(ga, NaN32)
@test all(gaNaN .=== [NaN32 7.0f0; 2.0f0 NaN32])
gaMi = replace_missing(ga, missing)
@test all(gaMi .=== ga)

@testset "boolmask" begin
    @test boolmask(Array, ga) == [false true; true false]
    @test boolmask(ga) == [false true; true false]
    @test boolmask(ga99) == [false true; true false]
    @test boolmask(gaNaN) == [false true; true false]
    @test dims(boolmask(ga)) == (X(Base.OneTo(2), mode=NoIndex()), Y(Base.OneTo(2), mode=NoIndex()))
end

@testset "missingmask" begin
    @test all(missingmask(A) .=== [missing true; true missing])
    @test all(missingmask(ga) .=== [missing true; true missing])
    @test all(missingmask(ga99) .=== [missing true; true missing])
    @test all(missingmask(gaNaN) .=== [missing true; true missing])
    @test dims(missingmask(ga)) == (X(Base.OneTo(2), mode=NoIndex()), Y(Base.OneTo(2), mode=NoIndex()))
end

@testset "mask" begin
    A1 = [missing 1; 2 3]
    ga1 = GeoArray(A1, (X, Y); missingval=missing)
    @test all(mask(ga1; to=ga) .=== [missing 1; 2 missing])
    ga2 = replace_missing(ga1 .* 1.0; missingval=NaN)
    @test all(mask(ga2; to=ga) .=== [NaN 1.0; 2.0 NaN])
    ga3 = replace_missing(ga1; missingval=-9999)
    mask!(ga3; to=ga)
    @test all(ga3 .=== [-9999 1; 2 -9999])
end

@testset "classify" begin
    A1 = [missing 1; 2 3]
    ga1 = GeoArray(A1, (X, Y); missingval=missing)
    @test all(classify(ga1, 1=>99, 2=>88, 3=>77) .=== [missing 99; 88 77])
    @test all(classify(ga1, 1=>99, 2=>88, 3=>77; others=0) .=== [missing 99; 88 77])
    @test all(classify(ga1, 1=>99, 2=>88; others=0) .=== [missing 99; 88 0])

    A2 = [1.0 2.5; 3.0 4.0]
    ga2 = GeoArray(A2, (X, Y); missingval=missing)
    @test classify(ga2, (2, 3)=>:x, >(3)=>:y) == [1.0 :x; 3.0 :y]
    @test classify(ga2, (>=(1), <(2))=>:x, >=(3)=>:y) == [:x 2.5; :y :y]
    classify!(ga2, (1, 2.5)=>0.0, >=(3)=>-1.0; lower=(>), upper=(<=))
    @test ga2 == [1.0 0.0; -1.0 -1.0]
    classify!(ga2, [1 2.5 0.0; 2.5 4.0 -1.0]; lower=(>), upper=(<=))
    @test ga2 == [1.0 0.0; -1.0 -1.0]
end

@testset "points" begin
    ga = GeoArray(A, (X(9.0:1.0:10.0), Y(0.1:0.1:0.2)); missingval=missing)
    @test all(collect(points(ga; order=(Y, X))) .=== [missing (0.2, 9.0); (0.1, 10.0) missing])
    @test all(collect(points(ga; order=(X, Y))) .=== [missing (9.0, 0.2); (10.0, 0.1) missing])
    @test all(points(ga; order=(X, Y), ignore_missing=true) .===
              [(9.0, 0.1) (9.0, 0.2); (10.0, 0.1) (10.0, 0.2)])
end

@testset "extract" begin
    A = [1 2; 3 4]
    ga = GeoArray(A, (X(9.0:1.0:10.0), Y(0.1:0.1:0.2)); name=:test, missingval=missing)
    @testset "points" begin
        @test all(extract(ga, [(9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]) .===
                  [(X=9.0, Y=0.1, test=1), (X=10.0, Y=0.2, test=4), (X=10.0, Y=0.3, test=missing)])
        @test all(extract(ga, [(0.1, 9.0), (0.2, 10.0), (0.3, 10.0)]; order=(Y, X)) .===
                  [(Y=0.1, X=9.0, test=1), (Y=0.2, X=10.0, test=4), (Y=0.3, X=10.0, test=missing)])
    end
    @testset "Tables.jl compatible" begin
        @test all(extract(ga, [(X=9.0, Y=0.1), (X=10.0, Y=0.2), (X=10.0, Y=0.3)]) .===
              [(X=9.0, Y=0.1, test=1), (X=10.0, Y=0.2, test=4), (X=10.0, Y=0.3, test=missing)])
        @test all(extract(ga, [(X=9.0, Y=0.1), (X=10.0, Y=0.2), (X=10.0, Y=0.3)]; order=(Y=>:Y, X=>:X)) .===
              [(Y=0.1, X=9.0, test=1), (Y=0.2, X=10.0, test=4), (Y=0.3, X=10.0, test=missing)])
    end
end

@testset "trim, crop, extend" begin
    A = [missing missing missing
         missing 2.0     0.5
         missing 1.0     missing]
    ga = GeoArray(A, (X(1.0:1.0:3.0), Y(1.0:1.0:3.0)); missingval=missing)
    # Test with missing on all sides
    ga_r = rot180(ga)
    trimmed = trim(ga)
    trimmed_r = trim(ga_r)
    @test all(trimmed .=== [2.0 0.5; 1.0 missing])
    @test all(trimmed_r .=== [missing 1.0; 0.5 2.0])
    cropped = crop(ga; to=trimmed)
    cropped_r = crop(ga_r; to=trimmed_r)
    @test all(cropped .=== trimmed)
    @test all(cropped_r .=== trimmed_r)
    extended = extend(cropped, ga)[1]
    extended_r = extend(cropped_r; to=ga_r)
    @test all(extended .=== ga)
    @test all(extended_r .=== ga_r)
end

@testset "mosaic" begin
    reg1 = GeoArray([0.1 0.2; 0.3 0.4], (X(2.0:1.0:3.0), Y(5.0:1.0:6.0)))
    reg2 = GeoArray([1.1 1.2; 1.3 1.4], (X(3.0:1.0:4.0), Y(6.0:1.0:7.0)))
    irreg1 = GeoArray([0.1 0.2; 0.3 0.4], (X([2.0, 3.0]), Y([5.0, 6.0])))
    irreg2 = GeoArray([1.1 1.2; 1.3 1.4], (X([3.0, 4.0]), Y([6.0, 7.0])))

    span_x1 = Explicit(vcat((1.5:1.0:2.5)', (2.5:1.0:3.5)'))
    span_x2 = Explicit(vcat((2.5:1.0:3.5)', (3.5:1.0:4.5)'))
    exp1 = GeoArray([0.1 0.2; 0.3 0.4], (X([2.0, 3.0]; mode=Sampled(span=span_x1)), Y([5.0, 6.0])))
    exp2 = GeoArray([1.1 1.2; 1.3 1.4], (X([3.0, 4.0]; mode=Sampled(span=span_x2)), Y([6.0, 7.0])))
    @test val(span(mosaic(first, exp1, exp2), X)) == [1.5 2.5 3.5; 2.5 3.5 4.5]
    @test all(mosaic(first, reg1, reg2) .===
              mosaic(first, irreg1, irreg2) .===
              mosaic(first, irreg1, irreg2) .===
              [0.1 0.2 missing;
               0.3 0.4 1.2;
               missing 1.3 1.4])
    @test all(mosaic(last, reg1, reg2) .===
              mosaic(last, irreg1, irreg2) .===
              mosaic(last, exp1, exp2) .=== [0.1 0.2 missing;
                                             0.3 1.1 1.2;
                                             missing 1.3 1.4])

    # 3 dimensions
    A1 = GeoArray(ones(2, 2, 2), (X(2.0:-1.0:1.0), Y(5.0:1.0:6.0), Ti(DateTime(2001):Year(1):DateTime(2002))))
    A2 = GeoArray(zeros(2, 2, 2), (X(3.0:-1.0:2.0), Y(4.0:1.0:5.0), Ti(DateTime(2002):Year(1):DateTime(2003))))
    @test all(mosaic(mean, A1, A2) |> parent .=== cat([missing missing missing
                                                 missing 1.0     1.0
                                                 missing 1.0     1.0    ],
                                                [0.0     0.0 missing
                                                 0.0     0.5     1.0
                                                 missing 1.0     1.0    ],
                                                [0.0     0.0     missing
                                                 0.0     0.0     missing
                                                 missing missing missing], dims=3))
end

@testset "resample" begin
    raster_path = maybedownload("https://download.osgeo.org/geotiff/samples/gdal_eg/cea.tif")

    output_res = 0.0027
    output_crs = EPSG(4326)
    resample_method = "near"

    ## Resample cea.tif manually with ArchGDAL
    wkt = convert(String, convert(WellKnownText, output_crs))
    AG_output = ArchGDAL.read(raster_path) do dataset
        ArchGDAL.gdalwarp([dataset], ["-t_srs", "$(wkt)",
                                "-tr", "$(output_res)", "$(output_res)",
                                "-r", "$(resample_method)"]) do warped
            ArchGDAL.read(ArchGDAL.getband(warped, 1))
        end
    end

    ## Resample cea.tif using resample
    cea = read(geoarray(raster_path))
    GD_output = resample(cea, output_res; crs=output_crs, method=resample_method)

    cea_permuted = permutedims(read(geoarray(raster_path)), (Y, X, Band))
    GD_permuted_output = resample(cea_permuted, output_res; crs=output_crs, method=resample_method)

    # Compare ArchGDAL, resample and permuted resample
    @test AG_output == GD_output[Band(1)] == permutedims(GD_permuted_output, (X, Y, Band))[Band(1)]
    @test abs(step(dims(GD_output, Y))) ≈ abs(step(dims(GD_output, X))) ≈
          abs(step(dims(GD_permuted_output, X))) ≈ output_res

    @testset "snapped size and dim index match" begin
        snaptarget = GD_output
        snapped = resample(cea; to=snaptarget)
        @test size(snapped) == size(snaptarget)
        @test isapprox(index(snaptarget, Y), index(snapped, Y))
        @test isapprox(index(snaptarget, X), index(snapped, X))
    end
end


text:
# TODO: implement readgexf

"""
savegexf(f::IO, g::SimpleGraph, gname::AbstractString)

Writes a graph `g` with name `gname`
to a file `f` in the
[Gexf](http://gexf.net/format/) format.

Returns 1 (number of graphs written).
"""
function savegexf(f::IO, g::SimpleGraph, gname::AbstractString)
    xdoc = XMLDocument()
    xroot = create_root(xdoc, "gexf")
    set_attribute(xroot,"xmlns","http://www.gexf.net/1.2draft")
    set_attribute(xroot,"version","1.2")
    set_attribute(xroot,"xmlns:xsi","http://www.w3.org/2001/XMLSchema-instance")
    set_attribute(xroot,"xsi:schemaLocation","http://www.gexf.net/1.2draft/gexf.xsd")

    xmeta = new_child(xroot, "meta")
    xdesc = new_child(xmeta, "description")
    add_text(xdesc, gname)
    xg = new_child(xroot, "graph")
    strdir = is_directed(g) ? "directed" : "undirected"
    set_attribute(xg,"defaultedgetype",strdir)

    xnodes = new_child(xg, "nodes")
    for i=1:nv(g)
        xv = new_child(xnodes, "node")
        set_attribute(xv,"id","$(i-1)")
    end

    xedges = new_child(xg, "edges")
    m=0
    for e in edges(g)
        xe = new_child(xedges, "edge")
        set_attribute(xe,"id","$m")
        set_attribute(xe,"source","$(src(e)-1)")
        set_attribute(xe,"target","$(dst(e)-1)")
        m+=1
    end

    show(f, xdoc)
    return 1
end

filemap[:gexf] = (NI, NI, savegexf, NI)


text:
#####
##### Gradient AD implementation using ForwardDiff
#####

import .ForwardDiff
using BenchmarkTools: @belapsed

struct ForwardDiffLogDensity{L, C} <: ADGradientWrapper
    ℓ::L
    gradientconfig::C
end

function Base.show(io::IO, ℓ::ForwardDiffLogDensity)
    print(io, "ForwardDiff AD wrapper for ", ℓ.ℓ,
          ", w/ chunk size ", length(ℓ.gradientconfig.seeds))
end

_default_chunk(ℓ) = ForwardDiff.Chunk(dimension(ℓ))

# defined to make the tag match
_logdensity_closure(ℓ) = x -> logdensity(ℓ, x)

function _default_gradientconfig(ℓ, chunk::ForwardDiff.Chunk)
    ForwardDiff.GradientConfig(_logdensity_closure(ℓ), zeros(dimension(ℓ)), chunk)
end

function _default_gradientconfig(ℓ, chunk::Integer)
    _default_gradientconfig(ℓ, ForwardDiff.Chunk(chunk))
end

"""
    ADgradient(:ForwardDiff, ℓ; chunk, gradientconfig)
    ADgradient(Val(:ForwardDiff), ℓ; chunk, gradientconfig)

Wrap a log density that supports evaluation of `Value` to handle `ValueGradient`, using
`ForwardDiff`.

Keywords are passed on to `ForwardDiff.GradientConfig` to customize the setup. In
particular, chunk size can be set with a `chunk` keyword argument (accepting an integer or a
`ForwardDiff.Chunk`).
"""
function ADgradient(::Val{:ForwardDiff}, ℓ;
                    chunk = _default_chunk(ℓ),
                    gradientconfig = _default_gradientconfig(ℓ, chunk))
    ForwardDiffLogDensity(ℓ, gradientconfig)
end

function logdensity_and_gradient(fℓ::ForwardDiffLogDensity, x::AbstractVector)
    @unpack ℓ, gradientconfig = fℓ
    buffer = _diffresults_buffer(ℓ, x)
    result = ForwardDiff.gradient!(buffer, _logdensity_closure(ℓ), x, gradientconfig)
    _diffresults_extract(result)
end

"""
$(SIGNATURES)

Default chunk sizes to try for benchmarking. Fewer than `M`, always contains `1` and `N`.
"""
function heuristic_chunks(N, M = 20)
    step = max(N ÷ M, 1)
    Ns = 1:step:N
    if N ∉ Ns
        Ns = vcat(Ns, N)
    end
    Ns
end

"""
$(SIGNATURES)

Benchmark a log density problem with various chunk sizes using ForwardDiff.

`chunks`, which defaults to all possible chunk sizes, determines the chunks that are tried.

The function returns `chunk => time` pairs, where `time` is the benchmarked runtime in
seconds, as determined by `BenchmarkTools.@belapsed`. The gradient is evaluated at `x`
(defaults to zeros).

*Runtime may be long* because of tuned benchmarks, so when `markprogress == true` (the
default), dots are printed to mark progress.

This function is not exported, but part of the API when `ForwardDiff` is imported.
"""
function benchmark_ForwardDiff_chunks(ℓ;
                                      chunks = heuristic_chunks(dimension(ℓ), 20),
                                      markprogress = true,
                                      x = zeros(dimension(ℓ)))
    map(chunks) do chunk
        ∇ℓ = ADgradient(Val(:ForwardDiff), ℓ; chunk = ForwardDiff.Chunk(chunk))
        markprogress && print(".")
        chunk => @belapsed logdensity_and_gradient($(∇ℓ), $(x))
    end
end


text:
@testset "Subpopulation age report" begin

@testset "function subpopulationAgeReport" begin
    subpop1 = Subpopulation( "Subpop1" )
    addSubpopulationCondition!( subpop1, MPcondition( "branch", !=, "A" ) )
    subpop2 = Subpopulation( "Subpop2" )
    addSubpopulationCondition!( subpop2, MPcondition( "time in node", >=, 84 ) )
    subpop3 = Subpopulation( "Subpop3" )
    addSubpopulationCondition!( subpop3, MPcondition( "started as", ==,
        "Reserve junior" ) )
    subpop4 = Subpopulation( "Subpop4" )
    addSubpopulationCondition!( subpop4, MPcondition( "tenure", <, 12 ) )
    subpop5 = Subpopulation( "Subpop5" )
    addSubpopulationCondition!( subpop5, MPcondition( "was", !=, "B junior" ) )
    subpop6 = Subpopulation( "Subpop6" )
    subpop7 = Subpopulation( "Subpop7" )
    addSubpopulationCondition!( subpop7, MPcondition( "age", >=, 240 ) )

    setSubpopulationSourceNode!.( [subpop1, subpop2, subpop3, subpop4,
        subpop5, subpop6, subpop7], ["", "", "", "Reserve junior", "Career",
        "Branch A", "Foo"] )

    @test isempty( subpopulationAgeReport( mpSim, [-12.0, -6.0], 12, :age,
        subpop1, subpop2, subpop3, subpop4, subpop5, subpop6, subpop7 ) )
    @test isempty( subpopulationAgeReport( mpSim, [12.0, 6.0], -12, :age,
        subpop1, subpop2, subpop3, subpop4, subpop5, subpop6, subpop7 ) )
    @test isempty( subpopulationAgeReport( mpSim, [12.0, 6.0], 12, :age,
        subpop7 ) )
    @test isempty( subpopulationAgeReport( mpSim, [12.0, 6.0], 12, :foo,
        subpop1, subpop2, subpop3, subpop4, subpop5, subpop6, subpop7 ) )
    @test isempty( subpopulationAgeReport( mpSim, -12, 12, :age, subpop1,
        subpop2, subpop3, subpop4, subpop5, subpop6, subpop7 ) )

    report = subpopulationAgeReport( mpSim, 12, 12, :age, subpop1, subpop2,
        subpop3, subpop4, subpop5, subpop6, subpop7 )
    @test all( haskey.( Ref( report ), string.( "Subpop", 1:6 ) ) ) &&
        !haskey( report, "Subpop7" )
    @test all( map( collect( keys( report ) ) ) do subpopName
        subpopReport = report[subpopName]
        return all( isa.( subpopReport[:, "min"], Missing ) .|
            ( 240 .<= subpopReport[:, "min"] .<= 360 ) )
    end )
    @test all( map( collect( keys( report ) ) ) do subpopName
        subpopReport = report[subpopName]
        return all( isa.( subpopReport[:, "max"], Missing ) .|
            ( 240 .<= subpopReport[:, "max"] .<= 360 ) )
    end )

    report = subpopulationAgeReport( mpSim, 12, 12, :tenure, subpop1, subpop2,
        subpop3, subpop4, subpop5, subpop6, subpop7 )
    @test all( map( collect( keys( report ) ) ) do subpopName
        subpopReport = report[subpopName]
        return all( isa.( subpopReport[:, "min"], Missing ) .|
            ( 0 .<= subpopReport[:, "min"] .<= 120 ) )
    end )
    @test all( map( collect( keys( report ) ) ) do subpopName
        subpopReport = report[subpopName]
        return all( isa.( subpopReport[:, "max"], Missing ) .|
            ( 0 .<= subpopReport[:, "max"] .<= 120 ) )
    end )

    report = subpopulationAgeReport( mpSim, 12, 12, :timeInNode, subpop1,
        subpop2, subpop3, subpop4, subpop5, subpop6, subpop7 )

    @test all( map( collect( keys( report ) ) ) do subpopName
        subpopReport = report[subpopName]
        return all( isa.( subpopReport[:, "min"], Missing ) .|
            ( 0 .<= subpopReport[:, "min"] .<= 120 ) )
    end )
    @test all( map( collect( keys( report ) ) ) do subpopName
        subpopReport = report[subpopName]
        return all( isa.( subpopReport[:, "max"], Missing ) .|
            ( 0 .<= subpopReport[:, "max"] .<= 120 ) )
    end )
end  # @testset "function subpopulationAgeReport"

end  # @testset "Subpopulation age report"


text:
"""
Overview: collections of routines for O(n) access to data
"""

"""
> find_index_global(combs,ndofs,p)
It computes the position in the cmap array of a given combination
- combs : total combinations
- ndofs : number of degrees of freedom
- p : order of the expansion
"""
function find_index_global(combs,ndofs,p)
  #
  # It computes the position in the cmap array of a given combination
  #
  pos = 1
  for i = 1:p
    pos += (combs[i]-1)*ndofs^(p-i)
  end
  return pos
end

"""
> find_index_map(j::Int64,k::Int64,s::Int64,p::Int64,ndofs::Int64,index_comb::Array{Int64})
It finds the position of a map for a certain index_comb during time-derivatives residual calculation.
- j : counter
- k : counter
- s : counter
- p : expansion order
- ndofs : number of degrees of freedom
- index_comb : reference monomials exponentials
"""
function find_index_map(j::Int64,k::Int64,s::Int64,p::Int64,
                        ndofs::Int64,index_comb::Array{Int64})
  pos = 1
  @inbounds @simd for i = 1:k
    pos += (index_comb[i]-1)*ndofs^(p-j-i+1)
  end
  @inbounds pos += (s-1)*ndofs^(p-j-k)
  @inbounds @simd for i = k+j+1:p
    pos += (index_comb[i]-1)*ndofs^(p-i)
  end
  return pos
end

"""
> find_index_rdyn(j::Int64,k::Int64,ndofs::Int64,index_comb::Array{Int64})
It finds the position of a reduced dynamics vector for a certain index_comb during time-derivatives residual calculation.
- j : counter
- k : counter
- ndofs : number of degrees of freedom
- index_comb : reference monomials exponentials
"""
function find_index_rdyn(j::Int64,k::Int64,ndofs::Int64,index_comb::Array{Int64})
  pos = 1
  @inbounds @simd for i = k+1:k+j
    pos += (index_comb[i]-1)*ndofs^(k+j-i)
  end
  return pos
end

"""
> find_position_g(j::Int64,p::Int64,ndofs::Int64,index_comb::Array{Int64})
It finds the maps position to compute their contribution to the quadratic nonlinearities.
- j : counter
- p : order
- ndofs : number of degrees of freedom
- index_comb : reference monomials exponentials
"""
function find_position_g(j::Int64,p::Int64,ndofs::Int64,
                          index_comb::Array{Int64})
  pos1 = 1
  pos2 = 1
  @inbounds @simd for i = 1:j
    pos1 += (index_comb[i]-1)*ndofs^(j-i)
  end
  @inbounds @simd for i = j+1:p
    pos2 += (index_comb[i]-1)*ndofs^(p-i)
  end
  return pos1,pos2
end

"""
> find_position_g(j::Int64,k::Int64,p::Int64,ndofs::Int64,index_comb::Array{Int64})
It finds the maps position to compute their contribution to the cubic nonlinearities.
- j : counter
- k : counter
- p : order
- ndofs : number of degrees of freedom
- index_comb : reference monomials exponentials
"""
function find_position_h(j::Int64,k::Int64,p::Int64,
                          ndofs::Int64,index_comb::Array{Int64})
  #
  # It finds the maps position to compute their contribution to the
  # cubic nonlinearities.
  #
  pos1 = 1
  pos2 = 1
  pos3 = 1
  @inbounds @simd for i = 1:j
    pos1 += (index_comb[i]-1)*ndofs^(j-i)
  end
  @inbounds @simd for i = 1:k
    pos2 += (index_comb[j+i]-1)*ndofs^(k-i)
  end
  @inbounds @simd for i = j+k+1:p
    pos3 += (index_comb[i]-1)*ndofs^(p-i)
  end
  return pos1,pos2,pos3
end


text:
module PracticePackage

include("add.jl")

export add

end


text:
# nlyse COVID-19 arrival time in countries
# this corresponds to Fig.2 in
# Blasius B. (2020) Power-law distribution in the number of confirmed COVID-19 cases. Chaos
# https://doi.org/10.1063/5.0013031

module A

using CSV, DataFrames, Polynomials, PyPlot, StatsBase
using LsqFit

include("read_data.jl")
include("tools.jl")

function time_exponents()
# compute time dependent exponents
   t,C_ts,D_ts,countries = read_all_data()

   (days, ncountries) = size(C_ts)
   println("days ", days)
   println("total _countries ", ncountries)

   # calculate exponent from log-likelihood
   expsC = zeros(days)
   expsD = zeros(days)
   for t=1:days
       C = C_ts[t,:]
       C = C[findall(x->x>0, C)]  # remove zeros
       mu, st = estimate(C)
       expsC[t] = mu

       D = D_ts[t,:]
       D = D[findall(x->x>0, D)]  # remove zeros
       mu, st = estimate(D)
       expsD[t] = mu
   end
   expsC, expsD
end

function prevalence_overtime(fig, fig_x, fig_y1,fig_width_y)
   t,I_ts,D_ts,countries = read_all_data()

   (days, ncountries) = size(I_ts)
   println("days ", days)
   println("total _countries ", ncountries)

   first_date = Date(2020,1,22)  # first entry in database: Jan 22, 2020

   # make sure here to include only the data until the time span
   # that should be investigated
   last_date = first_date + Day(days-1)
   println(last_date)

   cols = ["k", "r", "b", "g", "m", "y"]
   times = [0,14,28,42,56,70]   # vector of time delays
   times = times .+ (days - 76)
   dates = string.(last_date .- Day.(times))
   println(dates)

   for i = 1: length(times)
       t = times[i] # time difference
       c = cols[i]  # color

       I = I_ts[end .- t,:]
       ind = findall(x->x>0, I)
       I = sort(I[ind], rev=true)
       nI = length(I)
       println("countries, n ",nI)

       D = D_ts[end .- t,:]
       ind = findall(x->x>0, D)
       D = sort(D[ind], rev=true)
       nD = length(D)

       # log -bin Histogram
       ex1I,ey1I = log2histogram(I)
       x1I = log.(ex1I)
       y1I = log.(ey1I)
       ex1D,ey1D = log2histogram(D)
       x1D = log.(ex1D)
       y1D = log.(ey1D)

   ax = PyPlot.axes([fig_x, fig_y1, 0.8, fig_width_y])
       loglog(ex1I,ey1I,"-", color=c)
       scatter(ex1I,ey1I,s=50, color=c, label="distribution on "*dates[i])
       myfont = Dict("fontname"=>"Sans","fontsize"=>20)
       xticks(fontsize=14)
       yticks(fontsize=14)
       xlim(1,1e6)
       xlabel("cases, n",fontsize=16)
       ylabel(L"frequency, $P(n)$",fontsize=16, labelpad=10)
       legend(loc="lower left",fancybox="true", fontsize=12)
       ax1 = ax.transAxes
       text(-0.15, 0.95, "a)", fontsize=18, fontweight="bold", transform = ax1)

   ax = PyPlot.axes([0.67, 0.82, 0.26, 0.15]) # Inset: cumulative distribution
       loglog(I,1:nI,"-", color=c)
       ylim(1,200)
       xlim(1,5e5)
       xticks(fontsize=12)
       yticks(fontsize=12)
       myfont = Dict("fontname"=>"Sans","fontsize"=>10)
       xlabel("cases, n", fontsize=14)
       ylabel(L"$N \cdot C(n)$", fontsize=14)
   end

end

function arrival_histogram(fig, fig_x, fig_y2,fig_width_y)

  c,d,t = arrival_times()
  ind = sortperm(t)
  t = t[ind]; c=c[ind]; d = d[ind]
  # first entry Jan 22 2020
  # at this day 6 countries already invaded, China, Japan, South Korea, Thailand, Taiwan, US)
  t = t .- t[1]

  h = StatsBase.fit(Histogram, t,  nbins=100, closed=:left)
  h_weight = h.weights
  h_bins = h.edges[1][2:end]

  println(size(h_bins))
  println(size(h_weight))
  @. model(x, p) = p[1]*exp(x*p[2])
  p0 = [0.5, 0.5]
  expfit = curve_fit(model, collect(h_bins)[1:57], h_weight[1:57], p0)
  println("exp fit: ",coef(expfit))
  cf = coef(expfit)

  alphaC, alphaD =  time_exponents()
  ndays = length(alphaC)

  h_dates = Date(2020,1,21) + Day.(h_bins)

  ind_march22 = Dates.value(Date(2020,3,22) - Date(2020,1,22)) + 1

  show_days= Dates.value(Date(2020,4,05) - Date(2020,1,22)) + 1

   majorformatter = matplotlib.dates.DateFormatter("%d.%m")
   majorlocator = matplotlib.dates.DayLocator(interval=7)

   ax = PyPlot.axes([fig_x, fig_y2, 0.8, fig_width_y])
       bar(h_dates,h_weight)
       plot(h_dates[1:61], cf[1]*exp.(cf[2]*h_bins[1:61]),"k", linewidth=2,
       linestyle="dashed")
       plot([Date(2020,3,22),Date(2020,3,22)], [0,15],"-r", linewidth=2)
       ax.xaxis_date()
       ax.xaxis.set_major_formatter(majorformatter)
       ax.xaxis.set_major_locator(majorlocator)
       xticks(Date(2020,1,26) .+ Day.(0:7:10*7),fontsize=14)
       yticks(fontsize=14)
       xlim(Date(2020,1,21),Date(2020,4,05))
       ylim(0,15)
       ylabel("countries",fontsize=16, labelpad=10)
       xlabel("arrival date", fontsize=16)
       for label in ax.get_xticklabels()
           label.set_ha("right")
           label.set_rotation(30)
       end
       ax1 = ax.transAxes
       text(-0.15, 0.95, "b)", fontsize=18, fontweight="bold", transform = ax1)

   ax1 = PyPlot.axes([0.25, 0.3, 0.28, 0.15]) # Inset: cumulative distribution
       plot(1:ndays,alphaC, linewidth=2, color="black")
       # plot the exponent of deaths from day 26, the first day with at least
       # 5 countries with a confirmed death case
       plot(26:ndays,alphaD[26:ndays], linewidth=2, color="blue")
       plot([ind_march22, ind_march22],[1,2],"r-")
       xticks(fontsize=12)
       yticks(fontsize=12)
       xlim(0,show_days)
       ylim(1,2)
       ax1.yaxis.grid()
       xlabel("time (days)", fontsize=14)
       ylabel(L"exponent $\mu$", fontsize=14, labelpad=10)
end

function make_plot()
    fig_sep_y   = 0.1
    fig_width_y = 0.4
    fig_width_y = 0.39
    fig_x = 0.15
    fig_y2=0.1
    fig_y1=fig_y2 + fig_width_y + fig_sep_y

    fig_Iy = 0.25
    fig_Iwidth_y = 0.1
    fig_Iwidth_y = 0.12
    fig_Iwidth_x = 0.22
    fig_Iwidth_x = 0.26

    fig=figure("prevalence over time1", figsize=(8,12))
    clf()
    prevalence_overtime(fig, fig_x, fig_y1,fig_width_y)
    arrival_histogram(fig, fig_x, fig_y2,fig_width_y)

   #savefig("../figures/distribution_time.png")
end

make_plot()

end  #module


text:
module SimpleUpdate

using TensorOperations, TensorNetworkTensors, LinearAlgebra, Parameters
export SimpleUpdateProblem, simpleupdate,
       update!, update,
       extract_weight!, extract_weights!,
       energy

include("simpleupdateproblem.jl")

export tfisinghamiltonian, tfisingpropagator
include("ising.jl")
end # module


text:
nx = 256
os = 3
buf = getfield(BUF, :buf)
x = getfield(BUF, :int32_buf)
nx_add = 1400000
nx_new = 36000
pref = path * "/SampleFiles/"
cfile = pref * "Restricted/03_02_27_20140927.sjis.ch"
files = String[ "UW/00012502123W"               "uw"            "_"
                "SEGY/03.334.12.09.00.0362.1"   "passcal"       "pa-full"
                "SEGY/03.334.12.09.00.0362.1"   "passcal"       "passcal"
                "SAC/test_le.sac"               "sac"           "_"
                "SEED/test.mseed"               "mseed"         "_"
                "Restricted/2014092709*.cnt"    "win32"         "win"
                "ASCII/0215162000.c00"          "lennartz"      "_"
                "ASCII/geo-slist.csv"           "geocsv.slist"  "_"
                "Restricted/SHW.UW.mseed"       "mseed"         "lo-mem"
                "Restricted/test_rev_1.segy"    "segy"          "full"
                "Restricted/test_rev_1.segy"    "segy"          "_"
                "SAC/test_be.sac"               "sac"           "full"
                "SAC/test_be.sac"               "sac"           "_"
                "ASCII/geo-tspair.csv"          "geocsv"        "_"
                ]
checkbuf_8!(buf, 65536)
checkbuf_8!(buf, 4*(os + nx))
checkbuf!(x, os + nx)

# test fillx_i16_le!
y = rand(Int16, nx)
copyto!(buf, 1, reinterpret(UInt8, y), 1, 2*nx)
fillx_i16_le!(x, buf, nx, os)
@test x[1+os:nx+os] == y

# test fillx_i32_le!
y = rand(Int32, nx)
copyto!(buf, 1, reinterpret(UInt8, y), 1, 4*nx)
fillx_i32_le!(x, buf, nx, os)
@test x[1+os:nx+os] == y

fillx_i32_be!(x, buf, nx, os)
@test x[1+os:nx+os] == bswap.(y)

# Do not run on Appveyor; it can't access restricted files, so this breaks
printstyled("  read_data\n", color=:light_green)
nf = size(files,1)
for n = 1:nf
  fname = pref * files[n,1]
  if (occursin("Restricted", fname)==false || (has_restricted==true))
    fwild = fname[1:end-1] * "*"
    f_call = files[n,2]
    opt = files[n,3]
    printstyled(string("    ", n, "/", nf, " ", f_call, "\n"), color=:light_green)
    if opt == "pa-full"
      S = read_data("passcal", fname, full=true)
      S = read_data("passcal", fwild, full=true)
      S = read_data(fname)
      S = read_data(fname, full=true)
    elseif opt == "win"
      S = read_data(f_call, fwild, cf=cfile)
      S = read_data(fwild, cf=cfile)
    elseif opt == "slist"
      S = read_data("geocsv.slist", fname)
      S = read_data("geocsv.slist", fwild)
      S = read_data(fwild)
    elseif opt == "lo-mem"
      S = read_data(f_call, fname, nx_new=nx_new, nx_add=nx_add)
      S = read_data(f_call, fwild, nx_new=nx_new, nx_add=nx_add)
    elseif opt == "full"
      S = read_data(f_call, fname, full=true)
      S = read_data(f_call, fwild, full=true)
      S = read_data(fwild, full=true)
    else
      S = read_data(f_call, fname)
      if f_call == "uw"
        fwild = fname[1:end-3]*"*"*"W"
      end
      S = read_data(f_call, fwild)
      S = read_data(fwild)
    end
  end
end

@test_throws ErrorException read_data("deez", "nutz.sac")
resize!(x, 65535)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ResourceManagementErrorWithDetails <: SwaggerModel
    code::Any # spec type: Union{ Nothing, String } # spec name: code
    message::Any # spec type: Union{ Nothing, String } # spec name: message
    target::Any # spec type: Union{ Nothing, String } # spec name: target
    details::Any # spec type: Union{ Nothing, Vector{ResourceManagementErrorWithDetails} } # spec name: details

    function ResourceManagementErrorWithDetails(;code=nothing, message=nothing, target=nothing, details=nothing)
        o = new()
        validate_property(ResourceManagementErrorWithDetails, Symbol("code"), code)
        setfield!(o, Symbol("code"), code)
        validate_property(ResourceManagementErrorWithDetails, Symbol("message"), message)
        setfield!(o, Symbol("message"), message)
        validate_property(ResourceManagementErrorWithDetails, Symbol("target"), target)
        setfield!(o, Symbol("target"), target)
        validate_property(ResourceManagementErrorWithDetails, Symbol("details"), details)
        setfield!(o, Symbol("details"), details)
        o
    end
end # type ResourceManagementErrorWithDetails

const _property_map_ResourceManagementErrorWithDetails = Dict{Symbol,Symbol}(Symbol("code")=>Symbol("code"), Symbol("message")=>Symbol("message"), Symbol("target")=>Symbol("target"), Symbol("details")=>Symbol("details"))
const _property_types_ResourceManagementErrorWithDetails = Dict{Symbol,String}(Symbol("code")=>"String", Symbol("message")=>"String", Symbol("target")=>"String", Symbol("details")=>"Vector{ResourceManagementErrorWithDetails}")
Base.propertynames(::Type{ ResourceManagementErrorWithDetails }) = collect(keys(_property_map_ResourceManagementErrorWithDetails))
Swagger.property_type(::Type{ ResourceManagementErrorWithDetails }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ResourceManagementErrorWithDetails[name]))}
Swagger.field_name(::Type{ ResourceManagementErrorWithDetails }, property_name::Symbol) =  _property_map_ResourceManagementErrorWithDetails[property_name]

function check_required(o::ResourceManagementErrorWithDetails)
    true
end

function validate_property(::Type{ ResourceManagementErrorWithDetails }, name::Symbol, val)
end


text:
# Test using the Boston housing and neighborhood dataset
@testset "Bostonhsg" begin

    boston = sdataset("Bostonhsg")

    # ----------------
    #  Distance  with 3
    # ----------------
    Wdist3 = dnearneigh(boston, threshold = 3)

    @test nobs(Wdist3) == 506
    @test mean(Wdist3) ≈ 44.711462 atol = 1e-5
    @test median(Wdist3) ≈ 29.0
    @test minimum(Wdist3) == 0
    @test maximum(Wdist3) == 122
    @test nislands(Wdist3) == 4
    @test islands(Wdist3) == [55, 65, 286, 287]

    # Moran's I
    mboston = moran(boston.MEDV, Wdist3, permutations = 0)

    # ----------------
    #  Distance  with 4
    # ----------------
    Wdist4 = dnearneigh(boston, threshold = 4)

    @test nobs(Wdist4) == 506
    @test mean(Wdist4) ≈ 72.225296 atol = 1e-5
    @test median(Wdist4) ≈ 56.0
    @test minimum(Wdist4) == 1
    @test maximum(Wdist4) == 178
    @test nislands(Wdist4) == 0
    @test islands(Wdist4) == []

    # Moran's I with 999 permutations
    mboston = moran(boston.MEDV, Wdist4, rng = StableRNG(1234567), permutations = 999)

    @test score(mboston) ≈ 0.267418 atol = 1e-5
    @test pvalue(mboston) ≈ 0.001 atol = 1e-5
    @test std(mboston) ≈ 0.013811 atol = 1e-5
    @test zscore(mboston) ≈ 19.493602 atol = 1e-5

    # ----------------
    #  KNN with k = 10
    # ----------------
    Wknn10 = knearneigh(boston, k = 10)

    @test nobs(Wknn10) == 506
    @test mean(Wknn10) == 10
    @test median(Wknn10) ≈ 10.0
    @test minimum(Wknn10) == 10
    @test maximum(Wknn10) == 10
    @test nislands(Wknn10) == 0

    # Moran I with 999 permutations
    mboston = moran(boston.MEDV, Wknn10, rng = StableRNG(1234567), permutations = 999)

    @test score(mboston) ≈ 0.530743 atol = 1e-5
    @test std(mboston) ≈ 0.018956 atol = 1e-5
    @test zscore(mboston) ≈ 28.080397 atol = 1e-5
    @test pvalue(mboston) ≈ 0.001 atol = 1e-5

    # Moran's I with 9999 permutations
    mboston = moran(boston.MEDV, Wknn10, rng = StableRNG(1234567), permutations = 9999)

    @test score(mboston) ≈ 0.530743 atol = 1e-5
    @test std(mboston) ≈ 0.018541 atol = 1e-5
    @test zscore(mboston) ≈ 28.727399 atol = 1e-5
    @test pvalue(mboston) ≈ 0.0001 atol = 1e-5

    # Geary's c I with 9999 permutations
    cboston = geary(boston.MEDV, Wknn10, rng = StableRNG(1234567), permutations = 9999)

    @test score(cboston) ≈ 0.432582 atol = 1e-5
    @test std(cboston) ≈ 0.022467 atol = 1e-5
    @test zscore(cboston) ≈ -25.249835 atol = 1e-5
    @test pvalue(cboston) ≈ 0.0001 atol = 1e-5

    # Local Moran
    lmboston = localmoran(boston.MEDV, Wknn10, rng = StableRNG(1234567), permutations = 9999)

    # ----------------
    #  Plot Recipes
    # ----------------

    # Test Moran's I distribution
    @test_nowarn RecipesBase.apply_recipe(Dict{Symbol, Any}(), mboston)

    # Test Moran's Plot
    @test_nowarn RecipesBase.apply_recipe(Dict{Symbol, Any}(), boston.MEDV, Wknn10)
    @test_nowarn RecipesBase.apply_recipe(Dict{Symbol, Any}([(:standardize, false)]), boston.MEDV, Wknn10)

    # Test LISA Cluster Map with points
    @test_nowarn RecipesBase.apply_recipe(Dict{Symbol, Any}(), boston, lmboston)
    @test_nowarn RecipesBase.apply_recipe(Dict{Symbol, Any}(), boston.geometry, lmboston)

end