text:
# test_response_format.jl
#
# Author: Tom Price
# Date: Dec 2018
#
# This Julia script tests the JSON data structures
# that are returned in the body of responses to GET requests
#
# *     based on documentation of REST API from juliaapi_new.txt
#
# *     Julia listens on http://localhost:8081 via HttpServer.jl and dispatch.jl
#       and sends responses to the Rails app on http://localhost:3000
#
# *     the relevant controller in the Rails app is
#       chaipcr/web/app/controllers/experiments_controller.rb
#
# *     currently calls are POST-ed to
#       http://127.0.0.1:8081/experiments/#{experiment.id}/standard_curve
#       http://127.0.0.1:8081/experiments/#{experiment.id}/amplification
#       http://127.0.0.1:8081/experiments/#{experiment.id}/meltcurve
#       http://127.0.0.1:8081/experiments/#{experiment.id}/analyze

import JSON, DataStructures.OrderedDict
import FactCheck: facts, context, @fact, clear_results, exitstatus, less_than_or_equal

FactCheck.setstyle(:default)

# ********************************************************************************
#
# call: experiments/:experiment_id/standard_curve
#
# ********************************************************************************

function standard_curve_response_test()
    response=JSON.parse("""{
        "targets":
        [
            {"target_id": 1, "slope": 9.99, "offset": 9.99, "efficiency": 1.02, "r2": 0.99},
            {"target_id": 2, "slope": 9.99, "offset": 9.99, "efficiency": 0.98, "r2": 0.99},
            {"target_id": 3, "error": "xxxxx"},
            {"target_id": 4, "error": "xxxxx"}
        ]
    }"""; dicttype=OrderedDict)
    verify_response(
        Action_DICT["standard_curve"](),
        response
    )
end

# ********************************************************************************
#
# call: experiments/:experiment_id/amplification
#
# ********************************************************************************

# single channel

function singlechannel_amplification_response_test()
    response=JSON.parse(
        # after deconvolution (if dual channel) and adjusting well-to-well variation,
        # before baseline subtraction
    """{
        "rbbs_ary3": [
            [
                [1.0101, 1.0102,    1.0140],
                [1.0201, 1.0202,    1.0240],

                [1.1601, 1.1602,    1.1640]
            ]
        ],
    """*
        # after baseline subtraction
    """
        "blsub_fluos": [
            [
                [1.0101, 1.0102,    1.0140],
                [1.0201, 1.0202,    1.0240],

                [1.1601, 1.1602,    1.1640]
            ]
        ],
    """*
        # first derivative of the amp curve
    """
        "dr1_pred": [
            [
                [1.0101, 1.0102,    1.0140],
                [1.0201, 1.0202,    1.0240],

                [1.1601, 1.1602,    1.1640]
            ]
        ],
    """*
        # second derivative of the amp curve
    """
        "dr2_pred": [
            [
                [1.0101, 1.0102,    1.0140],
                [1.0201, 1.0202,    1.0240],

                [1.1601, 1.1602,    1.1640]
            ]
        ],

        "cq": [
            [1.01, 1.02,    1.16]
        ],
    """*
        # starting quantity from absolute quantification
    """
        "d0": [
            [1.01, 1.02,    1.16]
        ],
    """*
        # fluorescence threshold if Ct method is used:
        # [channel_1] for single channel,
        # [channel_1, channel_2] for dual channel,
        # [] empty for automatic detection from data
    """
        "ct_fluos": [1.0],
        "assignments_adj_labels_dict": {
            "rbbs_ary3":   ["aa", "aa",     "Aa"],
            "blsub_fluos": ["aA", "aa",     "aa"],
            "d0":          ["aa", "AA",     "Aa"],
            "cq":          ["aa", "Aa",     "aa"]
        },
        "valid": true
    }"""; dicttype=OrderedDict)
    verify_response(
        Action_DICT["amplification"](),
        response
    )
end

# dual channel

function dualchannel_amplification_response_test()
    response=JSON.parse(
        # after deconvolution (if dual channel) and adjusting well-to-well variation,
        # before baseline subtraction
    """{
        "rbbs_ary3": [
            [
                [1.0101, 1.0102,    1.0140],
                [1.0201, 1.0202,    1.0240],

                [1.1601, 1.1602,    1.1640]
            ],
            [
                [2.0101, 2.0102,    2.0140],
                [2.0201, 2.0202,    2.0240],

                [2.1601, 2.1602,    2.1640]
            ]
        ],
    """*
        # after baseline subtraction
    """
        "blsub_fluos": [
            [
                [1.0101, 1.0102,    1.0140],
                [1.0201, 1.0202,    1.0240],

                [1.1601, 1.1602,    1.1640]
            ],
            [
                [2.0101, 2.0102,    2.0140],
                [2.0201, 2.0202,    2.0240],

                [2.1601, 2.1602,    2.1640]
            ]
        ],
    """*
        # first derivative of the amp curve
    """
        "dr1_pred": [
            [
                [1.0101, 1.0102,    1.0140],
                [1.0201, 1.0202,    1.0240],

                [1.1601, 1.1602,    1.1640]
            ],
            [
                [2.0101, 2.0102,    2.0140],
                [2.0201, 2.0202,    2.0240],

                [2.1601, 2.1602,    2.1640]
            ]
        ],
    """*
        # second derivative of the amp curve
    """
        "dr2_pred": [
            [
                [1.0101, 1.0102,    1.0140],
                [1.0201, 1.0202,    1.0240],

                [1.1601, 1.1602,    1.1640]
            ],
            [
                [2.0101, 2.0102,    2.0140],
                [2.0201, 2.0202,    2.0240],

                [2.1601, 2.1602,    2.1640]
            ]
        ],

        "cq": [
            [1.01, 1.02,    1.16],
            [2.01, 2.02,    2.16]
        ],
    """*
        # starting quantity from absolute quantification
    """
        "d0": [
            [1.01, 1.02,    1.16],
            [2.01, 2.02,    2.16]
        ],
    """*
        # fluorescence threshold if Ct method is used:
        # [channel_1] for single channel,
        # [channel_1, channel_2] for dual channel,
        # [] empty for automatic detection from data
    """
        "ct_fluos": [1.0, 2.0],
        "assignments_adj_labels_dict": {
            "rbbs_ary3":   ["aa", "aa",     "Aa"],
            "blsub_fluos": ["aA", "aa",     "aa"],
            "d0":          ["aa", "AA",     "Aa"],
            "cq":          ["aa", "Aa",     "aa"]
        },
        "valid": true
    }"""; dicttype=OrderedDict)
    verify_response(
        Action_DICT["amplification"](),
        response
    )
end

# error response body:

function error_amplification_response_test()
    response=JSON.parse("""{
        "valid": false,
        "error": "xxxx"
    }"""; dicttype=OrderedDict)
    verify_response(
        Action_DICT["amplification"](),
        response
    )
end

# ********************************************************************************
#
# call: experiments/:experiment_id/meltcurve
#
# ********************************************************************************

# each matrix r is [ temperature_1 fluo_1 derivative_1;
#                    temperature_2 fluo_2 derivative_2;
#                    ...           ...    ...
#                    temperature_n fluo_n derivative_n]
# n_grid >> 0 is the same for each well
#
# in JSON
r="[[0.101, 0.102,     0.199]," * # temp
  " [0.201, 0.202,     0.299]," * # fluo
  " [0.301, 0.302,     0.399]]"   # slope
r1_01=r; r1_02=r; r1_16=r
r2_01=r; r2_02=r; r2_16=r

# each matrix s is [ temp_max_1 peak_area_1;
#                    temp_max_2 peak_area_2;
#                    ...           ...    ...
#                    temp_max_n peak_area_n;]
# n_peaks >=0 may be different for each well
#
# in JSON
s="[[0.101, 0.102,     0.104]," * # Tm
  " [0.201, 0.202,     0.204]]"   # area
s1_01=s; s1_02=s; s1_16=s
s2_01=s; s2_02=s; s2_16=s

# single channel data

function singlechannel_meltcurve_response_test()
    response=JSON.parse("""{
        "melt_curve_data":     [
            [ $(r1_01), $(r1_02),     $(r1_16) ]
        ],
        "melt_curve_analysis": [
            [ $(s1_01), $(s1_02),     $(s1_16) ]
        ],
        "valid": true
    }"""; dicttype=OrderedDict)
    println(response)
    println(JSON.json(response))
    verify_response(
        Action_DICT["meltcurve"](),
        response
    )
end

# dual channel data

function dualchannel_meltcurve_response_test()
    response=JSON.parse("""{
        "melt_curve_data": [
            [ $(r1_01), $(r1_02),     $(r1_16) ],
            [ $(r2_01), $(r2_02),     $(r2_16) ]
        ],
        "melt_curve_analysis": [
            [ $(s1_01), $(s1_02),     $(s1_16) ],
            [ $(s2_01), $(s2_02),     $(s2_16) ]
        ],
        "valid": true
    }"""; dicttype=OrderedDict)
    verify_response(
        Action_DICT["meltcurve"](),
        response
    )
end

# error response body:

function error_meltcurve_response_test()
    response=JSON.parse("""{
        "valid": false,
        "error": "xxxx"
    }"""; dicttype=OrderedDict)
    verify_response(
        Action_DICT["meltcurve"](),
        response
    )
end

# ********************************************************************************
#
# call: system/loadscript?script=path%2Fto%2Fanalyze.jl
#
# ********************************************************************************

# success response body:

function loadscript_response_test()
    response=JSON.parse("""{
        "script": "path/to/analyze.jl",
        "valid": true
    }"""; dicttype=OrderedDict)
    verify_response(
        Action_DICT["loadscript"](),
        response
    )
end

# error response body:

function error_loadscript_response_test()
    response=JSON.parse("""{
        "valid": false,
        "error": "xxxx"
    }"""; dicttype=OrderedDict)
    verify_response(
        Action_DICT["loadscript"](),
        response
    )
end

# ********************************************************************************
#
# call: experiments/:experiment_id/thermal_performance_diagnostic
#
# ********************************************************************************

# success response body (thermal_performance_diagnostic):

function thermal_performance_diagnostic_response_test()
    response=JSON.parse("""{
        "Heating":  {
            "AvgRampRate": [5.3743,true],
            "TotalTime": [8001,false],
            "MaxBlockDeltaT": [1.27,true]
        },
        "Cooling": {
            "AvgRampRate": [0,false],
            "TotalTime": [9000,false],
            "MaxBlockDeltaT": [0.94,true]
        },
        "Lid": {
            "HeatingRate": [1.3031,true],
            "TotalTime": [32999,false]
        },
        "valid": true
    }"""; dicttype=OrderedDict)
    verify_response(
        Action_DICT["thermal_performance_diagnostic"](),
        response
    )
end

# ********************************************************************************
#
# call: experiments/:experiment_id/thermal_consistency
#
# ********************************************************************************

# success response body (thermal_consistency):

function thermal_consistency_response_test()
    response=JSON.parse("""{
        "tm_check": [
            {
                "Tm": [79.0773,true],
                "area": 1361.5005
            },
            {
                "Tm": [78.5635,true],
                "area": 1763.7998
            },
            {
                "Tm": [78.4396,true],
                "area": 1434.3995
            },
            {
                "Tm": [78.2172,true],
                "area": 1672.9787
            },
            {
                "Tm": [77.8384,true],
                "area": 639.2474
            },
            {
                "Tm": [78.3816,true],
                "area": 1076.2143
            },
            {
                "Tm": [78.6829,true],
                "area": 895.6206
            },
            {
                "Tm": [78.6403,true],
                "area": 366.6082
            },
            {
                "Tm": [78.0791,true],
                "area": 255.8169
            },
            {
                "Tm": [77.5054,true],
                "area": 193.9588
            },
            {
                "Tm": [78.0182,true],
                "area": 1114.7039
            },
            {
                "Tm": [78.171,true],
                "area": 1324.3671
            },
            {
                "Tm": [78.1117,true],
                "area": 1219.4364
            },
            {
                "Tm": [77.8815,true],
                "area": 499.9462
            },
            {
                "Tm": [78.3167,true],
                "area": 1097.0781
            },
            {
                "Tm": [79.2487,true],
                "area": 1318.3026
            }
        ],
        "delta_Tm": [1.7434,true],
        "valid": true
    }"""; dicttype=OrderedDict)
    verify_response(
        Action_DICT["thermal_consistency"](),
        response
    )
end

# ********************************************************************************
#
# call: experiments/:experiment_id/optical_cal
#
# ********************************************************************************

# valid

function valid_optical_cal_response_test()
    response=JSON.parse("""{
        "valid": true
    }"""; dicttype=OrderedDict)
    verify_response(
        Action_DICT["optical_cal"](),
        response
    )
end

# invalid

function invalid_optical_cal_response_test()
    response=JSON.parse("""{
        "valid": false,
        "error": "xxxx"
    }"""; dicttype=OrderedDict)
    verify_response(
        Action_DICT["optical_cal"](),
        response
    )
end

# ********************************************************************************
#
# call: experiments/:experiment_id/optical_test_single_channel
#
# ********************************************************************************

function singlechannel_optical_response_test()
    response=JSON.parse("""{
        "optical_data": [{
            "baseline": 1704,
            "excitation": 84756,
            "valid": true
        },{
            "baseline": 1844,
            "excitation": 76751,
            "valid": true
        },{
            "baseline": 1729,
            "excitation": 75655,
            "valid": true
        },{
            "baseline": 1695,
            "excitation": 32242,
            "valid": true
        },{
            "baseline": 1683,
            "excitation": 71588,
            "valid": true
        },{
            "baseline": 1837,
            "excitation": 96184,
            "valid": true
        },{
            "baseline": 1696,
            "excitation": 90374,
            "valid": true
        },{
            "baseline": 1848,
            "excitation": 88519,
            "valid": true
        },{
            "baseline": 1675,
            "excitation": 105414,
            "valid": true
        },{
            "baseline": 1731,
            "excitation": 90122,
            "valid": true
        },{
            "baseline": 1718,
            "excitation": 94174,
            "valid": true
        },{
            "baseline": 1803,
            "excitation": 110185,
            "valid": true
        },{
            "baseline": 1797,
            "excitation": 92436,
            "valid": true
        },{
            "baseline": 1773,
            "excitation": 35866,
            "valid": true
        },{
            "baseline": 1789,
            "excitation": 101293,
            "valid": true
        },{
            "baseline": 1864,
            "excitation": 98680,
            "valid": true
        }],
        "valid": true
    }"""; dicttype=OrderedDict)
    verify_response(
        Action_DICT["optical_test_single_channel"](),
        response
    )
end

# ********************************************************************************
#
# call: experiments/:experiment_id/optical_test_dual_channel
#
# *******************************************************************************

function dualchannel_optical_response_test()
    response=JSON.parse("""{
        "optical_data": [
            {
                "baseline": [[7373,true],[1790,true]],
                "FAM": [[56997,false],[25204,false]],
                "HEX": [[10827,false],[24547,false]],
                "water": [[7373,true],[1790,true]]
            },
            {
                "baseline": [[16590,true],[1908,true]],
                "FAM": [[113931,false],[44197,false]],
                "HEX": [[21176,false],[43755,false]],
                "water": [[16590,true],[1908,true]]
            },
            {
                "baseline": [[10622,true],[1945,true]],
                 "FAM": [[97683,false],[40424,false]],
                 "HEX": [[14754,false],[39349,false]],
                "water": [[10622,true],[1945,true]]
            },
            {
                "baseline": [[3194,true],[1692,true]],
                "FAM": [[15820,false],[7494,false]],
                "HEX": [[4289,false],[7489,false]],
                "water": [[3194,true],[1692,true]]
            },
            {
                "baseline": [[7210,true],[1623,true]],
                "FAM": [[56260,false],[22039,false]],
                "HEX": [[10460,false],[22624,false]],
                "water": [[7210,true],[1623,true]]
            },
            {
                "baseline": [[9309,true],[1958,true]],
                "FAM": [[75858,false],[30542,false]],
                "HEX": [[14908,false],[31364,false]],
                "water": [[9309,true],[1958,true]]
            },
            {
                "baseline": [[9923,true],[1804,true]],
                "FAM": [[76166,false],[32064,false]],
                "HEX": [[19045,false],[32810,false]],
                "water": [[9923,true],[1804,true]]
            },
            {
                "baseline": [[17224,true],[2063,true]],
                "FAM": [[110970,false],[39879,false]],
                "HEX": [[33050,false],[41685,false]],
                "water": [[17224,true],[2063,true]]
            },
            {
                "baseline": [[6116,true],[2098,true]],
                "FAM": [[75856,false],[35634,false]],
                "HEX": [[9767,false],[33825,false]],
                "water": [[6116,true],[2098,true]]
            },
            {
                "baseline": [[11647,true],[1967,true]],
                "FAM": [[118257,false],[43666,false]],
                "HEX": [[19407,false],[42969,false]],
                "water": [[11647,true],[1967,true]]
            },
            {
                "baseline": [[10132,true],[2558,true]],
                "FAM": [[80505,false],[35494,false]],
                "HEX": [[14507,false],[34011,false]],
                "water": [[10132,true],[2558,true]]
            },
            {
                "baseline": [[7019,true],[2757,true]],
                "FAM": [[66396,false],[25384,false]],
                "HEX": [[10391,false],[24408,false]],
                "water": [[7019,true],[2757,true]]
            },
            {
                "baseline": [[12806,true],[1900,true]],
                "FAM": [[117630,false],[45047,false]],
                "HEX": [[21790,false],[42365,false]],
                "water": [[12806,true],[1900,true]]
            },
            {
                "baseline": [[9926,true],[1843,true]],
                "FAM": [[70441,false],[24076,false]],
                "HEX": [[17115,false],[22416,false]],
                "water": [[9926,true],[1843,true]]
            },
            {
                "baseline": [[6189,true],[1822,true]],
                "FAM": [[49142,false],[24322,false]],
                "HEX": [[11908,false],[22722,false]],
                "water": [[6189,true],[1822,true]]
            },
            {
                "baseline": [[10674,true],[2133,true]],
                "FAM": [[106136,false],[47080,false]],
                "HEX": [[22462,false],[42102,false]],
                "water": [[10674,true],[2133,true]]
            }
        ],

        "Ch1:Ch2": {
            "FAM": [1.406863,-0.425347,1.633333,4.803922,-9.868056,-0.235926,0.594856,-4.9,-3.811111,6.186869,1.020833,0.340278,0.704023,-0.888889,-0.106456,1.696728],
            "HEX": [0.870219,0.629768,3.175926,3.296024,1.361111,-0.297743,0.506897,-1.341241,-7.712963,0.291667,1.841503,0.680556,-1.681373,-1.852004,-0.397863,2.807292]
        },
        "valid": true
    }"""; dicttype=OrderedDict)
    verify_response(
        Action_DICT["optical_test_dual_channel"](),
        response
    )
end

# error response body:

function error_dualchannel_optical_response_test()
    response=JSON.parse("""{
        "valid": false,
        "error": "xxxx"
    }"""; dicttype=OrderedDict)
    verify_response(
        Action_DICT["optical_test_dual_channel"](),
        response
    )
end

# ********************************************************************************
#
#                               R U N   E X A M P L E S
#
# ********************************************************************************

function verify_response_examples()
    examples = [
        :standard_curve_response_test,
        :singlechannel_amplification_response_test,
        :dualchannel_amplification_response_test,
        :error_amplification_response_test,
        :singlechannel_meltcurve_response_test,
        :dualchannel_meltcurve_response_test,
        :error_meltcurve_response_test,
        :loadscript_response_test,
        :error_loadscript_response_test,
        :thermal_performance_diagnostic_response_test,
        :thermal_consistency_response_test,
        :valid_optical_cal_response_test,
        :invalid_optical_cal_response_test,
        :singlechannel_optical_response_test,
        :dualchannel_optical_response_test,
        :error_dualchannel_optical_response_test
    ]
    FactCheck.clear_results()
    OrderedDict(map(examples) do f
        f =>
            try
                !getfield(Main,f)()
            catch err
                err
            end

    end)
end

# Usage:
# verify_response_examples() # every test should return true

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
"""
    Coulomb{T<:AbstractFloat} <: EmpiricalPotential{NamedTuple{()},NamedTuple{(:rcutoff,)}}

The Coulomb, or electrical potential is a simple two-body intermolecular potential that describes the electrical potential energy between two atoms with charge ``q_{1}`` and ``q_{2}``, at a distance ``r``, is given by

```math
\\begin{equation}
V_{LJ}(r; \\epsilon, \\sigma, rcutoff, species) =
\\begin{cases}
    0 & r > rcutoff \\\\
    \\frac{k_{e} q_{1} q_{2}}{r} & r < rcutoff.
\\end{cases}
\\end{equation}
```
where ``k_{e}`` is known as Coulomb's constant.

Users must supply the electric charges of the atoms, ``q_{1}`` and ``q_{2}``, and the radial cutoff.
"""
struct Coulomb{T<:AbstractFloat} <: EmpiricalPotential{NamedTuple{()},NamedTuple{(:rcutoff,)}}
    q₁::T
    q₂::T
    rcutoff::T
    species::Tuple
    function Coulomb(q₁, q₂, rcutoff, species)
        q₁, q₂, rcutoff = promote(austrip(q₁), austrip(q₂), austrip(rcutoff))
        new{typeof(rcutoff)}(q₁, q₂, rcutoff, Tuple(species))
    end
    function Coulomb{T}(; q₁::T, q₂::T, rcutoff::T, species::Tuple) where {T}
        new{T}(q₁, q₂, rcutoff, species)
    end
end

get_rcutoff(c::Coulomb) = c.rcutoff
get_species(c::Coulomb) = c.species

potential_energy(R::AbstractFloat, c::Coulomb) = kₑ * c.q₁ * c.q₂ / R
force(R::AbstractFloat, r::SVector{3}, c::Coulomb) = (kₑ * c.q₁ * c.q₂ / R^3)r

ext: jl


text:
using SolverPeeker
using Test
using JuMP

"""
Initiate some models
"""

# Model with linear constraints and objective
model1 = Model()
@variable(model1, x >= 0, Int)
@variable(model1, y <= 5)
@variable(model1, z)
@objective(model1, Max, 3x+5y-2z)
@constraint(model1, x+2y <= 45)
@constraint(model1, z - y >= 4)
@constraint(model1, 0<= 2z <= 7)

# Model with nothging specified. Should have objective 0.
model2 = Model()

# Objective registered as 0 in model(?)
model_nl = Model()
@variable(model_nl, x, start = 0.0)
@variable(model_nl, y, start = 0.0)
@NLobjective(model_nl, Min, (1 - x)^2 + 100 * (y - x^2)^2)
@constraint(model_nl, x + y == 10)
#@NLobjective(model_nl, Min, sin(x))

@testset "recognize" begin

    """
    Constant objective test
    """
    """
    @test recognize(model1, ObjectiveIsConstant()) == false

    @test recognize(model2, ObjectiveIsConstant()) == true

    # Add variable and constant objective. Should still be true:
    @variable(model2, x)
    @objective(model2, Min, 5)
    @test recognize(model2, ObjectiveIsConstant()) == true
    """

end

max_stars_repo_name: henriean/solver-peeker

max_issues_repo_name: henriean/solver-peeker

max_forks_repo_name: henriean/solver-peeker


text:
module BeamPropagation

using
    StaticArrays,
    Unitful,
    LoopVectorization,
    StatsBase,
    StructArrays

macro params(fields_tuple)
    fields = fields_tuple.args
    esc(
        quote
            NamedTuple{Tuple($fields)}(($fields_tuple))
        end)
end
export @params

macro with_unit(arg1, arg2)
    arg2 = @eval @u_str $arg2
    return convert(Float64, upreferred(eval(arg1) .* arg2).val)
end
export @with_unit

function dtstep_euler!(particles, f, abstol, p, dt_min, dt_max)

    @inbounds for i in 1:size(particles, 1)
        particles.a[i] = f(particles.idx[i], particles[i].r, particles[i].v, p)
    end

    @turbo for i in 1:size(particles, 1)
        r = particles.r[i]
        v = particles.v[i]
        a = particles.a[i]
        dt = particles.dt[i]

        # Perform an Euler method step (first order)
        particles.v[i] = v + a * dt
        particles.r[i] = r + v * dt
    end
    return nothing
end

function dtstep_eulerrich!(particles, f, abstol, p, dt_min, dt_max)

    @inbounds for i in 1:size(particles, 1)
        particles.a[i] = f(particles.idx[i], particles.r[i], particles.v[i], p)
    end

    @turbo for i in 1:size(particles, 1)
        r = particles.r[i]
        v = particles.v[i]
        a = particles.a[i]
        dt = particles.dt[i]

        # Perform an Euler method step (first order)
        particles.v1[i] = v + a * dt
        particles.r1[i] = r + v * dt
    end

    @inbounds for i in 1:size(particles, 1)
        particles.a[i] = f(particles.idx[i], particles.r1[i], particles[i].v1, p)
    end

    @turbo for i in 1:size(particles, 1)
        r = particles.r[i]
        v = particles.v[i]
        r1 = particles.r1[i]
        v1 = particles.v1[i]
        a = particles.a[i]
        dt = particles.dt[i]

        # Perform a Heun method step (second order)
        v2 = 0.5v + 0.5 * (v1 + a * dt)
        particles.v2[i] = v2
        particles.r2[i] = 0.5r + 0.5 * (r1 + v2 * dt)
    end

    @inbounds for i in 1:size(particles, 1)
        r1 = particles.r1[i]
        r2 = particles.r2[i]
        v2 = particles.v2[i]
        dt = particles.dt[i]

        if particles.use_adaptive[i]
            error = sqrt(sum((r2 - r1)).^2)
            particles.error[i] = error

            if !iszero(error)
                dt_corr = dt * 0.9 * sqrt(abstol / (2 * error))
            else
                dt_corr = dt_max
            end
            new_dt = min(max(dt_corr, dt_min), dt_max)
            particles.dt[i] = new_dt

            if error < abstol
                particles.r[i] = r2
                particles.v[i] = v2
            end
        else
            particles.r[i] = r2
            particles.v[i] = v2
        end

    end

    return nothing
end

mutable struct Particle
    r::SVector{3, Float64}
    v::SVector{3, Float64}
    a::SVector{3, Float64}
    dt::Float64
    use_adaptive::Bool
    error::Float64
    idx::Int64
    r1::SVector{3, Float64}
    v1::SVector{3, Float64}
    r2::SVector{3, Float64}
    v2::SVector{3, Float64}
    dead::Bool
end
export Particle

function save!(save, idxs, rs, vs, as, states, s)
    @inbounds for i in 1:size(rs,1)
        save(idxs[i], rs[i], vs[i], as[i], states[i], s)
    end
    return nothing
end

function initialize_dists_particles!(r, v, a, i_start, particles, dt, use_adaptive)
    @inbounds for i in 1:size(particles, 1)
        particles.r[i] = particles.r1[i] = SVector.(rand(r[1]), rand(r[2]), rand(r[3]))
        particles.v[i] = particles.v1[i] = SVector.(rand(v[1]), rand(v[2]), rand(v[3]))
        particles.a[i] = SVector.(rand(a[1]), rand(a[2]), rand(a[3]))
        particles.dead[i] = false
        particles.idx[i] = i_start + i - 1
        particles.dt[i] = dt
        particles.use_adaptive[i] = use_adaptive
        particles.error[i] = 0.0
    end
    return nothing
end
export initialize_dists_particles!

function discard_particles!(particles, discard)
    @inbounds for i in 1:size(particles, 1)
        particles.dead[i] = discard(particles.r[i], particles.v[i])
    end
    StructArrays.foreachfield(x -> deleteat!(x, particles.dead), particles)
    return nothing
end

# function propagate_particles!(r, v, a, alg, particles, f::F1, save::F2, discard::F3, save_every, delete_every, max_steps, update, p, s, dt, use_adaptive, dt_min, dt_max, abstol) where {F1, F2, F3}
#
#     initialize_dists_particles!(r, v, a, 1, particles, dt, use_adaptive)
#
#     step = 0
#     while (step <= max_steps)
#
#         update(particles, p, s, dt)
#
#         if step % save_every == 0
#             save(particles, p, s)
#         end
#
#         if step % delete_every == 0
#             discard_particles!(particles, discard)
#         end
#
#         if alg == "euler"
#             dtstep_euler!(particles, f, abstol, p, dt_min, dt_max)
#         elseif alg == "rkf12"
#             dtstep_eulerrich!(particles, f, abstol, p, dt_min, dt_max)
#         end
#
#         step += 1
#     end
#
#     return nothing
# end
# export propagate_particles!

function propagate_particles!(r, v, a, alg, particles, f::F1, save::F2, discard::F3, save_every, delete_every, max_steps, update, p, s, dt, use_adaptive, dt_min, dt_max, abstol) where {F1, F2, F3}

    n = length(particles)
    n_threads = Threads.nthreads()
    chunk_size = ceil(Int64, n / n_threads)

    Threads.@threads for i in 1:n_threads

        start_idx   = (i-1)*chunk_size+1
        end_idx     = min(i*chunk_size, n)
        chunk_idxs  = start_idx:end_idx
        actual_chunk_size = length(chunk_idxs)

        particles_chunk = particles[chunk_idxs]
        initialize_dists_particles!(r, v, a, start_idx, particles_chunk, dt, use_adaptive)

        for step in 0:(max_steps - 1)

            update(particles_chunk, p, s, dt)

            if step % save_every == 0
                save(particles_chunk, p, s)
            end

            if step % delete_every == 0
                discard_particles!(particles_chunk, discard)
            end

            if alg == "euler"
                dtstep_euler!(particles_chunk, f, abstol, p, dt_min, dt_max)
            elseif alg == "rkf12"
                dtstep_eulerrich!(particles_chunk, f, abstol, p, dt_min, dt_max)
            end

            if iszero(length(particles_chunk))
                break
            end
        end
    end

    return nothing
end
export propagate_particles!

end


text:
@testset "utils" begin
    @testset "counts" begin
        path = BGEN.datadir("example.8bits.bgen")
        b = Bgen(path)
        for v in iterator(b)
            cnt = counts!(b, v)
            dose = ref_allele_dosage!(b, v)
            correct_cnt = zeros(Int, 512)
            for v in dose
                if !isnan(v)
                    ind = convert(Int, round(v * 255)) + 1
                    correct_cnt[ind] += 1
                end
            end
            correct_cnt[512] = count(isnan.(dose))
            @test all(cnt .== correct_cnt)
        end
    end
end


text:
using Base.Threads
using Formatting
using OrderedCollections
using Statistics
using Dates

@with_kw mutable struct FarmZone <: AgComponent
    name::String
    climate::Any
    fields::Array{FarmField}
    water_sources::Array{WaterSource}
    _irrigation_volume_by_source::DataFrame = DataFrame()
    manager = nothing
end

# function Base.show(io::IO, z::FarmZone)
#     compact = get(io, :compact, false)

#     name = z.name
#     fields = Tuple((f.name, f.total_area_ha) for f in z.fields)
#     zone_ws = Tuple(ws.name for ws in z.water_sources)
#     manager = z.manager.name

#     if compact
#         print(z, " ($name)")
#     else
#         println("""
#         FarmZone ($name):

#         Fields: $(length(fields))

#         Water Sources: $(zone_ws)

#         Manager: $(manager)
#         """)
#     end
# end

"""Use allocation volume from a particular water source.

Volumes in ML.
"""
function use_allocation!(ws::WaterSource, vol_ML::Float64)::Nothing
    ws.allocation -= vol_ML

    if ws.allocation < 0.0
        if isapprox(0.0, ws.allocation, atol=1e-8)
            ws.allocation = 0.0
            return
        end

        ws_name = ws.name
        msg = "Allocation cannot be below 0 ML! Currently: $(ws.allocation)ML\n"
        msg *= "Tried to use: $(vol_ML)ML\n"
        msg *= "From: $ws_name\n"
        throw(ArgumentError(msg))
    end

    return nothing
end

function Base.getproperty(z::FarmZone, v::Symbol)
    if v == :avail_allocation
        """Available water allocation in ML."""
        all_allocs = [ws.allocation for ws in z.water_sources]
        return round(sum(all_allocs), digits=4)
    elseif v == :total_area_ha
        return sum([f.total_area_ha for f in z.fields])
    elseif v == :total_irrigated_area
        return sum([!isnothing(f.irrigated_area) ? f.irrigated_area : 0.0 for f in z.fields])
    end

    return getfield(z, v)
end

"""Determine the possible irrigation area using water from each water source."""
function possible_area_by_allocation(zone::FarmZone, field::FarmField, req_water_ML::Float64)::NamedTuple
    @assert in(field.name, [f.name for f::FarmField in zone.fields]) "Field must be in zone"

    zone_ws::Tuple = Tuple(zone.water_sources)
    tmp = NamedTuple{Tuple(Symbol(ws.name) for ws::WaterSource in zone_ws)}(
        possible_irrigation_area(field, ws.allocation, req_water_ML) for ws::WaterSource in zone_ws
    )

    return tmp
end

"""The total area marked for irrigation."""
function irrigated_area(zone::FarmZone)::Float64
    fields::Array{FarmField} = zone.fields
    return sum([f.irrigated_area for f in fields])
end

# """Calculate ML per ha to apply."""
# function calc_irrigation_water(zone::FarmZone, field::FarmField)::Float64
#     req_water_ML::Float64 = required_water(field)

#     # Can only apply water that is available
#     ML_per_ha::Float64 = zone.avail_allocation / field.irrigated_area
#     if req_water_ML < ML_per_ha
#         irrig_water = req_water_ML
#     else
#         irrig_water = ML_per_ha
#     end

#     return irrig_water
# end

"""Apply irrigation water to field"""
function apply_irrigation!(field::CropField,
                          ws::WaterSource, water_to_apply_mm::Float64,
                          area_to_apply::Float64)::Nothing
    vol_ML_ha::Float64 = water_to_apply_mm / mm_to_ML
    vol_ML::Float64 = vol_ML_ha * area_to_apply
    use_allocation!(ws, vol_ML)

    apply::Float64 = (water_to_apply_mm * field.irrigation.efficiency)
    field.soil_SWD::Union{Float64, AgParameter} -= max(0.0, apply)::Float64
    field.irrigated_volume = (ws.name, vol_ML)

    return nothing
end

"""Apply rainfall and ET to influence soil water deficit."""
function apply_rainfall!(zone::FarmZone, dt::Date)::Nothing
    climate::Climate = zone.climate::Climate
    data::DataFrame = climate.data::DataFrame
    idx::BitArray = (climate.time_steps::Array .== dt)::BitArray
    subset::DataFrame = data[idx, :]

    z_name::String = zone.name::String
    @inbounds for f::FarmField in zone.fields
        # get rainfall and et for datetime
        zf_id::String = "$(z_name)_$(f.name)"
        rain_col::Symbol = Symbol(zf_id, "_rainfall")
        et_col::Symbol = Symbol(zf_id, "_ET")

        rainfall::Float64, et::Float64 = subset[1, rain_col], subset[1, et_col]

        update_SWD!(f, rainfall, et)
    end

    return nothing
end

"""Apply rainfall and ET to influence soil water deficit."""
function apply_rainfall!(zone::FarmZone, dt_idx::Int64)::Nothing
    data::DataFrame = zone.climate.data::DataFrame

    z_name::String = zone.name::String
    @inbounds for f::FarmField in zone.fields
        # get rainfall and et for datetime
        zf_id::String = "$(z_name)_$(f.name)"
        rain_col::Symbol = Symbol(zf_id, "_rainfall")
        et_col::Symbol = Symbol(zf_id, "_ET")

        rainfall::Float64, et::Float64 = data[dt_idx, rain_col], data[dt_idx, et_col]

        update_SWD!(f, rainfall, et)
    end

    return nothing
end

"""Update water allocations"""
function update_available_water!(zone::FarmZone, allocations::NamedTuple)::Nothing
    for (k, v) in pairs(allocations)
        for ws in zone.water_sources
            if ws.name == string(k)
                ws.allocation = v
            end
        end
    end
end

"""Log irrigation volumes from water sources"""
function log_irrigation_by_water_source(zone::FarmZone, f::FarmField, dt::Date)::Nothing

    # Construct log structure if needed
    if nrow(zone._irrigation_volume_by_source) == 0
        tmp_dict::OrderedDict = OrderedDict()
        tmp_dict[:Date] = Date[]
        for ws::WaterSource in zone.water_sources
            tmp_dict[Symbol(ws.name)] = Float64[]
        end

        zone._irrigation_volume_by_source = DataFrame(; tmp_dict...)
    end

    tmp::Array{Float64} = Float64[0.0 for ws in zone.water_sources]
    for (i::Int64, ws::WaterSource) in enumerate(zone.water_sources)
        try
            tmp[i] += f.irrigation_from_source[ws.name]
        catch e
            if isa(e, KeyError)
                continue
            end
            throw(e)
        end
    end

    push!(zone._irrigation_volume_by_source, [dt tmp...])

    return nothing
end

"""Collate logged values, summing on identical datetimes"""
function collate_field_logs(zone::FarmZone, sym::Symbol; last=false)::OrderedDict{Date, Float64}
    target_log::Dict{Date, Float64} = Dict{Date, Float64}()
    for f::FarmField in zone.fields
        tmp = getfield(f, sym)

        if last
            tmp = Dict(sort(collect(tmp))[end])
        end

        target_log = merge(+, target_log, tmp)
    end

    return OrderedDict(sort(collect(target_log)))
end

function aggregate_field_logs(field_logs::DataFrame)::DataFrame
    collated::DataFrame = aggregate(groupby(field_logs, :Date), sum)

    collated[:, Symbol("Dollar per ML")] = collated[:, :income_sum] ./ collated[:, :irrigated_volume_sum]
    collated[:, Symbol("ML per Irrigated Yield")] = collated[:, :irrigated_volume_sum] ./ collated[:, :irrigated_yield_sum]
    collated[:, Symbol("Dollar per Ha")] = collated[:, :income_sum] ./ (collated[:, :dryland_area_sum] .+ collated[:, :irrigated_area_sum])
    collated[:, Symbol("Mean Irrigated Yield")] = collated[:, :irrigated_yield_sum] ./ collated[:, :irrigated_area_sum]
    collated[:, Symbol("Mean Dryland Yield")] = collated[:, :dryland_yield_sum] ./ collated[:, :dryland_area_sum]

    collated[isnan.(collated[!,Symbol("Mean Dryland Yield")]), Symbol("Mean Dryland Yield")] .= 0
    collated[isnan.(collated[!,Symbol("Mean Irrigated Yield")]), Symbol("Mean Irrigated Yield")] .= 0

    return collated
end

"""Collate logged values, aggregating to the zonal level based on identical datetimes."""
function collate_field_logs(zone::FarmZone)::DataFrame
    tmp::DataFrame = reduce(vcat, [f._seasonal_log for f in zone.fields])
    collated::DataFrame = aggregate_field_logs(tmp)

    return collated
end

"""Collate logged values, aggregating to the zonal level based on identical datetimes."""
function collate_field_logs(seasonal_logs::Dict)::DataFrame
    s_logs = values(seasonal_logs)
    tmp::DataFrame = reduce(vcat, s_logs)
    collated::DataFrame = aggregate_field_logs(tmp)

    return collated
end

function water_used_by_source(zone::FarmZone)::DataFrame
    return water_used_by_source(zone, nothing)
end

function water_used_by_source(zone::FarmZone, dt)::DataFrame
    ws_irrig = zone._irrigation_volume_by_source

    if nrow(ws_irrig) == 0
        return DataFrame(Date=dt, surface_water_sum=0.0, groundwater_sum=0.0)
    end

    if !(isnothing(dt))
        ws_irrig = ws_irrig[ws_irrig[:Date] .== dt, :]
    end

    if nrow(ws_irrig) == 0
        if isnothing(dt)
            return aggregate(groupby(ws_irrig, :Date), sum)
        end

        # Catch empty subset
        return DataFrame(Date=dt, surface_water_sum=0.0, groundwater_sum=0.0)
    end

    return aggregate(groupby(ws_irrig, :Date), sum)
end

"""Collect model run results for a FarmZone"""
function collect_results(zone::FarmZone; last=false)::NamedTuple
    field_logs = Dict(
        f.name => f._seasonal_log
        for f in zone.fields
    )

    collated = collate_field_logs(field_logs)

    irrig_ws::DataFrame = water_used_by_source(zone)

    collated_res = hcat(collated, irrig_ws[:, setdiff(names(irrig_ws), [:Date])])

    return (zone_results=collated_res, field_results=field_logs)
end

function create(data::Dict{Symbol}, climate_data::Climate)::FarmZone
    spec = deepcopy(data)

    if haskey(spec, :component)
        _ = pop!(spec, :component)
    end

    name = spec[:name]

    if haskey(spec, :water_sources)
        water_specs::Dict{Symbol, Dict} = spec[:water_sources]
        pump_specs::Dict{Symbol, Dict} = spec[:pump_spec]
        water_sources::Array{WaterSource} = create(water_specs, pump_specs)
    else
        water_sources = []
    end

    # This will be used in future to provide list of irrigations/crops that could be considered
    # crop_specs::Dict{Symbol, Dict} = spec[:crop_spec]
    # irrig_spec::Dict{Symbol, Dict} = spec[:irrigation_spec]

    field_specs = deepcopy(spec[:fields])
    for (fk, f) in field_specs
        f[:irrigation] = create(collect(values(f[:irrigation_spec]))[1])
        f[:crop_rotation] = [create(c, climate_data.time_steps[1])
                             for c in collect(values(f[:crop_rotation_spec]))]
        f[:crop] = f[:crop_rotation][1]

        # Clean up unneeded specs
        delete!(f, :irrigation_spec)
        delete!(f, :crop_rotation_spec)
    end

    fields = [create(v) for (k,v) in field_specs]

    zone_spec::Dict{Symbol, Any} = Dict(
        :name => name,
        :climate => climate_data,
        :fields => fields,
        :water_sources => water_sources
    )

    zone = FarmZone(; zone_spec...)

    return zone
end

function reset!(z::FarmZone)::Nothing
    for f in z.fields
        reset!(f)

        initial_dt = z.climate.time_steps[1]

        f._next_crop_idx = 1
        set_next_crop!(z.manager, f, initial_dt)

    end

    return
end

max_stars_repo_name: ConnectedSystems/agtor-jl

max_issues_repo_name: ConnectedSystems/agtor-jl

max_forks_repo_name: ConnectedSystems/agtor-jl


text:
# 3 dimensionsal evaluation functions.
# Transforms code is scattered, and will be moved to better location with stable 3D Examples
# also worth it to consolidate if other affine transforms package is available

# using Euler angles for linear sampling in product of potentials
# Gaussian model for prior
function evalPotential(obs::PriorPose3, Xi::Array{Graphs.ExVertex,1}; N::Int64=200)
  mu = veeEuler(obs.Zi)
  return rand( MvNormal(mu, obs.Cov), N)
end

# Project all particles (columns) Xval with Z, that is for all  SE3(Xval[:,i])*Z
function projectParticles(Xval::Array{Float64,2}, Z::SE3, Cov::Array{Float64,2})
  # TODO optimize for more speed with threads and better memory management
  r,c = size(Xval)
  RES = zeros(r,c) #*cz

  ent, x = SE3(0), SE3(0)
  ENT = rand( MvNormal(zeros(6), Cov), c )

  j=1
  # for j in 1:cz
    for i in 1:c
      x.R, x.t = convert(SO3,Euler(Xval[4:6,i])), Xval[1:3,i]
      ent.R, ent.t = convert(SO3,so3(ENT[4:6,i])), ENT[1:3,i]
      newval = x*Z*ent
      RES[1:r,i*j] = veeEuler(newval)
    end
  # end
  #
  return RES
end

# Still limited to linear sampler, then reprojected onto ball -- TODO upgrade manifold sampler
function evalPotential(odom::Pose3Pose3, Xi::Array{Graphs.ExVertex,1}, Xid::Int64; N::Int64=200)
    # rz,cz = size(odom.Zij)
    Xval = Array{Float64,2}()
    # implicit equation portion -- bi-directional pairwise function made explicit here
    if Xid == Xi[1].index #odom.
        # reverse direction
        println(" reverse direction")
        Z = inverse(odom.Zij)
        Xval = getVal(Xi[2])
    elseif Xid == Xi[2].index
        # forward direction
        println(" forward direction")
        Z = odom.Zij
        Xval = getVal(Xi[1])
    else
        error("Bad evalPairwise Pose3Pose3")
    end

    return projectParticles(Xval, Z, odom.Cov)
end


text:
## variate_size ##############################################################

"""
    variate_size(distribution) :: Dims

Return the `size` of each sample from `distribution`, when it can be
determined in advance. This is used in particular by `Pack`.

# Examples
```julia-repl
julia> variate_size(Normal())
()

julia> variate_size(Fill(Normal(), 2, 3))
(2, 3)

julia> rand(SubSeq(1:9, 0.3), 2)
2-element Array{Array{Int64,1},1}:
 [4, 7]
 [5]

julia> variate_size(SubSeq(1:9, 0.3)) # no pre-determined size
ERROR: MethodError: no method matching variate_size(::SubSeq{Array{Int64,1},UnitRange{Int64}})
[...]
```
"""
function variate_size end

variate_size(d::Distribution{<:Number}) = ()

## Zip #######################################################################

# note: Zip(a, b...) is similar to RandomExtensions.make(Tuple, a, b...)

struct Zip{T,X} <: Distribution{T}
    xs::X

    function Zip(as...)
        T = Tuple{map(gentype, as)...}
        xs = map(a -> a isa DataType ? Uniform(a) : a, as)
        new{T,typeof(xs)}(xs)
    end
end

support(z::Zip) = Iterators.product((support(x) for x in z.xs)...)

function pmf(z::Zip, xs::Tuple)
    length(xs) == length(z.xs) || return 0.0
    prod(pmf(z.xs[i], xs[i]) for i in 1:length(xs))
end

function Base.show(io::IO, z::Zip)
    print(io, "Zip(")
    join(io, z.xs, ", ")
    print(io, ')')
end

### sampling

Sampler(::Type{RNG}, z::Zip, n::Repetition) where {RNG<:AbstractRNG} =
    SamplerTag{typeof(z)}(map(x -> sampler(RNG, x, n), z.xs))

function reset!(sp::SamplerTag{<:Zip}, n...)
    foreach(sp.data) do x
        reset!(x, n...)
    end
    sp
end

rand(rng::AbstractRNG, sp::SamplerTag{<:Zip{T}}) where {T} =
    map(x -> rand(rng, x), sp.data)::T

function rand!(rng::AbstractRNG, t::T,
               sp::SamplerTag{<:Zip{S}},
               ::Val{N}) where N where S where T<:Tuple
    N < 2 && throw(ArgumentError("can not mutate a tuple"))
    fieldcount(T) != fieldcount(S) &&
        throw(ArgumentError("tuples must have same size"))

    for (x, s) in zip(t, sp.data)
        rand!(rng, x, s, Val(N-1))
    end
    t
end

## Fill ######################################################################

struct Fill{X,T,N} <: Distribution{Array{T,N}}
    x::X
    dims::Dims{N}

    Fill(x::X, dims::Dims{N}) where {X,N} = new{X,gentype(x),N}(x, dims)
end

Fill(x, dims::Integer...) where {X} = Fill(x, Dims(dims))

Fill(::Type{X}, dims::Dims{N})    where {X,N} = Fill(Uniform(X), dims)
Fill(::Type{X}, dims::Integer...) where {X}   = Fill(Uniform(X), Dims(dims))

variate_size(f::Fill) = f.dims

Sampler(RNG::Type{<:AbstractRNG}, f::Fill, n::Repetition) =
    SamplerTag{typeof(f)}((x    = sampler(RNG, f.x, Val(Inf)),
                           dims = f.dims))

rand(rng::AbstractRNG, sp::SamplerTag{<:Fill}) =
    rand!(rng, gentype(sp)(undef, sp.data.dims), sp, Val(1))

function rand!(rng::AbstractRNG, a::AbstractArray, sp::SamplerTag{<:Fill}, ::Val{N}) where N
    dims = sp.data.dims
    if dims != size(a)
        if length(dims) == 1 && a isa AbstractVector
            resize!(a, dims[1])
        else
            throw(ArgumentError("can not resize destination array"))
        end
    end
    x = reset!(sp.data.x, prod(dims))
    for i in eachindex(a)
        if N == 1
            @inbounds a[i] = rand(rng, x)
        else
            @inbounds rand!(rng, a[i], x, Val(N-1))
        end
    end
    return a
end

# TODO: create a lazy fpmf dict which doesn't store all values, but only
# the pmf of f.x, it's then cheap to compute values and keys on the fly
function pmf(f::Fill; normalized=true)
    A = eltype(f)
    xf = pmf(f.x, normalized=normalized)
    fpmf = Dict{A,Float64}()
    n = prod(f.dims)
    # TODO: don't use Iterators.product, which would need "static" n
    for xs in Iterators.product(Iterators.repeated(xf, n)...)
        e = A(undef, f.dims)
        a = copyto!(e, first.(xs))
        @assert !haskey(fpmf, a)
        fpmf[a] = prod(last, xs)
    end
    PMF(f, fpmf, normalized=normalized; count=abs(xf.count)^n)
end

function Base.show(io::IO, f::Fill)
    print(io, "Fill(", f.x, ", ")
    join(io, f.dims, ", ")
    print(io, ')')
end

## Pack ######################################################################

struct Pack{X,T,N} <: Distribution{Array{T,N}}
    x::X
    dims::Dims{N} # (inner_dims..., outer_dims...)
    dim::Int # length of inner_dims

    function Pack(x::X, dims::Dims{N}) where {X,N}
        A = gentype(x)
        inner = variate_size(x)
        dims = (inner..., dims...)
        new{X,eltype(gentype(x)),length(dims)}(x, dims, length(inner))
    end
end

Pack(x, dims::Integer...) where {X} = Pack(x, Dims(dims))

Pack(::Type{X}, dims::Dims{N})    where {X,N} = Pack(Uniform(X), dims)
Pack(::Type{X}, dims::Integer...) where {X}   = Pack(Uniform(X), Dims(dims))

function Base.show(io::IO, p::Pack)
    len = length(variate_size(p.x))
    print(io, "Pack(", p.x, ", ")
    join(io, p.dims[len+1:end], ", ")
    print(io, ')')
end

Sampler(RNG::Type{<:AbstractRNG}, p::Pack, n::Repetition) =
    SamplerTag{typeof(p)}((x    = sampler(RNG, p.x, Val(Inf)),
                           dims = p.dims,
                           dim  = p.dim))

function rand!(rng::AbstractRNG, A::AbstractArray{T,N},
               sp::SamplerTag{Pack{X,T,N}}, ::Val{1}) where {X,T,N}

    size(A) == sp.data.dims || throw(DimensionMismatch(
    "size of destination array ($(size(A))) does not match" *
        "dimensions of Pack ($(sp.data.dims))"))

    inner_idxs = CartesianIndices(sp.data.dims[1:sp.data.dim])
    idxs = CartesianIndices(sp.data.dims[sp.data.dim+1:end])
    x = sp.data.x

    reset!(x, length(idxs))

    if gentype(x) <: AbstractArray
        for idx in idxs
            rand!(rng, view(A, inner_idxs, idx), x, Val(1))
        end
    else
        for idx in idxs # similar to Fill
            @inbounds A[idx] = rand(rng, x)
        end
    end
    A
end

rand(rng::AbstractRNG, sp::SamplerTag{Pack{X,T,N}}) where {X,T,N} =
    rand!(rng, Array{T,N}(undef, sp.data.dims), sp, Val(1))


text:
module SerendipityRefFEsTests

using Test
using InplaceArrays.Fields
using InplaceArrays.ReferenceFEs

reffe = SerendipityRefFE(Float64,QUAD,2)
@test reffe.data.dofs.nodes == Point{2,Float64}[
  (0.0, 0.0), (1.0, 0.0), (0.0, 1.0), (1.0, 1.0),
  (0.5, 0.0), (0.5, 1.0), (0.0, 0.5), (1.0, 0.5)]
@test reffe.face_own_nodeids == [[1], [2], [3], [4], [5], [6], [7], [8], Int[]]

reffe = SerendipityRefFE(Float64,QUAD,4)
@test reffe.face_own_nodeids == [[1], [2], [3], [4], [5, 6, 7], [8, 9, 10], [11, 12, 13], [14, 15, 16], [17]]

reffe = SerendipityRefFE(Float64,HEX,2)
@test reffe.face_own_nodeids == [
  [1], [2], [3], [4], [5], [6], [7], [8],
  [9], [10], [11], [12], [13], [14], [15], [16],
  [17], [18], [19], [20], Int[], Int[], Int[], Int[], Int[], Int[], Int[]]

end # module


text:
function makefigure4(dist::DataFrame;)

  ###Common Data and Parameters

  local sdist::SubDataFrame = view(dist, dist.primary, :)

  F4_STYLE_IDX::Dict = Dict(:textsize => 10.0,
      :linesize => 0.1,
      :shapesize => 0.25,
      :xdomain => [-12, 12],
      :width => 400,
      :height => 300,
      :xaxisname=>"Payment Evt Day",
      :legendtextsize => 9.0,
      :yaxisname => "#",
      :ydomain => [0.,15.],
      :seriescolor => :black)

  ###Panel C
  F4_FX = :day
  F4_FOCAL_FIELDS = [:lnumtrd, :ldolvol]

  figure4(sdist,
    Fx = F4_FX,
    focalfields = F4_FOCAL_FIELDS,
    styleidx = F4_STYLE_IDX)

end

function figure4(dist::AbstractDataFrame;
  focalfields::Vector{Symbol} = error("focalfields is required"),
  Fx::Symbol = error("Fx is required"),
  styleidx::Dict = error("style idx is required"))

  figure4fullpath::String = "$FIGURE_PATH\\figure4-$(REPLICATION_TYPE[])-$(DIST_TYPE[])-$(OUT_SUFFIX[]).pdf"

  agg::DataFrame = aggregatebyday(v->(mean(skipmissing(v))), dist, selectcols=[Fx; focalfields])
  long::DataFrame = melt(agg, Fx, focalfields)

  #maybe add both the top dividend rate and the marginal rate?
  figure4 = makefigure4graph(long,
    Fy = :value,
    Fx = Fx,
    Fgroup = :variable,
    styleidx=styleidx)
  save(figure4fullpath, figure4)
end

function makefigure4graph(df::DataFrame;
  Fy::Symbol = error("Fy is required"),
  Fx::Symbol = error("Fx is required"),
  Fgroup = error("Fgroup is required"),
  styleidx::Dict = error("styleidx is required"))

  ###NOTE: Begin the JSON for the plot
  #some plot parameters are below

  series = @vlplot() +
    @vlplot( #first part contains the line and circle dots
      data = df,
      width=styleidx[:width],
      height=styleidx[:height],
      title=nothing,
      mark={:line,
        clip=true,
        size=styleidx[:linesize]},
      x={Fx,
        axis={title=nothing,#styleidx[:xaxisname],
          formatType = :number,
          format= ".4",
          grid=false,
          titleFontSize=styleidx[:textsize],
          tickCount=12,
          titleFontWeight=:normal,
          labels=true,
          labelOpacity = {
            condition = {test = "datum.$(Fx) % 10 == 0",
              value = 0.0},
            value = 1.0}},
        type=:quantitative,
        scale={domain=styleidx[:xdomain]}},
      y={Fy,
        axis={title=styleidx[:yaxisname],
            grid=false,
            tickCount=4,
            titleFontSize=styleidx[:textsize],
            titleFontWeight=:normal},
          type=:quantitative,
          scale={zero=false,
            domain=styleidx[:ydomain]}
        },
      color={value=styleidx[:seriescolor]},
      detail={Fgroup,
        type=:nominal}
    ) +
    @vlplot(
      data = df,
      mark={:point,
        clip=true,
        size=styleidx[:shapesize]},
      y={Fy,
        type=:quantitative},
      x={Fx,
        type=:quantitative,
        scale={domain=styleidx[:xdomain]}},
      shape={Fgroup,
        type=:nominal,
        legend=nothing},
      color={value=styleidx[:seriescolor]},
    ) + @vlplot( #handles horizontal text
      data={values=[
        {xloc=-11, tx="Log Dollar Volume of Trades", yloc=3},
        {xloc=2, tx="Log Number of Trades", yloc=12}]},
      mark={:text,
        size=styleidx[:textsize],
        align=:left},
      x={:xloc,
        type=:quantitative},
      y={:yloc,
        type=:quantitative},
      color = {value=:black},
      text = {:tx, type=:nominal}
    )

  return series
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
@testset "ReLU Last Layer - Inactive" begin

    small_nnet_file = "$(@__DIR__)/../examples/networks/small_nnet_id.nnet"

    small_nnet = read_nnet(small_nnet_file, last_layer_activation = ReLU())

    # The input set is always [-0.9:0.9]
    in_hyper  = Hyperrectangle(low = [-0.9], high = [0.9])
    in_hpoly  = convert(HPolytope, in_hyper)

    # superset of the output
    out_superset    = Hyperrectangle(low = [-10.0], high = [10.0])
    # includes points in the output region
    out_overlapping = Hyperrectangle(low = [-10.0], high = [-5.0])

    @testset "Group 1" begin
        problem_holds    = Problem(small_nnet, in_hpoly, convert(HPolytope, out_superset))
        problem_violated = Problem(small_nnet, in_hpoly, convert(HPolytope, out_overlapping))

        for solver in [MaxSens(resolution = 0.6), ExactReach(), Ai2(), Ai2h(), Box()]
            holds    = solve(solver, problem_holds)
            violated = solve(solver, problem_violated)

            @testset "$(typeof(solver))" begin
                @test holds.status    ∈ (:holds, :unknown)
                @test violated.status ∈ (:violated, :unknown)
            end
        end

    end

    @testset "Group 2, 3, 4, 6" begin

        problem_holds    = Problem(small_nnet, in_hyper, HPolytope([HalfSpace([1.], 10.)]))     # y < 10.0
        problem_violated = Problem(small_nnet, in_hyper, HPolytope([HalfSpace([-1.], -20.)]))     # y > 20.0

        group2 = [NSVerify(), MIPVerify(), ILP()]
        group3 = [ConvDual(), Duality()]
        group4 = [FastLin(), FastLip()]
        group6 = [Reluplex(), Planet()]

        for solver in [group2; group3; group4; group6]
            holds    = solve(solver, problem_holds)
            violated = solve(solver, problem_violated)

            @testset "$(typeof(solver))" begin
                @test holds.status    ∈ (:holds, :unknown)
                @test violated.status ∈ (:violated, :unknown)
            end
        end
    end

    @testset "Group 5" begin
        problem_holds    = Problem(small_nnet, in_hyper, out_superset)
        problem_violated = Problem(small_nnet, in_hyper, out_overlapping)

        for solver in [ReluVal(max_iter = 10), DLV(), Sherlock(ϵ = 0.5), BaB(), Neurify()]
            holds    = solve(solver, problem_holds)
            violated = solve(solver, problem_violated)

            @testset "$(typeof(solver))" begin
                @test holds.status    ∈ (:holds, :unknown)
                @test violated.status ∈ (:violated, :unknown)
            end
        end
    end

    @testset "Certify" begin
        ### Certify - only works for single hidden layer
        tiny_nnet = read_nnet("$(@__DIR__)/../examples/networks/tiny_nnet.nnet")
        solver_certify = Certify()
        in_set  = Hyperrectangle([2.0], [.5])
        out_set = HPolytope(ones(1,1), [2.5])
        problem_certify = Problem(tiny_nnet, in_set, out_set)
        @test @no_error solve(solver_certify, problem_certify)
    end
end


text:
isdefined(Base, :__precompile__) && __precompile__()

module Dubins

using Memento
using Compat

import Compat: @__MODULE__

const LOGGER = getlogger(@__MODULE__)
setlevel!(LOGGER, "info")

__init__() = Memento.register(LOGGER)

include("typedefs.jl")
include("paths.jl")
include("path_fcns.jl")

end


text:
module Generate

import Downloads
import JSON3
import StructTypes
import ZipFile

export download_fhir_json_schema
export output_fhir_types

include("types.jl")

include(joinpath(dirname(dirname(@__DIR__)), "src", "fhir-to-julia.jl"))
include("download.jl")
include("generate-fhir-types.jl")
include("map-fhir-types.jl")

end # end module Generate


text:
struct AwaitingInput end
struct CodeFinished end
const IntcodeOutput = Union{Int64, AwaitingInput, CodeFinished}

mutable struct Tape
    tape::Array{Int64,1}
    pos::Int64
    relbase::Int64
    inch::Channel{Int64}
    outch::Channel{IntcodeOutput}
end

function Tape(tape, pos, relbase)
    Tape(tape, pos, relbase, Channel{Int64}(1), Channel{IntcodeOutput}(1))
end

function Tape(tape)
    Tape(tape, 0, 0)
end

function Base.copy(tape)
    Tape(copy(tape.tape), tape.pos, tape.relbase)
end

function read(tape::Tape, pos)
    return tape.tape[pos + 1]
end

function write!(tape::Tape, pos, value)
    tape.tape[pos + 1] = value
end

const POSITION_MODE = 0
const IMMEDIATE_MODE = 1
const RELATIVE_MODE = 2

function get_op_mode(tape::Tape, offset)
    instruction = read(tape, tape.pos)
    param_mode = div(instruction, 10^(1+offset)) % 10
end

function get_arg(tape::Tape, offset)
    mode = get_op_mode(tape, offset)
    if mode == POSITION_MODE
        pos_arg = read(tape, tape.pos + offset)
        arg = read(tape, pos_arg)
    elseif mode == IMMEDIATE_MODE
        arg = read(tape, tape.pos + offset)
    elseif mode == RELATIVE_MODE
        pos_arg = read(tape, tape.pos + offset)
        arg = read(tape, tape.relbase + pos_arg)
    else
        error("invalid mode: " * mode * " !")
    end
    arg
end

function set_arg!(tape::Tape, offset, value)
    mode = get_op_mode(tape, offset)
    pos = read(tape, tape.pos + offset)
    if mode == RELATIVE_MODE
        pos += tape.relbase
    end
    write!(tape, pos, value)
end

function op_add!(tape)
    arg1 = get_arg(tape, 1)
    arg2 = get_arg(tape, 2)
    result = arg1 + arg2
    set_arg!(tape, 3, result)
    tape.pos += 4
end

function op_mul!(tape)
    arg1 = get_arg(tape, 1)
    arg2 = get_arg(tape, 2)
    result = arg1 * arg2
    set_arg!(tape, 3, result)
    tape.pos += 4
end

function op_input!(tape)
    put!(tape.outch, AwaitingInput())
    input = take!(tape.inch)
    set_arg!(tape, 1, input)
    tape.pos += 2
end

function op_output!(tape)
    output = get_arg(tape, 1)
    put!(tape.outch, output)
    tape.pos += 2
end

function op_jumptrue!(tape)
    arg1 = get_arg(tape, 1)
    arg2 = get_arg(tape, 2)
    if arg1 != 0
        tape.pos = arg2
    else
        tape.pos += 3
    end
end

function op_jumpfalse!(tape)
    arg1 = get_arg(tape, 1)
    arg2 = get_arg(tape, 2)
    if arg1 == 0
        tape.pos = arg2
    else
        tape.pos += 3
    end
end

function op_lessthan!(tape)
    arg1 = get_arg(tape, 1)
    arg2 = get_arg(tape, 2)
    set_arg!(tape, 3, arg1 < arg2 ? 1 : 0)
    tape.pos += 4
end

function op_equal!(tape)
    arg1 = get_arg(tape, 1)
    arg2 = get_arg(tape, 2)
    set_arg!(tape, 3, arg1 == arg2 ? 1 : 0)
    tape.pos += 4
end

function op_adjustrelbase!(tape)
    arg = get_arg(tape, 1)
    tape.relbase += arg
    tape.pos += 2
end

function run_opcode!(tape)
    instruction = read(tape, tape.pos)
    opcode = instruction % 100
    if opcode == 99     return opcode
    elseif opcode == 1  op_add!(tape)
    elseif opcode == 2  op_mul!(tape)
    elseif opcode == 3  op_input!(tape)
    elseif opcode == 4  op_output!(tape)
    elseif opcode == 5  op_jumptrue!(tape)
    elseif opcode == 6  op_jumpfalse!(tape)
    elseif opcode == 7  op_lessthan!(tape)
    elseif opcode == 8  op_equal!(tape)
    elseif opcode == 9  op_adjustrelbase!(tape)
    else
        error("unknown opcode " * string(opcode) * "!")
    end
    opcode
end

function run_program!(tape)
    while true
        opcode = run_opcode!(tape)
        if opcode == 99
            put!(tape.outch, CodeFinished())
            return tape
        end
    end
end

function set_memory!(tape, size)
    memory = zeros(Int64, size)
    memory[1:length(tape.tape)] = tape.tape
    tape.tape = memory
    tape
end

function interact!(tape)
    while true
        out = take!(tape.outch)
        if isa(out, AwaitingInput)
            print("INPUT: ")
            input = readline()
            put!(tape.inch, parse(Int64, input))
        elseif isa(out, CodeFinished)
            println("FINISHED")
            return
        else
            println("OUTPUT: ", out)
        end
    end
end

function parse_data(program_str)
    [parse(Int64, v) for v in split(program_str, ',')]
end

max_stars_repo_name: st--/advent-of-code

max_issues_repo_name: st--/advent-of-code

max_forks_repo_name: st--/advent-of-code


text:
module TestBincount

using Test
import TreeParzen: Bincounts

@test_throws DimensionMismatch Bincounts.bincount([1], [1.0, 2.0], 3)

@test Bincounts.bincount([1], 3) == [1.0, 0.0, 0.0]

end
true


text:
# This file is part of Kpax3. License is MIT.

# suppose we
# a) split a cluster with 15 units into two clusters of 8 and 7 units
#    respectively, moving from k = 5 to k = 6
# b) move a unit from a cluster with 15 units to a new cluster (split)
# c) move a unit from a cluster with 15 units to another cluster with 7 units

function test_mcmc_partition_ratios()
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  settings = Kpax3.KSettings(ifile, ofile)

  x = Kpax3.AminoAcidData(settings)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(x.data, settings.γ, settings.r)

  state = Kpax3.AminoAcidState(x.data, [1; 1; 1; 1; 1; 1], priorR, priorC, settings)
  support = Kpax3.MCMCSupport(state, priorC)

  support.vi = 8
  support.vj = 7

  for (α, θ) in ((0.4, -0.3), (0.4, 0.0), (0.4, 2.1), (0.0, 2.1), (-2.4, 10))
    ep = Kpax3.EwensPitman(α, θ)

    lr1 = Kpax3.logdPriorRow(50, 6, [22;  8; 7; 5; 1; 7], ep) - Kpax3.logdPriorRow(50, 5, [22; 15; 7; 5; 1], ep)

    Kpax3.logratiopriorrowsplit!(6, ep, support)
    @test isapprox(support.lograR, lr1, atol=ε)

    Kpax3.logratiopriorrowmerge!(5, ep, support)
    @test isapprox(support.lograR, -lr1, atol=ε)

    lr2 = Kpax3.logdPriorRow(50, 6, [22; 14; 7; 5; 1; 1], ep) - Kpax3.logdPriorRow(50, 5, [22; 15; 7; 5; 1], ep)

    @test isapprox(Kpax3.logratiopriorrowsplit(6, 15, ep), lr2, atol=ε)
    @test isapprox(Kpax3.logratiopriorrowmerge(5, 14, ep), -lr2, atol=ε)

    lr3 = Kpax3.logdPriorRow(50, 5, [22; 14; 8; 5; 1], ep) - Kpax3.logdPriorRow(50, 5, [22; 15; 7; 5; 1], ep)

    @test isapprox(Kpax3.logratiopriorrowmove(15, 7, ep), lr3, atol=ε)
  end

  nothing
end

test_mcmc_partition_ratios()


text:
read_data() = readlines("day08.txt")

# read_data() = split("""
# nop +0
# acc +1
# jmp +4
# acc +3
# jmp -3
# acc -99
# acc +1
# jmp -4
# acc +6""", "\n")

function parse_instructions()
    lines = read_data()
    instructions = split.(lines, " ")
    return map(x -> (Symbol(x[1]), parse(Int, x[2])), instructions)
end

function execute(instructions)
    val = 0
    len = length(instructions)
    visited = falses(len)
    i = 1
    @inbounds while true
        # needed by Part2
        i > len && return (val, :terminated)

        # Next instruction
        instr = instructions[i]
        visited[i] && return (val, :infinite)

        # Remember this instruction
        # push!(visited, i)
        visited[i] = true

        # Execute!
        if instr[1] === :acc
            val += instr[2]
        elseif instr[1] === :jmp
            i += instr[2] - 1
        end
        i += 1
    end
end

part1() = execute(parse_instructions())

# part 2 brute force

function part2()
    prog = parse_instructions()
    sz = length(prog)
    jmps = findall(x -> x[1] == :jmp, prog)
    for j in jmps
        np = modified_program(prog, j, :nop)
        result = execute(np)
        result[2] == :terminated && return result[1]
    end
    # NOTE:
    # I didn't have to complete the code by replacing nop with jmp because I have
    # already reached an answer. Otherwise, it would be simialr code as follows.
    # nops = [i for i in 1:sz if prog[i][1] == :nop && prog[i][2] != 0]
    return "not found"
end

function modified_program(prog, i, chg)
    prog = copy(prog)
    prog[i] = (chg, prog[i][2])
    return prog
end


text:
function loopwinFixsize(nSNP, nReads, rawdata, snpID, bin2geno, nPer,
                        SNPSize, flagAnnotate, snpWtType, keepIdx,
                        nPerKeep, kernel, kinship, test,
                        XtNullBasis, WPreSim, nBlockAscent, nMMmax,
                        nNullSimPts, tolX, pvalueComputing, device,
                        nNullSimNewtonIter, rankQPhi, evalPhiAdj,
                        XPhitNullBasis, yShift, y, KPhiAdj,
                        weightedW, QPhi, PrePartialSumW,
                        PreTotalSumW, windowSize, nPreRank, X, ynew)

  QRes = Array(Float64, nPerKeep, rankQPhi);
  nSNPredVec = fill(0, nReads);
  snpSqrtWts = Float64[];
  tmpvec = similar(yShift);
  tmpvecQRes = Array(Float64, rankQPhi);
  yWork = similar(evalPhiAdj);
  partialSumWConst = Array(Float64, nNullSimPts);
  totalSumWConst = Array(Float64, nNullSimPts);
  subXPhitSV = Array(Float64, size(XPhitNullBasis, 2), rankQPhi);
  pSubXPhitSV = pointer(subXPhitSV);
  offset = 0;
  oneConst = ones(nPerKeep);
  geno = Array(Float64, nPer, SNPSize);
  genoOri = Array(Float64, nPer, SNPSize);
  snpIDred = similar(snpID);
  tmpMat = Array(Float64, windowSize, size(XPhitNullBasis, 2));
  VWorkSqrt = Array(Float64, length(evalPhiAdj), windowSize);
  VWorkSqrt2 = Array(Float64, length(evalPhiAdj), windowSize);
  mafOri = Array(Float64, SNPSize);
  maf = Array(Float64, SNPSize);
  sumIsnan = Array(Float64, 1, SNPSize);
  nChrObs = Array(Float64, 1, SNPSize);
  tmpvecKernel = Array(Float64, windowSize);
  if pvalueComputing == "chi2"
    nPtsChi2 = 300;
    partialSumW = Array(Float64, nPtsChi2);
    totalSumW = Array(Float64, nPtsChi2);
    lambda = Array(Float64, 1, nPtsChi2);
    W = Array(Float64, windowSize, nPtsChi2);
    tmpmat0 = Array(Float64, windowSize, nPtsChi2);
    tmpmat1 = Array(Float64, windowSize, nPtsChi2);
    tmpmat2 = Array(Float64, windowSize, nPtsChi2);
    tmpmat3 = Array(Float64, windowSize, nPtsChi2);
    tmpmat4 = Array(Float64, windowSize, nPtsChi2);
    tmpmat5 = Array(Float64, windowSize, nPtsChi2);
    denomvec = Array(Float64, nPtsChi2);
    d1f = Array(Float64, nPtsChi2);
    d2f = Array(Float64, nPtsChi2);
    if test == "eScore"
      simnull = Array(Float64, nPtsChi2);
    else
      simnull = Float64[];
    end
  else
    nPtsChi2 = 300;
    partialSumW = Array(Float64, nNullSimPts);
    totalSumW = Array(Float64, nNullSimPts);
    lambda = Array(Float64, 1, nNullSimPts);
    W = Array(Float64, windowSize, nNullSimPts);
    tmpmat0 = Array(Float64, windowSize, nNullSimPts);
    tmpmat1 = Array(Float64, windowSize, nNullSimPts);
    tmpmat2 = Array(Float64, windowSize, nNullSimPts);
    tmpmat3 = Array(Float64, windowSize, nNullSimPts);
    tmpmat4 = Array(Float64, windowSize, nNullSimPts);
    tmpmat5 = Array(Float64, windowSize, nNullSimPts);
    denomvec = Array(Float64, nNullSimPts);
    d1f = Array(Float64, nNullSimPts);
    d2f = Array(Float64, nNullSimPts);
    if test == "eScore"
      simnull = Array(Float64, nNullSimPts);
    else
      simnull = Float64[];
    end
  end

  results = Array(Any, nReads * iceil(SNPSize / windowSize), 3);
  idx = 0;

  for idxRead = 1 : nReads

    # read the binary data
    if idxRead == nReads
      curSNPSize = nSNP - (idxRead - 1) * SNPSize;
    else
      curSNPSize = SNPSize;
    end

    readgeno!(genoOri, curSNPSize, nPer, SNPSize, idxRead, bin2geno, rawdata,
              offset, flagAnnotate);

    tmpMatIsnan = isnan(genoOri);
    sum!(sumIsnan, tmpMatIsnan);
    for i = 1 : curSNPSize
      nChrObs[1, i] = 2 * (nPer - sumIsnan[1, i]);
    end

    # minor allele frequencies and get rid of mono-allelic SNPs
    counter = 0;
    offset = (idxRead - 1) * SNPSize;
    for i = 1 : curSNPSize
      mafOri[i] = 0.0;
      for j = 1 : nPer
        if !tmpMatIsnan[j, i]
          mafOri[i] += genoOri[j, i];
        end
      end
      mafOri[i] = mafOri[i] / nChrObs[1, i];
      if mafOri[i] != 0
        counter += 1;
        maf[counter] = mafOri[i];
        snpIDred[counter] = snpID[i + offset];
        pGenoOri = pointer(genoOri) + (i - 1) * nPer * sizeof(Float64);
        pGeno = pointer(geno) + (counter - 1) * nPer * sizeof(Float64);
        BLAS.blascopy!(nPer, pGenoOri, 1, pGeno, 1);
      end
    end
    nSNPredVec[idxRead] = curSNPSize - counter;

    if counter == 0
      #println("No valid SNPs at the group which starting from row ",
      #        offset + 1);
      continue;
    end

    # impute missing genotype by expected minor allele counts
    for j = 1 : counter
      for i = 1 : nPer
        if isnan(geno[i, j])
          geno[i, j] = 2 * maf[j];
        end
      end
    end

    # boundary of windows
    if idxRead == nReads
      nSNPUpdate = nSNP - (idxRead - 1) * SNPSize - nSNPredVec[idxRead];
    else
      nSNPUpdate = SNPSize - nSNPredVec[idxRead];
    end
    gLB = [1:windowSize:nSNPUpdate];
    gUB = [gLB[2:end]-1; nSNPUpdate];
    nGrps = length(gLB);
    vc0List = zeros(nGrps);
    vc1List = zeros(nGrps);
    pvalList = zeros(nGrps);

    # looping over windows
    for g = 1 : nGrps

      # get SNP weights
      grpSize = gUB[g] - gLB[g] + 1;
      if isempty(snpWtType)
        # constant weights
        if g > 1 && g < nGrps
          snpSqrtWts = fill!(snpSqrtWts, 1);
        else
          snpSqrtWts = ones(grpSize);
        end
      elseif snpWtType == "beta"
        # weights by beta density evaluated at MAF
        snpSqrtWts = sqrt(Distributions.pdf(Distributions.Beta(1, 25), maf[gLB[g] : gUB[g]]));
      elseif snpWtType == "invvar"
        # weights by inverse variance maf*(1-maf)
        snpSqrtWts = 1 ./ sqrt( maf[gLB[g] : gUB[g]] .*
                               (1 - maf[gLB[g] : gUB[g]]) );
      end

      # retrieve group genotypes
      if kernel == "GRM"
        gMAF = maf[gLB[g] : gUB[g]];
        gSNP = geno[keepIdx, gLB[g] : gUB[g]];
        # center and scale genotype matrix by MAF
        BLAS.ger!(-2.0, oneConst, gMAF, gSNP);
        if g == nGrps && grpSize != windowSize
          tmpvecKernel = tmpvecKernel[1 : grpSize];
        end
        for i = 1 : grpSize
          tmpvecKernel[i] = snpSqrtWts[i] / sqrt(2 * gMAF[i] * (1 - gMAF[i]));
        end
        scale!(gSNP, tmpvecKernel);
        #gSNP = gSNP .* (snpSqrtWts ./ sqrt(2 * gMAF .* (1 - gMAF)))';
        if g == nGrps && grpSize != windowSize
          tmpvecKernel = [tmpvecKernel; Array(Float64, windowSize - grpSize)];
        end
      elseif kernel == "IBS1"
        # Table 3 in notes
        gSNP = hcat(geno[keepIdx, gLB[g] : gUB[g]],
                    2 - geno[keepIdx, gLB[g] : gUB[g]]);
        scale!(gSNP, [snpSqrtWts; snpSqrtWts]);
      elseif kernel == "IBS2"
        # Table 2 in notes
        gSNP = hcat(geno[keepIdx, gLB[g] : gUB[g]] .>= 1,
                    geno[keepIdx, gLB[g] : gUB[g]] .<= 1);
        scale!(gSNP, [snpSqrtWts; snpSqrtWts]);
      elseif kernel == "IBS3"
        # Table 1 in notes
        gSNP = hcat(-2 * (geno[keepIdx, gLB[g] : gUB[g]] .<= 1) + 1,
                    2 * (geno[keepIdx, gLB[g] : gUB[g]] .>= 1) - 1,
                    fill(sqrt(2), nPerKeep, gUB[g] - gLB[g] + 1));
        scale!(gSNP, [snpSqrtWts; snpSqrtWts; snpSqrtWts]);
      end

      # testing

      if kinship == "none"

        if test == "eRLRT"
          (b, vc0List[g], vc1List[g], pvalList[g]) =
            vctest(ynew, [], XtNullBasis * gSNP, WPreSim = WPreSim, tests = test,
                   Vform = "half", nBlockAscent = nBlockAscent, nMMmax = nMMmax,
                   nNullSimPts = nNullSimPts, pvalueComputings = pvalueComputing,
                   nNullSimNewtonIter = nNullSimNewtonIter, tolX = tolX,
                   devices = device, PrePartialSumW = PrePartialSumW,
                   PreTotalSumW = PreTotalSumW, partialSumWConst = partialSumWConst,
                   totalSumWConst = totalSumWConst, windowSize = windowSize,
                   partialSumW = partialSumW, totalSumW = totalSumW,
                   lambda = lambda, W = W, nPreRank = nPreRank,
                   tmpmat0 = tmpmat0, tmpmat1 = tmpmat1, tmpmat2 = tmpmat2,
                   tmpmat3 = tmpmat3, tmpmat4 = tmpmat4, tmpmat5 = tmpmat5,
                   denomvec = denomvec, d1f = d1f, d2f = d2f, offset = offset,
                   nPtsChi2 = nPtsChi2, simnull = simnull);
        else
          (b, vc0List[g], vc1List[g], pvalList[g]) =
            vctest(y, X, gSNP, WPreSim = WPreSim, tests = test,
                   Vform = "half", nBlockAscent = nBlockAscent, nMMmax = nMMmax,
                   nNullSimPts = nNullSimPts, pvalueComputings = pvalueComputing,
                   nNullSimNewtonIter = nNullSimNewtonIter, tolX = tolX,
                   devices = device, PrePartialSumW = PrePartialSumW,
                   PreTotalSumW = PreTotalSumW, partialSumWConst = partialSumWConst,
                   totalSumWConst = totalSumWConst, windowSize = windowSize,
                   partialSumW = partialSumW, totalSumW = totalSumW,
                   lambda = lambda, W = W, nPreRank = nPreRank,
                   tmpmat0 = tmpmat0, tmpmat1 = tmpmat1, tmpmat2 = tmpmat2,
                   tmpmat3 = tmpmat3, tmpmat4 = tmpmat4, tmpmat5 = tmpmat5,
                   denomvec = denomvec, d1f = d1f, d2f = d2f, offset = offset,
                   nPtsChi2 = nPtsChi2, simnull = simnull);
        end

      else

        # obtain some orthonormal vectors of space R^n - [QX,QPHI,S]
        # See Golub and Van Loan (1996) Algorithm 12.4.2 on p602
        if g == nGrps && grpSize != windowSize
          tmpMat = tmpMat[1 : grpSize, :];
        end
        BLAS.gemm!('T', 'N', 1.0, gSNP, XPhitNullBasis, 0.0, tmpMat);
        (UXPhitS, svalXPhitS, VXPhitS) = svd(tmpMat, thin = false);
        if g == nGrps && grpSize != windowSize
          tmpMat = [tmpMat; Array(Float64, windowSize - grpSize,
                                  size(XPhitNullBasis, 2))];
        end
        rankXPhitS =
          countnz(svalXPhitS .>
                  max(size(XPhitNullBasis, 2),
                      size(gSNP, 2) * eps(svalXPhitS[1])));
        pXPhitSV = pointer(VXPhitS) +
          rankXPhitS * size(XPhitNullBasis, 2) * sizeof(Float64);
        BLAS.blascopy!(size(XPhitNullBasis, 2) * rankQPhi,
                       pXPhitSV, 1, pSubXPhitSV, 1);
        BLAS.gemm!('N', 'N', 1.0, XPhitNullBasis, subXPhitSV, 0.0, QRes);

        # working Y and working non-trivial variance component
        #yWork = (PhiAdjsvd[:U]' * (yShift + KPhiAdj * (QRes' * y))) ./ sqrt(evalPhiAdj);
        BLAS.blascopy!(length(yShift), yShift, 1, tmpvec, 1);
        BLAS.gemv!('T', 1.0, QRes, y, 0.0, tmpvecQRes);
        BLAS.gemv!('N', 1.0, KPhiAdj, tmpvecQRes, 1.0, tmpvec);
        #yWork = BLAS.gemv('T', 1.0, PhiAdjsvd[:U], tmpvec) ./ sqrt(evalPhiAdj);
        BLAS.gemv!('T', 1.0, weightedW, tmpvec, 0.0, yWork);
        #for i = 1 : length(yWork)
        #  yWork[i] = yWork[i] * InvSqrtevalPhiAdj[i];
        #end
        #VWorkSqrt = PhiAdjsvd[:U]' * (QPhi' * gSNP) ./ sqrt(evalPhiAdj);
        if g == nGrps && grpSize != windowSize
          VWorkSqrt = VWorkSqrt[:, 1 : grpSize];
          VWorkSqrt2 = VWorkSqrt2[:, 1 : grpSize];
        end
        BLAS.gemm!('T', 'N', 1.0, QPhi, gSNP, 0.0, VWorkSqrt2);
        BLAS.gemm!('T', 'N', 1.0, weightedW, VWorkSqrt2, 0.0, VWorkSqrt);
        #scale!(InvSqrtevalPhiAdj, VWorkSqrt);

        # obtain LRT and p-value
        (b, vc0List[g], vc1List[g], pvalList[g]) =
          vctest(yWork, [], VWorkSqrt, WPreSim = WPreSim, tests = test,
                 Vform = "half", nBlockAscent = nBlockAscent, nMMmax = nMMmax,
                 nNullSimPts = nNullSimPts, pvalueComputings = pvalueComputing,
                 nNullSimNewtonIter = nNullSimNewtonIter, tolX = tolX,
                 devices = device, PrePartialSumW = PrePartialSumW,
                 PreTotalSumW = PreTotalSumW, partialSumWConst = partialSumWConst,
                 totalSumWConst = totalSumWConst, windowSize = windowSize,
                 partialSumW = partialSumW, totalSumW = totalSumW,
                 lambda = lambda, W = W, nPreRank = nPreRank,
                 tmpmat0 = tmpmat0, tmpmat1 = tmpmat1, tmpmat2 = tmpmat2,
                 tmpmat3 = tmpmat3, tmpmat4 = tmpmat4, tmpmat5 = tmpmat5,
                 denomvec = denomvec, d1f = d1f, d2f = d2f, offset = offset,
                 nPtsChi2 = nPtsChi2, simnull = simnull);

        # recover VWorkSqrt
        if g == nGrps && grpSize != windowSize
          VWorkSqrt = hcat(VWorkSqrt,
                           Array(Float64, length(evalPhiAdj),
                                 windowSize - grpSize));
          VWorkSqrt2 = hcat(VWorkSqrt2,
                            Array(Float64, length(evalPhiAdj),
                                  windowSize - grpSize));
        end

      end

      idx += 1;
      results[idx, :] = hcat(snpIDred[gLB[g]], snpIDred[gUB[g]], pvalList[g]);

    end

  end

  return results[1:idx, :];

end


text:
##
using Revise
using GLMakie
using Pkg
Pkg.activate(".")
using ProteinStructureReconstruction
using Images
using HDF5

##
module ZStackIO

	import Images
    using RecursiveArrayTools: VectorOfArray
	export find_zstack_folder_by_number, read_STORM_folder, BoundingBox, process_STORM_data
	unravel_array(vec_arr) = convert(Array, VectorOfArray(vec_arr))

	function find_zstack_folder_by_number(dir::AbstractString = "./"; nb::Integer)
		name_regex = Regex("^Stack $(nb)")
		for path in readdir(dir)
			if occursin(name_regex, path) && isdir(dir)
				return joinpath(dir, path)
			end
		end
		error("The z stack with the given number does not exist!")
	end

	function read_STORM_folder(dir::AbstractString = "./")
		root, folders, files = first(walkdir(dir))
        print(folders)
		names = map(target_name_from_foldername, folders)
		image_stacks = map(x -> read_subfolder(root, x), folders)

		if "overlay.tif" in files
			overlay = Images.load(joinpath(root, "overlay.tif"))
		else
			error("overlay image does not exist!")
		end

		return overlay, names, image_stacks
	end

	function read_subfolder(root, subdir::AbstractString)
		files = readdir(joinpath(root, subdir), join = true)
		sort!(files, by = slice_nb_from_filename)
		images = map(Images.load, files)
		return images
	end

	function slice_nb_from_filename(filename::AbstractString)
		m = match(r"z(?<slice_nb>\d+) \d+.tif$", filename)
		return parse(Int, m[:slice_nb])
	end

	function target_name_from_foldername(foldername::AbstractString)
		re_test = match(r"^(?<exp_name>[A-Z, a-z]*)_[\w, \s]*", foldername)
		return re_test[:exp_name]
	end

	struct BoundingBox{T<:Integer}
		lower_left::Tuple{T, T}
		upper_right::Tuple{T, T}
	end
	BoundingBox(diagonal_vertices::AbstractVector) = BoundingBox(diagonal_vertices...)

	Base.size(bb::BoundingBox) = bb.upper_right .- bb.lower_left .+ 1
	Base.size(bb::BoundingBox, dim::Integer) = Base.size(bb)[dim]

	function Base.to_indices(A::AbstractMatrix, inds, I::Tuple{BoundingBox})
		x_min, y_min = I[1].lower_left
		x_max, y_max = I[1].upper_right
		return Base.to_indices(A, inds, (x_min:x_max, y_min:y_max))
	end

	function find_scale_bar_length(overlay_grayed::AbstractMatrix)
		iswhite = isequal.(overlay_grayed, 1.0)
		white_components = Images.label_components(iswhite)
		bboxes = BoundingBox.(Images.component_boxes(white_components))

		# first bbox is for the null label (label = 0)
		bbox_sizes = size.(bboxes[2:end])

		# scale bar has smallest x (vertical) size
		scale_bar_index = argmin(first.(bbox_sizes))

		scale_bar_length = bbox_sizes[scale_bar_index][2]
		return scale_bar_length
	end

	function find_image_bbox(overlay_grayed::AbstractMatrix)
	# know apriori that the non image features (scale, colorbar) can be cropped out with this choice of indexing
	overlay_cropped = overlay_grayed[:, 1:1500]

	# then find the bounding box for the relevant image region
	isnonzero = Int.(overlay_cropped .> 0.0)
	bboxes = BoundingBox(Images.component_boxes(isnonzero)[end])
	end

	function process_STORM_data(overlay::AbstractMatrix, names, image_stacks)
		overlay_grayed = Images.Gray.(overlay)

		scale_bar_len = find_scale_bar_length(overlay_grayed)
		bbox = find_image_bbox(overlay_grayed)

		cropped_stacks = unravel_array(map(x -> crop_image_stack(x, bbox), image_stacks))
        return cropped_stacks, scale_bar_len, names
	end

	function crop_image_stack(image_stack::AbstractVector{T}, bbox::BoundingBox) where {T <: AbstractMatrix}
		cropped_stack = map(img -> getindex(img, bbox), image_stack)
		return unravel_array(cropped_stack)
	end
end
##
stack, scale_bar_len, names = ZStackIO.find_zstack_folder_by_number("./Data/STORM/", nb = 2) |> ZStackIO.read_STORM_folder |> args -> ZStackIO.process_STORM_data(args...);

##
scale = Float32.([5.0*10^-6/scale_bar_len, 5.0*10^-6/scale_bar_len, 2.0*10^-7])
##
rawview.(stack)
##
filename = "./Data/STORM_stack2.h5"
fid = h5open(filename, "w")


text: # Parametric active contours


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule QDLDL_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("QDLDL")
JLLWrappers.@generate_main_file("QDLDL", UUID("9ae34a36-e1b6-54e9-a33d-8bba9913f854"))
end  # module QDLDL_jll


text:
import Base: ==, convert

import ACEbase: read_dict, write_dict

export Dat, eval_obs, observations, vec_obs, devec_obs,
       set_weights!

# -----------------------------------------------------------------
# registry of observation types

# TODO: not sure yet that we really need this - leave it for now but revisit.

const _obstypes = Dict{String, Any}()

function register_obstype!(key::String, ObsType; force = false)
   if !forces && haskey(_obstypes, key)
      error("""This observation type key is already registered. To overwrite it,
               use `register_obstype!(...; force=true)`.""")
   end
   _obstype[key] = ObsType
   return nothing
end

function empty_obstypes!()
   empty!(_obstypes)
   return nothing
end

# -----------------------------------------------------------------
# basic data structure to manage training data

"""
`Dat`: store one configuration (input, e.g., structure, state, ...)  that can
have multiple observations attached to it. Fields:
- `config::Any` : the structure
- `configtype::String` : Each `dat::Dat` belongs to a group identified by a string `dat.configtype` to allow filtering, and grouping.
- `obs::Vector{Any}`  : list of observations
- `meta::Dict{String, Any}` : any additional meta information that we may want to attach to this data point; this needs to be raw json.
"""
mutable struct Dat
   config                             # configuration
   configtype::String                 # group identifier
   obs::Vector{Any}                  # list of observations
   meta::Dict{String, Any}           # anything else ...
end

Dat(config, configtype="", obs=Any[]) =
      Dat(config, configtype, obs, Dict{String, Any}())

==(d1::Dat, d2::Dat) = (
      (d1.config == d2.config) &&
      (d1.configtype == d2.configtype) &&
      (d1.obs == d2.obs) &&
      (d1.meta == d2.meta)
   )

write_dict(d::Dat) =
   Dict("__id__" => "ACEfit_Dat",
         "config" => write_dict(d.config),
         "configtype" => d.configtype,
         "obs" => write_dict.(d.obs),
         "meta" => d.meta)

Dat(D::Dict) = Dat( read_dict(D["config"]),
                    D["configtype"],
                    Vector{Any}(D["obs"]),
                    Dict{String, Any}(D["meta"]) )

read_dict(::Val{:ACEfit_Dat}, D::Dict) = Dat(D)

observations(d::Dat) = d.obs

cost(d::Dat) = 1

# -----------------------------------------------------------------
# Abstract Observations interface
# at this point it is not clear that we need an abstract Observations
# type hierarchy. Better to try and live without for now.

# prototypes

"""
convert some real data, in some generic format, into a vector to be stored
in a `Dat` or Lsq system. E.g.,
```julia
V = virial(...)::Matrix{Float64}
obsV = ObsVirial(V)
vec_obs(obsV::ObsVirial) = obsV.V[ [1,2,3,5,6,9] ]  # NB: V is symmetric
```
"""
function vec_obs end

"""
convert a Vector{T} to some real data, e.g.,
```julia
x::Vector{Float64}
devec_obs(::Type{ObsVirial}, x) = [ x[1] x[2] x[3];
                                    x[2] x[4] x[5];
                                    x[3] x[5] x[6] ]
```
"""
function devec_obs end

"""
Evaluate a specific observation type: Given an observation `obs`,
a model `model` and  a configuration `cfg = dat.config`, the call
```julia
eval_obs(obstype, model, cfg)
```
must return the corresponding observation. For example, if
`obs::ObsPotentialEnergy` and `cfg = at::Atoms`, and `model` is an interatomic
potential, then
```julia
eval_obs(::Type{TOBS}, model, cfg) where {TOBS <: ObsPotentialEnergy} =
      TOBS( energy(model, cfg) )
```
"""
function eval_obs end

# the next prototype is again experimental- renaming this to basis_obs
# instead of dispatching means that we don't have to worry about an abstract
# type hierarchy, we should discuss whether it obscures or clarifies

"""
`function basis_obs end` - evaluate an observation on a basis and a config;
used for llsq assembly,
```julia
basis_obs(obstype, basis, cfg)
```
For examle, for a potential energy the implementatin might look like this:
```julia
basis_obs(::Type{TOBS}, basis, at) where {TOBS <: ObsPotentialEnergy} =
   TOBS.( energy(basis, at) )
```
Note how the `TOBS` operation is now broadcast since `energy(basis, at)`
returns a vector of observations (one for each  basis functions) whereas
for a *model* we would use `eval_obs`.
"""
function basis_obs end

# a first attempt at an interface for weights...
# this feels like we should start dispatching and
# requiring o::AbstractObservation

function set_weight!(o, w)
   o.weight = w
   return o
end

function get_weight(o)
   if !hasproperty(o, :weight)
      return 1
   end
   return o.weight
end

function _compute_weight(OBSTYPE, cfg, weights::Dict, weighthooks::Dict)
   if haskey(weights, OBSTYPE)
      w = weights[OBSTYPE]
   elseif haskey(weights, "default")
      w = weights["default"]
   else
      w = 1.0
   end
   if haskey(weighthooks, OBSTYPE)
      w = weighthooks[OBSTYPE](w, cfg)
   elseif haskey(weighthooks, "default")
      w = weighthooks["default"](w, cfg)
   end
end

function set_weights!(dat::Dat, weights::Dict, weighthooks::Dict)
   key = dat.configtype
   if !haskey(weights, key)
      key = "default"
   end

   key_h = dat.configtype
   if !haskey(weighthooks, key_h)
      key_h = "default"
   end

   for o in dat.obs
      w = _compute_weight(typeof(o), dat.config,
                          weights[key], weighthooks[key_h])
      set_weight!(o, w)
   end
   return dat
end

_id_weighhooks = Dict("default" => Dict("default" => (w, cfg) -> w))

"""
TODO: needs documentation
"""
function set_weights!(data::AbstractVector{<: Dat},
                      weights::Dict,
                      weighthooks::Dict = _id_weighhooks)
   for dat in data
      set_weights!(dat, weights, weighthooks)
   end
   return data
end


text:
using Zygote: @adjoint

"""
    dispatch_diff_apply!(reg::ArrayReg, circuit::AbstractBlock, diff_params::Vector) -> ArrayReg

Dispatch parameters to differentiable nodes, and evaluate the output.
"""
function dispatch_diff_apply!(reg::ArrayReg, circuit::AbstractBlock, diff_params::Vector)
    dispatch_to_diff!(circuit, diff_params)
    apply!(reg, circuit)
end

@adjoint function dispatch_diff_apply!(reg::ArrayReg{B,T,AT}, circuit::AbstractBlock, diff_params::Vector) where {B,T,AT}
    out = dispatch_diff_apply!(reg, circuit, diff_params)
    saved_out = copy(out)
    out, function (adjout)
        dpsi, dparams = classical_autodiff!(circuit, saved_out, ArrayReg{B,T,AT}(adjout.x.state/2); output_eltype=Any)
        return (dpsi, nothing, dparams)
    end
end

@adjoint Base.copy(reg::ArrayReg) = Base.copy(reg), adjreg->(adjreg,)


text:
export ReinforcePolicy

using Flux: softmax
using LinearAlgebra: dot
using StatsBase

"""
    ReinforcePolicy(approximator, α, γ)
This is the very basic reinforce algorithm.
TODO: implement some other variants.
"""
Base.@kwdef struct ReinforcePolicy{A<:AbstractApproximator} <: AbstractPolicy
    approximator::A
    α::Float64
    γ::Float64
end

(π::ReinforcePolicy)(obs) =
    obs |> get_state |> π.approximator |> softmax |> x -> Weights(x, 1.0) |> sample

RLBase.get_prob(π::ReinforcePolicy, s) = s |> π.approximator |> softmax
RLBase.get_prob(π::ReinforcePolicy, s, a) = get_prob(π, s)[a]

# TODO: handle neural network q approximator
function RLBase.update!(π::ReinforcePolicy{<:LinearQApproximator}, t::AbstractTrajectory)
    if length(t) > 0 && get_trace(t, :terminal)[end]
        states, actions, rewards = get_trace(t, :state), get_trace(t, :action), get_trace(t, :reward)
        Q, α, γ = π.approximator, π.α, π.γ
        gains = discount_rewards(rewards, γ)
        γₜ = 1.0

        for (i, (s, a, g)) in enumerate(zip(states, actions, gains))
            # !!! we will multiply `Q.feature_func(s, a)` in the `LinearQApproximator` again!
            update!(
                Q,
                (s, a) => α * γₜ * g *
                          (Q.feature_func(s, a) .-
                           sum(x -> get_prob(π, s, x) .* Q.feature_func(s, x), Q.actions)),
            )
            γₜ *= γ
        end
    end
end


text:
# Copyright (c) 2013: Clp.jl contributors
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

using Test

@testset "MathOptInterface" begin
    include("MOI_wrapper.jl")
end


text:
const tutorialpath = normpath(joinpath(@__DIR__, "..", "examples", "tutorial"))
const tworegionpath = joinpath(tutorialpath, "02-two-region-model")
const oneregionpath = joinpath(tutorialpath, "01-one-region-model")

function run_oneregion()
    include(joinpath(oneregionpath, "one-region-model.jl"))
end

function run_tworegion()
    include(joinpath(tworegionpath, "main.jl"))
end


text:
using AdFem
using Test
using PyPlot
using PyCall
using LinearAlgebra

include("MFEM/MUtils.jl")
# sess = Session();

# include("invkernel.jl")
# include("mfem.jl")