text:
using Dates
using Statistics
using LinearAlgebra
import TOML

const SENTINEL_WAVELENGTH = 5.5465763  # cm
const PHASE_TO_CM = SENTINEL_WAVELENGTH / ( 4 * π )
const P2MM = 365 * 10 * PHASE_TO_CM   # mm / year

const DateOrNone = Union{Date,Nothing}

_stack_size_mb(filename, dset) =
    prod(size(filename, dset)) * sizeof(eltype(filename, dset)) / 1e6
_strnothing(aa) = isnothing(aa) ? "nothing" : aa   # Since we can't string(nothing)

function run_inversion(config_file::Dict{AbstractString,Any})
    conf_dict = TOML.parsefile(config_file)
    # Convert to symbols so it works to pass as kwargs
    symbol_dict = Dict(Symbol(k) => v for (k, v) in conf_dict)
    outfile, outdset = run_inversion(symbol_dict...)
    h5writeattr(outfile, outdset, conf_dict)
end

# overwrites in the dset fail now... is that fine? force user to delete?
function run_inversion(;
    unw_stack_file::String=UNW_FILENAME,
    input_dset::String=STACK_FLAT_SHIFTED_DSET,
    outfile::String="",
    # outgroup::String="velos",
    stack_average::Bool=false,
    constant_velocity::Bool=true,
    ignore_geo_file::String="",
    max_temporal_baseline::Int=500,
    prune_outliers=true,
    sigma=4,
    gap=1,
    min_date::DateOrNone=nothing,
    max_date::DateOrNone=nothing,
    alpha::Real=0.0,
    L1::Bool=false,
    reference_station=nothing, # Note: these only work for stackavg...
    ref_row=nothing,
    ref_col=nothing,
    use_distributed=true,
)
    if isempty(outfile)
        outfile = _default_outfile()
    end

    # averaging or linear means output will is 3D array (not just map of velocities)
    is_3d = !(stack_average || constant_velocity)
    outdset = is_3d ? "stack" : "velos"

    # Check if this file/dset already exists
    isfile(outfile) &&
    outdset in names(outfile) &&
    error("$outdset already exists in $outfile")

    # the valid igram indices is out of all possible layers in the stack
    geolist, intlist, valid_igram_indices = load_geolist_intlist(
        unw_stack_file,
        ignore_geo_file,
        max_temporal_baseline,
        min_date=min_date,
        max_date=max_date,
    )

    println("$outdset geolist range: : $(extrema(geolist))")
    println("Using $input_dset as input")
    # Now: can we load the input stack into memory? or do we need distributed?
    stack_file_size = _stack_size_mb(unw_stack_file, input_dset)
    can_fit_mem = (stack_file_size * 8) < getmemavail()  # Rough padding for total memory check

    # Note: as of version 1.3.0rc2, my threaded version on preloaded matrix is slower than
    # distributed reading/writing... not sure why but oh well
    # TODO: figure if it's worth having non-process version
    println("Stack size: $stack_file_size, avail RAM: $(getmemavail()), fitting in ram: $can_fit_mem")
    println("Using Distributed to solve: $use_distributed")

    if stack_average
        # # Dont need shift for avg
        # input_dset = STACK_FLAT_DSET
        println("Averaging stack for solution")
        velo_file_out = run_stackavg(
            unw_stack_file,
            input_dset,
            outfile,
            outdset,
            geolist,
            intlist;
            reference_station=reference_station,
            ref_row=ref_row,
            ref_col=ref_col,
        )
    else
        println("Running SBAS for solution")

        velo_file_out = run_sbas(
            unw_stack_file,
            input_dset,
            outfile,
            outdset,
            geolist,
            intlist,
            valid_igram_indices,
            constant_velocity,
            alpha,
            L1,
            prune_outliers,
            sigma,
        )
        # if use_distributed
        # else
        #     # If we want to read the whole stack in at once:
        #     @time unw_stack = h5read(unw_stack_file, input_dset)[:, :, valid_igram_indices]
        #     velo_file_out = run_sbas(unw_stack, outfile, outdset, geolist, intlist,
        #                              constant_velocity, alpha, L1, prune)
        # end

    end

    ####################

    dem_rsc = Sario.load_dem_from_h5(unw_stack_file)
    Sario.save_dem_to_h5(outfile, dem_rsc, overwrite=true)
    Sario.save_geolist_to_h5(outfile, outdset, geolist, overwrite=true)
    save_reference(outfile, unw_stack_file, outdset, input_dset)

    return outfile, outdset
end

function load_geolist_intlist(
    unw_stack_file,
    ignore_geo_file,
    max_temporal_baseline;
    min_date::DateOrNone=nothing,
    max_date::DateOrNone=nothing,
)
    @show unw_stack_file
    geolist = Sario.load_geolist_from_h5(unw_stack_file)
    intlist = Sario.load_intlist_from_h5(unw_stack_file)

    # If we are ignoreing some indices, remove them from for all pixels
    geo_idxs, igram_idxs = find_valid_indices(
        geolist,
        intlist,
        min_date,
        max_date,
        ignore_geo_file,
        max_temporal_baseline,
    )
    return geolist[geo_idxs], intlist[igram_idxs], igram_idxs
end

"""Cut down the full list of interferograms and geo dates

- Cuts date ranges (e.g. piecewise linear solution) with  `min_date` and `max_date`
- Can ignore whole dates by reading `ignore_geo_file`
- Prunes long time baseline interferograms with `max_temporal_baseline`
"""
function find_valid_indices(
    geo_date_list::Array{Date,1},
    igram_date_list::Array{Igram,1},
    min_date::DateOrNone=nothing,
    max_date=DateOrNone = nothing,
    ignore_geo_file::String="",
    max_temporal_baseline::Int=500,
)

    ig1 = length(igram_date_list)  # For logging purposes, what do we start with
    if isempty(ignore_geo_file)
        println("Not ignoring any .geo dates")
        ignore_geos = []
    else
        ignore_geos = sort(Sario.find_geos(filename=ignore_geo_file, parse=true))
        println("Ignoring the following .geo dates:")
        println(ignore_geos)
    end

    # TODO: do I want to be able to pass an array of dates to ignore?

    # First filter by remove igrams with either date in `ignore_geo_file`
    valid_geos = [g for g in geo_date_list if !(g in ignore_geos)]
    valid_igrams = [ig for ig in igram_date_list if !(ig in ignore_geos)]
    println("Ignoring $(ig1 - length(valid_igrams)) igrams listed in $ignore_geo_file")

    # Remove geos and igrams outside of min/max range
    if !isnothing(min_date)
        println("Keeping data after min_date: $min_date")
        valid_geos = [g for g in valid_geos if g > min_date]
        valid_igrams = [ig for ig in valid_igrams if (ig[1] > min_date && ig[2] > min_date)]
    end
    if !isnothing(max_date)
        println("Keeping data only before max_date: $max_date")
        valid_geos = [g for g in valid_geos if g < max_date]
        valid_igrams = [ig for ig in valid_igrams if (ig[1] < max_date && ig[2] < max_date)]
    end

    # This is just for logging purposes:
    too_long_igrams =
        filter(ig -> temporal_baseline(ig) > max_temporal_baseline, valid_igrams)
    println("Ignoring $(length(too_long_igrams)) igrams with longer baseline than $max_temporal_baseline days")

    ### Remove long time baseline igrams ###
    valid_igrams =
        filter(ig -> temporal_baseline(ig) <= max_temporal_baseline, valid_igrams)

    ### Collect remaining geo dates and igrams
    valid_geo_indices = indexin(valid_geos, geo_date_list)
    valid_igram_indices = indexin(valid_igrams, igram_date_list)

    println("Ignoring $(ig1 - length(valid_igrams)) igrams total")

    return valid_geo_indices, valid_igram_indices
end

"""Get min,max pairs to split up interval of dates
gap=1 means form adjacent pairs
gap=2 skips a date in between, makes overlapping pairs to get larger, redundant estimates
"""
function _get_pairs(dates, gap=1)
    isempty(dates) && return (nothing, nothing)
    length(dates) < gap && error("Need at least $gap dates if gap=$gap")
    dates_pad = [nothing; dates; nothing]
    return collect(zip(dates_pad[1:end - gap], dates_pad[(gap + 1):end]))
end

function save_reference(
    h5file,
    unw_stack_file,
    dset_name,
    stack_flat_shifted_dset,
    overwrite=true,
)
    # Delete if exists
    if overwrite
        h5open(h5file) do f
            if exists(attrs(f[dset_name]), REFERENCE_ATTR)
                a_delete(f[dset_name], REFERENCE_ATTR)
            end
            if exists(attrs(f[dset_name]), REFERENCE_STATION_ATTR)
                a_delete(f[dset_name], REFERENCE_STATION_ATTR)
            end
        end
    end
    # Read ref. unfo from unw file, save what was used to deformation result file
    reference = get(h5readattr(unw_stack_file, stack_flat_shifted_dset), REFERENCE_ATTR, "")
    h5writeattr(h5file, dset_name, Dict(REFERENCE_ATTR => reference))
    reference_station =
        get(h5readattr(unw_stack_file, stack_flat_shifted_dset), REFERENCE_STATION_ATTR, "")
    h5writeattr(h5file, dset_name, Dict(REFERENCE_STATION_ATTR => reference_station))
end


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges
const MOIBC = MathOptInterface.Bridges.Constraint

include("../utilities.jl")

mock = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()))
config = MOIT.TestConfig()

@testset "ZeroOne" begin
    bridged_mock = MOIBC.ZeroOne{Float64}(mock)

    bridge_type = MOIBC.ZeroOneBridge{Float64}
    @test MOI.supports_constraint(bridge_type, MOI.SingleVariable, MOI.ZeroOne)
    @test MOIBC.concrete_bridge_type(bridge_type,
                                     MOI.SingleVariable,
                                     MOI.ZeroOne) == bridge_type

    @test MOI.supports(bridged_mock, MOI.ConstraintPrimalStart(), bridge_type)

    MOIT.basic_constraint_tests(bridged_mock, config, include = [(MOI.SingleVariable, MOI.ZeroOne)])

    MOIU.set_mock_optimize!(mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1, 0, 0, 1, 1]))
    MOIT.knapsacktest(bridged_mock,config)

    ci = first(MOI.get(bridged_mock,
                       MOI.ListOfConstraintIndices{MOI.SingleVariable, MOI.ZeroOne}()))

    @test MOI.get(bridged_mock, MOI.ConstraintPrimal(), ci) == 1

    test_delete_bridge(bridged_mock, ci, 5,
                       ((MOI.SingleVariable, MOI.Integer, 0),
                        (MOI.SingleVariable, MOI.Interval{Float64}, 0),),
                       num_bridged = 5)

    MOIU.set_mock_optimize!(mock,
        (mock::MOIU.MockOptimizer) -> begin
            MOI.set(mock, MOI.ObjectiveBound(), 20.0)
            MOIU.mock_optimize!(mock, [4, 5, 1])
        end)
    MOIT.int1test(bridged_mock,config)

    MOIU.set_mock_optimize!(mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0; zeros(10)]))
    MOIT.int3test(bridged_mock,config)

    ci = first(MOI.get(bridged_mock,
                       MOI.ListOfConstraintIndices{MOI.SingleVariable, MOI.ZeroOne}()))

    @testset "$attr" for attr in [MOI.ConstraintPrimalStart(),]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        value = 1.0
        MOI.set(bridged_mock, attr, ci, value)
        @test MOI.get(bridged_mock, attr, ci) ≈ value
    end

    s = """
    variables: x, y
    cy: y == 1.0
    cx: x in ZeroOne()
    minobjective: x
    """
    model = MOIU.Model{Float64}()
    MOIU.loadfromstring!(model, s)
    sb = """
    variables: x, y
    cy: y == 1.0
    intg: x in Integer()
    intv: x in Interval(0.0,1.0)
    minobjective: x
    """
    modelb = MOIU.Model{Float64}()
    MOIU.loadfromstring!(modelb, sb)

    MOI.empty!(bridged_mock)
    @test MOI.is_empty(bridged_mock)
    MOIU.loadfromstring!(bridged_mock, s)
    MOIU.test_models_equal(bridged_mock, model, ["x", "y"], ["cy", "cx"])

    # setting names on mock
    ci = first(MOI.get(mock, MOI.ListOfConstraintIndices{
        MOI.SingleVariable, MOI.Integer}()))
    MOI.set(mock, MOI.ConstraintName(), ci, "intg")
    ci = first(MOI.get(mock, MOI.ListOfConstraintIndices{
        MOI.SingleVariable, MOI.Interval{Float64}}()))
    MOI.set(mock, MOI.ConstraintName(), ci, "intv")

    MOIU.test_models_equal(mock, modelb, ["x", "y"], ["cy", "intg", "intv"])

end

ext: jl


text:
using CSV, DataFrames, StatsPlots, StatsBase

categories = String[
    "food", "animal", "clothing", "tool", "sports equipment", "vegetable", "vehicle", "musical instrument",
    "fruit", "body part", "dessert", "toy", "container", "part of car", "weapon", "bird", "furniture",
    "kitchen tool", "office supply", "clothing accessory", "kitchen appliance", "plant", "insect",
    "home décor", "medical equipment", "electronic device", "drink"
]

function parse_categories(datafile::String)
    df = DataFrame(CSV.File(datafile, header = true))

    # check that all categories are defined
    for i in skipmissing(df.top_down)
        if i ∉ categories
            error("There is no category called \"$(i)\".")
        end
    end

    frequencies = countmap(skipmissing(df.category))

    theme(:solarized)

    plt = bar(
        categories,
        Int[get(frequencies, c, 0) for c in categories],
        legend = false,
        xticks = (0.5:length(categories), categories),
        xrotation = 45,
        fontfamily = font("Times"),
        xlabel = "Object Category",
        ylabel = "Frequency",
        title = "Frequency of Object Categories in chosen ~200 Distractors"
    )

    # plt = bar(
    #     collect(keys(frequencies)),
    #     collect(values(frequencies));
    #     legend = false,
    #     xticks = (0.5:length(keys(frequencies)), collect(keys(frequencies))),
    #     xrotation = 45,
    #     fontfamily = font("Times"),
    #     xlabel = "Category",
    #     ylabel = "Frequency",
    #     title = "Frequency of Categories in chosen 200 Distractors"
    # )

    savefig(plt, "/Users/jakeireland/Desktop/categories.pdf")
end

parse_categories("/Users/jakeireland/projects/GorillaSCExperiment/data/distractor_categories.csv")

max_stars_repo_name: jakewilliami/GorillaSCExperiment

max_issues_repo_name: jakewilliami/GorillaSCExperiment

max_forks_repo_name: jakewilliami/GorillaSCExperiment


text:
begin
    import Pkg
    Pkg.activate(".")

    using MCDTS
    using DelayEmbeddings
    using DynamicalSystemsBase
    using DelimitedFiles
    using ChaosTools
    using Random

    using PyPlot
    pygui(true)

    lo = Systems.lorenz()
    tr = trajectory(lo, 500; dt = 0.01, Ttr = 10)

    # Lyapunov exponent and time
    λ = ChaosTools.lyapunov(lo, 100000, dt=0.01; Ttr=1000)
    lyap_time = Int(floor((1/λ) / 0.01))

    x = tr[:,1]
    y = tr[:,2]
    w1 = DelayEmbeddings.estimate_delay(tr[:,1],"mi_min")
    w2 = DelayEmbeddings.estimate_delay(tr[:,2],"mi_min")
    w3 = DelayEmbeddings.estimate_delay(tr[:,3],"mi_min")
    theiler = maximum([w1,w2,w3])
end

## Predictions based on embedding
begin
    T_steps = 8*lyap_time
    x1 = tr[1:end-T_steps,1]
    x2 = tr[end-T_steps+1:end,1]
    y1 = tr[1:end-T_steps,2]
    y2 = tr[end-T_steps+1:end,2]

    dmax = 10
    # cao
    𝒟, τ_tde1, _ = optimal_traditional_de(x1, "afnn"; dmax = dmax, w = w1)
    optimal_d_tde1 = size(𝒟, 2)
    τ_cao = [(i-1)*τ_tde1 for i = 1:optimal_d_tde1]
    Y_cao = genembed(x1, τ_cao.*(-1))

    # kennel's method
    𝒟, τ_tde2, _ = optimal_traditional_de(x1, "fnn"; dmax = dmax, w = w1)
    optimal_d_tde2 = size(𝒟, 2)
    τ_kennel = [(i-1)*τ_tde2 for i = 1:optimal_d_tde2]
    Y_kennel = genembed(x1, τ_kennel.*(-1))

    # hegger's method
    𝒟, τ_tde3, _ = optimal_traditional_de(x1, "ifnn"; dmax = dmax, w = w1)
    optimal_d_tde3 = size(𝒟, 2)
    τ_hegger = [(i-1)*τ_tde3 for i = 1:optimal_d_tde3]
    Y_hegger = genembed(x1, τ_hegger.*(-1))
end

# mcdts
begin
    using MCDTS
    Random.seed!(1234)
    taus = 0:150
    trials2 = 30
    trials = 10
    Tw_in = 50 # prediction horizon in
    max_depth = 10
    Tw_out = 150 # prediction horizon out
    KNN = 5 # nearest neighbors for pred method
    K = KNN
    samplesize = 0.2
    PredMeth = MCDTS.local_model("zeroth", KNN, Tw_out, Tw_in, trials2)
    PredLoss = MCDTS.PredictionLoss(4)
    PredType = MCDTS.MCDTSpredictionType(PredLoss, PredMeth)
    optmodel = MCDTS.MCDTSOptimGoal(MCDTS.Prediction_error(PredType,0,samplesize), MCDTS.Continuity_function(13,samplesize))
end
@time tree = mcdts_embedding(Dataset(x1), optmodel, w1, taus, trials; max_depth = max_depth, choose_func = (L->(MCDTS.softmaxL(L,β=1.))), verbose=true)
best_node = MCDTS.best_embedding(tree)
τ_mcdts = best_node.τs
# τ_mcdts = [0, 18, 123, 82, 66, 114, 29, 73]
Y_mcdts = genembed(x1, τ_mcdts.*(-1))

begin
    using MCDTS
    Random.seed!(1234)
    taus = 0:150
    trials2 = 30
    trials = 10
    max_depth = 10
    Tw = 100 # prediction horizon out
    KNN = 5 # nearest neighbors for pred method
    K = KNN
    samplesize = 0.2
    PredMeth = MCDTS.local_model("linear", KNN, Tw, 1, trials2)
    PredLoss = MCDTS.PredictionLoss(4)
    PredType = MCDTS.MCDTSpredictionType(PredLoss, PredMeth)
    optmodel = MCDTS.MCDTSOptimGoal(MCDTS.Prediction_error(PredType,0,samplesize), MCDTS.Continuity_function(13,samplesize))
end
@time tree = mcdts_embedding(Dataset(x1), optmodel, w1, taus, trials; max_depth = max_depth, verbose=true)
best_node = MCDTS.best_embedding(tree)
τ_mcdts2 = best_node.τs
#τ_mcdts2 = [0,18,9]
Y_mcdts2 = genembed(x1, τ_mcdts.*(-1))

## make predictions based on these embedding models
offset = 300
prediction_cao = deepcopy(Y_cao[1:end-offset])
prediction_kennel = deepcopy(Y_kennel[1:end-offset])
prediction_hegger = deepcopy(Y_hegger[1:end-offset])
prediction_mcdts = deepcopy(Y_mcdts[1:end-offset])

K1 = 5
#T_steps = 483
for T = 1:T_steps
    println(T)
    # iterated one step
    predicted_cao, _ = MCDTS.local_zeroth_prediction(prediction_cao, K1; theiler = w1)
    push!(prediction_cao, predicted_cao)
    predicted_kennel, _ = MCDTS.local_zeroth_prediction(prediction_kennel, K1; theiler = w1)
    push!(prediction_kennel, predicted_kennel)
    predicted_hegger, _ = MCDTS.local_zeroth_prediction(prediction_hegger, K1; theiler = w1)
    push!(prediction_hegger, predicted_hegger)
    predicted_mcdts, _ = MCDTS.local_zeroth_prediction(prediction_mcdts, K1; theiler = w1)
    push!(prediction_mcdts, predicted_mcdts)
    # linear
    # predicted_cao, _ = MCDTS.local_linear_prediction(prediction_cao, K1; theiler = w1)
    # push!(prediction_cao, predicted_cao)
    # predicted_kennel, _ = MCDTS.local_linear_prediction(prediction_kennel, K1; theiler = w1)
    # push!(prediction_kennel, predicted_kennel)
    # predicted_hegger, _ = MCDTS.local_linear_prediction(prediction_hegger, K1; theiler = w1)
    # push!(prediction_hegger, predicted_hegger)
    # predicted_mcdts, _ = MCDTS.local_linear_prediction(prediction_mcdts, K1; theiler = w1)
    # push!(prediction_mcdts, predicted_mcdts)
end

## Plot predictions
begin
    sp = length(tr[1:end-offset])-T_steps
    t2 = (-sp+1:T_steps) ./ lyap_time

    figure(figsize=(20,10))
    subplot(4,1,1)
    plot(t2[8000:end], tr[8000:end-offset,1], ".-", label="true data")
    plot(t2[end-T_steps:end], prediction_cao[end-T_steps:end,1], ".-", label="Cao")
    title("x-component (iterated one-step)")
    xlim(-5, 12)
    ylim(-20,20)
    xlabel("Lyapunov time units")
    legend()
    grid()

    subplot(4,1,2)
    plot(t2[8000:end], tr[8000:end-offset,1], ".-", label="true data")
    plot(t2[end-T_steps:end], prediction_kennel[end-T_steps:end,1], ".-", label="Kennel")
    title("x-component (iterated one-step)")
    xlim(-5, 12)
    ylim(-20,20)
    xlabel("Lyapunov time units")
    legend()
    grid()

    subplot(4,1,3)
    plot(t2[8000:end], tr[8000:end-offset,1], ".-", label="true data")
    plot(t2[end-T_steps:end], prediction_hegger[end-T_steps:end,1], ".-", label="Hegger")
    title("x-component (iterated one-step)")
    xlim(-5, 12)
    ylim(-20,20)
    xlabel("Lyapunov time units")
    legend()
    grid()

    subplot(4,1,4)
    plot(t2[8000:end], tr[8000:end-offset,1], ".-", label="true data")
    plot(t2[end-T_steps:end], prediction_mcdts[end-T_steps:end,1], ".-", label="MCDTS")
    title("x-component (iterated one-step)")
    xlim(-5, 12)
    ylim(-20,20)
    xlabel("Lyapunov time units")
    legend()
    grid()

end


text:
# Contains structures and methods related to transient types

struct TcspcTransient
    data::AbstractVector{UInt32}
    t::TcspcTimeAxis
    # Inner constructor for main validation
    function TcspcTransient(tax::TcspcTimeAxis,
        data::AbstractVector{<:Integer})
        # Ensure that number of bins match
        nbins = length(data)
        if nbins != tax.nbins
            error("Number of elements ($nbins) must match length of axis ($(tax.nbins))")
        end
        new(data,tax)
    end
end

function TcspcTransient(tcal::AbstractFloat,
    data::AbstractVector{<:Integer},
    ψ::AbstractFloat)
    # Determine number of bins from data
    nbins = length(data)
    if nbins != nextpow(2,nbins)
        error("nbins ($nbins) must be power of 2")
    end
    # Build time axis
    t = TcspcTimeAxis(tcal,nbins,ψ)
    TcspcTransient(t,data)
end

# Multiple dispatch for padded input
function TcspcTransient(start::AbstractFloat,
    stop::AbstractFloat,data::AbstractVector{<:Integer},
    ψ::AbstractFloat)
    # Construct time axis from start, stop,
    # number of data elements
    nbins = length(data)
    tax = TcspcTimeAxis(start,stop,nbins,ψ)
    # Determine whether we need to pad zeros
    if tax.nbins != nbins
        # Determine where to start padding on both sides
        leftPadEnd = time2chan(tax,start)
        rightPadStart = time2chan(tax,stop)
        # Construct a new data vector
        temp = zeros(eltype(data),size(tax))
        # Fill new array with data
        temp[leftPadEnd:(rightPadStart - 1)] = data
        data = temp
    end
    TcspcTransient(tax.tcal,data)
end

# Helper methods
function poisson_weights(trans::TcspcTransient)
    # Returns 1 / data for all counts above 15
    data = zeros(size(trans.data))
    data[findall(<=(15),trans.data)] .= 1/15
    data[findall(==(0),trans.data)] .= 1
    dind = findall(>(15),trans.data)
    data[dind] = 1 ./ trans.data[dind]
    return data
end

function guess_start_stop(trans::TcspcTransient)
    # Guesses where the start of the curve is,
    # as well as where to stop fitting data
    maxtup = findmax(trans.data)
    # Find point of highest rise in curve
    x = 1:trans.t.nbins
    y = trans.data
    itp = interpolate((x,),y,Gridded(Linear()))
    deriv = only.(Interpolations.gradient.(Ref(itp),x))
    dmaxtup = findmax(deriv)
    # Good approx: 2 * dmaxtup[2] - maxtup[2]
    fit_start = 2 * dmaxtup[2] - maxtup[2]
    # Find first element less than 50
    fit_end = findfirst(<(50),y)
    (fit_start,fit_end)
end

function estimate_irf(trans::TcspcTransient)
    # Estimates instrument response function for curve
    maxtup = findmax(trans.data)
    chantup = guess_start_stop(trans)
    # Offset ≈ tcal * fit_start + (maxtup[2] - fit_start) / 2
    offschan = round(Int,chantup[1] + (maxtup[2] - chantup[1]) / 2)
    # Determine time
    offs = chan2time(trans.t,offschan)
    # FWHM ≈ tcal * (maxtup[2] - fit_start)
    fwhm = chan2time(trans.t,round(Int,maxtup[2] - chantup[1]))
    # N = sum(data)
    N = sum(trans.data)
    # Simulate IRF, then rescale it
    irf = simulate_irf(trans.t,N,fwhm,offs)
    mxirftup = findmax(irf)
    # Scale = mxIRF / mxData
    sc = mxirftup[1] / maxtup[1]
    round(Int,sc * irf)
end


text:
using Test
using ClimaCoupler
using ClimateMachine.GenericCallbacks, ClimateMachine.ODESolvers
using Dates, Unitful, LinearAlgebra

#= Coupled Oscillators
https://scholar.harvard.edu/files/schwartz/files/lecture3-coupled-oscillators.pdf

Two coupled oscillators start from rest but with at least
one perturbed from equilibrium (hence, motion). In our simulation set up,
each individual mass is only aware of its own position, not that of the
other. Through the coupler, the force each mass exerts on the other is passed
to the other mass. A coupler is certainly not necessary to solve the problem,
we could set up a simple matrix system and also have an analytical solution,
however this is a good demonstration of the spectrum of coupling set ups:
from tight monolithic coupling (matrix system) to the black box (our coupler).

| k   κ   k |
|~~~◯~~~◯~~~|
|   m   m   |
|   1   2   |

F11 = -kx1 - κx1
F21 = κx2
F11 + F21 = m x1''

F22 = -kx2 - κx2
F12 = κx1
F22 + F12 = m x2''

dx/dt = v; dv/dt = a = F/m
=#
@testset "Coupled Timestepping" begin
    FT = Float64

    function exact_sol(t, k, κ, m, x1_0, x2_0)
        ωs, ωf = sqrt(k / m), sqrt((k + 2κ) / m)
        As, Af = x1_0 + x2_0, x1_0 - x2_0
        x1 = 0.5 * (As * cos.(ωs .* t) + Af * cos.(ωf .* t))
        x2 = 0.5 * (As * cos.(ωs .* t) - Af * cos.(ωf .* t))
        return x1, x2
    end

    mutable struct RecordState
        u::Vector{Any}
        counter::Int
    end
    RecordState() = RecordState(Any[], 0)

    rs1, rs2 = RecordState(), RecordState()

    mutable struct SpringRHS{T}
        v::T
        F_in::T
    end

    function (f!::SpringRHS)(du, u, p, t, args...; increment = false)
        f!.v .= u[1]
        du[1] = u[2] .+ increment * du[1] # u[1] is position
        du[2] = (-(p.k + p.κ) * f!.v[1] + f!.F_in[1]) / p.m .+ increment * du[2] # u[2] is velocity
    end

    p = (k = FT(4), κ = FT(2), m = FT(1))
    cdt = FT(0.01)
    endtime = 20
    t = 0.0:(cdt * 10):endtime
    nsteps1, nsteps2 = 2, 3

    timestepper(f, state, nsteps) = LSRK54CarpenterKennedy(f, state, dt = cdt / nsteps, t0 = 0.0)

    function main(::Type{FT}) where {FT}
        x1, x2 = [FT(1), FT(0)], [FT(0), FT(0)] # initial position & velocity
        f1! = SpringRHS([FT(0)], [FT(0)])
        f2! = SpringRHS([FT(0)], [FT(0)])

        # cb1, cb2: record total iterations for each model
        cb1 = EveryXSimulationSteps(1) do
            rs1.counter += 1
        end
        cb2 = EveryXSimulationSteps(1) do
            rs2.counter += 1
        end

        mass1 = CplModel(nothing, f1!, nothing, x1, timestepper(f1!, x1, nsteps1), nsteps1, cb1)
        mass2 = CplModel(nothing, f2!, nothing, x2, timestepper(f2!, x2, nsteps2), nsteps2, cb2)

        coupler = CplState()
        coupler_register!(coupler, :F12, [p.κ * x1[1]], nothing, DateTime(0), u"N")
        coupler_register!(coupler, :F21, [p.κ * x2[1]], nothing, DateTime(0), u"N")

        comp1 = (pre_step = prestep1, component_model = mass1, post_step = poststep1)
        comp2 = (pre_step = prestep2, component_model = mass2, post_step = poststep2)
        components = (mass1 = comp1, mass2 = comp2)
        cpl_solver = CplSolver(component_list = components, coupler = coupler, coupling_dt = cdt, t0 = 0.0)

        # set up coupler callbacks
        # cb3, cb4: record state info with time.
        push!(rs1.u, deepcopy(mass1.state[1]))
        cb3 = EveryXSimulationSteps(10) do
            push!(rs1.u, deepcopy(mass1.state[1]))
        end
        push!(rs2.u, deepcopy(mass2.state[1]))
        cb4 = EveryXSimulationSteps(10) do
            push!(rs2.u, deepcopy(mass2.state[1]))
        end
        cbvector = (cb3, cb4)

        return cpl_solver, cbvector
    end

    # presteps: recv force from other mass
    function prestep1(csolver)
        mass1 = csolver.component_list.mass1.component_model
        mass1.discretization.F_in .= coupler_get(csolver.coupler, :F21, nothing, DateTime(0), u"N")
    end

    function prestep2(csolver)
        mass2 = csolver.component_list.mass2.component_model
        mass2.discretization.F_in .= coupler_get(csolver.coupler, :F12, nothing, DateTime(0), u"N")
    end

    # poststeps: send force upon other mass
    function poststep1(csolver)
        mass1 = csolver.component_list.mass1.component_model
        coupler_put!(csolver.coupler, :F12, p.κ * mass1.state[1], nothing, DateTime(0), u"N")
    end

    function poststep2(csolver)
        mass2 = csolver.component_list.mass2.component_model
        coupler_put!(csolver.coupler, :F21, p.κ * mass2.state[1], nothing, DateTime(0), u"N")
    end

    function run(cpl_solver, params, numberofsteps, cbvector)
        solve!(nothing, cpl_solver, params; numberofsteps = numberofsteps, callbacks = cbvector)
    end

    simulation, cbvector = main(FT)
    nsteps = Int(endtime / cdt)
    run(simulation, p, nsteps, cbvector)

    # subiteration execution
    @test rs2.counter == (nsteps2 / nsteps1) * rs1.counter
    # we solve the problem well enough (toggle plots on to check)
    @test norm(exact_sol(t, p..., 1.0, 0.0)[1] .- rs1.u) <= 1e-1
    @test norm(exact_sol(t, p..., 1.0, 0.0)[2] .- rs2.u) <= 1e-1
end

# Quality Assurance Plotting
# using Plots
# p1 = plot(t, exact_sol(t, p..., 1.0, 0.0)[1], label="exact",
#         xlabel="t", ylabel="x1", title="mass 1")
# plot!(t, rs1.u, label="sim", ls=:dash)
# p2 = plot(t, exact_sol(t, p..., 1.0, 0.0)[2], label="exact",
#         xlabel="t", ylabel="x2", title="mass 2")
# plot!(t, rs2.u, label="sim", ls=:dash)
# plot(p1, p2, layout = (2, 1), legend = true, lw=3)

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
#module ExClark
include(joinpath("..",  "src","Lyap.jl"))
include(joinpath("..",  "src","LinProc.jl"))
using Lyap
using LinProc

include(joinpath("..",  "src","misc.jl"))

srand(4)
#srand(5)

si = 0.1; include("excoeff2.jl")

#T = 1.2
T = 4
N = 1201 #design points
#
N = 2201 #design points
K = 50000 #samples
subsample = 100
 # if v=re() is distributed with log density le()
# and yy is a path with drift LinProc.Bcirc(T, v, b, sigma, B, beta, lambda) and diffusion sigma
# (lambda has to be computed for B(T, v) and a(T,v)!)
# then the expectation of llikelicirc(... yy, ...) is 1/pbar
# where pbar = exp(lp(..., B, beta, lambda) - le(v0)))

v0max = [NaN, NaN]

function raremh(K, N, E, sizeE, prop, q, v0,  B, A, bundlesize = 100)
	global res = zeros(int(K/subsample)+1, 11)
	println("Compute P(X in E), K=$K")

	llmin = Inf
	llmax = -Inf
	lplambdamax = -Inf
	lplambdamin = Inf

	L = L2 = 0.0
	V = 0.0
	V2 = 0.0
 	# avoid cancellation
	V0 = BigFloat(0.0)
	V02 = BigFloat(0.0)

	Dt = diff(linspace(0., T, N))
	dt = Dt[1]
	Smax = LinProc.taui(T-dt,T)
	S = linspace(0.,Smax, N)
	Ds = diff(S)
	ds = Ds[1]
	global t = zeros(2)

	lambdav0 = Lyap.lyap(B', -A)
	lpv0 =  LinProc.lp(T, u, v0, B, beta, lambdav0)
	Cpv0 = exp(lpv0)
	yy = zeros(length(u), N)
	y = zeros(length(u))

	for k in bundlesize:bundlesize:K

		if (OS_NAME != :Windows) print("$k $V\r") end

		t[1] += @elapsed for n in 1:bundlesize
	 		DW = randn(2, N-1) .* sqrt(dt)
		 	yy = LinProc.eulerv(0.0, u, b, sigma, Dt, DW)
			v = yy[1:2, N]
			V += 1.*E(v)
			V2 += E(v).^2
		end

		t[2] += @elapsed begin
			# propose new choice of
		 	v0prop = prop(v0)
			lpv0prop = LinProc.lp(T, u, v0prop, B, beta, lambdav0)

		 	if rand() < min(1.0, exp(lpv0prop - lpv0))
		 		v0 = v0prop
		 		lpv0 = lpv0prop
		 	end
		#	println("v0 (", v0[1], ",", v0[2], ")", lpv0, " ", Cpv0/k)

			Cpv0 += exp(lpv0)*bundlesize

	 		# find lambda at the endpoint
			lambda = Lyap.lyap(B', -a(T,v0))
			lplambda = LinProc.lp(T, u, v0, B, beta, lambda)

	 		#yy = LinProc.eulerv(0.0, u, v0, LinProc.Bcirc(T, v0, b, sigma, B, beta, lambda), sigma, Dt, DW)
		 	#ll = LinProc.llikelixcirc(0, T, yy, b, a, B, beta, lambda)
		 	WN = randn(2, bundlesize*(N-1))

			ll = LinProc.transd(S, 0., T, u, v0, b, sigma, a, B, beta, lambda, WN, bundlesize)

			for n in 1:bundlesize
				l =   exp(lplambda + ll[n] - lpv0)*Cpv0/k *sizeE
				L += l
				L2 += l^2
				assert(E(v0) == 1.)
				V0 += E(v0) * l
				V02 += (E(v0) * l).^2

				# running mean and sum of squares
		 		if (ll[n] > llmax) vmax = v0 end
				llmin = min(llmin, ll[n])
				llmax = max(llmax, ll[n])
				lplambdamax =  max(lplambdamax, lplambda)
				lplambdamin = min(lplambdamin, lplambda)
			end
		end

	#println("$L $V0 $V")
		if (0 == k % subsample)
			print("$k:")

			p =  mc2(k, L, L2)
			println(" v ", strmc2(k, V, V2)," v0 ", strmc2(k,float64(V0), float64(V02)), " < min max's ll ", round(llmin, 1), " ", round(llmax, 1), " lp max ", round(lplambdamax,1),">")
			#," v0 ", mc2(k,float64(V0), float64(V02)), " < p $p max's ll ", roundv(llmax, 1)..., " lp ", roundv(lplambdamax,1)...,", v0max ", roundv(vmax,2)...," >" ))
			z1 = mc2(k, V, V2)
			z2 = mc2(k,float64(V0), float64(V02))

			res[k/subsample, :] = [k z1[1] z1[2] z2[1] z2[2] round(llmin, 1)  round(llmax, 1)  round(lplambdamin,1) round(lplambdamax,1) t[1] t[2] ]

		end
	end
	println("Speed ratio: ", t[2]/t[1])
end

function rare(K, N, E, re, pe,  B, A)

	println("Compute P(X in E), K=$K")

	llmin = Inf
	llmax = -Inf
	lplambdamax = -Inf
	lplambdamin = Inf

	L = L2 = 0.0
	V = 0.0
	V2 = 0.0
 	# avoid cancellation
	V0 = BigFloat(0.0)
	V02 = BigFloat(0.0)

	Dt = diff(linspace(0., T, N))
	dt = Dt[1]
	Smax = LinProc.taui(T-dt,T)
	S = linspace(0.,Smax, N)
	Ds = diff(S)
	ds = Ds[1]

	for k in 1:K

		if (OS_NAME != :Windows) print("$k $V\r") end

 		DW = randn(2, N-1) .* sqrt(dt)
	 	yy = LinProc.eulerv(0.0, u, b, sigma, Dt, DW)

		v = yy[1:2, N]

		v0 = re()

 		# find lambda at the endpoint
		lambda = Lyap.lyap(B', -a(T,v0))
		lplambda = LinProc.lp(T, u, v0, B, beta, lambda)

 		#yy = LinProc.eulerv(0.0, u, v0, LinProc.Bcirc(T, v0, b, sigma, B, beta, lambda), sigma, Dt, DW)
	 	#ll = LinProc.llikelixcirc(0, T, yy, b, a, B, beta, lambda)
		DW = randn(2, N-1) .* sqrt(ds)
		u0 =  LinProc.UofX(0,u,  T, v0,  B, beta)

		yy = LinProc.eulerv(0.0, u0, LinProc.bU(T, v0, b, a, B, beta, lambda),  (s,x) -> sqrt(T)*sigma(LinProc.ddd(s,x, 0.0, T, v0,  B, beta)...), Ds, DW)

		ll = LinProc.llikeliU(S, yy, T, v0, b, a,  B, beta, lambda)

		l =  exp(lplambda + ll)/pe(v0)

 		# running mean and sum of squares
 		if (ll > llmax) vmax = v0 end
		llmin = min(llmin, ll)
		llmax = max(llmax, ll)
		lplambdamax =  max(lplambdamax, lplambda)
		lplambdamin = min(lplambdamin, lplambda)

		L += l
		L2 += l^2
		assert(E(v0) == 1.)
		V0 += E(v0) * l
		V02 += (E(v0) * l).^2
		V += 1.*E(v)
		V2 += E(v).^2
		#println("$L $V0 $V")
		if (0 == k % subsample)
			print("$k:")

			p =  mc2(k, L, L2)
			println(" v ", strmc2(k, V, V2)," v0 ", strmc2(k,float64(V0), float64(V02)), " < min max's ll ", round(llmin, 1), " ", round(llmax, 1), " lp max ", round(lplambdamax,1),">")
			#," v0 ", mc2(k,float64(V0), float64(V02)), " < p $p max's ll ", roundv(llmax, 1)..., " lp ", roundv(lplambdamax,1)...,", v0max ", roundv(vmax,2)...," >" ))
			z1 = mc2(k, V, V2)
			z2 = mc2(k,float64(V0), float64(V02))

			res[k/subsample, :] = [k z1[1] z1[2] z2[1] z2[2] round(llmin, 1)  round(llmax, 1)  round(lplambdamin,1) round(lplambdamax,1)  ]

		end
	end
end

function dens(K, N, v0, t, T, B, A)

	println("Compute P(XT in dv), K=$K, T=$T, t=$t")
	println("E ~p(t, Xt)  =  ~p(0,u) * E exp(D(X°)(t)) ")

	llmax = [Inf,-Inf]
	llxmax = [Inf,-Inf]

	# avoid cancellation
	Lx = BigFloat(0.0)
	Lx2 = BigFloat(0.0)
	Lo = BigFloat(0.0)
	Lo2 = BigFloat(0.0)

	Dt = diff(linspace(0., t, N))
	dt = Dt[1]
	Smax = LinProc.taui(t,T)
	S = linspace(0.,Smax, N)
	Ds = diff(S)
	ds = Ds[1]
	lambda = Lyap.lyap(B', -A)
	lplambda = LinProc.lp(t, u, v0, B, beta, lambda) #
	gamma = inv(A)
	for k in 1:K

		if (OS_NAME != :Windows) print("$k \r") end

		DW = randn(2, N-1) .* sqrt(dt)
	 	X = LinProc.eulerv(0.0, u, b, sigma, Dt, DW)
		x = X[1:2, N]
		if (T-t) > 0.01 #if T-t too small, lp becomes too unstable for
			llx = LinProc.lp(T-t, x, v0, B, beta, lambda)
		else
			llx = LinProc.lp0(T-t, x, v0, B*v0 + beta, gamma)
		end

 		DW = randn(2, N-1) .* sqrt(ds)
		u0 =  LinProc.UofX(0,u,  T, v0,  B, beta)
		yy = LinProc.eulerv(0.0, u0, LinProc.bU(T, v0, b, a, B, beta, lambda),  (s,z) -> sqrt(T)*sigma(LinProc.ddd(s,z, 0.0, T, v0,  B, beta)...),  Ds, DW)

		ll = LinProc.llikeliU(S, yy, T, v0, b, a,  B, beta, lambda)

		lo =  exp(lplambda + ll)

		# running mean and sum of squares
		llmax = [min(llmax[1], ll), max(llmax[2], ll)]
		llxmax = [min(llxmax[1], llx), max(llxmax[2], llx)]

		Lo += lo
		Lo2 += lo^2
		Lx += exp(llx)
		Lx2 += exp(2*llx)

		if (0 == k % subsample)
			print("$k:")
			res[k/subsample,:] = [k mc2(k, float64(Lx), float64(Lx2))... mc2(k,float64(Lo), float64(Lo2))... roundv(llmax, 3)... roundv(llxmax,3)...]

			println(" ", strmc2(k, float64(Lx), float64(Lx2)), " ~ ", strmc2(k,float64(Lo), float64(Lo2)), " < max's llo ", roundv(llmax, 1), " llx ", roundv(llxmax,1)," >" )

		end
	end
	res

end

#v = 1.4*LinProc.mu(T, u, B, beta) # point not to far away from the mean

v= [-0.323933,  0.494032] #same as A->B in Example 2

##### proposal distribution to explore rare event E

#test, whether X(T) in E
#ep = 0.05
ep = 0.1
function E(x)
 max(abs(x - v)) <= ep
end
sizeE = (2*ep)^2

# a density, supported on E
function pe(x)
 E(x) / (2*ep)^2
end

# sampling from pe
function re()
  v + 2*ep*rand(2)-ep
end

#rare(K, N, E, re, pe)
#
println("try ExClark.D1()\n\t (= dens(K, N, v, t, T, B, a(T,v)))")

function D1()
 dens(K, N, v, 0.999*T, T, exp(-0.2*T)*bB, a(T,v))
end
function D2()
 rare(K, N, E, re, pe, exp(-0.2*T)*bB, a(T,v) )
end

function D3()
 raremh(K, N, E, sizeE,  v -> re(), null, v, exp(-0.2*T)*bB,a(T,v), 10)
end
#end


text:
# test_parallel_derivatives.jl

function test_parallel_derivatives()

  @testset "--- Checking Complete Derivative w.r.t angle of attack for a 2D airfoil ---" begin

    fname = "./input_vals_airfoil_parallel.jl"

    # Get the adjoint vector
    mesh, sbp, eqn, opts = solvePDE(fname)
    lift = createFunctional(mesh, sbp, eqn, opts, 1)
    val = EulerEquationMod.evalFunctional(mesh, sbp, eqn, opts, lift)

    @testset "Checking functional J in parallel" begin

      serial_vals = readdlm("BoundaryForceData_values.dat")
      @test isapprox( val, serial_vals[1]) atol= 1e-14
    end # End  testset("Checking ∂J/∂aoa")

    h = 1e-20
    pert = Complex128(0, h)

    pc, lo = getNewtonPCandLO(mesh, sbp, eqn, opts)
    ls = StandardLinearSolver(pc, lo, eqn.comm, opts)
    adjoint_vec = zeros(Complex128, mesh.numDof)
    calcAdjoint(mesh, sbp, eqn, opts, ls, lift, adjoint_vec, recalc_jac=true, recalc_pc=true)

    # compute dR/daoa and dJ/daoa
    eqn.params.aoa += pert
    evalResidual(mesh, sbp, eqn, opts)
    array3DTo1D(mesh, sbp, eqn, opts, eqn.res, eqn.res_vec)
    dRdaoa = imag(eqn.res_vec)/h

    val = evalFunctional(mesh, sbp, eqn, opts, lift)
    dJdaoa = imag(val)/h
    eqn.params.aoa -= pert

    term2 = dot(adjoint_vec, dRdaoa)
    term2 = MPI.Allreduce(term2, MPI.SUM, eqn.comm)
    DJDaoa = dJdaoa + term2

    # Check complete derivatives w.r.t alpha using finite difference
    pert = 1e-6
    eqn.params.aoa += pert
    solvePDE(mesh, sbp, eqn, opts, mesh)
    val2 = evalFunctional(mesh, sbp, eqn, opts, lift)
    eqn.params.aoa -= pert

    DJDaoa_fd = (val2 - real(val))/pert
    err_val = norm(DJDaoa - DJDaoa_fd, 2)

    @test isapprox( err_val, 0.0) atol= 1e-6

#=
    lift = lift.lift_val
    pc, lo = getNewtonPCandLO(mesh, sbp, eqn, opts)
    ls = StandardLinearSolver(pc, lo, eqn.comm, opts)
    adjoint_vec = zeros(Complex128, mesh.numDof)

    calcAdjoint(mesh, sbp, eqn, opts, ls, lift, adjoint_vec, recalc_jac=true, recalc_pc=true)
    dJdaoa = EulerEquationMod.eval_dJdaoa(mesh, sbp, eqn, opts, lift, "lift", adjoint_vec)

    # Check complete derivatives w.r.t alpha using finite difference
    pert = 1e-6
    eqn.params.aoa += pert
    solvePDE(mesh, sbp, eqn, opts, mesh)
    EulerEquationMod.evalFunctional(mesh, sbp, eqn, opts, lift)
    lift_pert = lift.lift_val
    drag_pert = lift.drag_val
    eqn.params.aoa -= pert

    dJdaoa_fd = (lift_pert - lift)/pert

    @testset "Check complete derivative dJ/daoa" begin
      deriv_err = norm(dJdaoa_fd - dJdaoa, 2)
      @test isapprox( deriv_err, 0.0) atol=1e-6
    end
=#
  end # End Facts Checking Complete Derivative w.r.t angle of attack for a 2D airfoil

  return nothing
end

add_func1!(EulerTests, test_parallel_derivatives, [TAG_PARALLEL_DERIVATIVES, TAG_LONGTEST])


text:
include(joinpath(Pkg.dir("MathProgBase"),"test","nlp.jl"))
convexnlptest(KnitroSolver())
rosenbrocktest(KnitroSolver())
nlptest(KnitroSolver())
nlptest_nohessian(KnitroSolver())


text:
function build_arch(name::String,A)
    num_fifos = 2

    f = open(name * ".json")
    s = JSON.parse(f)

    arch = TopLevel{A,2}(s["name"])
    row = s["dimensions"]["row"]
    col = s["dimensions"]["column"]
    processor_add = s["processor_tile"]["add"]
    processor_remove = s["processor_tile"]["remove"]
    one_port = s["memory"]["1_port"]
    two_port = s["memory"]["2_port"]
    num_links = s["configurations"]["num_links"]
    input_handler = s["input_handler"]
    output_handler = s["output_handler"]
    special_mem_connection = s["configurations"]["special_mem_connection"]

    # memories connected in a non-asap way are dealt with below
    offset = Array{Tuple{Int64,Int64},1}()
    if(special_mem_connection)
      tmp = s["configurations"]["offsets"]
      for off in tmp
        push!(offset, stringtotuple(off))
      end
    end

    close(f)

    # preallocate arrays to store Tuple addresses
    in_addr = Array{Tuple{Int64,Int64},1}()
    out_addr = Array{Tuple{Int64,Int64},1}()
    proc = Array{Tuple{Int64,Int64},1}()

    # push the tuples into their respective arrays
    for handler in input_handler
      push!(in_addr,stringtotuple(handler))
    end
    for handler in output_handler
      push!(out_addr,stringtotuple(handler))
    end

    # get all the processors
    # add/remove processors to and from the processor array
    for r in 0:row-1, c in 0:col-1
      match = false
      for rem in processor_remove
        if(rem == "($r,$c)")
          match = true
        end
      end
      match && continue
      push!(proc,(r,c))
    end
    for p in processor_add
      push!(proc,stringtotuple(p))
    end

    one_port_mem = Array{Tuple{Int64,Int64},1}()
    two_port_mem = Array{Tuple{Int64,Int64},1}()
    memory_1port = Array{AsapMapper.MemoryLocation{2},1}()
    memory_2port = Array{AsapMapper.MemoryLocation{2},1}()
    extra_two_port_mem = Array{Tuple{Int64,Int64},1}()

    # get all mems
    count = 1
    for mem in one_port
      addr = stringtotuple(mem)
      push!(one_port_mem,addr)
      if(special_mem_connection)
        push!(memory_1port, MemoryLocation(addr,[offset[count]]))
        count += 1
      else
        push!(memory_1port, MemoryLocation(addr,[(-1,0)]))
      end
    end
    for mem in two_port
      addr = stringtotuple(mem)
      push!(two_port_mem,addr)
      if(special_mem_connection)
        second_address = other_address(addr,offset,count)
        push!(memory_2port, MemoryLocation(addr,second_address,[offset[count],offset[count+1]]))
        push!(extra_two_port_mem,second_address)
        count += 2
      else
        push!(memory_2port, MemoryLocation(addr,[(-1,0),(-1,1)]))
      end
    end
    mem_neighbor_addrs = vcat(mem_neighbor(memory_2port),mem_neighbor(memory_1port))

    ####################
    # Normal Processor #
    ####################
    # Get a processor tile and instantiate it.
    processor = build_processor_tile(num_links)
    for p in proc
      overlap = false
      one_match = false
      two_match = false
      ex_two_match = false
      for mem_neighbor_addr in mem_neighbor_addrs
        if(mem_neighbor_addr == CartesianIndex(p[1],p[2]))
          overlap = true
        end
      end
      for mem in one_port_mem
        if(mem == (p[1],p[2]))
          one_match = true
        end
      end
      for mem in two_port_mem
        if(mem == (p[1],p[2]))
          two_match = true
        end
      end
      for mem in extra_two_port_mem
        if(mem == (p[1],p[2]))
          ex_two_match = true
        end
      end
      (overlap || one_match || two_match || ex_two_match) && continue
      add_child(arch, processor, CartesianIndex(p[1],p[2]))
    end

    ####################
    # Memory Processor #
    ####################
    # Instantiate a memory processor at every address neighboring a memory.
    memory_processor = build_processor_tile(num_links, include_memory = true)
    for mem_neighbor_addr in mem_neighbor_addrs
      add_child(arch, memory_processor, mem_neighbor_addr)
    end

    #################
    # 2 Port Memory #
    #################
    memory_twoport = build_memory(2)
    for mem in two_port_mem
          add_child(arch, memory_twoport, CartesianIndex(mem[1],mem[2]))
    end

    #################
    # 1 Port Memory #
    #################
    memory_oneport = build_memory(2)
    for mem in one_port_mem
          add_child(arch, memory_oneport, CartesianIndex(mem[1],mem[2]))
    end

    #################
    # Input Handler #
    #################
    input_handler = build_input_handler(1)
    for handler in in_addr
      add_child(arch, input_handler, CartesianIndex(handler[1],handler[2]))
    end

    ##################
    # Output Handler #
    ##################
    output_handler = build_output_handler(1)
    for handler in out_addr
      add_child(arch, output_handler, CartesianIndex(handler[1],handler[2]))
    end

    connect_processors(arch, num_links)
    connect_io(arch, num_links)
    if (special_mem_connection)
      if(one_port_mem != nothing)
        connect_memories_cluster(arch, memory_1port)
      end
      if(two_port_mem != nothing)
        connect_memories_cluster(arch, memory_2port)
      end
    else
      connect_memories(arch)
    end

    return arch
end

function stringtotuple(string)

  row,column = split(string,",")
  row = parse(Int64,split(row,"(")[2])
  column = parse(Int64,split(column,")")[1])

  return (row,column)
end

function other_address(addr, offset, count)
  if(offset[count][1] != offset[count+1][1])
    second_addr = (addr[1]+offset[count+1][1],addr[2])
  elseif(offset[count][2] != offset[count+1][2])
    second_addr =  (addr[1],addr[2]+offset[count+1][2])
  end
  return second_addr
end


text:
"""
Computes mean resultant vector length for circular data.

	1. α: sample of angles in radians
    - w: number of incidences in case of binned angle data
    - d: spacing of bin centers for binned data, used to correct for bias in estimation of r, in radians
    - dims: compute along this dimension(default=1)

	return: mean resultant length
"""
function circ_r(α; w = ones(size(α)), d = 0, dims = 1)
  # compute weighted sum of cos and sin of angles
  r = sum(w .* cis.(α); dims)

  # obtain length
  r = abs.(r) ./ sum(w; dims)

  # for data with known spacing, apply correction factor to correct for bias in the estimation of r (see Zar, p. 601, equ. 26.16)
  if d != 0
    c = d / 2 / sin(d / 2)
    r *= c
  end
  length(r) == 1 ? r[1] : r
end

"""
Computes the mean direction for circular data.

	1. α: sample of angles in radians
    - w: weightings in case of binned angle data
    - dims: compute along this dimension(default=1)

	return:
	- μ: mean direction
	- ul: upper 95% confidence limit
	- ll: lower 95% confidence limit
"""
function circ_mean(α; w = ones(size(α)), dims = 1)
  # compute weighted sum of cos and sin of angles
  r = sum(w .* cis.(α); dims)

  # obtain mean by
  μ = angle.(r)
  μ = length(μ) == 1 ? μ[1] : μ

  # confidence limits
  t = circ_confmean(α; xi = 0.05, w, d = 0, dims)
  ul = μ .+ t
  ll = μ .- t

  return (; μ, ul, ll)
end

"""
Transforms p-axial data to a common scale.

  1. α: sample of angles in radians
  2. p: number of modes(default=1)

  return: transformed data
"""
circ_axial(α, p = 1) = mod2pi.(α * p)

"""
Computes the median direction for circular data.

  1. α: sample of angles in radians

  return: median direction
"""
function circ_median(α::AbstractVector)
  α = mod2pi.(α)
  n = length(α)

  dd = circ_dist2(α)
  m1 = dropdims(count(>=(0), dd; dims = 1), dims = 1)
  m2 = dropdims(count(<=(0), dd; dims = 1), dims = 1)

  dm = abs.(m1 .- m2)
  m = minimum(dm)
  md = circ_mean(α[dm.==m]).μ
  αμ = circ_mean(α).μ
  if abs(circ_dist(αμ, md)) > abs(circ_dist(αμ, md + π))
    md = mod2pi(md + π)
  end
  return md
end

"""
Computes circular standard deviation for circular data (equ. 26.20, Zar).

	1. α: sample of angles in radians
   - w: weightings in case of binned angle data
   - d: spacing of bin centers for binned data, used to correct for bias in estimation of r, in radians
   - dims: compute along this dimension(default=1)

	return:
  - s: angular deviation
  - s0: circular standard deviation
"""
function circ_std(α; w = ones(size(α)), d = 0, dims = 1)
  # compute mean resultant vector length
  r = circ_r(α; w, d, dims)

  s = sqrt.(2 .* (1 .- r))      # 26.20
  s0 = sqrt.(-2 .* log.(r))   # 26.21
  return (; s, s0)
end

"""
Computes circular variance for circular data (equ. 26.17/18, Zar).

	1. α: sample of angles in radians
  - w: number of incidences in case of binned angle data
  - d: spacing of bin centers for binned data, used to correct for bias in estimation of r, in radians
	- dims: compute along this dimension(default=1)

	return:
	- S: circular variance 1-r
	- s: angular variance 2(1-r)
"""
function circ_var(α; w = ones(size(α)), d = 0, dims = 1)
  # compute mean resultant vector length
  r = circ_r(α; w, d, dims)

  # apply transformation to var
  S = 1 .- r
  s = 2 * S
  return (; S, s)
end

"""
Computes the confidence limits on the mean for circular data.

	1. α: sample of angles in radians
   - xi: (1-xi) confidence limits are computed, default=0.05
   - w: number of incidences in case of binned angle data
   - d: spacing of bin centers for binned data, used to correct for bias in estimation of r, in radians
   - dims: compute along this dimension(default=1)

	return: mean ± d yields upper/lower (1-xi)% confidence limit
"""
function circ_confmean(α; xi = 0.05, w = ones(size(α)), d = 0, dims = 1)
  # compute ingredients for conf. lim.
  r = circ_r(α; w, d, dims)
  n = sum(w; dims)
  R = n .* r
  c2 = quantile(Chisq(1), 1 - xi)

  # check for resultant vector length and select appropriate formula
  t = zeros(size(r))

  for i = 1:length(r)
    if r[i] < 0.9 && r[i] > sqrt(c2 / 2 / n[i])
      t[i] = sqrt((2 * n[i] * (2 * R[i]^2 - n[i] * c2)) / (4 * n[i] - c2))  # equ. 26.24
    elseif r[i] >= 0.9
      t[i] = sqrt(n[i]^2 - (n[i]^2 - R[i]^2) * exp(c2 / n[i]))      # equ. 26.25
    else
      t[i] = NaN
      @warn "Requirements for confidence levels not met."
    end
  end

  # apply final transform
  t = acos.(clamp.(t ./ R, -1, 1))
  t = length(t) == 1 ? t[1] : t
end

"""
Calculates a measure of angular skewness.

	1. α: sample of angles in radians
    - w: weightings in case of binned angle data
    - dims: compute along this dimension(default=1)

	return:
	- b: skewness (from Pewsey)
	- b0: alternative skewness measure (from Fisher)
"""
function circ_skewness(α; w = ones(size(α)), dims = 1)
  # compute neccessary values
  R = circ_r(α; w, d = 0, dims)
  θ = circ_mean(α; w, dims).μ
  _, ρ₂, μ₂ = circ_moment(α; w, p = 2, cent = false, dims)

  # compute skewness
  if ndims(θ) > 0
    θ₂ = repeat(θ, outer = Int.(size(α) ./ size(θ)))
  else
    θ₂ = θ
  end
  b = sum(w .* sin.(2 * circ_dist(α, θ₂)); dims) ./ sum(w; dims)
  b0 = ρ₂ .* sin.(circ_dist(μ₂, 2θ)) ./ (1 .- R) .^ (3 / 2)   # (formula 2.29)
  b = length(b) == 1 ? b[1] : b
  b0 = length(b0) == 1 ? b0[1] : b0
  return (; b, b0)
end

"""
Calculates a measure of angular kurtosis.

  1. α: sample of angles in radians
    - w: weightings in case of binned angle data
    - dims: compute along this dimension(default=1)

  return:
  - k: kurtosis (from Pewsey)
  - k0: kurtosis (from Fisher)
"""
function circ_kurtosis(α; w = ones(size(α)), dims = 1)
  # compute mean direction
  R = circ_r(α; w, d = 0, dims)
  θ = circ_mean(α; w, dims).μ
  _, ρ₂, _ = circ_moment(α; w, p = 2, cent = true, dims)
  _, _, μ₂ = circ_moment(α; w, p = 2, cent = false, dims)

  # compute skewness
  if ndims(θ) > 0
    θ₂ = repeat(θ, outer = Int.(size(α) ./ size(θ)))
  else
    θ₂ = θ
  end
  k = sum(w .* cos.(2 * circ_dist(α, θ₂)); dims) ./ sum(w; dims)
  k0 = (ρ₂ .* cos.(circ_dist(μ₂, 2θ)) .- R .^ 4) ./ (1 .- R) .^ 2    # (formula 2.30)
  k = length(k) == 1 ? k[1] : k
  k0 = length(k0) == 1 ? k0[1] : k0
  return (; k, k0)
end

"""
Calculates the complex p-th centred or non-centred moment of the angular data in angle.

  1. α: sample of angles in radians
   - w: number of incidences in case of binned angle data
   - p: p-th moment to be computed, default=1
   - cent: if true, central moments are computed, default = false
   - dims: compute along this dimension(default=1)

  return:
  - mp: complex p-th moment
  - ρp: magnitude of the p-th moment
  - μp: angle of th p-th moment
"""
function circ_moment(α; w = ones(size(α)), p = 1, cent = false, dims = 1)
  if cent
    θ = circ_mean(α; w, dims).μ
    if ndims(θ) > 0
      v = Int.(size(α) ./ size(θ))
      θ = repeat(θ, outer = v)
    end
    α = circ_dist(α, θ)
  end

  n = size(α, dims)
  c̄ = sum(cos.(p * α) .* w; dims) / n
  s̄ = sum(sin.(p * α) .* w; dims) / n
  mp = c̄ .+ im * s̄
  ρp = abs.(mp)
  μp = angle.(mp)
  mp = length(mp) == 1 ? mp[1] : mp
  ρp = length(ρp) == 1 ? ρp[1] : ρp
  μp = length(μp) == 1 ? μp[1] : μp
  return (; mp, ρp, μp)
end

"""
Pairwise difference α-β around the circle computed efficiently.

  1. α: sample of linear random variable
  2. β: sample of linear random variable or one single angle

  return: matrix with differences
"""
circ_dist(α, β) = angle.(cis.(α) ./ cis.(β))

"""
All pairwise difference α-β around the circle computed efficiently.

  1. α: sample of linear random variable
  2. β: sample of linear random variable

  return: matrix with pairwise differences
"""
circ_dist2(α) = circ_dist2(α, α)
circ_dist2(α, β) = angle.(cis.(α) ./ cis.(β'))

"""
Computes descriptive statistics for circular data.

  1. α: sample of angles in radians
    - w: weightings in case of binned angle data
    - d: spacing of bin centers for binned data, used to correct for bias in estimation of r, in radians

  return: descriptive statistics
"""
function circ_stats(α::AbstractVector; w = ones(size(α)), d = 0)
  # mean
  mean = circ_mean(α; w).μ
  # median
  median = circ_median(α)
  # variance
  var = circ_var(α; w, d)
  # standard deviation
  std, std0 = circ_std(α; w, d)
  # skewness
  skewness, skewness0 = circ_skewness(α; w)
  # kurtosis
  kurtosis, kurtosis0 = circ_kurtosis(α; w)

  return (; mean, median, var,std,std0, skewness,skewness0,kurtosis,kurtosis0)
end

"""
Computes Rayleigh test for non-uniformity of circular data.

	H0: the population is uniformly distributed around the circle
  HA: the populatoin is not distributed uniformly around the circle
	Assumption: the distribution has maximally one mode and the data is
							sampled from a von Mises distribution!

	1. α: sample of angles in radians
    - w: number of incidences in case of binned angle data
    - d: spacing of bin centers for binned data, used to correct for bias in estimation of r, in radians

	return:
	- p: p-value of Rayleigh's test
	- z: value of the z-statistic
"""
function circ_rtest(α; w = ones(size(α)), d = 0)
  r = circ_r(α; w, d)
  n = sum(w)

  # compute Rayleigh's R (equ. 27.1)
  R = n * r

  # compute Rayleigh's z (equ. 27.2)
  z = R^2 / n

  # compute p value using approxation in Zar, p. 617
  p = exp(sqrt(1 + 4n + 4(n^2 - R^2)) - (1 + 2n))
  return (; p, z)
end

"""
Computes Omnibus or Hodges-Ajne test for non-uniformity of circular data.

	H0: the population is uniformly distributed around the circle
	HA: the population is not distributed uniformly around the circle
	Alternative to the Rayleigh and Rao's test. Works well for unimodal,
	bimodal or multimodal data. If requirements of the Rayleigh test are met, the latter is more powerful.

	1. α: sample of angles in radians
    - sz: step size for evaluating distribution, default 1 degree
    - w: number of incidences in case of binned angle data

	return:
	- p: p-value
	- m: minimum number of samples falling in one half of the circle
"""
function circ_otest(α; sz = deg2rad(1), w = ones(size(α)))
  α = mod2pi.(α)
  n = sum(w)
  dg = 0:sz:π

  m1 = zeros(size(dg))
  m2 = zeros(size(dg))
  for i = 1:length(dg)
    m1[i] = sum((α .> dg[i]) .& (α .< π + dg[i]) .* w)
    m2[i] = n - m1[i]
  end
  m = min(minimum(m1), minimum(m2))

  if n > 50
    # approximation by Ajne (1968)
    A = π * sqrt(n) / 2 / (n - 2m)
    p = sqrt(2π) / A * exp(-π^2 / 8 / A^2)
  else
    # exact formula by Hodges (1955)
    # p = 2^(1-n) * (n-2m) * nchoosek(n,m)  # revised below for numerical stability
    p = exp((1 - n) * log(2) + log(n - 2m) + loggamma(n + 1) - loggamma(m + 1) - loggamma(n - m + 1))
  end
  return (; p, m)
end

"""
Calculates Rao's spacing test by comparing distances between points on a circle to those expected from a uniform distribution.
  H0: Data is distributed uniformly around the circle.
  H1: Data is not uniformly distributed around the circle.

  Alternative to the Rayleigh test and the Omnibus test. Less powerful
  than the Rayleigh test when the distribution is unimodal on a global scale but uniform locally.

  Due to the complexity of the distributioin of the test statistic, we resort to the tables published by
  Russell, Gerald S. and Levitin, Daniel J.(1995) An expanded table of probability values for rao's spacing test, Communications in Statistics - Simulation and Computation
  Therefore the reported p-value is the smallest α level at which the test would still be significant.
  If the test is not significant at the α=0.1 level, we return the critical value for α = 0.05 and p = 0.5.

  1. α: sample of angles in radians

  return:
  - p: smallest p-value at which test would be significant
  - u: computed value of the test-statistic u
  - UC: critical value of the test statistic at sig-level
"""
function circ_raotest(α)
  # for the purpose of the test, convert to angles
  α = sort(rad2deg.(α))
  n = length(α)

  # compute test statistic
  u = 0
  λ = 360 / n
  for j = 1:n-1
    ti = α[j+1] - α[j]
    u += abs(ti - λ)
  end

  tn = 360 - α[n] + α[1]
  u = 0.5(u + abs(tn - λ))

  getVal = (N, u) -> begin
      # Table II from Russel and Levitin, 1995
      α = [0.001, 0.01, 0.05, 0.10]
      table = [
        4 247.32 221.14 186.45 168.02
        5 245.19 211.93 183.44 168.66
        6 236.81 206.79 180.65 166.30
        7 229.46 202.55 177.83 165.05
        8 224.41 198.46 175.68 163.56
        9 219.52 195.27 173.68 162.36
        10 215.44 192.37 171.98 161.23
        11 211.87 189.88 170.45 160.24
        12 208.69 187.66 169.09 159.33
        13 205.87 185.68 167.87 158.50
        14 203.33 183.90 166.76 157.75
        15 201.04 182.28 165.75 157.06
        16 198.96 180.81 164.83 156.43
        17 197.05 179.46 163.98 155.84
        18 195.29 178.22 163.20 155.29
        19 193.67 177.08 162.47 154.78
        20 192.17 176.01 161.79 154.31
        21 190.78 175.02 161.16 153.86
        22 189.47 174.10 160.56 153.44
        23 188.25 173.23 160.01 153.05
        24 187.11 172.41 159.48 152.68
        25 186.03 171.64 158.99 152.32
        26 185.01 170.92 158.52 151.99
        27 184.05 170.23 158.07 151.67
        28 183.14 169.58 157.65 151.37
        29 182.28 168.96 157.25 151.08
        30 181.45 168.38 156.87 150.80
        35 177.88 165.81 155.19 149.59
        40 174.99 163.73 153.82 148.60
        45 172.58 162.00 152.68 147.76
        50 170.54 160.53 151.70 147.05
        75 163.60 155.49 148.34 144.56
        100 159.45 152.46 146.29 143.03
        150 154.51 148.84 143.83 141.18
        200 151.56 146.67 142.35 140.06
        300 148.06 144.09 140.57 138.71
        400 145.96 142.54 139.50 137.89
        500 144.54 141.48 138.77 137.33
        600 143.48 140.70 138.23 136.91
        700 142.66 140.09 137.80 136.59
        800 142.00 139.60 137.46 136.33
        900 141.45 139.19 137.18 136.11
        1000 140.99 138.84 136.94 135.92]

      ridx = findfirst(table[:, 1] .>= N)
      cidx = findfirst(table[ridx, 2:end] .< u)

      if isnothing(cidx)
        UC = table[ridx, end-1]
        p = 0.5
      else
        UC = table[ridx, cidx+1]
        p = α[cidx]
      end
      return p, UC
    end

  # get critical value from table
  p, UC = getVal(n, u)

  return (; p, u, UC)
end

"""
Computes V test for non-uniformity of circular data with a specified mean direction.

  H0: the population is uniformly distributed around the circle
  HA: the population is not distributed uniformly around the circle but has a mean of `m`.

  !!!note
    Not rejecting H0 may mean that the population is uniformly distributed around the circle OR
    that it has a mode but that this mode is not centered at `m`.

  The V test has more power than the Rayleigh test and is preferred if there is reason to believe in a specific mean direction.

  1. α: sample of angles in radians
  - m: suspected mean direction, default=0
  - w: number of incidences in case of binned angle data
  - d: spacing of bin centers for binned data, used to correct for bias in estimation of r, in radians

  return:
  - p: p-value of V test
  - v: value of the V statistic
"""
function circ_vtest(α; m = 0, w = ones(size(α)), d = 0)
  # compute some ingredients
  r = circ_r(α; w, d)
  μ, = circ_mean(α; w)
  n = sum(w)

  # compute Rayleigh's R (equ. 27.1)
  R = n * r

  # compute the V statistic (equ. 27.5)
  v = R * cos(μ - m)

  # compute u (equ. 27.6)
  u = v * sqrt(2 / n)

  # compute p-value from one tailed normal approximation
  p = 1 - cdf(Normal(), u)

  return (; p, v)
end

"""
Tests for significance of the median.

  H0: the population has median angle `md`
  HA: the population has not median angle `md`

  1. α: sample of angles in radians
  - md: median to test, default=0

  return: p-value
"""
function circ_medtest(α; md = 0)
  n = length(α)

  # compute deviations from median
  d = circ_dist(α, md)

  n1 = sum(d .< 0)
  n2 = sum(d .> 0)

  # compute p-value with binomial test
  sum(pdf.(Binomial(n, 0.5), [0:min(n1, n2); max(n1, n2):n]))
end

"""
One-Sample test for the mean angle.
  H0: the population has mean `m`.
  HA: the population has not mean `m`.

  !!!note: This is the equvivalent to a one-sample t-test with specified mean direction.

  1. α: sample of angles in radians
  - m: assumed mean direction, default=0
  - w: number of incidences in case of binned angle data
  - d: spacing of bin centers for binned data, used to correct for bias in estimation of r, in radians
  - xi: alpha level of the test

  return:
  - h: false if H0 can not be rejected, true otherwise
  - μ: mean
  - ul: upper (1-xi) confidence level
  - ll: lower (1-xi) confidence level
"""
function circ_mtest(α; xi = 0.05, m = 0, w = ones(size(α)), d = 0)
  # compute ingredients
  μ, = circ_mean(α; w)
  t = circ_confmean(α; xi, w, d)
  ul = μ + t
  ll = μ - t

  # compute test via confidence limits (example 27.3)
  h = abs(circ_dist2(m, μ)) > t

  return (; h, μ, ul, ll)
end

"""
Parametric Watson-Williams multi-sample test for equal means. Can be used as a one-way ANOVA test for circular data.

  H0: the s populations have equal means
  HA: the s populations have unequal means

  !!! note
  Use with binned data is only advisable if binning is finer than 10 deg.
  In this case, α is assumed to correspond to bin centers.

  The Watson-Williams two-sample test assumes underlying von-Mises distributrions.
  All groups are assumed to have a common concentration parameter k.

  1. α: angles in radians
  2. idx: indicates which population the respective angle in α comes from, 1:s
    - w: number of incidences in case of binned angle data

  return:
  - p: p-value of the Watson-Williams multi-sample test. Discard H0 if p is small.
  - F: F statistics
"""
function circ_wwtest(α, idx; w = ones(size(α)))
  # number of groups
  u = unique(idx)
  s = length(u)

  # number of samples
  n = sum(w)

  # compute relevant quantitites
  pn = zeros(s)
  pr = zeros(s)
  for t = 1:s
    pidx = idx .== u[t]
    pn[t] = sum(pidx .* w)
    pr[t] = circ_r(α[pidx]; w = w[pidx])
  end

  r = circ_r(α; w)
  rw = sum(pn .* pr) / n

  # make sure assumptions are satisfied
  checkAssumption =
    (rw, n) -> begin
      if n >= 11 && rw < 0.45
        @warn "Test not applicable. Average resultant vector length < 0.45."
      elseif n < 11 && n >= 7 && rw < 0.5
        @warn "Test not applicable. Average number of samples per population 6 < x < 11 and average resultant vector length < 0.5."
      elseif n >= 5 && n < 7 && rw < 0.55
        @warn "Test not applicable. Average number of samples per population 4 < x < 7 and average resultant vector length < 0.55."
      elseif n < 5
        @warn "Test not applicable. Average number of samples per population < 5."
      end
    end
  checkAssumption(rw, mean(pn))

  # test statistic
  kk = circ_kappa(rw)
  β = 1 + 3 / (8 * kk)    # correction factor
  A = sum(pr .* pn) - r * n
  B = n - sum(pr .* pn)

  F = β * (n - s) * A / (s - 1) / B
  # p = 1 - fcdf(F, s - 1, n - s)
  p = 1 - cdf(FDist(s - 1, n - s), F)

  # fprintf('\nANALYSIS OF VARIANCE TABLE (WATSON-WILLIAMS TEST)\n\n');
  # fprintf('%s\t\t\t\t%s\t%s\t\t%s\t\t%s\t\t\t%s\n', ' ' ,'d.f.', 'SS', 'MS', 'F', 'P-Value');
  # fprintf('--------------------------------------------------------------------\n');
  # fprintf('%s\t\t\t%u\t\t%.2f\t%.2f\t%.2f\t\t%.4f\n', 'Columns', s-1 , A, A/(s-1), F, pval);
  # fprintf('%s\t\t%u\t\t%.2f\t%.2f\n', 'Residual ', n-s, B, B/(n-s));
  # fprintf('--------------------------------------------------------------------\n');
  # fprintf('%s\t\t%u\t\t%.2f', 'Total   ',n-1,A+B);
  # fprintf('\n\n')

  # table = ["Source",'d.f.','SS','MS','F','P-Value';
  #          'Columns', s-1 , A, A/(s-1), F, pval;
  #          'Residual ', n-s, B, B/(n-s), [], [];
  #          'Total',n-1,A+B,[],[],[]]

  return (; p, F)
end

"""
Parametric two-way ANOVA for circular data with interations.

  !!!note
  The test assumes underlying von-Mises distributrions.
  All groups are assumed to have a common concentration parameter k, between 0 and 2.

  1. α: angles in radians
  2. idp: indicates the level of factor 1 (1:p)
  3. idq: indicates the level of factor 2 (1:q)
    - inter: whether to include effect of interaction
    - fn: string array containing names of the factors

  return:
  - p: pvalues for factors and interaction
  - s: statistic of each p value
"""
function circ_hktest(α, idp, idq; inter = true, fn = ['A', 'B'])
  # number of groups for every factor
  pu = unique(idp)
  p = length(pu)
  qu = unique(idq)
  q = length(qu)

  # number of samples
  n = length(α)

  # compute important sums for the test statistics
  cn = zeros(p, q)
  cr = zeros(p, q)
  pm = zeros(p)
  pr = zeros(p)
  pn = zeros(p)
  qm = zeros(q)
  qr = zeros(q)
  qn = zeros(q)
  for pp = 1:p
    p_id = idp .== pu[pp] # indices of factor1 = pp
    for qq = 1:q
      q_id = idq .== qu[qq] # indices of factor2 = qq
      idx = p_id .& q_id
      cn[pp, qq] = sum(idx)     # number of items in cell
      cr[pp, qq] = cn[pp, qq] * circ_r(α[idx]) # R of cell
    end
    # R and mean angle for factor 1
    pr[pp] = sum(p_id) * circ_r(α[p_id])
    pm[pp] = circ_mean(α[p_id]).μ
    pn[pp] = sum(p_id)
  end

  # R and mean angle for factor 2
  for qq = 1:q
    q_id = idq .== qu[qq]
    qr[qq] = sum(q_id) * circ_r(α[q_id])
    qm[qq] = circ_mean(α[q_id]).μ
    qn[qq] = sum(q_id)
  end

  # R and mean angle for whole sample (total)
  tr = n * circ_r(α)

  # estimate kappa
  kk = circ_kappa(tr / n)

  # different formulas for different width of the distribution
  if kk > 2
    # large kappa

    # effect of factor 1
    eff_1 = sum(pr .^ 2 ./ dropdims(sum(cn, dims = 2), dims = 2)) - tr .^ 2 / n
    df_1 = p - 1
    ms_1 = eff_1 / df_1

    # effect of factor 2
    eff_2 = sum(qr .^ 2 ./ dropdims(sum(cn, dims = 1), dims = 1)) - tr .^ 2 / n
    df_2 = q - 1
    ms_2 = eff_2 / df_2

    # total effect
    eff_t = n - tr .^ 2 / n
    df_t = n - 1

    m = mean(cn)

    if inter
      # correction factor for improved F statistic
      β = 1 / (1 - 1 / (5kk) - 1 / (10kk^2))

      # residual effects
      eff_r = n - sum(cr .^ 2 ./ cn)
      df_r = p * q * (m - 1)
      ms_r = eff_r / df_r

      # interaction effects
      eff_i =
        sum(cr .^ 2 ./ cn) - sum(qr .^ 2 ./ qn) - sum(pr .^ 2 ./ pn) +
        tr .^ 2 / n
      df_i = (p - 1) * (q - 1)
      ms_i = eff_i / df_i

      # interaction test statistic
      FI = ms_i / ms_r
      pI = 1 - cdf(FDist(df_i, df_r), FI)

    else
      # residual effect
      eff_r = n - sum(qr .^ 2 ./ qn) - sum(pr .^ 2 ./ pn) + tr .^ 2 / n
      df_r = (p - 1) * (q - 1)
      ms_r = eff_r / df_r

      # interaction effects
      eff_i = []
      df_i = []
      ms_i = []

      # interaction test statistic
      FI = []
      pI = NaN
      β = 1
    end

    # compute all test statistics as F = β * MS(A) / MS(R)
    F1 = β * ms_1 / ms_r
    p1 = 1 - cdf(FDist(df_1, df_r), F1)

    F2 = β * ms_2 / ms_r
    p2 = 1 - cdf(FDist(df_2, df_r), F2)
    s = [F1, F2, FI]
  else
    # small kappa

    # correction factor
    rr = besseli(1, kk) / besseli(0, kk)
    f = 2 / (1 - rr^2)

    chi1 = f * (sum(pr .^ 2 ./ pn) - tr .^ 2 / n)
    df_1 = 2 * (p - 1)
    p1 = 1 - cdf(Chisq(df_1), chi1)

    chi2 = f * (sum(qr .^ 2 ./ qn) - tr .^ 2 / n)
    df_2 = 2 * (q - 1)
    p2 = 1 - cdf(Chisq(df_2), chi2)

    chiI =
      f * (
        sum(cr .^ 2 ./ cn) - sum(pr .^ 2 ./ pn) - sum(qr .^ 2 ./ qn) +
        tr .^ 2 / n
      )
    df_i = (p - 1) * (q - 1)
    pI = 1 - cdf(Chisq(df_i), chiI)
    s = [chi1, chi2, chiI]
  end
  p = [p1, p2, pI]
  #
  # if kk>2
  #
  #   fprintf('\nANALYSIS OF VARIANCE TABLE (HIGH KAPPA MODE)\n\n');
  #
  #   fprintf('%s\t\t\t\t%s\t%s\t\t%s\t\t%s\t\t\t%s\n', ' ' ,'d.f.', 'SS', 'MS', 'F', 'P-Value');
  #   fprintf('--------------------------------------------------------------------\n');
  #   fprintf('%s\t\t\t\t%u\t\t%.2f\t%.2f\t%.2f\t\t%.4f\n', fn{1}, df_1 , eff_1, ms_1, F1, p1);
  #   fprintf('%s\t\t\t\t%u\t\t%.2f\t%.2f\t%.2f\t\t%.4f\n', fn{2}, df_2 , eff_2, ms_2, F2, p2);
  #   if (inter)
  #       fprintf('%s\t\t%u\t\t%.2f\t%.2f\t%.2f\t\t%.4f\n', 'Interaction', df_i , eff_i, ms_i, FI, pI);
  #   end
  #   fprintf('%s\t\t%u\t\t%.2f\t%.2f\n', 'Residual ', df_r, eff_r, ms_r);
  #   fprintf('--------------------------------------------------------------------\n');
  #   fprintf('%s\t\t%u\t\t%.2f', 'Total   ',df_t,eff_t);
  #   fprintf('\n\n')
  # else
  #   fprintf('\nANALYSIS OF VARIANCE TABLE (LOW KAPPA MODE)\n\n');
  #
  #   fprintf('%s\t\t\t\t%s\t%s\t\t\t%s\n', ' ' ,'d.f.', 'CHI2', 'P-Value');
  #   fprintf('--------------------------------------------------------------------\n');
  #   fprintf('%s\t\t\t\t%u\t\t%.2f\t\t\t%.4f\n', fn{1}, df_1 , chi1, p1);
  #   fprintf('%s\t\t\t\t%u\t\t%.2f\t\t\t%.4f\n', fn{2}, df_2 , chi2, p2);
  #   if (inter)
  #       fprintf('%s\t\t%u\t\t%.2f\t\t\t%.4f\n', 'Interaction', df_i , chiI, pI);
  #   end
  #   fprintf('--------------------------------------------------------------------\n');
  #   fprintf('\n\n')
  #
  # end
  return (; p, s)
end

"""
A parametric two-sample test to determine whether two concentration parameters are different.

  H0: The two concentration parameters are equal.
  HA: The two concentration parameters are different.

  Assumptions: both samples are drawn from von Mises type distributions and their joint resultant vector length should be > .7

  1. α₁: sample of angles in radians
  2. α₂: sample of angles in radians

  return:
  - p: p-value that samples have different concentrations
  - f: f-statistic calculated
"""
function circ_ktest(α₁, α₂)
  n1 = length(α₁)
  n2 = length(α₂)

  R1 = n1 * circ_r(α₁)
  R2 = n2 * circ_r(α₂)

  # make sure that r̄ > .7
  r̄ = (R1 + R2) / (n1 + n2)

  if r̄ < 0.7
    @warn "Resultant vector length should be > 0.7"
  end

  # calculate test statistic
  f = ((n2 - 1) * (n1 - R1)) / ((n1 - 1) * (n2 - R2))
  if f > 1
    p = 2 * (1 - cdf(FDist(n1, n2), f))
  else
    f = 1 / f
    p = 2 * (1 - cdf(FDist(n2, n1), f))
  end
  return (; p, f)
end

"""
Tests for symmetry about the median.
  H0: the population is symmetrical around the median
  HA: the population is not symmetrical around the median

  1. α: sample of angles in radians

  return: p-value
"""
function circ_symtest(α)
  # compute median
  md = circ_median(α)

  # compute deviations from median
  d = circ_dist(α, md)

  # compute wilcoxon sign rank test
  pvalue(SignedRankTest(d))
end

const kuipertable= [5	1.45800000000000	1.56500000000000	1.68200000000000	1.76300000000000	1.83800000000000	1.92000000000000	1.97000000000000;
                    6	1.47100000000000	1.58200000000000	1.71100000000000	1.79300000000000	1.86700000000000	1.95700000000000	2.02000000000000;
                    7	1.48300000000000	1.59800000000000	1.72700000000000	1.81400000000000	1.89400000000000	1.98700000000000	2.05100000000000;
                    8	1.49300000000000	1.60800000000000	1.74100000000000	1.83000000000000	1.91100000000000	2.00900000000000	2.07700000000000;
                    9	1.50000000000000	1.61800000000000	1.75200000000000	1.84300000000000	1.92600000000000	2.02700000000000	2.09700000000000;
                    10	1.50700000000000	1.62500000000000	1.76100000000000	1.85400000000000	1.93800000000000	2.04100000000000	2.11300000000000;
                    11	1.51300000000000	1.63100000000000	1.76900000000000	1.86200000000000	1.94800000000000	2.05300000000000	2.12600000000000;
                    12	1.51700000000000	1.63700000000000	1.77600000000000	1.87000000000000	1.95700000000000	2.06200000000000	2.13700000000000;
                    13	1.52200000000000	1.64200000000000	1.78200000000000	1.87600000000000	1.96400000000000	2.07100000000000	2.15400000000000;
                    14	1.52500000000000	1.64600000000000	1.78700000000000	1.88200000000000	1.97000000000000	2.07800000000000	2.15400000000000;
                    15	1.52900000000000	1.65000000000000	1.79100000000000	1.88700000000000	1.97600000000000	2.08500000000000	2.16100000000000;
                    16	1.53200000000000	1.65300000000000	1.79500000000000	1.89200000000000	1.98100000000000	2.09000000000000	2.16800000000000;
                    17	1.53400000000000	1.65700000000000	1.79900000000000	1.89600000000000	1.98600000000000	2.09600000000000	2.17300000000000;
                    18	1.53700000000000	1.65900000000000	1.80200000000000	1.89900000000000	1.99000000000000	2.10000000000000	2.17800000000000;
                    19	1.53900000000000	1.66200000000000	1.80500000000000	1.90300000000000	1.99300000000000	2.10400000000000	2.18300000000000;
                    20	1.54100000000000	1.66400000000000	1.80800000000000	1.90600000000000	1.99700000000000	2.10800000000000	2.18700000000000;
                    21	1.54300000000000	1.66700000000000	1.81000000000000	1.90800000000000	2	                2.11200000000000	2.19100000000000;
                    22	1.54500000000000	1.66900000000000	1.81300000000000	1.91100000000000	2.00300000000000	2.11500000000000	2.19400000000000;
                    23	1.54700000000000	1.67000000000000	1.81500000000000	1.91300000000000	2.00500000000000	2.11800000000000	2.19800000000000;
                    24	1.54900000000000	1.67200000000000	1.81700000000000	1.91600000000000	2.00800000000000	2.12100000000000	2.20100000000000;
                    25	1.55000000000000	1.67400000000000	1.81900000000000	1.91800000000000	2.01000000000000	2.12300000000000	2.20300000000000;
                    30	1.55600000000000	1.68100000000000	1.82600000000000	1.92600000000000	2.01900000000000	2.13400000000000	2.21500000000000;
                    35	1.56100000000000	1.68600000000000	1.83200000000000	1.93300000000000	2.02600000000000	2.14100000000000	2.22300000000000;
                    40	1.56500000000000	1.69000000000000	1.83700000000000	1.93800000000000	2.03200000000000	2.14800000000000	2.23000000000000;
                    45	1.56800000000000	1.69400000000000	1.84100000000000	1.94200000000000	2.03600000000000	2.15200000000000	2.23500000000000;
                    50	1.57100000000000	1.69700000000000	1.84400000000000	1.94600000000000	2.04000000000000	2.15700000000000	2.23900000000000;
                    100	1.58800000000000	1.71400000000000	1.86200000000000	1.96500000000000	2.06000000000000	2.17800000000000	2.26200000000000;
                    200	1.60000000000000	1.72600000000000	1.87600000000000	1.97900000000000	2.07500000000000	2.19400000000000	2.27900000000000;
                    500	1.61000000000000	1.73700000000000	1.88700000000000	1.99000000000000	2.08700000000000	2.20700000000000	2.29200000000000;
                    501	1.62000000000000	1.74700000000000	1.89800000000000	2.00100000000000	2.09800000000000	2.21800000000000	2.30300000000000]

"""
The Kuiper two-sample test tests whether the two samples differ significantly.
The difference can be in any property, such as mean location and dispersion.
It is a circular analogue of the Kolmogorov-Smirnov test.

  H0: The two distributions are identical.
  HA: The two distributions are different.

  1. α₁: sample of angles in radians
  2. α₂: sample of angles in radians
    - n: resolution at which the cdf are evaluated

  return:
  - p: p-value; the smallest of .10, .05, .02, .01, .005, .002, .001,
        for which the test statistic is still higher than the respective critical value.
        this is due to the use of tabulated values. if p>.1, pval is set to 1.
  - k: test statistic
  - K: critical value
"""
function circ_kuipertest(α₁, α₂; n = 100)
  n = length(α₁)
  m = length(α₂)

  # create cdfs of both samples
  ϕ₁, cdf1 = circ_samplecdf(α₁; n)
  _, cdf2 = circ_samplecdf(α₂; n)

  # maximal difference between sample cdfs
  dplus, gdpi = findmax([0; cdf1 .- cdf2])
  dminus, gdmi = findmax([0; cdf2 .- cdf1])

  # calculate k-statistic
  k = n * m * (dplus + dminus)

  # find p-value
  p, K = kuiperlookup(min(n, m), k / sqrt(n * m * (n + m)))
  K *= sqrt(n * m * (n + m))

  return (; p, k, K)
end

function kuiperlookup(n, k)
  α = [0.10, 0.05, 0.02, 0.01, 0.005, 0.002, 0.001]
  nn = kuipertable[:, 1]

  # find correct row of the table
  row = findfirst(n.==nn)
  if isnothing(row)
    # find closest value if no entry is present
    row = length(nn) - sum(n .< nn)
    if row == 0
      error("n too small.")
    else
      @warn "n=$n not found in table, using closest n=$(nn[row]) present."
    end
  end

  # find minimal p-value and test-statistic
  idx = findlast(kuipertable[row, 2:end] .< k)
  if isnothing(idx)
    p = 1
    K = NaN
  else
    p = α[idx]
    K = kuipertable[row, idx+1]
  end

  return (; p, K)
end

"""
Circular correlation coefficient for two circular random variables.

  1. α₁: sample of angles in radians
  2. α₂: sample of angles in radians

  return:
  - ρ: correlation coefficient
  - p: p-value
"""
function circ_corrcc(α₁, α₂)
  # compute mean directions
  n = length(α₁)
  ᾱ₁ = circ_mean(α₁).μ
  ᾱ₂ = circ_mean(α₂).μ

  # compute correlation coeffcient from p. 176
  num = sum(sin.(α₁ .- ᾱ₁) .* sin.(α₂ .- ᾱ₂))
  den = sqrt(sum(sin.(α₁ .- ᾱ₁) .^ 2) .* sum(sin.(α₂ .- ᾱ₂) .^ 2))
  ρ = num / den

  # compute pvalue
  l20 = mean(sin.(α₁ .- ᾱ₁) .^ 2)
  l02 = mean(sin.(α₂ .- ᾱ₂) .^ 2)
  l22 = mean((sin.(α₁ .- ᾱ₁) .^ 2) .* (sin.(α₂ .- ᾱ₂) .^ 2))

  ts = sqrt((n * l20 * l02) / l22) * ρ
  p = 2 * (1 - cdf(Normal(), abs(ts)))
  return (; ρ, p)
end

"""
Correlation coefficient between one circular and one linear random variable.

  1. α: sample of angles in radians
  2. x: sample of linear random variable

  return:
  - ρ: correlation coefficient
  - p: p-value
"""
function circ_corrcl(α, x)
  n = length(α)

  # compute correlation coefficent for sin and cos independently
  rxs = cor(x, sin.(α))
  rxc = cor(x, cos.(α))
  rcs = cor(sin.(α), cos.(α))

  # compute angular-linear correlation (equ. 27.47)
  ρ = sqrt((rxc^2 + rxs^2 - 2 * rxc * rxs * rcs) / (1 - rcs^2))

  # compute pvalue
  p = 1 - cdf(Chisq(2), n * ρ^2)
  return (; ρ, p)
end

"""
Computes an approximation to the ML estimate of the concentration parameter kappa of the von Mises distribution.

  1. α: angles in radians OR α is length resultant
    - w: number of incidences in case of binned angle data

  return: estimated value of kappa
"""
function circ_kappa(α; w = ones(size(α)))
  N = length(α)

  if N > 1
    R = circ_r(α; w)
  else
    R = α
  end

  if R < 0.53
    κ = 2R + R^3 + 5R^5 / 6
  elseif R >= 0.53 && R < 0.85
    κ = -0.4 + 1.39R + 0.43 / (1 - R)
  else
    κ = 1 / (R^3 - 4R^2 + 3R)
  end

  if N < 15 && N > 1
    if κ < 2
      κ = max(κ - 2 * (N * κ)^-1, 0)
    else
      κ = (N - 1)^3 * κ / (N^3 + N)
    end
  end
  κ
end

"""
Performs a simple agglomerative clustering of angular data.

1. α: sample of angles in radians
  - k: number of clusters desired, default=2

return:
- cid: cluster id for each entry of α
- α: sorted angles, matched with cid
- μ: mean direction of angles in each cluster
"""
function circ_clust(α; k = 2)

  n = length(α)
  n < k && error("Not enough data for clusters.")

  # prepare data
  cid = 1:n

  # main clustering loop
  num_unique = length(unique(cid))

  while (num_unique > k)

    # find centroid means...

    # calculate the means for each putative cluster
    μ = fill(NaN, n)
    for j = 1:n
      if sum(cid .== j) > 0
        μ[j], = circ_mean(α(cid .== j))
      end
    end

    # find distance between centroids...
    μdist = abs.(circ_dist2(μ))

    # find closest pair of clusters/datapoints
    mindist = minimum(μdist[tril(ones(size(μdist), size(μdist)), -1)==1])
    row, col = findall(μdist .== mindist)

    # update cluster id's
    cid[cid.==maximum(row)] = minimum(col)

    # update stop criteria
    num_unique = length(unique(cid))

  end

  # renumber cluster ids (so cids [1 3 7 10] => [1 2 3 4])
  cid2 = deepcopy(cid)
  uniquecids = unique(cid)
  for j = 1:length(uniquecids)
    cid[cid2.==uniquecids[j]] = j
  end

  # compute final cluster means
  μ = fill(NaN, num_unique)
  for j = 1:num_unique
    if sum(cid .== j) > 0
      μ[j], = circ_mean(α[cid.==j]')
    end
  end

  return (; cid, α, μ)
end

"""
Helper function for circ_kuipertest. Evaluates CDF of sample.

  1. α: sample of angles in radians
    - n: resolution at which the cdf are evaluated

  return:
  - ϕ: angles at which CDF are evaluated
  - c: CDF values at ϕ
"""
function circ_samplecdf(α; n = 100)
  ϕ = range(0, 2π, length = n + 1)

  # ensure all points in α are on interval [0, 2pi)
  α = sort(mod2pi.(α))

  dp = 1 / length(α) # incremental change in probability
  c = cumsum(map((l, r) -> dp * sum(l .<= α .< r), ϕ[1:end-1], ϕ[2:end]))
  ϕ = ϕ[1:end-1]

  return (; ϕ, c)
end


text:
#
# DON'T USE THIS (OR FIX IT FOR ME..)
# This is a broken version of reweighter that uses
# julia's standard NLsolve package. I'll return to it when
# I have a free day.
#
using NLsolve

Buffer = Dict{Symbol,Any}

"
   Implements the micro data weighting procedures from:
   * Creedy 2003 http://www.treasury.govt.nz/publications/research-policy/wp/2003/03-17/twp03-17.pdf
   * Creedy 2003  http://www.business.curtin.edu.au/files/creedy2.pdf
   * Jean-Claude Deville and Carl-Erik Sarndal http://www.jstor.org/stable/2290268
"
@enum DistanceFunctionType chi_square d_and_s_type_a d_and_s_type_b constrained_chi_square d_and_s_constrained

"
Make a weights vector which weights the matrix `data`
so when summed the col totals sum to `target_populations`
See the Creedy Paper for `function_type`
If using one of the constrained types,
the output weights should be no more than ru*the initial weight,
no less than rl
Returns a Dict with :=>weights and some extra info on convergence.
data : KxJ matrix where k is num observations and J is num constraints;
see:
Microdata Adjustment by the Minimum Information Loss Principle Joachim Merz; FFB Discussion Paper No. 10 July 1994
for a good discussion on how to lay out the dataset

intial_weights, new_weights : K length vector
target_populations - J length vector;

tolx, tolf, max_iterations : see Solve_Non_Linear_Equation_System in the parent
ru/rl max/min acceptable values of ratio of final_weight/initial_weight (for constrained distance functions)

note: chi-square is just there for checking purposes; use Do_Basic_Reweighting if that's all you need.

"
function doreweighting2(
    ;
    data               :: AbstractArray{ <:Real, 2 },
    initial_weights    :: AbstractArray{ <:Real, 1 }, # a column
    target_populations :: AbstractArray{ <:Real, 1 }, # a row
    functiontype       :: DistanceFunctionType,
    upper_multiple     :: Real = 0.0,
    lower_multiple     :: Real = 0.0,
    tolx               :: Real = 0.000001,
    tolf               :: Real = 0.000001 ) :: Buffer

    nrows, ncols = size( data )
    # ncols = size( data )[2]
    @assert ncols == size( target_populations )[1]
    @assert nrows == size( initial_weights )[1]
    a = target_populations - (initial_weights'*data)'
    ru = upper_multiple # shorthand
    rl = lower_multiple

    #
    # instead of minimising the loss function (G) we just need to find the zeros
    # of the 1st derivatives. This computes the derivatives and the J
    # of the G function, hence the slightly confusing names (F=>1st deriv of G)
    # See table 3 and eqns 25 and 26 in Creedy 2003 http://ideas.repec.org/p/nzt/nztwps/03-17.html
    # This version is in the 'fj' form needed for the NLsolve package if we want to
    # generate the derivative and J in one go, as we do since it's a potentially huge loop
    #
    function computefandj!(
        F       :: Union{Vector,Nothing},
        J       :: Union{Matrix,Nothing},
        lamdas  :: Vector )
        if J !== nothing
            J = zeros( ncols, ncols )
        end
        z = zeros( ncols )
        for row in 1:nrows
            rv = data[row,:]
            u = (rv' * lamdas)[1]
            println( "u=$u lamdas = $lamdas")
            d_g_m1 = 0.0
            g_m1 = 0.0
            if functiontype == chi_square
                d_g_m1 = 1.0;
                g_m1 = 1.0 + u;
            elseif functiontype == constrained_chi_square
                if( u < ( rl - 1.0 ))
                   g_m1 = rl
                   d_g_m1 = 0.0
                elseif( u > ( ru - 1.0 ))
                   g_m1 = ru
                   d_g_m1 = 0.0
                else
                   g_m1 = 1.0 + u
                   d_g_m1 = 1.0
                end
            elseif functiontype == d_and_s_type_a
               g_m1 = ( 1.0 -  u/2.0 ) ^ ( -2 )
               d_g_m1 = ( 1.0 - u/2.0 ) ^ ( -3 )
            elseif functiontype == d_and_s_type_b
               g_m1 = ( 1.0- u ) ^ (-1 )
               d_g_m1 = ( 1.0 - u ) ^ ( -2 )
            elseif functiontype == d_and_s_constrained
               alpha = ( ru - rl ) / (( 1.0 - rl )*( ru - 1.0 ))
               g_m1 = rl*(ru-1.0)+ru*(1.0-rl)*exp( alpha*u )/((ru-1.0)+(1.0-rl)*(exp( alpha*u )))
               d_g_m1 = g_m1 * ( ru - g_m1 ) *
                 ((( 1.0 - rl )*alpha*exp( alpha*u )) /
                  (( ru - 1.0 ) + (( 1.0 - rl ) * exp( alpha*u ))))
           end # function cases
           for col in 1:ncols
               z[col] += initial_weights[row]*data[row,col]*(g_m1-1.0)
               ## the J
               if J != nothing
                   # println( "HESS CALL")
                   for c2 in 1:ncols
                       zz = initial_weights[row]*data[row,col]*data[row,c2]
                       J[col,c2] += zz*d_g_m1
                   end
               end
           end
        end # obs loop
        if F !== nothing
            F = a - z
        end
        println( "F $F")
        println( "J $J" )
        println( "lamdas $lamdas" )
    end # nested function

    J = zeros( ncols, ncols )
    F = zeros( ncols )
    initial_x = zeros( ncols )
    # initial_x = [0.214763; -0.296946; -0.0128647; 0.078507]
    objective = only_fj!( computefandj! ) # this means: we're providing one function 'fj' which
                                   # computes both the first derivatives value and the J
    rc = nlsolve(
        only_fj!( computefandj! ),
        F,
        J,
        initial_x,
        inplace        = false,
        method         = :newton,
        autoscale      = true,
        extended_trace = true )
        # xtol           = tolx,
        # ftol           = tolf )
    new_weights = copy(initial_weights)
    x = ( rc.zero )

    # construct the new weights from the x
    if converged( rc )
        for r in 1:nrows
            row = data[r,:]
            u = (row'*x)[1]
            g_m1 = 0.0
            if functiontype == chi_square
                g_m1 = 1.0 + u;
            elseif functiontype == constrained_chi_square
                if( u < ( rl - 1.0 ))
                   g_m1 = rl
                elsif( u > ( ru - 1.0 ))
                   g_m1 = ru
                else
                   g_m1 = 1.0 + u
                end
            elseif functiontype == d_and_s_type_a
               g_m1 = ( 1.0 -  u/2.0 ) ^ ( -2 )
            elseif functiontype == d_and_s_type_b
               g_m1 = ( 1.0- u ) ^ (-1 )
            elseif functiontype == d_and_s_constrained
               alpha = ( ru - rl ) / (( 1.0 - rl )*( ru - 1.0 ))
               g_m1 = rl*(ru-1.0)+ru*(1.0-rl)*exp( alpha*u )/((ru-1.0)+(1.0-rl)*(exp( alpha*u )))
           end # function cases
            #
            # Creedy wp 03/17 table 3
            #
            new_weights[r] = initial_weights[r]*g_m1
        end
    end # converged
    return Buffer(:x=>x, :rc => rc, :weights=>new_weights ) # , :converged => converge )
end # do reweighting

"
This is a route-1 approach to Chi-square reweighting.
The iterative main method should produce identical results
when method=chi_square. This is kept here mainly for testing.
Note the weights can be negative.
See the Creedy Papers.
"
function dochisquarereweighting(
    data               :: AbstractArray{ <:Real, 2 },
    initial_weights    :: AbstractArray{ <:Real, 1 }, # a row
    target_populations :: AbstractArray{ <:Real, 1 } ) :: Array{ <:Real }

    nrows = size( data )[1]
    ncols = size( data )[2]

    row = zeros( ncols )
    populations = zeros( ncols )
    x  = zeros( ncols )
    weights = zeros( nrows )
    m = zeros( ncols, ncols )
    for r in 1:nrows
        row = data[r,:]
        m += initial_weights[r]*(row*row')
        for c in 1:ncols
            populations[c] += (row[c]*initial_weights[r])'
        end
    end
    x = (m^-1)*(target_populations-populations)
    for r in 1:nrows
        row = data[r,:]
        weights[r] = initial_weights[r]*(1.0 + (row'*x)[1])
    end
    return weights;
end


text:
using FourierCore, FourierCore.Grid, FourierCore.Domain
using FFTW, LinearAlgebra, BenchmarkTools, Random, JLD2
using GLMakie
rng = MersenneTwister(1234)
# Random.seed!(123456789)
Random.seed!(12)
# jld_name = "high_order_timestep_spatial_tracer_"
jld_name = "blocky"
include("transform.jl")
include("random_phase_kernel.jl")
# using GLMakie
using CUDA
arraytype = CuArray
Ω = S¹(4π)^2
N = 2^7 # number of gridpoints
grid = FourierGrid(N, Ω, arraytype=arraytype)
nodes, wavenumbers = grid.nodes, grid.wavenumbers

x = nodes[1]
y = nodes[2]
kˣ = wavenumbers[1]
kʸ = wavenumbers[2]
# construct filter
kxmax = maximum(kˣ)
kymax = maximum(kˣ)
filter = @. (kˣ)^2 + (kʸ)^2 ≤ ((kxmax / 2)^2 + (kymax / 2)^2)

# now define the random field
phase_speed = 0.2 # default is 1.0
wavemax = 3
𝓀 = arraytype(collect(-wavemax:0.5:wavemax))
𝓀ˣ = reshape(𝓀, (length(𝓀), 1))
𝓀ʸ = reshape(𝓀, (1, length(𝓀)))
inertial_exponent = -3.0
stream_function_exponent = (inertial_exponent - 1) / 4;
A = @. 0.2 * (𝓀ˣ * 𝓀ˣ + 𝓀ʸ * 𝓀ʸ)^(stream_function_exponent)
A[A.==Inf] .= 0.0
φ = arraytype(2π * rand(size(A)...))
field = arraytype(zeros(N, N))

##
# Fields
# velocity
ψ = arraytype(zeros(ComplexF64, N, N))
u = similar(ψ)
v = similar(ψ)

# theta
θ = similar(ψ)
∂ˣθ = similar(ψ)
∂ʸθ = similar(ψ)
κΔθ = similar(ψ)
θ̇ = similar(ψ)
s = similar(ψ)
θ̅ = similar(ψ)
k₁ = similar(ψ)
k₂ = similar(ψ)
k₃ = similar(ψ)
k₄ = similar(ψ)
θ̃ = similar(ψ)

# source
s = similar(ψ)
@. s = cos(kˣ[5] * x)

# phase
φ̇ = similar(A)

# operators
∂x = im * kˣ
∂y = im * kʸ
Δ = @. ∂x^2 + ∂y^2

# plan ffts
P = plan_fft!(ψ)
P⁻¹ = plan_ifft!(ψ)

# number of gridpoints in transition is about λ * N / 2
bump(x; λ=10 / N, width=π / 2) = 0.5 * (tanh((x + width / 2) / λ) - tanh((x - width / 2) / λ))
bumps(x; λ=20 / N, width=1.0) = 0.25 * (bump(x, λ=λ, width=width) + bump(x, λ=λ, width=2.0 * width) + bump(x, λ=λ, width=3.0 * width) + bump(x, λ=λ, width=4.0 * width))

##
Δx = x[2] - x[1]
κ = 1.0 * Δx^2 # 1.0 / N * 2^(2)  # roughly 1/N for this flow
# κ = 2 / 2^8 # fixed diffusivity
# κ = 2e-4
Δt = (Δx) / (4π) * 5

# take the initial condition as negative of the source
tic = Base.time()

# save some snapshots
ψ_save = typeof(real.(Array(ψ)))[]
θ_save = typeof(real.(Array(ψ)))[]

r_A = Array(@. sqrt((x - 2π)^2 + (y - 2π)^2))
θ_A = [bumps(r_A[i, j]) for i in 1:N, j in 1:N]
θ .= CuArray(θ_A)
# @. θ = bump(sqrt(x^2 + y^2)) # scaling so that source is order 1
θclims = extrema(Array(real.(θ))[:])
P * θ # in place fft
@. κΔθ = κ * Δ * θ
P⁻¹ * κΔθ # in place fft
s .= -κΔθ * 0.0
P⁻¹ * θ # in place fft
θ̅ .= 0.0

t = [0.0]
tend = 100 # 5000

iend = ceil(Int, tend / Δt)

tlist = collect(0:tend)
indlist = [ceil(Int, t / Δt) for t in tlist]
indlist[1] = 1
lagrangian_list = zeros(length(indlist))
eulerian_list = copy(lagrangian_list)

params = (; ψ, A, 𝓀ˣ, 𝓀ʸ, x, y, φ, u, v, ∂ˣθ, ∂ʸθ, s, P, P⁻¹, filter)

size_of_A = size(A)

realizations = 100
tmpA = []
for j in 1:realizations

    # new realization of flow
    rand!(rng, φ)
    φ .*= 2π
    event = stream_function!(ψ, A, 𝓀ˣ, 𝓀ʸ, x, y, φ)
    wait(event)

    # get horizontal velocity
    P * ψ
    u0 = ∂y .* ψ
    P⁻¹ * ψ
    P⁻¹ * u0

    # initialize tracer
    θ .= u0
    if (j == 1) | (j == realizations)
        push!(tmpA, real.(Array(u)))
    end

    t[1] = 0.0

    ii = 1
    for i = 1:iend
        # fourth order runge kutta on deterministic part
        # keep ψ frozen is the correct way to do it here

        # the below assumes that φ is just a function of time
        θ_rhs_new!(k₁, θ, params)
        @. θ̃ = θ + Δt * k₁ * 0.5

        φ_rhs_normal!(φ̇, φ, rng)
        @. φ += phase_speed * sqrt(Δt / 2) * φ̇

        θ_rhs_new!(k₂, θ̃, params)
        @. θ̃ = θ + Δt * k₂ * 0.5
        θ_rhs_new!(k₃, θ̃, params)
        @. θ̃ = θ + Δt * k₃

        φ_rhs_normal!(φ̇, φ, rng)
        @. φ += phase_speed * sqrt(Δt / 2) * φ̇

        θ_rhs_new!(k₄, θ̃, params)
        @. θ += Δt / 6 * (k₁ + 2 * k₂ + 2 * k₃ + k₄)

        t[1] += Δt

        if i in indlist
            # get horizontal velocity
            P * ψ
            @. u = ∂y * ψ
            P⁻¹ * ψ
            P⁻¹ * u

            lagrangian_list[ii] += real(mean(θ .* u0)) / realizations
            eulerian_list[ii] += real(mean(u .* u0)) / realizations

            ii += 1
            # println("saveing at ", i)
        end

    end
    println("finished realization ", j)
    @. θ̅ += θ / realizations
end

toc = Base.time()
println("the time for the simulation was ", toc - tic, " seconds")

fig = Figure()
ax = Axis(fig[1,1]; xlabel = "log10(time)", ylabel = "autocorrelation", xlabelsize =30, ylabelsize =30)

logtlist = log10.(tlist .+1)
ln1 = lines!(ax, logtlist, lagrangian_list, color=:blue, label = "Lagrangian")
ln2 = lines!(ax, logtlist, eulerian_list, color=:orange, label = "Eulerian")
axislegend(ax, position=:rc)
display(fig)

#=
x_A = Array(x)[:] .- 2π
θ_A = tmpA[1]
θ_F = tmpA[2]
θ̅_F = Array(real.(θ̅))

fig = Figure(resolution=(2048, 512))
ax1 = Axis(fig[1, 1], title="t = 0")
ax2 = Axis(fig[1, 2], title="instantaneous t = " * string(tend))
ax3 = Axis(fig[1, 4], title="ensemble average t = " * string(tend))
println("the extrema of the end field is ", extrema(θ_F))
println("the extrema of the ensemble average is ", extrema(θ̅_F))
colormap = :bone_1
colormap = :nipy_spectral
heatmap!(ax1, x_A, x_A, θ_A, colormap=colormap, colorrange=(0.0, 1.0), interpolate=true)
hm = heatmap!(ax2, x_A, x_A, θ_F, colormap=colormap, colorrange=(0.0, 1.0), interpolate=true)
hm_e = heatmap!(ax3, x_A, x_A, θ̅_F, colormap=colormap, colorrange=(0.0, 0.2), interpolate=true)
Colorbar(fig[1, 3], hm, height=Relative(3 / 4), width=25, ticklabelsize=30, labelsize=30, ticksize=25, tickalign=1,)
Colorbar(fig[1, 5], hm_e, height=Relative(3 / 4), width=25, ticklabelsize=30, labelsize=30, ticksize=25, tickalign=1,)
display(fig)
=#

max_stars_repo_name: sandreza/FourierCore

max_issues_repo_name: sandreza/FourierCore

max_forks_repo_name: sandreza/FourierCore