text:
function init_pseudo_observable_mappings!(m::Model1002)

    pseudo_names = [:y_t, :y_f_t, :NaturalRate, :π_t, :OutputGap, :ExAnteRealRate, :LongRunInflation,
                    :MarginalCost, :Wages, :FlexibleWages, :Hours, :FlexibleHours, :z_t,
                    :Expected10YearRateGap, :NominalFFR, :Expected10YearRate,
                    :Expected10YearNaturalRate,
                    :ExpectedNominalNaturalRate, :NominalRateGap, :LaborProductivityGrowth, :u_t]
                    # :i_f_t, :R_t, :c_f_t, :qk_f_t, :k_f_t, :kbar_f_t, :u_f_t, :rk_f_t, :w_f_t, :L_f_t, :rktil_f_t, :n_f_t, :c_t,
                    # :b_t, :r_f_t]

    if subspec(m) == "ss12"
        to_add = [:g_t, :b_t, :μ_t, :λ_f_t, :λ_w_t, :rm_t, :σ_ω_t, :μ_e_t,
                  :γ_t, :π_star_t, :e_lr_t, :e_tfp_t, :e_gdpdef_t, :e_corepce_t, :e_gdp_t, :e_gdi_t]
        pseudo_names = vcat(pseudo_names, to_add)
    end

    if subspec(m) in ["ss13", "ss14", "ss15", "ss16", "ss17", "ss18", "ss19"]
        push!(pseudo_names, :Sinf_t, :Sinf_w_coef_t, :ι_p, :πtil_t, :πtil_t1, :e_tfp_t)
        if subspec(m) in ["ss14", "ss15", "ss16", "ss18", "ss19"]
            push!(pseudo_names, :e_tfp_t1)
        end
    end

    if haskey(get_settings(m), :add_covid_pseudoobs)
        if get_setting(m, :add_covid_pseudoobs) && subspec(m) in ["ss59", "ss60", "ss61", "ss62", "ss63", "ss64", "ss65", "ss66", "ss67", "ss68", "ss69", "ss70", "ss71", "ss72", "ss73", "ss74", "ss75", "ss76", "ss77", "ss78", "ss79", "ss80", "ss81", "ss82", "ss83", "ss84", "ss85"]
            push!(pseudo_names, :ziid, :varphiiid, :biidc)
        end
    end

    if haskey(get_settings(m), :add_pseudo_gdp)
        if get_setting(m, :add_pseudo_gdp) && subspec(m) in ["ss59", "ss60", "ss61", "ss62", "ss63", "ss64", "ss65", "ss66", "ss67", "ss68", "ss69", "ss70", "ss71", "ss72", "ss73", "ss74", "ss75", "ss76", "ss77", "ss78", "ss79", "ss80", "ss81", "ss82", "ss83", "ss84", "ss85"]
            push!(pseudo_names, :PseudoGDP)
        end
    end

    if haskey(get_settings(m), :add_pseudo_corepce)
        if get_setting(m, :add_pseudo_corepce) && subspec(m) in ["ss59", "ss60", "ss61", "ss62", "ss63", "ss64", "ss65", "ss66", "ss67", "ss68", "ss69", "ss70", "ss71", "ss72", "ss73", "ss74", "ss75", "ss76", "ss77", "ss78", "ss79", "ss80", "ss81", "ss82", "ss83", "ss84", "ss85"]
            push!(pseudo_names, :PseudoCorePCE)
        end
    end

    if haskey(get_settings(m), :add_NominalWageGrowth)
        if get_setting(m, :add_NominalWageGrowth)
            push!(pseudo_names, :NominalWageGrowth)
        end
    end
    if haskey(get_settings(m), :add_ztil)
        if get_setting(m, :add_ztil)
            push!(pseudo_names, :ztil)
        end
    end
    if haskey(get_settings(m), :add_zp)
        if get_setting(m, :add_zp)
            push!(pseudo_names, :zp)
        end
    end
    if haskey(get_settings(m), :add_pgap)
        if get_setting(m, :add_pgap)
            push!(pseudo_names, :pgap)
        end
    end
    if haskey(get_settings(m), :add_ygap)
        if get_setting(m, :add_ygap)
            push!(pseudo_names, :ygap)
        end
    end
    if haskey(get_settings(m), :add_altpolicy_pgap)
        if get_setting(m, :add_altpolicy_pgap)
            push!(pseudo_names, :pgap)
        end
    end
    if haskey(get_settings(m), :add_altpolicy_ygap)
        if get_setting(m, :add_altpolicy_ygap)
            push!(pseudo_names, :ygap)
        end
    end

    if haskey(get_settings(m), :add_urhat)
        if get_setting(m, :add_urhat)
            push!(pseudo_names, :urhat)
        end
    end

    if haskey(get_settings(m), :add_rw)
        if get_setting(m, :add_rw)
            push!(pseudo_names, :rw)
            push!(pseudo_names, :Rref)
        end
    end

    if haskey(get_settings(m), :add_laborshare_measurement)
        if get_setting(m, :add_laborshare_measurement)
            push!(pseudo_names, :laborshare_t)
        end
    end

    if haskey(get_settings(m), :add_laborproductivity_measurement)
        if get_setting(m, :add_laborproductivity_measurement)
            push!(pseudo_names, :laborproductivity)
        end
    end

    if haskey(get_settings(m), :add_nominalgdp_level)
        if get_setting(m, :add_nominalgdp_level)
            push!(pseudo_names, :NominalGDPLevel)
        end
    end

    if haskey(get_settings(m), :add_nominalgdp_growth)
        if get_setting(m, :add_nominalgdp_growth)
            push!(pseudo_names, :NominalGDPGrowth)
        end
    end

    if haskey(get_settings(m), :add_flexible_price_growth)
        if get_setting(m, :add_flexible_price_growth)
            push!(pseudo_names, :FlexibleGDPGrowth, :FlexibleConsumptionGrowth, :FlexibleInvestmentGrowth)
        end
    end

    if haskey(get_settings(m), :add_cumulative)
        if get_setting(m, :add_cumulative)
            push!(pseudo_names, :AccumOutputGap, :GDPLevel,
                  :TechnologyLevel, :FlexibleGDPLevel, :ConsumptionLevel,
                  :FlexibleConsumptionLevel, :InvestmentLevel,
                  :FlexibleInvestmentLevel)
        end
    end

    if haskey(get_settings(m), :add_laborproductivitygrowth_nome_measurement)
        if get_setting(m, :add_laborproductivitygrowth_nome_measurement)
            push!(pseudo_names, :LaborProductivityGrowthNoME)
        end
    end

    if haskey(get_settings(m), :add_Epi_t_measurement)
        if get_setting(m, :add_Epi_t_measurement)
            push!(pseudo_names, :Epi_t)
        end
    end

    # Create PseudoObservable objects
    pseudo = OrderedDict{Symbol,PseudoObservable}()
    for k in pseudo_names
        pseudo[k] = PseudoObservable(k)
    end

    # Fill in names and reverse transforms
    pseudo[:y_t].name = "Output Growth"
    pseudo[:y_t].longname = "Output Growth Per Capita"

    pseudo[:y_f_t].name = "Flexible Output Growth"
    pseudo[:y_f_t].longname = "Output that would prevail in a flexible-price economy."

    pseudo[:NaturalRate].name = "Real Natural Rate"
    pseudo[:NaturalRate].longname = "The real interest rate that would prevail in a flexible-price economy."
    pseudo[:NaturalRate].rev_transform = quartertoannual

    pseudo[:π_t].name = "Inflation"
    pseudo[:π_t].longname = "Inflation"
    pseudo[:π_t].rev_transform = quartertoannual

    pseudo[:OutputGap].name = "Output Gap"
    pseudo[:OutputGap].longname = "Output Gap"

    pseudo[:ExAnteRealRate].name = "Ex Ante Real Rate"
    pseudo[:ExAnteRealRate].longname = "Ex Ante Real Rate"
    pseudo[:ExAnteRealRate].rev_transform = quartertoannual

    pseudo[:LongRunInflation].name = "Long Run Inflation"
    pseudo[:LongRunInflation].longname = "Long Run Inflation"
    pseudo[:LongRunInflation].rev_transform = quartertoannual

    pseudo[:MarginalCost].name = "Marginal Cost"
    pseudo[:MarginalCost].longname = "Marginal Cost"

    pseudo[:Wages].name = "Wages"
    pseudo[:Wages].longname = "Wages"

    pseudo[:FlexibleWages].name = "Flexible Wages"
    pseudo[:FlexibleWages].longname = "Wages that would prevail in a flexible-wage economy"

    pseudo[:Hours].name = "Hours"
    pseudo[:Hours].longname = "Hours"

    pseudo[:FlexibleHours].name     = "Flexible Hours"
    pseudo[:FlexibleHours].longname = "Flexible Hours"

    pseudo[:z_t].name     = "z_t (Technology Growth minus Steady State Growth)"
    pseudo[:z_t].longname = "z_t (Technology Growth minus Steady State Growth)"

    pseudo[:Expected10YearRateGap].name     = "Expected 10-Year Rate Gap"
    pseudo[:Expected10YearRateGap].longname = "Expected 10-Year Rate Gap"
    pseudo[:Expected10YearRateGap].rev_transform = quartertoannual

    pseudo[:NominalFFR].name     = "Nominal FFR"
    pseudo[:NominalFFR].longname = "Nominal FFR at an annual rate"
    pseudo[:NominalFFR].rev_transform = quartertoannual

    pseudo[:Expected10YearRate].name     = "Expected 10-Year Rate"
    pseudo[:Expected10YearRate].longname = "Expected 10-Year Interest Rate"
    pseudo[:Expected10YearRate].rev_transform = quartertoannual

    pseudo[:Expected10YearNaturalRate].name     = "Expected 10-Year Natural Rate"
    pseudo[:Expected10YearNaturalRate].longname = "Expected 10-Year Natural Rate of Interest"
    pseudo[:Expected10YearNaturalRate].rev_transform = quartertoannual

    pseudo[:ExpectedNominalNaturalRate].name     = "Expected Nominal Natural Rate"
    pseudo[:ExpectedNominalNaturalRate].longname = "Natural Rate + Expected Inflation"
    pseudo[:ExpectedNominalNaturalRate].rev_transform = quartertoannual

    pseudo[:NominalRateGap].name     = "Nominal Rate Gap"
    pseudo[:NominalRateGap].longname = "Nominal FFR - Nominal Natural Rate"
    pseudo[:NominalRateGap].rev_transform = quartertoannual

    pseudo[:LaborProductivityGrowth].name     = "Labor Productivity Growth"
    pseudo[:LaborProductivityGrowth].longname = "Labor Productivity Growth Rate"
    pseudo[:LaborProductivityGrowth].rev_transform = quartertoannual

    pseudo[:u_t].name     = "u_t"
    pseudo[:u_t].longname = "u_t"

    # pseudo[:i_f_t].name     = "i_f_t"
    # pseudo[:i_f_t].longname = "i_f_t"

    # pseudo[:R_t].name     = "R_t"
    # pseudo[:R_t].longname = "R_t"

    # pseudo[:c_f_t].name     = "c_f_t"
    # pseudo[:c_f_t].longname = "c_f_t"

    # pseudo[:r_f_t].name     = "r_f_t"
    # pseudo[:r_f_t].longname = "r_f_t"

    # pseudo[:qk_f_t].name     = "qk_f_t"
    # pseudo[:qk_f_t].longname = "qk_f_t"

    # pseudo[:k_f_t].name     = "k_f_t"
    # pseudo[:k_f_t].longname = "k_f_t"

    # pseudo[:kbar_f_t].name     = "kbar_f_t"
    # pseudo[:kbar_f_t].longname = "kbar_f_t"

    # pseudo[:u_f_t].name     = "u_f_t"
    # pseudo[:u_f_t].longname = "u_f_t"

    # pseudo[:rk_f_t].name     = "rk_f_t"
    # pseudo[:rk_f_t].longname = "rk_f_t"

    # pseudo[:w_f_t].name     = "w_f_t"
    # pseudo[:w_f_t].longname = "w_f_t"

    # pseudo[:L_f_t].name     = "L_f_t"
    # pseudo[:L_f_t].longname = "L_f_t"

    # pseudo[:rktil_f_t].name     = "rktil_f_t"
    # pseudo[:rktil_f_t].longname = "rktil_f_t"

    # pseudo[:b_t].name     = "b_t"
    # pseudo[:b_t].longname = "b_t"

    # pseudo[:n_f_t].name     = "n_f_t"
    # pseudo[:n_f_t].longname = "n_f_t"

    if subspec(m) in ["ss13", "ss14", "ss15", "ss16", "ss17", "ss18", "ss19"]
        pseudo[:Sinf_t].name     = "Sinf_t"
        pseudo[:Sinf_t].longname = "Sinf_t, PDV of Emc_t"
        pseudo[:Sinf_w_coef_t].name     = "Sinf_w_coef_t"
        pseudo[:Sinf_w_coef_t].longname = "Sinf_w_coef_t, PDV of Emc_t multiplied by coefficient"
        pseudo[:ι_p].name = "iota_p"
        pseudo[:ι_p].longname = "iota_p"
        pseudo[:πtil_t].name     = "pitil_t"
        pseudo[:πtil_t].longname = "Fundamental Inflation"
        pseudo[:πtil_t1].name     = "pitil_t1"
        pseudo[:πtil_t1].longname = "Fundamental Inflation Lag 1"
        pseudo[:e_tfp_t].name = "e_tfp_t"
        pseudo[:e_tfp_t].longname = "e_tfp_t"
        if subspec(m) in ["ss14", "ss15", "ss16", "ss18", "ss19"]
            pseudo[:e_tfp_t1].name = "e_tfp_t1"
            pseudo[:e_tfp_t1].longname = "e_tfp_t1"
        end
    end

    # Other exogenous processes
    if subspec(m) == "ss12"
        for i in to_add
            pseudo[i].name = string(DSGE.detexify(i))
            pseudo[i].longname = string(DSGE.detexify(i))
        end
    end

    if haskey(get_settings(m), :add_NominalWageGrowth)
        if get_setting(m, :add_NominalWageGrowth)
            pseudo[:NominalWageGrowth].name = "Nominal Wage Growth"
            pseudo[:NominalWageGrowth].longname = "Nominal Wage Growth"
        end
    end

    if haskey(get_settings(m), :add_laborshare_measurement)
        if get_setting(m, :add_laborshare_measurement)
            pseudo[:laborshare_t].name     = "Log Labor Share"
            pseudo[:laborshare_t].longname = "Log Labor Share"
        end
    end

    if haskey(get_settings(m), :add_nominalgdp_level)
        if get_setting(m, :add_nominalgdp_level)
            pseudo[:NominalGDPLevel].name     = "Nominal GDP Level"
            pseudo[:NominalGDPLevel].longname = "Nominal GDP Level"
        end
    end

    if haskey(get_settings(m), :add_nominalgdp_growth)
        if get_setting(m, :add_nominalgdp_growth)
            pseudo[:NominalGDPGrowth].name     = "Nominal GDP Growth"
            pseudo[:NominalGDPGrowth].longname = "Nominal GDP Growth"
        end
    end

    if haskey(get_settings(m), :add_cumulative)
        if get_setting(m, :add_cumulative)
            pseudo[:AccumOutputGap].name     = "Accumulated Output Gap"
            pseudo[:AccumOutputGap].longname = "Accumulated Output Gap"
            pseudo[:GDPLevel].name     = "GDP Level"
            pseudo[:GDPLevel].longname = "GDP Level"
            pseudo[:TechnologyLevel].name     = "Technology (z) Level"
            pseudo[:TechnologyLevel].longname = "Technology (z) Level"
            pseudo[:FlexibleGDPLevel].name     = "Flexible GDP Level"
            pseudo[:FlexibleGDPLevel].longname = "Flexible GDP Level"
            pseudo[:ConsumptionLevel].name     = "Consumption Level"
            pseudo[:ConsumptionLevel].longname = "Consumption Level"
            pseudo[:FlexibleConsumptionLevel].name     = "Flexible Consumption Level"
            pseudo[:FlexibleConsumptionLevel].longname = "Flexible Consumption Level"
            pseudo[:InvestmentLevel].name     = "Investment Level"
            pseudo[:InvestmentLevel].longname = "Investment Level"
            pseudo[:FlexibleInvestmentLevel].name     = "Flexible Investment Level"
            pseudo[:FlexibleInvestmentLevel].longname = "Flexible Investment Level"
        end
    end

    if haskey(get_settings(m), :add_flexible_price_growth)
        if get_setting(m, :add_flexible_price_growth)
            pseudo[:FlexibleGDPGrowth].name     = "Flexible GDP Level"
            pseudo[:FlexibleGDPGrowth].longname = "Flexible GDP Growth"
            pseudo[:FlexibleConsumptionGrowth].name     = "Flexible Consumption Growth"
            pseudo[:FlexibleConsumptionGrowth].longname = "Flexible Consumption Growth"
            pseudo[:FlexibleInvestmentGrowth].name     = "Flexible Investment Growth"
            pseudo[:FlexibleInvestmentGrowth].longname = "Flexible Investment Growth"
        end
    end

    if haskey(get_settings(m), :add_ztil)
        if get_setting(m, :add_ztil)
            pseudo[:ztil].name     = "ztil"
            pseudo[:ztil].longname = "ztil"
        end
    end

    if haskey(get_settings(m), :add_zp)
        if get_setting(m, :add_zp)
            pseudo[:zp].name     = "zp"
            pseudo[:zp].longname = "zp"
        end
    end

    if haskey(get_settings(m), :add_altpolicy_pgap)
        if get_setting(m, :add_altpolicy_pgap)
            pseudo[:pgap].name     = "pgap"
            pseudo[:pgap].longname = "pgap"
        end
    end

    if haskey(get_settings(m), :add_altpolicy_ygap)
        if get_setting(m, :add_altpolicy_ygap)
            pseudo[:ygap].name     = "ygap"
            pseudo[:ygap].longname = "ygap"
        end
    end

    if haskey(get_settings(m), :add_urhat)
        if get_setting(m, :add_urhat)
            pseudo[:urhat].name     = "urhat"
            pseudo[:urhat].longname = "urhat"
        end
    end

    if haskey(get_settings(m), :add_rw)
        if get_setting(m, :add_rw)
            pseudo[:rw].name     = "rw"
            pseudo[:rw].longname = "rw"
            pseudo[:Rref].name     = "Rref"
            pseudo[:Rref].longname = "Rref"
        end
    end

    if haskey(get_settings(m), :add_covid_pseudoobs)
        if get_setting(m, :add_covid_pseudoobs) && subspec(m) in ["ss59", "ss60", "ss61", "ss62", "ss63", "ss64", "ss65", "ss66", "ss67", "ss68", "ss69", "ss70", "ss71", "ss72", "ss73", "ss74", "ss75", "ss76", "ss77", "ss78", "ss79", "ss80", "ss81", "ss82", "ss83", "ss84", "ss85"]
            pseudo[:ziid].name     = "ziid"
            pseudo[:ziid].longname = "ziid"
            pseudo[:biidc].name     = "biidc"
            pseudo[:biidc].longname = "biidc"
            pseudo[:varphiiid].name     = "varphiiid"
            pseudo[:varphiiid].longname = "varphiiid"
        end
    end

    if haskey(get_settings(m), :add_pseudo_gdp)
        if get_setting(m, :add_pseudo_gdp) && subspec(m) in ["ss59", "ss60", "ss61", "ss62", "ss63", "ss64", "ss65", "ss66", "ss67", "ss68", "ss69", "ss70", "ss71", "ss72", "ss73", "ss74", "ss75", "ss76", "ss77", "ss78", "ss79", "ss80", "ss81", "ss82", "ss83", "ss84", "ss85"]
            pseudo[:PseudoGDP].name     = "GDP Growth Pseudo-observable"
            pseudo[:PseudoGDP].longname = "GDP Growth Pseudo-observable"
            pseudo[:PseudoGDP].rev_transform = loggrowthtopct_annualized_percapita
        end
    end

    if haskey(get_settings(m), :add_pseudo_corepce)
        if get_setting(m, :add_pseudo_corepce) && subspec(m) in ["ss59", "ss60", "ss61", "ss62", "ss63", "ss64", "ss65", "ss66", "ss67", "ss68", "ss69", "ss70", "ss71", "ss72", "ss73", "ss74", "ss75", "ss76", "ss77", "ss78", "ss79", "ss80", "ss81", "ss82", "ss83", "ss84", "ss85"]
            pseudo[:PseudoCorePCE].name     = "Core PCE Pseudo-observable"
            pseudo[:PseudoCorePCE].longname = "Core PCE Pseudo-observable"
            pseudo[:PseudoCorePCE].rev_transform = loggrowthtopct_annualized
        end
    end

    if haskey(get_settings(m), :add_laborproductivity_measurement)
        if get_setting(m, :add_laborproductivity_measurement)
            pseudo[:laborproductivity].name     = "Log Labor Productivity"
            pseudo[:laborproductivity].longname = "Log Labor Productivity"
        end
    end

    if haskey(get_settings(m), :add_laborproductivitygrowth_nome_measurement)
        if get_setting(m, :add_laborproductivitygrowth_nome_measurement)
            pseudo[:LaborProductivityGrowthNoME].name     = "Labor Productivity Growth (No ME)"
            pseudo[:LaborProductivityGrowthNoME].longname = "Labor Productivity Growth (No ME)"
        end
    end

    if haskey(get_settings(m), :add_Epi_t_measurement)
        if get_setting(m, :add_Epi_t_measurement)
            pseudo[:Epi_t].name     = "Short-term Inflation Expectations"
            pseudo[:Epi_t].longname = "Short-term Inflation Expectations"
        end
    end

    # Needed to implement pseudo-measurement equation correctly
    m <= Setting(:forward_looking_pseudo_observables, [:Expected10YearRateGap, :Expected10YearRate, :Expected10YearNaturalRate])

    # Add to model object
    m.pseudo_observable_mappings = pseudo
end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright � 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

function read_map(s)
  data = [strip(line) for line in s]
  height = length(data)
  width = length(data[1])
  start = ones(Int64, 2)

  map = zeros(Int64, width, height)
  for i=1:height
    for j=1:width
      if data[i][j] != STARTPOS
        map[i, j] = convert(Int64, data[i][j])
      else
        start = Any[i, j]
      end
    end
  end

  return (map, start)
end

function define_fitness(trailfile::AbstractString)
  file = joinpath(dirname(@__FILE__), trailfile)
  data = open(file)
  (map, start) = read_map(readlines(data))
  world = World(map, start)

  ex = quote
    function get_fitness(code)
      get_fitness(code, $world, $W_LEN)
    end
  end
  eval(ex)
end

function get_fitness(code, world::World, w_len::Float64)
  f = to_function(code)
  # copy world
  world = copy(world)

  # create ant
  ant = Ant(world.start[1], world.start[2])

  f(world, ant)

  eaten = convert(Float64, ant.eaten)
  return -eaten + W_LEN * length(string(code))
end

ext: jl


text:
using Test

using MarketData
using TimeSeries

using MarketTechnicals

@testset "Momentum" begin

@testset "rsi" begin
    # TTR: RSI(x["Close"])
    # TTR: 45.95363 50.50249 48.97830 48.83191 ...
    @test values(rsi(cl, wilder=true))[1:4] ≈ [45.95363, 50.50249, 48.97830, 48.83191] atol=.01
    # TTR value is 55.84922
    @test values(rsi(cl))[end] ≈ 55.849 atol=.01
    # TTR value is 55.95932
    @test values(rsi(cl, wilder=true))[end] ≈ 55.959 atol=.01
    # TTR value is 56.21947
    @test values(rsi(cl, 10))[end] ≈ 56.219 atol=.01
    @test values(rsi(ohlc))[end, :] ≈ [68.030, 61.872, 70.062, 55.849] atol=.01

    # seed is included (TODO: debate amongst yourselves)
    @test timestamp(rsi(cl))[1]              == Date(2000,1,24)
    @test timestamp(rsi(cl, wilder=true))[1] == Date(2000,1,24)
    @test timestamp(rsi(cl, 10))[1]          == Date(2000,1,18)
    @test timestamp(rsi(ohlc))[1]            == Date(2000,1,24)

    @test timestamp(rsi(cl))[end]              == Date(2001,12,31)
    @test timestamp(rsi(cl, wilder=true))[end] == Date(2001,12,31)
    @test timestamp(rsi(cl, 10))[end]          == Date(2001,12,31)
    @test timestamp(rsi(ohlc))[end]            == Date(2001,12,31)
end

@testset "macd" begin
    @test isapprox(values(macd(cl))[end, 1], -0.020, atol=.01)
    @test isapprox(values(macd(cl))[end, 2], 0.421, atol=.01) # TTR value with percent=FALSE is 0.421175152
    @test isapprox(values(macd(cl))[end, 3], 0.441, atol=.01) # TTR value with percent=FALSE is 0.4414275
    @test timestamp(macd(cl))[end] == Date(2001,12,31)
end

@testset "macd multi-column TimeArray" begin
    # multi-column TimeArray
    # TTR: MACD(..., maType="EMA", percent=0)
    ta = macd(ohlc[:Open, :Close])
    @test colnames(ta)[1:2]  == [:Open_macd, :Close_macd]
    @test timestamp(ta)[end] == Date(2001, 12, 31)
    @test isapprox(values(ta)[end, 3], 0.44254569, atol=.01)    # Open_dif
    @test isapprox(values(ta)[end, 5], 4.536854e-01, atol=.01)  # Open_signal
    @test isapprox(values(ta)[end, 4], 0.421175152, atol=.01)   # Close_dif
    @test isapprox(values(ta)[end, 6], 4.414275e-01, atol=.01)  # Close_signal
end

@testset "chaikinoscillator" begin
    """
    Quote from TTR

    > EMA(adl, 3) - EMA(adl, 10)
        [1]            NA            NA            NA            NA            NA
        [6]            NA            NA            NA            NA  -6851466.867
    [11]  -5508824.158  -4145747.583  -8413321.022 -10025864.828 -10655603.670
    [16]  -8761003.151  -8025814.774  -6995436.368  -6906221.064  -4303730.935
    [21]  -3660985.797  -3451576.257  -2336638.695  -1217362.729    546925.363
    [26]   1788691.016   1349030.229   1738675.716   1239751.751   1946952.206
    [31]   2893048.890   2658896.120   2715328.628   1937437.546   1910085.631
    [36]   2145262.005   1978248.295   1075784.464   1048400.357    814497.124
    [41]   2914695.065   2944188.520   3555908.642   3085073.413   2268377.529
    [46]   1942893.354   1778694.035   1810347.568   1241765.961   -261885.061
    [51]  -1237656.884   -576205.604    613370.430    649457.675   1738353.505
    [56]   3591576.201   2828408.458   1949186.773   1159350.740    626861.789
    [61]    183217.498  -1132617.797   -712700.508   -691386.657    -57610.324
    [66]    774998.562    561961.806   1084391.707    634202.420   -523959.226
    [71]  -1968845.762  -2496129.556  -3167709.694  -2000329.427   -230334.824
    [76]   -577397.260  -1747313.599   -806782.684    617588.662    443152.062
    [81]   1214832.531   1341044.310   1592063.283    936678.383    186056.580
    [86]  -1209996.220  -1516911.058  -1995687.520  -2776094.760  -4090715.405
    """

    ta = chaikinoscillator(ohlcv)
    @test colnames(ta)     == [:chaikinoscillator]
    @test meta(ta)         == meta(ta)
    @test timestamp(ta)[1] == timestamp(ohlcv)[10]
    @test isapprox(values(ta)[1], -6851466.867, atol=.01)
    @test isapprox(values(ta)[2], -5508824.158, atol=.01)
    @test isapprox(values(ta)[3], -4145747.583, atol=.01)
end

@testset "cci" begin
    # TTR::CCI value is -38.931614
    @test isapprox(values(cci(ohlc))[1]  , -38.931614, atol=.01)
    # TTR::CCI value is 46.3511339
    @test isapprox(values(cci(ohlc))[end], 46.3511339, atol=.01)
    @test timestamp(cci(ohlc))[end] == Date(2001, 12, 31)

    ta = cci(TimeArray(collect(Date(2011, 1, 1):Day(1):Date(2011, 1, 30)),
                       fill(42, (30, 4)), [:Open, :High, :Low, :Close]))
    @test all(values(ta) .== 0)
end

@testset "aroon" begin
    """
    Quote from TTR

    > aroon(x[c('High', 'Low')], 25)
            aroonUp aroonDn oscillator
    [24,]      NA      NA         NA
    [25,]      NA      NA         NA
    [26,]      48      28         20
    [27,]      44      24         20
    [28,]      40      20         20
    [29,]      36      16         20
    [30,]      32      12         20
    """

    ta = aroon(ohlc)
    @test colnames(ta)  == [:up, :down, :osc]
    @test timestamp(ta) == timestamp(ohlc[25:end])
    @test isapprox(values(ta)[2, 1], 48)
    @test isapprox(values(ta)[2, 2], 28)
    @test isapprox(values(ta)[2, 3], 20)
end

@testset "roc" begin
    ta = roc(cl, 3)
    @test colnames(ta) == [:Close_roc_3]
    @test isapprox(values(ta)[1], -0.15133107021618722, atol=.01)
    @test isapprox(values(ta)[2], -0.02926829268292683, atol=.01)
    @test isapprox(values(ta)[3], -0.06009615384615385, atol=.01)
end

@testset "adx" begin
    ta = adx(ohlc)
    @test colnames(ta) == [:adx, :dx, :di_plus, :di_minus]
    @test isapprox(values(ta)[1, 1], 10.5998, atol=.01)
    @test isapprox(values(ta)[1, 2], 0.3916,  atol=.01)
    @test isapprox(values(ta)[1, 3], 23.6226, atol=.01)
    @test isapprox(values(ta)[1, 4], 23.4383, atol=.01)
    @test timestamp(ta)[1] == Date(2000, 2, 10)
end

@testset "stochasticoscillator" begin
    """
    Quote from TTR
    > stoch(x[, c("High", "Low", "Close")], maType=SMA)
                fastK      fastD      slowD
    [18,] 0.671428571 0.69466667 0.67441270
    [19,] 0.432000000 0.59342857 0.64901587
    [20,] 0.492857143 0.53209524 0.60673016
    [21,] 0.392857143 0.43923810 0.52158730
    [22,] 0.217586207 0.36776683 0.44636672
    [23,] 0.326296296 0.31224655 0.37308382
    [24,] 0.500000000 0.34796083 0.34265807
    [25,] 0.724444444 0.51691358 0.39237365
    [26,] 0.754814815 0.65975309 0.50820917
    [27,] 0.801061008 0.76010676 0.64559114
    [28,] 0.839964633 0.79861349 0.73949111
    """

    ta = stochasticoscillator(ohlc)
    @test colnames(ta)  == [:fast_k, :fast_d, :slow_d]
    @test timestamp(ta) == timestamp(ohlc[18:end])
    @test isapprox(values(ta)[1, 1], 67.142857, atol=.01)
    @test isapprox(values(ta)[1, 2], 69.466667, atol=.01)
    @test isapprox(values(ta)[1, 3], 67.441270, atol=.01)
end

end  # @testset "Momentum"


text:
import Statistics: mean, var, cov

"""
    Product <: MultivariateDistribution

An N dimensional `MultivariateDistribution` constructed from a vector of N independent
`UnivariateDistribution`s.

```julia
Product(Uniform.(rand(10), 1)) # A 10-dimensional Product from 10 independent `Uniform` distributions.
```
"""
struct Product{
    S<:ValueSupport,
    T<:UnivariateDistribution{S},
    V<:AbstractVector{T},
} <: MultivariateDistribution{S}
    v::V
    function Product(v::V) where
        V<:AbstractVector{T} where
        T<:UnivariateDistribution{S} where
        S<:ValueSupport
        return new{S, T, V}(v)
    end
end

length(d::Product) = length(d.v)
function Base.eltype(::Type{<:Product{S,T}}) where {S<:ValueSupport,
                                                    T<:UnivariateDistribution{S}}
    return eltype(T)
end

_rand!(rng::AbstractRNG, d::Product, x::AbstractVector{<:Real}) =
    broadcast!(dn->rand(rng, dn), x, d.v)
_logpdf(d::Product, x::AbstractVector{<:Real}) =
    sum(n->logpdf(d.v[n], x[n]), 1:length(d))

mean(d::Product) = mean.(d.v)
var(d::Product) = var.(d.v)
cov(d::Product) = Diagonal(var(d))
entropy(d::Product) = sum(entropy, d.v)
insupport(d::Product, x::AbstractVector) = all(insupport.(d.v, x))

"""
    product_distribution(dists::AbstractVector{<:UnivariateDistribution})

Creates a multivariate product distribution `P` from a vector of univariate distributions.
Fallback is the `Product constructor`, but specialized methods can be defined
for distributions with a special multivariate product.
"""
function product_distribution(dists::AbstractVector{<:UnivariateDistribution})
    return Product(dists)
end

"""
    product_distribution(dists::AbstractVector{<:Normal})

Computes the multivariate Normal distribution obtained by stacking the univariate
normal distributions. The result is a multivariate Gaussian with a diagonal
covariance matrix.
"""
function product_distribution(dists::AbstractVector{<:Normal})
    µ = mean.(dists)
    σ = std.(dists)
    return MvNormal(µ, σ)
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   ModiaMath.SimulationEngine(ModiaMath/SimulationEngine/_module.jl)
#

mutable struct IntegratorData
    model::ModiaMath.AbstractSimulationModel    # Model specific data structure
    simulationState::ModiaMath.SimulationState
    statistics::ModiaMath.SimulationStatistics   # Simulation statistics
    z::Vector{Float64}          # Vector of event indicators (zero crossings). If one of z[i] passes
                                # zero, that is beforeEvent(z[i])*z[i] < 0, an event is triggered
                                # (allocated during instanciation according to nz).
    zDir::Vector{Cint}          # zDir[i] =  0: Root is reported for both crossing directions
                                #         =  1: Root is reported when crossing from negative to positive direction
                                #         = -1: Root is reported when crossing from positive to negative direction

    eweight::Sundials.N_Vector   # Weight vector
    last_t::Float64              # Time instant of the last call to idasol_f (residue function)
    last_norm_y::Float64         # Norm of y of the last call to idasol_f (residue function)
    hcur::MVector{1,Cdouble}     # Current step size argument for IDAGetCurrentStep
    order::MVector{1,Cint}       # Current order for IDAGetCurrentOrder
    r::Vector{Sundials.realtype} # Residues vector used in idasol_f function
    ida_mem::Ptr{Nothing}        # IDA pointer (to access all IDAgetXXX functions)

    y::Vector{Sundials.realtype}    # Julia vector of y wrapping y_N_Vector vector
    yp::Vector{Sundials.realtype}   # Julia vector of yp wrapping yp_N_Vector vector
    fulljac::Union{Matrix{Sundials.realtype}, Nothing}  # Julia matrix of fulljac wrapping fulljac_N_xxx

    function IntegratorData(model::ModiaMath.AbstractSimulationModel, simulationState::ModiaMath.SimulationState)
        ny = simulationState.nx
        nz = simulationState.nz
        new(model, simulationState, simulationState.statistics,
          ones(nz), fill(0, 0), Sundials.N_VNew_Serial(ny), typemin(Float64), 0.0,
          MVector{1,Cdouble}(0.0), MVector{1,Cint}(0), zeros(ny))
    end
end

include("sundials_ida_UserData.jl")

function updateIDAstatistics!(mem, stat::ModiaMath.SimulationStatistics)
    info = Clong[0]
    Sundials.IDAGetNumSteps(mem, info)
    stat.nSteps        += info[1]

    stat.sparseSolver ? Sundials.IDASlsGetNumJacEvals(mem, info) : Sundials.IDADlsGetNumJacEvals(mem, info)
    stat.nJac          += info[1]

    Sundials.IDAGetNumErrTestFails(mem, info)
    stat.nErrTestFails += info[1]

    h0 = [0.0]
    Sundials.IDAGetActualInitStep(mem, h0)
    stat.h0 = min(stat.h0, h0[1])
end

"""
    result = closeTimePoints(t1,t2,epsilon)

The function returns true if the two time points t1 and t2 are close together (with respect to epsilon)
"""
closeTimePoints(t1::Float64, t2::Float64, epsilon::Float64) = abs(t1 - t2) / max(abs(t2), 10 * epsilon) <= 100 * epsilon

"""
    ModiaMath.simulate!(simulationModel; log=false, startTime=NaN, stopTime=NaN,
                                         tolerance=NaN, interval=NaN,
                                         maxStepSize=NaN, maxNumberOfSteps=missing)

Simulates a DAE `simulationModel` that is defined with package `Modia`, package `Modia3D` or
with the `ModiaMath.@component` macro. The DAE is mathematically described as
implicit, index 1 DAE with events (containing an ODE or a semi-explicit
index 1 DAE as special cases). For details see [`ModiaMath.StructureOfDAE`](@ref).

During continuous integration, the DAE is solved with the
[Sundials](https://computation.llnl.gov/projects/sundials) IDA solver
(accessed via the Julia [Sundials](https://github.com/JuliaDiffEq/Sundials.jl) interface package).

Input arguments

- `simulationModel::ModiaMath.AbstractSimulationModel`: Model struct (generated with Modia, Modia3D or `ModiaMath.@component`).
- `log::Bool`: = true, if logging is enabled, otherwise it is disabled.
- `startTime::Float64`: Start time of the simulation in [s].
                        If startTime=NaN, the default startTime is used that is defined by the `simulationModel`.
- `stopTime::Float64`: Stop time of the simulation in [s].
                       If stopTime=NaN, the default stopTime is used that is defined by the `simulationModel`.
- `tolerance::Float64`: The relative tolerance for the integration.
  The absolute tolerance is computed as `0.1*tolerance*nominal(variable)` where
  `nominal(variable)` is the nominal value of the variable.
   If tolerance=NaN, the default tolerance is used that is defined by the `simulationModel`.
- `interval::Float64`: Output interval for results in [s]. If events occur, the event time instants are
  additionally added to the result.
  If interval=NaN, the default interval is used that is defined by the `simulationModel`.
- `maxStepSize::Float64`: Maximum absolute value of the step size.
  If `maxStepSize=NaN`, the default maximum step size of the integrator is used.
- `maxNumberOfSteps::Union{Int,Missing}`: Maximum number of
  steps to be taken by the solver in its attempt to reach the next output time.
  If `maxNumberOfSteps=missing`, the default value of the integrator is used (= 500).
"""
function simulate!(model::ModiaMath.AbstractSimulationModel;
                   tolerance=NaN,
                   startTime=NaN,
                   stopTime=NaN,
                   interval=NaN,
                   maxStepSize=NaN,
                   maxNumberOfSteps=missing,
                   log::Bool=false,
                   KLUorderingChoice::Int=1)

    use_fulljac = false

    sim           = model.simulationState
    sim.model     = model
    sim.tolerance = isnan(tolerance) ? sim.defaultTolerance : convert(Float64, tolerance)
    sim.startTime = isnan(startTime) ? sim.defaultStartTime : convert(Float64, startTime)
    sim.stopTime  = isnan(stopTime)  ? sim.defaultStopTime  : convert(Float64, stopTime)
    sim.interval  = isnan(interval)  ? (isnan(sim.defaultInterval) ? (sim.stopTime - sim.startTime)/500.0 : sim.defaultInterval)  : convert(Float64, interval)
    sim.time      = sim.startTime
    @assert(sim.stopTime >= sim.startTime)
    @assert(sim.interval >= 0.0)
    @assert(sim.tolerance > 0.0)

    stopTime2::Float64  = sim.stopTime
    startTime2::Float64 = sim.startTime
    tolerance::Float64  = sim.tolerance
    interval::Float64   = sim.interval
    logger              = sim.logger
    ModiaMath.setLog!(logger, log)

    if ModiaMath.isLogInfos(logger)
        println("... ModiaMath.simulate! (version ", ModiaMath.Version, " ", ModiaMath.Date, ") to simulate model: ", sim.name)
    end

    # Start timing measure
    cpuStart::UInt64 = time_ns()
    cpuLast::UInt64  = cpuStart
    cpuStartIntegration::UInt64 = cpuStart

    # Init variables
    nt::Int           = round(Int64, 1 + (stopTime2 - startTime2) / interval)
    t0::Float64       = startTime2
    tEnd::Float64     = stopTime2
    tReached::Float64 = t0
    epsilon::Float64  = eps()

    # Initialize simulation model and store result after initialization
    statistics = sim.statistics
    ModiaMath.reInitializeStatistics!(statistics, t0, stopTime2, interval, tolerance)

    if ModiaMath.isLogInfos(logger)
        println("      Initialization at time = ", t0, " s")
    end

    init = DAE.initialize!(model, sim, t0, nt, tolerance)
    @assert(length(init.y0)  == sim.nx)
    @assert(length(init.yp0) == sim.nx)
    @assert(tolerance > 0.0)

    if init.terminate || tEnd <= t0
        statistics.h0   = 0.0
        statistics.hMin = 0.0
        @goto TerminateSimulation
    end

    if ModiaMath.isLogInfos(logger)
        println("      Simulation started")
    end

    maxTime::Float64          = init.maxTime
    nextEventTime::Float64    = init.nextEventTime
    integrateToEvent::Bool = init.integrateToEvent

    if nextEventTime < t0
        nextEventTime = typemax(Float64)
    end

    # Initialize auxiliary variables
    simModel = IntegratorData(model, sim)
    stateEvent::Bool = false
    timeEvent::Bool  = false
    flag::Cint = 0
    ny = sim.nx
    nz = sim.nz

    y  = copy(init.y0)
    yp = copy(init.yp0)
    simModel.y  = y
    simModel.yp = yp
    if use_fulljac
        fulljac = zeros(sim.nx, sim.nx)
        simModel.fulljac = fulljac
    else
        simModel.fulljac = nothing
    end

    eventInfo = DAE.EventInfo()
    tret    = [0.0]
    jc::Int = 0  # communication point counter
    if ModiaMath.isLogInfos(logger) && nz > 0
        rootInfo = fill(Cint(0), nz)
    end

    # Create IDA
    mem = Sundials.IDACreate()
    if mem == C_NULL
        error("ModiaMath.simulate: Failed to allocate IDA solver object")
    end
    simModel.ida_mem = mem

    # Allocate N_Vector storage for y and yp
    y_N_Vector  = Sundials.N_VMake_Serial(ny, pointer(y))
    yp_N_Vector = Sundials.N_VMake_Serial(ny, pointer(yp))

    # Run IDA
    try
        # Set error handler function
        Sundials.IDASetErrHandlerFn(mem, old_cfunction(idasol_ErrHandlerFn, Nothing, Tuple{Cint,Cstring,Cstring,Cstring,Ref{typeof(IntegratorData)}}), pointer_from_objref(simModel))

        # Initialize IDA
        Sundials.IDAInit(mem, old_cfunction(idasol_f, Cint, Tuple{Sundials.realtype,Sundials.N_Vector,Sundials.N_Vector,Sundials.N_Vector,Ref{typeof(IntegratorData)}}), t0, y_N_Vector, yp_N_Vector)
        if use_fulljac
            # IDADlsJacFn
            Sundials.IDADlsSetJacFn(mem, old_cfunction(idasol_fulljac, Cint, Tuple{Sundials.realtype,Sundials.realtype,Sundials.N_Vector,Sundials.N_Vector,
                                                       Sundials.N_Vector,Sundials.SUNMatrix, Ref{typeof(IntegratorData)},Sundials.N_Vector,Sundials.N_Vector,Sundials.N_Vector}))
        end

        #IDASStolerances(mem, tolerance, 0.1*tolerance)
        tolAbs = 0.1 * tolerance * init.y_nominal
        for i in 1:ny
            if !init.y_errorControl[i]
                tolAbs[i] = 1e5 * init.y_nominal[i]   # switch tolerance control off
            end
        end

        Sundials.IDASVtolerances(mem, tolerance, tolAbs)
        Sundials.IDASetUserData(mem, simModel)

        if !isnan(maxStepSize)
            Sundials.IDASetMaxStep(mem, maxStepSize)
            sim.maxStepSize = maxStepSize
        elseif !isnan(sim.defaultMaxStepSize)
            Sundials.IDASetMaxStep(mem, sim.defaultMaxStepSize)
            sim.maxStepSize = sim.defaultMaxStepSize
        end

        if typeof(maxNumberOfSteps) != Missing
            Sundials.IDASetMaxNumSteps(mem, maxNumberOfSteps)
            sim.maxNumberOfSteps = maxNumberOfSteps
        elseif typeof(sim.defaultMaxNumberOfSteps) != Missing
            Sundials.IDASetMaxNumSteps(mem, sim.defaultMaxNumberOfSteps)
            sim.maxNumberOfSteps = sim.defaultMaxNumberOfSteps
        end

        # Sundials.IDASetId(mem, y_states_nvector)
        #if sim.sparse
        #   nnz_jac = nnz(sim.jac)
        #   Sundials.IDAKLU(mem, ny, nnz_jac)
        #   Sundials.IDAKLUSetOrdering(mem, KLUorderingChoice)
        #   Sundials.IDASlsSetSparseJacFn(mem);
        #else
        A  = Sundials.SUNDenseMatrix(ny, ny)
        LS = Sundials.SUNDenseLinearSolver(y, A)
        Sundials.IDADlsSetLinearSolver(mem, LS, A)
        #end

        # Initialize zero crossing function, if required
        hasZeroCross = nz > 0
        if hasZeroCross
            Sundials.IDARootInit(mem, nz, old_cfunction(idasol_g, Cint, Tuple{Sundials.realtype,Sundials.N_Vector,Sundials.N_Vector,Ptr{Sundials.realtype},Ref{typeof(IntegratorData)}}))
            Sundials.IDASetRootDirection(mem, sim.zDir)
        end

        # Initialize event variables
        tStop::Float64 = tEnd
        tNext = t0
        tc    = t0
        nTimeEvents = 0
        restart     = ModiaMath.NoRestart

        # Check for event at time = t0. If yes, trigger event and store result after event
        if closeTimePoints(t0, nextEventTime, epsilon)
            tReached = t0
            if ModiaMath.isLogEvents(logger)
                println("\n      Time event at time = $tReached s")
            end

            DAE.reset!(eventInfo)
            DAE.processEvent!(model, sim, tReached, y, yp, eventInfo)
            restart          = eventInfo.restart
            maxTime          = eventInfo.maxTime
            nextEventTime    = eventInfo.nextEventTime
            integrateToEvent = eventInfo.integrateToEvent

            statistics.nTimeEvents += 1
            if ModiaMath.isLogEvents(logger)
                println("        restart = ", restart)
            end

            # Reinitialize IDA if required by model
            if restart == ModiaMath.Restart
                statistics.nRestartEvents += 1
                Sundials.__IDAReInit(mem, tReached, y_N_Vector, yp_N_Vector);
            end
        end

        # Loop over communication points
        cpuStartIntegration = time_ns()

        while tReached < tEnd && restart != ModiaMath.Terminate  #------------------ while --------------
            if tReached >= tc
                # Last communication point reached -> increment communication point
                jc += 1
                tc = startTime2 + jc * interval
            end

            if closeTimePoints(nextEventTime, tc, epsilon)
                nextEventTime = tc
            end

            if nextEventTime <= tc # next tStop is a time event
                event = true
                tNext = nextEventTime
                tStop = integrateToEvent ? nextEventTime : min(tEnd, maxTime)
            else        # next tStop is a communication point (no event)
                event = false
                tNext = tc
                tStop = min(tEnd, maxTime)
            end

            # Integrate in direction of tStop
            #println("... tReached = ", tReached, ", tNext = ", tNext, ", tStop = ", tStop,
            #        ", y = ", y[1], ", yp = ", yp[1])
            Sundials.IDASetStopTime(mem, tStop)
            #flag = Sundials.IDASolve(mem, tNext, tret, simModel.y, simModel.yp, Sundials.IDA_NORMAL)
            flag = Sundials.__IDASolve(mem, tNext, tret, y_N_Vector, yp_N_Vector, Sundials.IDA_NORMAL)
            tReached = tret[1]
            #println("    tReached = ", tReached)

            stateEvent = flag == Sundials.IDA_ROOT_RETURN
            timeEvent  = event && closeTimePoints(tNext, tReached, epsilon)
            isEvent    = timeEvent || stateEvent

            # Store result point at tReached
            DAE.computeAndStoreResult!(model, sim, tReached, y, yp)

            # If event at tReached, handle the event
            if isEvent
                if ModiaMath.isLogEvents(logger)
                    if timeEvent && stateEvent
                        print("\n      Time and state (zero-crossing) event at time = $tReached s")
                    elseif timeEvent
                        println("\n      Time event at time = $tReached s")
                    elseif stateEvent
                        print("\n      State event (zero-crossing) at time = $tReached s")
                    end

                    if stateEvent
                       # Print information about the root
                        Sundials.IDAGetRootInfo(mem, rootInfo)
                        print(" (")
                        firstRoot = true
                        for iroot in eachindex(rootInfo)
                            if rootInfo[iroot] != 0
                                if firstRoot
                                    firstRoot = false
                                else
                                    print(", ")
                                end
                                print("z[", iroot, "]", rootInfo[iroot] > 0 ? " > 0" : " < 0")
                            end
                        end
                        print(")\n")
                    end
                end

                # Trigger event
                DAE.reset!(eventInfo)
                DAE.processEvent!(model, sim, tReached, y, yp, eventInfo)
                restart          = eventInfo.restart
                maxTime          = eventInfo.maxTime
                nextEventTime    = eventInfo.nextEventTime
                integrateToEvent = eventInfo.integrateToEvent

                if timeEvent
                    statistics.nTimeEvents += 1
                end

                if stateEvent
                    statistics.nStateEvents += 1
                end

                if restart == ModiaMath.Restart || restart == ModiaMath.FullRestart
                    statistics.nRestartEvents += 1
                end

                if ModiaMath.isLogEvents(logger)
                    println("        restart = ", restart)
                end

                # Handle restart flag
                if restart == ModiaMath.Restart
                    # Restart (dimensions do not change)
                    updateIDAstatistics!(mem, statistics)
                    Sundials.__IDAReInit(mem, tReached, y_N_Vector, yp_N_Vector);
                elseif restart == ModiaMath.FullRestart
                    error("FullRestart not yet supported")
                end
            end

            # print every 5 s the time in order that the user sees the progress
            if ModiaMath.isLogProgress(logger)
                cpuNew = time_ns()
                if (cpuNew - cpuLast) * 1e-9 > 5.0
                    cpuLast = cpuNew
                    @printf("      progress: integrated up to time = %.2g s\n", tReached)
                end
            end
        end #---------------------- end while ------------------------------------

        # Finalize statistics for simulation run
        updateIDAstatistics!(mem, statistics)
    finally
        # Free allocated memory
        Sundials.__N_VDestroy_Serial(y_N_Vector)
        Sundials.__N_VDestroy_Serial(yp_N_Vector)
        Sundials.IDAFree([mem])
    end

    @label TerminateSimulation
    if ModiaMath.isLogInfos(logger)
        println("\n      Simulation is terminated at time = ", tEnd, " s")
    end

    result = DAE.terminate!(model, sim, tReached, y, yp)
    statistics.cpuTimeInitialization = (cpuStartIntegration - cpuStart) * 1e-9
    statistics.cpuTimeIntegration    = (time_ns() - cpuStartIntegration) * 1e-9

    if ModiaMath.isLogStatistics(logger)
        println("\n      Statistics (get help with ?ModiaMath.SimulationStatistics):")
        display(statistics)
    end
    return result
end


text:
############################################################
###################### QMDP Solver #########################
############################################################

#=
POMDP Model Requirements:
    create_transition_distribution(pomdp)
    states(pomdp)
    actions(pomdp)
    iterator(space)
    transition(pomdp, s, a, dist)
    reward(pomdp, s, a, sp)
    pdf(dist, sp)
    discount(pomdp)
    stateindex(pomdp, sp)

    The user must implement the above functions to use QMDP.
=#

mutable struct QMDPSolver <: Solver
    max_iterations::Int64
    tolerance::Float64
    verbose::Bool
end
function QMDPSolver(;max_iterations::Int64=100, tolerance::Float64=1e-3, verbose=false)
    return QMDPSolver(max_iterations, tolerance, verbose)
end

@POMDP_require solve(solver::QMDPSolver, pomdp::POMDP) begin
    vi_solver = ValueIterationSolver(max_iterations=solver.max_iterations, belres=solver.tolerance)
    @subreq solve(vi_solver, pomdp)
end

function solve(solver::QMDPSolver, pomdp::POMDP; kwargs...)
    # deprecation warning - can be removed when Julia 1.0 is adopted
    if !isempty(kwargs)
        @warn "Keyword args for solve(::QMDPSolver, ::POMDP) are no longer supported. For verbose output, use the verbose option in the ValueIterationSolver"
    end
    vi_solver = ValueIterationSolver(max_iterations=solver.max_iterations, belres=solver.tolerance, verbose=solver.verbose, include_Q=true)
    vi_policy = solve(vi_solver, pomdp)

    return AlphaVectorPolicy(pomdp, vi_policy.qmat)
end


text:
module TensorAD
export DiffTensor

using OMEinsum, LinearAlgebra

const ADTypes = Union{Float32, Float64, ComplexF64, ComplexF32}

struct PartialBack
    input_ids::Union{UInt,Tuple}
    output_ids::Union{UInt,Tuple}
    back
end
struct Instruction
    info::String
    backs::NTuple{N,PartialBack} where N
end
struct Tape
    instructs::Vector{Instruction}
end
const GLOBAL_TAPE = Tape(Instruction[])

# to make objectid work consistently
mutable struct DiffTensor{T<:ADTypes,N,AT<:AbstractArray{T,N}} <: AbstractArray{T,N}
    data::AT
    requires_grad::Bool
    function DiffTensor(data::AT, requires_grad::Bool=true) where {T,N,AT<:AbstractArray{T,N}}
        if AT <: DiffTensor
            error("DiffTensor in DiffTensor is forbidden to prevent errors.")
        end
        new{T,N,AT}(data, requires_grad)
    end
end

function debug_info(f, args...; kwargs...)
    args = join(map(arg->"::$(typeof(arg))", args), ", ")
    if length(kwargs) != 0
        kwargs = join(["$k=$v" for (k,v) in kwargs], "")
        "∂$f($args; $kwargs)"
    else
        "∂$f($args)"
    end
end
Instruction(info::String, y::Union{DiffTensor, Tuple}, backs::Pair{<:Union{DiffTensor, NTuple{N,DiffTensor} where N}}...) = Instruction(info, PartialBack.(Ref(y), backs))
function PartialBack(y::Union{DiffTensor, NTuple}, x::Pair)
    PartialBack(getid(x.first), getid(y), x.second)
end

getid(x::DiffTensor) = objectid(x)
getid(x::Tuple) = getid.(x)
Base.zero(x::DiffTensor) = DiffTensor(zero(x.data), false)
Base.size(x::DiffTensor, indices...) = Base.size(x.data, indices...)
getdata(x::DiffTensor) = x.data
getgrad(d::AbstractDict, x::DiffTensor) = getgrad(d, getid(x), ()->zero(x))
function getgrad(d::AbstractDict, key::UInt, default)
    if haskey(d, key)
        return d[key]
    else
        return default()
    end
end
getgrad(d::AbstractDict, key::Tuple, default) = getgrad(Ref(d), key, Ref(default))
getgrad(d::AbstractDict, key::Tuple) = getgrad.(Ref(d), key)
function accumulate_gradient!(grad_storage, key::UInt, g)
    # accumulate or create
    if !haskey(grad_storage, key)
        grad_storage[key] = g
    else
        grad_storage[key] = grad_storage[key] + g
    end
end

function accumulate_gradient!(grad_storage, keys::Tuple, grads)
    for (t, g) in zip(keys, grads)
        accumulate_gradient!(grad_storage, t, g)
    end
end

Base.show(io::IO, ::MIME"text/plain", x::DiffTensor) = Base.show(io, x)
function Base.show(io::IO, x::DiffTensor)
    sz = join(string.(size(x)), "×")
    s = "$(typeof(x))[$sz]"
    if x.requires_grad
        s *= " ✓"
    end
    print(io, s)
end

Base.show(io::IO, ::MIME"text/plain", tape::Tape) = Base.show(io, tape)
function Base.show(io::IO, tape::Tape)
    for instruct in tape.instructs
        println(io, instruct.info)
    end
end

OMEinsum.asarray(x::Number, y::DiffTensor) = DiffTensor(asarray(x,y.data), y.requires_grad)
OMEinsum.asarray(x::DiffTensor, y::DiffTensor) = x  # fix ambiguity error

# does not return gradients in kwargs
# always return a tuple
function gradient(f, args...; kwargs...)
    empty!(GLOBAL_TAPE.instructs)
    y = f(args...; kwargs...)
    if ndims(y) != 0 || !(eltype(y) <: Real)
        @warn "differentiating a tensor not a scalar real number, got eltype `$(eltype(y))` and rank `$(ndims(y))`"
    end
    grad_storage = Dict{UInt,Any}()
    accumulate_gradient!(grad_storage, getid(y), DiffTensor(ones(eltype(y), size(y.data)...), false))
    back!(GLOBAL_TAPE, grad_storage)
    return getgrad(grad_storage, args)
end

function back!(y::Tape, grad_storage::AbstractDict)
    for i=length(y.instructs):-1:1
        instruct = y.instructs[i]
        debug_info, backs = instruct.info, instruct.backs
        @debug "$debug_info"
        for pb in backs
            dy = getgrad(grad_storage, pb.output_ids, ()->nothing)
            if dy !== nothing
                grads = pb.back(dy)
                accumulate_gradient!(grad_storage, pb.input_ids, grads)
            else
                @debug "can not find requested gradient for $pb"
            end
        end
    end
end

function hessian(f, x::DiffTensor)
    return jacobian(x->gradient(f, x)[1], x)
end

function jacobian(f, x::DiffTensor{T}) where T
    empty!(GLOBAL_TAPE.instructs)
    slices = DiffTensor{T,1}[]
    y = f(x)
    @debug GLOBAL_TAPE
    for i=1:length(y)
        @debug "jacobian, row $i"
        grad_storage = Dict{UInt,Any}()
        gy = zero(y)
        gy.data[i] = one(T)
        accumulate_gradient!(grad_storage, getid(y), gy)
        back!(GLOBAL_TAPE, grad_storage)
        push!(slices, vec(getgrad(grad_storage, x)))
    end
    return transpose(cat(slices...; dims=2))
end

include("rules.jl")
end


text:
#-----------------------------------------------------------------------------------------------------
# Measurement
#-----------------------------------------------------------------------------------------------------
export measure
function measure(
    o::Operation,
    v::AbstractVecOrMat
)
    dot(v, o * v)
end
#-----------------------------------------------------------------------------------------------------
function measure(
    o1::Operation,
    o2::Operation,
    v::AbstractVecOrMat
)
    dot(o1 * v, o2 * v)
end
#-----------------------------------------------------------------------------------------------------
export covmat
function covmat(
    ol::Vector{<:Operation},
    v::AbstractVector
)
    n = length(ol)
    am = Vector{Float64}(undef, n)
    cm = Matrix{Float64}(undef, n, n)
    for i=1:n
        am[i] = real(measure(ol[i], v))
    end
    for i=1:n
        for j=i:n
            cm[i,j] = real(measure(ol[i], ol[j], v)) - am[i] * am[j]
        end
    end
    Hermitian(cm)
end
#-----------------------------------------------------------------------------------------------------
function covmat(
    ol::Vector{<:Operation},
    v::AbstractMatrix
)
    n = length(ol)
    am = Vector{Float64}(undef, n)
    cm = Matrix{Float64}(undef, n, n)
    num = size(v, 2)
    for i=1:n
        am[i] = real(measure(ol[i], v))/num
    end
    for i=1:n
        for j=i:n
            cm[i,j] = real(measure(ol[i], ol[j], v))/num - am[i] * am[j]
        end
    end
    Hermitian(cm)
end
#-----------------------------------------------------------------------------------------------------
export covmatc
function covmatc(
    ol::Vector{<:Operation},
    v::AbstractVector
)
    n = length(ol)
    am = Vector{ComplexF64}(undef, n)
    cm = Matrix{ComplexF64}(undef, n, n)
    for i=1:n
        am[i] = measure(ol[i], v)
    end
    for i=1:n
        for j=i:n
            cm[i,j] = measure(ol[i], ol[j], v) - conj(am[i]) * am[j]
        end
    end
    Hermitian(cm)
end
#-----------------------------------------------------------------------------------------------------
function covmatc(
    ol::Vector{<:Operation},
    v::AbstractMatrix
)
    n = length(ol)
    am = Vector{ComplexF64}(undef, n)
    cm = Matrix{ComplexF64}(undef, n, n)
    num = size(v, 2)
    for i=1:n
        am[i] = measure(ol[i], v)/num
    end
    for i=1:n
        for j=i:n
            cm[i,j] = measure(ol[i], ol[j], v)/num - conj(am[i]) * am[j]
        end
    end
    Hermitian(cm)
end
#-----------------------------------------------------------------------------------------------------
# Entropy
#-----------------------------------------------------------------------------------------------------
export entropy
function entropy(
    v::AbstractVector,
    l::Integer,
    i::Integer;
    base=0,
    cutoff=1e-7
)
    base = base==0 ? round(Int64, length(v)^(1/l)) : base
    m = reshape(v, base^i, :)
    s = svdvals(m) .^ 2
    - sum(si * log(si) for si in s if si > cutoff)
end


text:
using StringEncodings

include("build_ReadStat.v0.1.1-dev.jl")


text:
using Documenter, BooSTjl

makedocs()

max_stars_repo_name: gabrielrvsc/BooSTjl

max_issues_repo_name: gabrielrvsc/BooSTjl

max_forks_repo_name: gabrielrvsc/BooSTjl


text:
"""
    encode(val, source_type)

Encodes `val::T` that comes from
a struct with declared type `source_type`.

Most of the time, `T == source_type`.
But fields with abstract or union types will
display a concrete `T` type and an abstract `source_type`.
When this happens, it is necessary to pass
type information to the encoded value.
"""
function encode end

"""
    decode(val, target_type, module)

Decodes `val::T` to a field with declared type `target_type`.

`val` is the value stored in the BSON.
"""
function decode end

"""
    encode_type(T)

Reports the type used to encode a type `T` to BSON.
"""
function encode_type end

#
# Native types supported by BSON
#

for tt in (String, Int32, Int64, DateTime, Float64, Bool, BSONObjectId)
    @eval begin
            function encode(val::$tt, ::Type{$tt})
                val
            end

            function decode(val::$tt, ::Type{$tt})
                val
            end

            function encode_type(::Type{$tt})
                $tt
            end
    end
end

#
# Integer numbers smaller than 32bits are encoded as Int32
#

for tt in (UInt8, UInt16, Int8, Int16, UInt32)
    @eval begin
            function encode(val::$tt, ::Type{$tt})
                Int32(val)
            end

            function decode(val::Int32, ::Type{$tt})
                ($tt)(val)
            end

            function encode_type(::Type{$tt})
                Int32
            end
    end
end

#
# UInt64 is reinterpreted into Int64
#
function encode(val::UInt64, ::Type{UInt64})
    reinterpret(Int64, val)
end

function decode(val::Int64, ::Type{UInt64})
    reinterpret(UInt64, val)
end

function encode_type(::Type{UInt64})
    Int64
end

#
# Integers can be decoded to Float64
#

function decode(val::Integer, ::Type{Float64})
    Float64(val)
end

#
# Int32 can be decoded as Int64
#

function decode(val::Int32, ::Type{Int64})
    Int64(val)
end

#
# Date is encoded as DateTime with zeroed Time
#

encode(val::Date, ::Type{Date}) = DateTime(val)

decode(val::DateTime, ::Type{Date}) = Date(val)

# accept String with standard format "yyyy-mm-dd"
decode(val::String, ::Type{Date}) = Date(val)

encode_type(::Type{Date}) = DateTime

#
# Vectors are encoded as BSON vectors with encoded values
#

function encode(val::Vector{T}, ::Type{Vector{T}}) where {T}
    [ encode(x, T) for x in val ]
end

function decode(val::Array, ::Type{Vector{T}}) where {T}
    T[ decode(x, T) for x in val ]
end

# in case type information is not available, assume Vector
function decode(val::Array, ::Type{Array})
    T = eltype(val)
    T[ decode(x, T) for x in val ]
end

function encode_type(val::Vector{T}) where {T}
    Vector{encode_type(T)}
end

#
# Symbols are encoded as strings
#

encode(val::Symbol, ::Type{Symbol}) = String(val)
decode(val::String, ::Type{Symbol}) = Symbol(val)
encode_type(::Type{Symbol}) = String

#
# Char is encoded as String
#
encode(val::Char, ::Type{Char}) = string(val)
decode(val::String, ::Type{Char}) = val[1]
encode_type(::Type{Char}) = String

#
# Nothing / Missing
#
#=
function decode(val::T, ::Type{Union{Nothing, T}}) where {T}
    println("called decode on union{Nothing, T}")
    decode(val, T)
end
=#

#
# DatePeriod and TimePeriod are encoded as Int
#

for tt in union(InteractiveUtils.subtypes(DatePeriod), InteractiveUtils.subtypes(TimePeriod))
    @eval begin
            function encode(val::$tt, ::Type{$tt})
                val.value
            end

            function decode(val::Int64, ::Type{$tt})
                ($tt)(val)
            end

            function decode(val::Int32, ::Type{$tt})
                ($tt)(val)
            end

            function encode_type(::Type{$tt})
                Int
            end
    end
end

#
# Dictionaries with String as symbols are natively supported by BSON.
# Edge case might fail: when V is an abstract type.
#

encode_dict_key(key::String) = key
decode_dict_key(key::String, ::Type{String}) = key

encode_dict_key(key::Symbol) = String(key)
decode_dict_key(key::String, ::Type{Symbol}) = Symbol(key)

encode_dict_key(key::Integer) = string(Int(key))
function decode_dict_key(key::String, ::Type{T}) where {T<:Integer}
    T(parse(Int, key))
end

function encode(val::Dict{K,V}, ::Type{Dict{K,V}}) where {K,V}
    encoded_dict = Dict{String, encode_type(V)}()
    for (k,v) in val
        encoded_dict[encode_dict_key(k)] = encode(v, V)
    end
    return encoded_dict
end

function decode(val::Dict, ::Type{Dict{K,V}}) where {K,V}
    decoded_dict = Dict{K,V}()
    for (k,v) in val
        decoded_dict[decode_dict_key(k, K)] = decode(val[k], V)
    end
    return decoded_dict
end

# encode call comming from struct field with abstract type
function encode(val::T, ::Type{A}) where {T, A}
    @assert T != A "encode method for type $T was not provided."
    return BSON("type" => typepathref(T), "value" => encode(val, T))
end

function decode(val::T, ::Type{A}) where {T<:Union{BSONSerializer.BSON, Dict}, A}
    if haskey(val, "value")
        return decode(val["value"], resolve_typepath(val["type"]))
    elseif haskey(val, "args")
        return decode(val, resolve_typepath(val["type"]))
    elseif T == A
        error("decode method for type $T was not provided.")
    else
        error("Can't decode from value of type $T to target type $A.")
    end
end

#
# encode Functions as singletons
#

function encode(f::Function, ::Type{T}) where {T<:Function}
    return typepathref(typeof(f))
end

function decode(val::Vector, ::Type{T}) where {T<:Function}
    datatype = resolve_typepath(val)
    @assert isdefined(datatype, :instance)
    return datatype.instance
end

function encode_type(::Type{T}) where {T<:Function}
    Vector{Any}
end


text:
include("./Find.jl")
main()


text:
using Test
using NeXLCore
using DataFrames

@testset "KRatio" begin
    kr = KRatio(
        [n"Fe K-L2", n"Fe K-L3"],
        Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
        Dict(:BeamEnergy => 15.0e3, :TakeOffAngle => deg2rad(38.0)),
        mat"Fe2O3",
        uv(0.77, 0.01),
    )
    @test n"Fe K-L2" in kr.xrays
    @test !(n"Ca K-L2" in kr.xrays)
    @test kr.unkProps[:BeamEnergy] == 20.0e3
    @test kr.unkProps[:TakeOffAngle] == deg2rad(40.0)
    @test kr.stdProps[:BeamEnergy] == 15.0e3
    @test kr.stdProps[:TakeOffAngle] == deg2rad(38.0)
    @test isapprox(kr.standard[n"Fe"], 0.6994315143211048, atol = 1e-8)
    @test value(kr.kratio) == 0.77
    @test σ(kr.kratio) == 0.01
    @test value(nonnegk(kr)) == 0.77
    @test σ(nonnegk(kr)) == 0.01
    kr = KRatio(
        [n"Fe K-L2", n"Fe K-L3"],
        Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
        Dict(:BeamEnergy => 15.0e3, :TakeOffAngle => deg2rad(38.0)),
        mat"Fe2O3",
        uv(-0.01, 0.01),
    )
    @test value(nonnegk(kr)) == 0.0
    @test σ(nonnegk(kr)) == 0.01
    @test element(kr) == n"Fe"
    kr2 = KRatio(
        [n"O K-L3"],
        Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
        Dict(:BeamEnergy => 15.0e3, :TakeOffAngle => deg2rad(38.0)),
        mat"Fe2O3",
        uv(0.6, 0.01),
    )
    krs = [kr, kr2]
    @test length(elms(krs)) == 2
    @test n"Fe" in elms(krs)
    @test n"O" in elms(krs)
    @test !(n"C" in elms(krs))

    df = asa(DataFrame, krs)
    @test df[1, "Standard"]=="Fe2O3"
    @test isapprox(df[1, "C[std]"],0.699432,atol=1.0e-6)
    @test isapprox(df[2, "C[std]"],1.0-0.699432,atol=1.0e-6)

    @test_throws ErrorException KRatio(
        [n"O K-L3", n"F K-L3"],
        Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
        Dict(:BeamEnergy => 15.0e3, :TakeOffAngle => deg2rad(38.0)),
        mat"Fe2O3",
        uv(0.6, 0.01),
    )
    @test_throws ErrorException KRatio(
        [n"F K-L3"],
        Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
        Dict(:BeamEnergy => 15.0e3, :TakeOffAngle => deg2rad(38.0)),
        mat"Fe2O3",
        uv(0.6, 0.01),
    )
end


text:
using OpenQuantumBase, Test

x = range(0,stop=10,length=100)
y1 = Array(x) + 1.0im*Array(x)
y2 = (10.0+10.0im) .- Array(x)
inter_complex = construct_interpolations(x, y1)
inter_real = construct_interpolations(x, 10 .- collect(x))
gridded_2_range_inter = construct_interpolations(collect(x), y1)

test_x = range(0,stop=10,length=50)
exp_complex = collect(test_x) + 1.0im*collect(test_x)
res_complex = inter_complex.(test_x)
exp_real = (10.0) .- collect(test_x)
res_real = inter_real.(test_x)

@test res_complex ≈ exp_complex
@test res_real ≈ exp_real
@test isa(inter_complex(-1), Complex)
@test inter_complex(-1) ≈ 0
@test inter_real(11) ≈ 0
@test gridded_2_range_inter.(test_x) ≈ exp_complex
@test gradient(inter_real, 2.3) ≈ -1
@test gradient(inter_real, [0.13, 0.21]) ≈ [-1, -1]

# Test for multi-demension array
y_array = transpose(hcat(y1, y2))
y_array_itp = construct_interpolations(x, y_array, order=1)
@test exp_complex ≈ y_array_itp(1, test_x)
@test y_array_itp([1, 2], -1) ≈ [0, 0]

# Test for extrapolation
x = range(1.0,stop=10.0)
y = 10 .- collect(x)
y_c = x + 1.0im*y

eitp = construct_interpolations(x, Array(x); extrapolation="line")
@test isapprox(eitp(0.0),0.0, atol=1e-8)
eitp = construct_interpolations(x, y_c; extrapolation="line")
@test imag(eitp(0)) ≈ 10
eitp = construct_interpolations(x, y_c; extrapolation="flat")
@test eitp(0) == 1.0 + 9.0im

x = [1.0, 3, 4, 5, 6, 7, 8, 9, 10]
@test_logs (:warn,"The grid is not uniform. Using grided linear interpolation.") construct_interpolations(x, 10 .- x, method="bspline")


text:
function hist_param(param,p1,p2,nhist)
np = length(param)
param_grid = zeros(2*nhist)
for i=1:nhist
  param_grid[2*i-1] = p1 + (p2-p1) * (i-1.0)/nhist
  param_grid[2*i  ] = p1 + (p2-p1) * i/nhist
end
param_hist = zeros(2*nhist)
nmiss = 0
for i=1:np
  ibin = floor(Int64,(param[i]-p1)/(p2-p1)*nhist)+1
  if ibin >= 1 && ibin <= nhist
    param_hist[2*ibin-1] += 1
    param_hist[2*ibin  ] += 1
  else
    nmiss +=1
  end
end
return [p1;param_grid;p2],[0;param_hist;0],nmiss
end

function rho_pearson(p1,p2)
@assert(length(p1)==length(p2))
np = length(p1)
rho = 0.0
m1 = mean(p1)
m2 = mean(p2)
s1 = std(p1)
s2 = std(p2)
for i=1:np
  rho += (p1[i]-m1)*(p2[i]-m2)
end
return rho/s1/s2/(np-1)
end

function acf_par(pvec,nstep_corr)
# Computes ACF:
nstep = length(pvec)
mpar = mean(pvec)
spar = std(pvec)
# Now compute ACF:
acf = zeros(nstep_corr)
for k=1:nstep_corr
  for m= 1:nstep-k
    acf[k] += (pvec[m]-mpar)*(pvec[m+k]-mpar)
  end
  acf[k] = acf[k]/spar^2/(nstep-k)
end
return acf
end

# Make a binned light curve:
include("model_tyc3559_v2.jl")
include("transit_orb.jl")
include("occultquad.jl")
include("kepler.jl")

using JLD
using PyPlot

#@load "tyc3559_mcmc_norv_fixedkernel_v03.jld"
@load "tyc3559_mcmc_norv_fixedkernel_v03_thin.jld"

#par_mcmc[:,:,4]=par_mcmc[:,:,4].*(1.0-par_mcmc[:,:,15])
#par_mcmc[:,:,13]=par_mcmc[:,:,13].*(1.0-par_mcmc[:,:,15])
#pname[4] = "D*(1-f3)"
#pname[13] = "a2c*(1-f3)"
# Derived parameters:

# 1). Density of star B:

using CGS
period = (par_mcmc[:,:,2]-par_mcmc[:,:,1])/nlast
tdur = par_mcmc[:,:,5]
b = par_mcmc[:,:,3]
rhoB = 3.0/pi^2 .*period./tdur.^3/(24.*3600.)^2/GRAV.*sqrt(1.0-b.^2)
mrhoB = mean(rhoB); srhoB = std(rhoB)
println("Density of star B: ",mrhoB," +- ",srhoB)

# make histogram:
nhist = 80
p1 = mrhoB-3*srhoB; p2 = mrhoB+3*srhoB
param_grid,param_hist,nmiss = hist_param(rhoB,p1,p2,nhist)

clf()
plot(param_grid,param_hist/maximum(param_hist))
plot([mrhoB-srhoB,mrhoB-srhoB],[0,1])
plot([mrhoB+srhoB,mrhoB+srhoB],[0,1])
plot([mrhoB,mrhoB],[0,1])
xlabel("Density of star B [g/cc]")
ylabel("Relative probability")

read(STDIN,Char)

# 2). Doppler mass of brown dwarf:

teff = 6500.0
x = PLANCK * C / (6000e-8)/BOLTZMANN/teff
# Approximate beaming with Blackbody [NEED TO REDO WITH ATMOSPHERE & KEPLER BANDPASS]:
#alpha_beam = (exp(x)*(3.0-x)-3.0)/(exp(x)-1.0)
alpha_beam = -0.898
# Compute radial velocity in km/s:
vr = par_mcmc[:,:,12]/(3.0-alpha_beam)*C/1e5
# Convert to mass of brown dwarf, assuming mass of star:
mB = 1.76
semi = (GRAV * mB * MSUN)^(1//3).* (period.*24*3600/2/pi).^(2//3)
aonr = period./tdur.*sqrt(1.0-b.^2)/pi
inc = acos(b./aonr)
vorb = 2*pi.*semi./(period.*24.*3600.)./1e5
# Mass of brown dwarf:
mBb1 = vr./vorb.*mB.*MSUN./MJUPITER./sin(inc)
mmBb1=mean(mBb1)
smBb1=std(mBb1)
println("Mass of Bb from Doppler: ",mmBb1," +- ",smBb1," M_Jupiter")

# 3a). Tidal distortion mass of brown dwarf:
# Equation 2 from Faigler & Mazeh:
q1 = par_mcmc[:,:,6]; q2 = par_mcmc[:,:,7]
u1 = 2.*sqrt(q1).*q2
u2 = sqrt(q1).*(1-2.*q2)
#g = 0.3+rand(size(q1))*0.7
g = 0.65
alpha_ellip = 0.15.*(15.0+u1).*(1.0+g)./(3.0-u1)
mBb2 = mB * MSUN .* aonr.^3./alpha_ellip.* abs(par_mcmc[:,:,13])./MJUPITER./sin(inc).^2
mmBb2 = mean(mBb2)
smBb2 = std(mBb2)
println("Mass of Bb from tides: ",mmBb2," +- ",smBb2," M_Jupiter")

# 3b). Compute the mass from the RV:
krv0 =  5.985117177227862
sig_krv = 0.43164531668541845
krv = krv0 + sig_krv*randn(nwalkers,nsteps)

mBb3 = (period.*24.*3600.*(krv.*1e5).^3.*(mB*MSUN)^2/(2pi*GRAV)).^(1//3)./sin(inc)./MJUPITER
mmBb3 = mean(mBb3)
smBb3 = std(mBb3)
println("Mass of Bb from RV: ",mmBb3," +- ",smBb3," M_Jupiter")

p1 = minimum([mmBb2-3smBb2,mmBb1-3smBb1,mmBb3-3smBb3])
p2 = maximum([mmBb2+3smBb2,mmBb1+3smBb1,mmBb3+3smBb3])
param_grid,param_hist1,nmiss = hist_param(mBb1,p1,p2,nhist)
param_grid,param_hist2,nmiss = hist_param(mBb2,p1,p2,nhist)
param_grid,param_hist3,nmiss = hist_param(mBb3,p1,p2,nhist)
clf()
#plot(param_grid,param_hist1)
#plot(param_grid,param_hist2)
#plot(param_grid,param_hist3)
#xlabel("Mass of brown dwarf [M_Jupiter]")
#ylabel("Relative probability")

fig, ax = subplots()
ax[:plot](param_grid,param_hist1./maximum(param_hist3),"r-",alpha=0.6,label="Beaming",linewidth=3)
ax[:plot](param_grid,param_hist2./maximum(param_hist3),"b-",alpha=0.6,label="Ellipsoidal",linewidth=3)
ax[:plot](param_grid,param_hist3./maximum(param_hist3),"g-",alpha=0.6,label="Doppler RV",linewidth=3)
ax[:set_xlabel]("Mass of brown dwarf [M_Jupiter]")
ax[:set_ylabel]("Relative probability")
ax[:axis]([0,200,0.,1.1])
ax[:legend](loc ="upper right")
PyPlot.savefig("brown_dwarf_mass_comparison.png",bbox_inches="tight")

read(STDIN,Char)

#
# 4). Brightness of brown dwarf from phase-curve & secondary eclipse:
fBb_phase = 2.*abs(par_mcmc[:,:,11])*1e6
fBb_occult = abs(par_mcmc[:,:,16])*1e6
mfBb_phase = mean(fBb_phase)
sfBb_phase = std(fBb_phase)
mfBb_occult = mean(fBb_occult)
sfBb_occult = std(fBb_occult)
println("Flux ratio of Bb to B+Bb, phase:  ",mfBb_phase," +- ",sfBb_phase)
println("Flux ratio of Bb to B+Bb, occult: ",mfBb_occult," +- ",sfBb_occult)
p1 = minimum([mfBb_phase-3*sfBb_phase,mfBb_occult-3*sfBb_occult])
p2 = maximum([mfBb_phase+3*sfBb_phase,mfBb_occult+3*sfBb_occult])
param_grid,param_hist1,nmiss = hist_param(fBb_phase,p1,p2,nhist)
param_grid,param_hist2,nmiss = hist_param(fBb_occult,p1,p2,nhist)
clf()
plot(param_grid,param_hist1)
plot(param_grid,param_hist2)
xlabel("fBb/(fB+fBb) [ppm]")
ylabel("Relative probability")
read(STDIN,Char)

# 5). Radius ratio of brown dwarf to star B:
depth = par_mcmc[:,:,4]
x = 1.0-sqrt(1.0-b.^2)
r21 = sqrt(depth.*(1.-u1./3-u2./6)./(1.0-u1.*x-u2.*x.^2))
mr21 = mean(r21); sr21=std(r21)
println("R_Bb/R_B: ",mr21," +- ",sr21)

# 6). Geometric albedo of the brown dwarf from occultation:
ag = fBb_occult.*1e-6.*(aonr./r21).^2
println("Geometric albedo: ",mean(ag)," +- ",std(ag))

# 7). Radius of brown dwarf:
rBb = r21.*(mB*MSUN./rhoB.*3/4/pi).^(1//3)./RJUPITER
mrBb = mean(rBb); srBb = std(rBb)
println("R_Bb: ",mrBb," +- ",srBb," [R_Jupiter]")
read(STDIN,Char)

# 8). Radius of Star B:

rB = (mB*MSUN./rhoB.*3/4/pi).^(1//3)./RSUN
mrB = mean(rB); srB = std(rB)
println("R_Bb: ",mrB," +- ",srB," [R_Jupiter]")

# 9).  Compute the density from RV:
rho_Bb = mBb3.*MJUPITER./(rBb.*RJUPITER).^3.*3./(4pi)
mrho_Bb = mean(rho_Bb); srho_Bb = std(rho_Bb)
println("Density of brown dwarf [g/cc]: ",mrho_Bb," +- ",srho_Bb," [g/cc]")

# # Compute autocorrelation functions:
# nstep_corr = 2000
# acorrlength = zeros(nvary,nwalkers)
# for i=1:nvary
#   clf()
#   acf = zeros(nwalkers,nstep_corr)
#   for j=1:nwalkers
# # First estimate mean & variance:
#     pvec = vec(par_mcmc[j,:,ivary[i]])
#     acf[j,:]=acf_par(pvec,nstep_corr)
#     acorrlength[i,j] = 0.0
#     k=1
#     while (acf[j,k] > 0.5) && (k < nstep_corr-1)
#       k+=1
#     end
#     acorrlength[i,j] = k
#     plot(vec(acf[j,:]))
#   end
#   println("Autocorrelation length for parameter ",pname[ivary[i]]," is ",maximum(acorrlength[i,:]))
# #  read(STDIN,Char)
# end
#
# println("Autocorrelation length for chains is: ",maximum(acorrlength))
# println("Effective length of MCMC: ",nwalkers*nsteps/maximum(acorrlength))
#
# read(STDIN,Char)
clf()

# Compute inter-parameter correlation
for i=1:nvary-1
  for j=i+1:nvary
    p1=vec(par_mcmc[:,:,ivary[i]])
    p2=vec(par_mcmc[:,:,ivary[j]])
    rho = rho_pearson(p1,p2)
    if abs(rho) > 0.5
      println("Parameters ",pname[ivary[i]]," and ",pname[ivary[j]]," correlate strongly  as ",rho)
      scatter(p1,p2)
      read(STDIN,Char)
      clf()
    end
  end
end

clf()

for i=1:nvary
  mp = zeros(nsteps)
  sp = zeros(nsteps)
  for j=1:nwalkers
    plot(vec(par_mcmc[j,1:nsteps,ivary[i]]))
  end
  for j=1:nsteps
    mp[j]=mean(par_mcmc[:,j,ivary[i]])
    sp[j]=std(par_mcmc[:,j,ivary[i]])
  end
  plot(mp+sp)
  plot(mp-sp)
  xlabel("MCMC step")
  ylabel(pname[ivary[i]])
  println("Hit return to continue")
  read(STDIN,Char)
  clf()
end

log_like_best =  -1e100
ibest = 0
jbest = 0
param_best = zeros(nparam)
for i=1:nsteps
  for j=1:nwalkers
    if log_like_mcmc[j,i] > log_like_best
      ibest = i
      jbest = j
      log_like_best = log_like_mcmc[j,i]
      param_best = vec(par_mcmc[j,i,:])
    end
  end
end

model = model_tyc3559_v2(time_phot,param_best)

# Compute residuals:

y = flux - model

period = (param_best[2]-param_best[1])/nlast

t0 = param_best[1]

nbin = 128
phasebin_resid = zeros(nbin)
phasebin_resid_med = zeros(nbin)
phasebin_raw = zeros(nbin)
phasebin_raw_med = zeros(nbin)
phasebin_model = zeros(nbin)
phase0 = mod((time_phot-t0)+period/4,period)/period
is = sortperm(phase0)
j=1
for ibin=1:nbin
  i=j
  while phase0[is[j+1]] < (ibin/nbin) && (j+1) < ntime
    j=j+1
  end
  phasebin_resid[ibin] = mean(y[is[i:j]])
  phasebin_resid_med[ibin] = median(y[is[i:j]])
  phasebin_raw[ibin] = mean(flux[is[i:j]])
  phasebin_raw_med[ibin] = median(flux[is[i:j]])
  phasebin_model[ibin] = mean(model[is[i:j]])
end

phase = zeros(nbin)
for i=1:nbin
  phase[i] = (i-0.5)/nbin
end

clf()
#scatter(phase,phasebin_resid)
scatter(phase,phasebin_raw)
scatter(phase,phasebin_raw_med)
plot(phase,phasebin_model)
plot(phase,phasebin_model)
model_reflect=((1+param_best[11]*cos((phase-0.25)*2pi))*(1.0-param_best[15])+param_best[15])*param_best[8]
model_beam   =((1+param_best[12]*sin((phase-0.25)*2pi))*(1.0-param_best[15])+param_best[15])*param_best[8]
model_ellip  =((1+param_best[13]*cos((phase-0.25)*4pi))*(1.0-param_best[15])+param_best[15])*param_best[8]
#  model[i]=((model_phot1*(1.0-param[16])+model_phot2*param[16])* (1.0-param[15])+param[15])*f0

plot(phase,model_reflect)
plot(phase,model_beam)
plot(phase,model_ellip)


text:
@doc """ check whether there are two gene whose location intersect
""" ->
function check_geneloc()
    chr_sted_gene = load(chr_sted_gene_fl, "chr_sted_gene_dict")
end

@doc """ all the chrsteds can't intersect with each other
""" ->
function check_chrsteds()

end


text:
# Activation functions
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

export ReLU, ReLUgrad
export LeakyReLU, LeakyReLUinv, LeakyReLUgrad
export Sigmoid, SigmoidInv, SigmoidGrad
export GaLU, GaLUgrad

###############################################################################
# Rectified linear unit (ReLU) (not invertible)

"""
    y = ReLU(x)

 Rectified linear unit (not invertible).

 See also: [`ReLUgrad`](@ref)
"""
function ReLU(x)
    return max.(0f0, x)
end

"""
    Δx = ReLUgrad(Δy, x)

 Backpropagate data residual through ReLU function.

 *Input*:

 - `Δy`: data residual

 - `x`: original input (since not invertible)

 *Output*:

 - `Δx`: backpropagated residual

 See also: [`ReLU`](@ref)
"""
function ReLUgrad(Δy, x)
    return Δy.*(sign.(x) .+ 1)/2
end

###############################################################################
# Leaky ReLU (invertible)

"""
    y = LeakyReLU(x; slope=0.01f0)

 Leaky rectified linear unit.

 See also: [`LeakyReLUinv`](@ref), [`LeakyReLUgrad`](@ref)
"""
function LeakyReLU(x; slope=0.01f0)
    return max.(0f0, x) + slope*min.(0f0, x)
end

"""
    x = LeakyReLUinv(y; slope=0.01f0)

 Inverse of leaky ReLU.

 See also: [`LeakyReLU`](@ref), [`LeakyReLUgrad`](@ref)
"""
function LeakyReLUinv(y; slope=0.01f0)
    return max.(0f0, y) + (1f0/slope)*min.(0f0, y)
end

"""
    Δx = ReLUgrad(Δy, y; slope=0.01f0)

 Backpropagate data residual through leaky ReLU function.

 *Input*:

 - `Δy`: residual

 - `y`: original output

 - `slope`: slope of non-active part of ReLU

 *Output*:

 - `Δx`: backpropagated residual

 See also: [`LeakyReLU`](@ref), [`LeakyReLUinv`](@ref)
"""
function LeakyReLUgrad(Δy, y; slope=0.01f0)
    x = LeakyReLUinv(y; slope=slope)  # recompute forward state
    p_mask = (sign.(x) .+ 1f0)/2f0
    return Δy.*p_mask + slope*Δy.*(1f0 .- p_mask)
end

###############################################################################
# Sigmoid (invertible if ouput nonzero)

"""
    y = Sigmoid(x)

 Sigmoid activation function.

 See also: [`SigmoidInv`](@ref), [`SigmoidGrad`](@ref)
"""
function Sigmoid(x)
    y = 1f0 ./ (1f0 .+ exp.(-x))
    return y
end

"""
    x = SigmoidInv(y)

 Inverse of Sigmoid function.

 See also: [`Sigmoid`](@ref), [`SigmoidGrad`](@ref)
"""
function SigmoidInv(y)
    if sum(isapprox.(y, 0f-6)) == 0
        x = -log.((1f0 .- y) ./ y)
    else
        throw("Input contains zeros.")
    end
    return x
end

"""
    Δx = SigmoidGrad(Δy, y; x=nothing)

 Backpropagate data residual through Sigmoid function.

 *Input*:

 - `Δy`: residual

 - `y`: original output

 - `x`: original input, if y not available (in this case, set y=nothing)

 *Output*:

 - `Δx`: backpropagated residual

 See also: [`Sigmoid`](@ref), [`SigmoidInv`](@ref)
"""
function SigmoidGrad(Δy, y; x=nothing)
    if ~isnothing(y) && isnothing(x)
        x = SigmoidInv(y)  # recompute forward state
    end
    Δx = Δy .* exp.(-x) ./ (1f0 .+ exp.(-x)).^2f0
    return Δx
end

###############################################################################
# Gated linear unit (GaLU) (not invertible)
# Adapted from Dauphin et al. (2017)

"""
    y = GaLU(x)

 Gated linear activation unit (not invertible).

 See also: [`GaLUgrad`](@ref)
"""
function GaLU(x::AbstractArray{Float32, 4})
    k = Int(size(x, 3) / 2)
    y = x[:, :, 1:k, :] .* Sigmoid(x[:, :, k+1:end, :])
    return y
end

function GaLU(x::AbstractArray{Float32, 5})
    k = Int(size(x, 4) / 2)
    y = x[:, :, :, 1:k, :] .* Sigmoid(x[:, :, :, k+1:end, :])
    return y
end

"""
    Δx = GaLUgrad(Δy, x)

 Backpropagate data residual through GaLU activation.

 *Input*:

 - `Δy`: residual

 - `x`: original input (since not invertible)

 *Output*:

 - `Δx`: backpropagated residual

 See also: [`GaLU`](@ref)
"""
function GaLUgrad(Δy::AbstractArray{Float32, 4}, x::AbstractArray{Float32, 4})
    k = Int(size(x, 3) / 2)
    x1 = x[:, :, 1:k, :]
    x2 = x[:, :, k+1:end, :]
    Δx = 0f0.*x
    Δx[:, :, 1:k, :] = Sigmoid(x2) .* Δy
    Δx[:, :, k+1:end, :] = SigmoidGrad(Δy, nothing; x=x2) .* x1
    return Δx
end

function GaLUgrad(Δy::AbstractArray{Float32, 5}, x::AbstractArray{Float32, 5})
    k = Int(size(x, 4) / 2)
    x1 = x[:, :, :, 1:k, :]
    x2 = x[:, :, :, k+1:end, :]
    Δx = 0f0.*x
    Δx[:, :, :, 1:k, :] = Sigmoid(x2) .* Δy
    Δx[:, :, :, k+1:end, :] = SigmoidGrad(Δy, nothing; x=x2) .* x1
    return Δx
end


text:
# Thi file is a part of the ExamplePackage.

"""
    sort_numbers(v)

Sort numbers so that numbers smaller than 5 are sorted into veca vector and
other numbers are sorted into vecb vector.
"""

function sort_numbers(v)
    veca = []
    vecb = []
    for i in v
        if i < 5
          push!(veca, i)
        else
          push!(vecb, i)
        end
    end
    return veca, vecb
end


text:
# These functions are used to hold "expensive" values in case they are
# needed again. Note that if the graph is modified, the cache should be
# wiped clean.

export cache_clear, cache_on, cache_off, cache_recall, cache_check, cache_save
# but we don't export cache_save

# We try to match the symbol's name with the function it saves.

"""
`cache_clear(G)` clears all items in `G`'s cache.

`cache_clear(G,item)` clears just that item.
"""
function cache_clear(G::SimpleGraph)
    if length(G.cache) > 0
        G.cache = Dict{Symbol,Any}()
    end
    nothing
end

cache_clear(G::SimpleGraph, item::Symbol) = delete!(G.cache, item)

"""
`cache_check(G,item)` checks if the symbol `item` is a valid key.
"""
cache_check(G::SimpleGraph, item::Symbol)::Bool = G.cache_flag && haskey(G.cache, item)

"""
`cache_recall(G,item)` retreives the value associated with `item`.

**WARNING**: No check is done to see if this value is defined. Be
sure to use `cache_check` first!
"""
cache_recall(G::SimpleGraph, item::Symbol) = deepcopy(G.cache[item])

"""
`cache_save(G,item,value)` saves `value` associated with
the symbol (key) `item` in the cache for this graph.
"""
function cache_save(G::SimpleGraph, item::Symbol, value)
    if G.cache_flag
        G.cache[item] = deepcopy(value)
    end
    nothing
end

"""
`cache_on(G)` activates results caching for this graph. See also: `cache_off`.
"""
cache_on(G::SimpleGraph) = G.cache_flag = true

"""
`cache_off(G)` deactivates cache checking. You may also wish to call
`cache_clear` to recover storage.

See also: `cache_on`.
"""
function cache_off(G::SimpleGraph)
    G.cache_flag = false
end


text:
module TilesetSplitter

using ..Ahorn
using Cairo

const tilesetCache = WeakKeyDict{Ahorn.Sprite, Array{Ahorn.Sprite, 2}}()

function generateTilesetSpriteMatrix(sprite::Ahorn.Sprite, tileWidth=8, tileHeight=8)
    if sprite.width % tileWidth != 0.0 sprite.height % tileHeight != 0.0
        error("Width and height must be divisible by tileWidth and tileHeight respectively.")
    end

    chunksX, chunksY = floor(Int, sprite.width / tileWidth), floor(Int, sprite.height / tileHeight)

    res = Array{Ahorn.Sprite, 2}(undef, chunksY, chunksX)

    for x in 0:1:chunksX - 1, y in 0:1:chunksY - 1
        surface = CairoARGBSurface(tileWidth, tileHeight)
        ctx = Ahorn.getSurfaceContext(surface)

        Ahorn.drawImage(ctx, sprite.surface, 0, 0, sprite.x + x * tileWidth, sprite.y + y * tileHeight, tileWidth, tileHeight, alpha=1.0)

        res[y + 1, x + 1] = Ahorn.Sprite(
            0,
            0,

            tileWidth,
            tileHeight,

            sprite.offsetX,
            sprite.offsetY,

            tileWidth,
            tileHeight,

            surface,
            sprite.filename
        )
    end

    return res
end

function getTilesetSpriteMatrix(sprite::Ahorn.Sprite, tileWidth=8, tileHeight=8)
    get!(tilesetCache, sprite) do
        return generateTilesetSpriteMatrix(sprite, tileWidth, tileHeight)
    end
end

end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn