text:
abstract type AbstractModelSet end

"""
    WindFarmModelSet(wakedeficitmodel, wake_deflection_model, wake_combination_model, local_ti_model)

Container for objects defining models to use in wind farm calculations

# Arguments
- `wake_defiict_model::AbstractWakeDeficitModel`: contains a struct defining the desired wake deficit model
- `wake_deflection_model::AbstractWakeDeflectionModel`: contains a struct defining the desired wake deflection model
- `wake_combination_model::AbstractWakeCombinationModel`: contains a struct defining the desired wake combination model
- `local_ti_model::AbstractTurbulenceIntensityModel`: contains a struct defining the desired turbulence intensity model
"""
struct WindFarmModelSet{DTM,DNM,CM,TIM} <: AbstractModelSet

    wake_deficit_model::DTM
    wake_deflection_model::DNM
    wake_combination_model::CM
    local_ti_model::TIM

end

"""
    rotate_to_wind_direction(xlocs, ylocs, wind_direction_met)

Rotates wind farm coordinates to be in wind direction reference where wind direction is to
the positive x.

# Arguments
- `xlocs::Array`: contains turbine east-west locations in the global reference frame
- `ylocs::Array`: contains turbine north-south locations in the global reference frame
- `wind_direction_met::Array`: contains wind direction in radians in meteorological standard
    system (N=0 rad, proceeds CW, wind from direction given)
"""
function rotate_to_wind_direction(xlocs, ylocs, wind_direction_met)
    # use radians

    # convert from meteorological polar system (CW, 0 rad.=N) to standard polar system (CCW, 0 rad.=E)
    wind_direction_cart = (3.0*pi/2.0 - wind_direction_met)

    if wind_direction_cart < 0.0
        wind_direction_cart += 2.0*pi
    end

    cos_wdr = cos(-wind_direction_cart)
    sin_wdr = sin(-wind_direction_cart)

    # convert to cartesian coordinates with wind to positive x
    x_cart = xlocs*cos_wdr - ylocs*sin_wdr
    y_cart = xlocs*sin_wdr + ylocs*cos_wdr

    return x_cart, y_cart
end

"""
    point_velocity(loc, turbine_x, turbine_y, turbine_z, turbine_yaw, turbine_ct, turbine_ai,
    rotor_diameter, hub_height, turbine_local_ti, sorted_turbine_index, wtvelocities,
    wind_resource, model_set::AbstractModelSet;
    wind_farm_state_id=1, downwind_turbine_id=0)

Calculates the wind speed at a given point for a given state

# Arguments
- `loc::Array{TF,3}`: Location of interest
- `turbine_x::Array{TF,nTurbines}`: turbine east-west locations in the global
    reference frame
- `turbine_y::Array{TF,nTurbines}`: turbine north-south locations in the global
    reference frame
- `turbine_z::Array{TF,nTurbines}`: turbine base height in the global reference frame
- `turbine_yaw::Array{TF,nTurbines}`: turbine yaw for the given wind direction in
    radians
- `turbine_ct::Array{TF,nTurbines}`: turbine north-south locations in the global
    reference frame
- `turbine_ai::Array{TF,nTurbines}`: turbine axial induction for the given state
- `rotor_diameter::Array{TF,nTurbines}`: turbine rotor diameters
- `hub_height::Array{TF,nTurbines}`: turbine hub heights
- `turbine_local_ti::Array{TF,nTurbines}`: turbine local turbulence intensity for
    the given state
- `sorted_turbine_index::Array{TF,nTurbines}`: array containing indices of wind turbines
    from most upwind to most downwind turbine in the given state
- `wtvelocities::Array{TF,nTurbines}`: effective inflow wind speed for given state
- `wind_resource::DiscretizedWindResource`: contains wind resource discreption (directions,
    speeds, frequencies, etc)
- `wind_farm_state_id::Int`: index to correct state to use from wind resource provided.
    Defaults to 1
- `downwind_turbine_id::Int`: index of wind turbine of interest (if any). If not a point for
    calculating effective wind speed of a turbine, then provide 0 (default)
"""
function point_velocity(locx, locy, locz, turbine_x, turbine_y, turbine_z, turbine_yaw, turbine_ct, turbine_ai,
                    rotor_diameter, hub_height, turbine_local_ti, sorted_turbine_index, wtvelocities,
                    wind_resource, model_set::AbstractModelSet;
                    wind_farm_state_id=1, downwind_turbine_id=0, eps=1E-6)

    wakedeficitmodel = model_set.wake_deficit_model
    wakedeflectionmodel = model_set.wake_deflection_model
    wakecombinationmodel = model_set.wake_combination_model

    # extract flow information
    wind_speed = wind_resource.wind_speeds[wind_farm_state_id]
    reference_height = wind_resource.measurement_heights[wind_farm_state_id]
    wind_shear_model = wind_resource.wind_shear_model
    shear_exponent = wind_shear_model.shear_exponent

    # set ground height
    ground_height = turbine_z[1]    # TODO: allow topology to be given

    # get number of turbines
    nturbines = length(turbine_x)

    # initialize deficit summation term to zero
    deficit_sum = 0.0

    # initialize point velocity with shear to zero
    point_velocity_with_shear = 0.0

    # loop through all turbines
    for u=1:nturbines

        # get index of upstream turbine
        upwind_turb_id = Int(sorted_turbine_index[u])

        # skip this loop if it would include a turbine's impact on itself)
        if upwind_turb_id==downwind_turbine_id; continue; end

        # downstream distance between upstream turbine and point
        x = locx - turbine_x[upwind_turb_id]

        if x > 0.0#-eps*rotor_diameter[upwind_turb_id]

            # check turbine relative locations and use a spline to smooth the deficit jump
            # at the rotor hub location
            if x < 0.0

                # set horizontal and vertical deflection of upstream "wake"
                # calculate wake deflection of the current wake at epsilon
                horizontal_deflection = wake_deflection_model(turbine_x[upwind_turb_id] + eps*rotor_diameter[upwind_turb_id], locy, locz, turbine_x, turbine_yaw, turbine_ct,
                                upwind_turb_id, rotor_diameter, turbine_local_ti, wakedeflectionmodel)
                vertical_deflection = 0.0

                # get deficit at downstream end of spline
                deltav_downstream = wake_deficit_model(turbine_x[upwind_turb_id] + eps*rotor_diameter[upwind_turb_id], locy, locz, turbine_x, turbine_y, turbine_z, horizontal_deflection, vertical_deflection,
                                            upwind_turb_id, downwind_turbine_id, hub_height, rotor_diameter, turbine_ai,
                                            turbine_local_ti, turbine_ct, turbine_yaw, wakedeficitmodel)

                # set deficit derivative at downstream end of spline
                deltav_deriv_down_stream = 0.0

                # set deficit at upstream end of spline
                deltav_upstream = 0.0

                # set deficit derivative at upstream end of spline
                deltav_deriv_up_stream = 0.0

                # calculate deltav using a hermite spline across the rotor hub region
                deltav = hermite_spline(locx, turbine_x[upwind_turb_id]-eps*rotor_diameter[upwind_turb_id], turbine_x[upwind_turb_id]+eps*rotor_diameter[upwind_turb_id], deltav_upstream, deltav_deriv_up_stream, deltav_downstream, deltav_deriv_down_stream)

            else

                # calculate wake deflection of the current wake at the point of interest
                horizontal_deflection = wake_deflection_model(locx, locy, locz, turbine_x, turbine_yaw, turbine_ct,
                                upwind_turb_id, rotor_diameter, turbine_local_ti, wakedeflectionmodel)

                vertical_deflection = 0.0

                # velocity difference in the wake
                deltav = wake_deficit_model(locx, locy, locz, turbine_x, turbine_y, turbine_z, horizontal_deflection, vertical_deflection,
                                upwind_turb_id, downwind_turbine_id, hub_height, rotor_diameter, turbine_ai,
                                turbine_local_ti, turbine_ct, turbine_yaw, wakedeficitmodel)
                # println("horizontal_deflection: ", horizontal_deflection)
            end

            # combine deficits according to selected wake combination method
            turb_inflow = wtvelocities[upwind_turb_id]
            deficit_sum = wake_combination_model(deltav, wind_speed, wtvelocities[upwind_turb_id], deficit_sum, wakecombinationmodel)
        end
    end

    # find velocity at point without shear
    point_velocity_without_shear = wind_speed - deficit_sum

    # adjust sample point velocity for shear
    point_velocity_with_shear = adjust_for_wind_shear(locz, point_velocity_without_shear, reference_height, ground_height, wind_shear_model)

    return point_velocity_with_shear

end

"""
    point_velocity(turbine_x, turbine_y, turbine_z, rotor_diameter, hub_height, turbine_yaw,
    sorted_turbine_index, ct_model, rotor_sample_points_y, rotor_sample_points_z, wind_resource,
    model_set::AbstractModelSet; wind_farm_state_id=1)

Calculates the wind speed at a given point for a given state

# Arguments
- `turbine_x::Array{TF,nTurbines}`: turbine east-west locations in the global
    reference frame
- `turbine_y::Array{TF,nTurbines}`: turbine north-south locations in the global
    reference frame
- `turbine_z::Array{TF,nTurbines}`: turbine base height in the global reference frame
- `rotor_diameter::Array{TF,nTurbines}`: turbine rotor diameters
- `hub_height::Array{TF,nTurbines}`: turbine hub heights
- `turbine_yaw::Array{TF,nTurbines}`: turbine yaw for the given wind direction in
    radians
- `sorted_turbine_index::Array{TF,nTurbines}`: turbine north-south locations in the
    global reference frame
- `ct_model::AbstractThrustCoefficientModel`: defines how the thrust coefficient changes
    with state etc
- rotor_sample_points_y::Array{TF,N}`: horizontal wind location of points to sample across
    the rotor swept area when calculating the effective wind speed for the wind turbine.
    Points are centered at the hub (0,0) and scaled by the radius (1=tip of blades)
- rotor_sample_points_z::Array{TF,N}`: vertical wind location of points to sample across the
    rotor swept area when calculating the effective wind speed for the wind turbine. Points
    are centered at the hub (0,0) and scaled by the radius (1=tip of blades)
- `wind_resource::DiscretizedWindResource`: wind resource discreption (directions, speeds,
    frequencies, etc)
- `model_set::AbstractModelSet`: defines wake-realated models to be used in analysis
- `wind_farm_state_id::Int`: index to correct state to use from wind resource provided.
    Defaults to 1
"""
function turbine_velocities_one_direction(turbine_x, turbine_y, turbine_z, rotor_diameter, hub_height, turbine_yaw,
                    sorted_turbine_index, ct_model, rotor_sample_points_y, rotor_sample_points_z, wind_resource,
                    model_set::AbstractModelSet; wind_farm_state_id=1, velocity_only=true)

    # get number of turbines and rotor sample point
    n_turbines = length(turbine_x)
    n_rotor_sample_points = length(rotor_sample_points_y)

    arr_type = promote_type(typeof(turbine_x[1]),typeof(turbine_y[1]),typeof(turbine_z[1]),typeof(rotor_diameter[1]),
                            typeof(hub_height[1]),typeof(turbine_yaw[1]))

    turbine_velocities = zeros(arr_type, n_turbines)
    turbine_ct = zeros(arr_type, n_turbines)
    turbine_ai = zeros(arr_type, n_turbines)
    turbine_local_ti = zeros(arr_type, n_turbines)

    for d=1:n_turbines

        # get index of downstream turbine
        downwind_turbine_id = Int(sorted_turbine_index[d])

        # initialize downstream wind turbine velocity to zero
        # println("start array: ", turbine_velocities[downwind_turbine_id])
        # wind_turbine_velocity = typeof(turbine_velocities[downwind_turbine_id])(0.0)
        wind_turbine_velocity = 0.0
        # turbine_velocities[downwind_turbine_id] = 0.0

        for p=1:n_rotor_sample_points

            # scale rotor sample point coordinate by rotor diameter (in rotor hub ref. frame)
            local_rotor_sample_point_y = rotor_sample_points_y[p]*0.5*rotor_diameter[downwind_turbine_id]
            local_rotor_sample_point_z = rotor_sample_points_z[p]*0.5*rotor_diameter[downwind_turbine_id]

            locx = turbine_x[downwind_turbine_id] .+ local_rotor_sample_point_y*sin(turbine_yaw[downwind_turbine_id])
            locy = turbine_y[downwind_turbine_id] .+ local_rotor_sample_point_y*cos(turbine_yaw[downwind_turbine_id])
            locz = turbine_z[downwind_turbine_id] .+ hub_height[downwind_turbine_id] + local_rotor_sample_point_z

            # calculate the velocity at given point
            point_velocity_with_shear = point_velocity(locx, locy, locz, turbine_x, turbine_y, turbine_z, turbine_yaw, turbine_ct, turbine_ai,
                                    rotor_diameter, hub_height, turbine_local_ti, sorted_turbine_index, turbine_velocities,
                                    wind_resource, model_set,
                                    wind_farm_state_id=wind_farm_state_id, downwind_turbine_id=downwind_turbine_id)

            # add sample point velocity to turbine velocity to be averaged later
            wind_turbine_velocity += point_velocity_with_shear

        end

        # final velocity calculation for downstream turbine (average equally across all points)
        wind_turbine_velocity /= n_rotor_sample_points

        turbine_velocities[downwind_turbine_id] = deepcopy(wind_turbine_velocity)

        # update thrust coefficient for downstream turbine
        turbine_ct[downwind_turbine_id] = calculate_ct(turbine_velocities[downwind_turbine_id], ct_model[downwind_turbine_id])

        # update axial induction for downstream turbine
        turbine_ai[downwind_turbine_id] = _ct_to_axial_ind_func(turbine_ct[downwind_turbine_id])

        # update local turbulence intensity for downstream turbine
        ambient_ti = wind_resource.ambient_tis[wind_farm_state_id]
        turbine_local_ti[downwind_turbine_id] = calculate_local_ti(turbine_x, turbine_y, ambient_ti, rotor_diameter, hub_height, turbine_yaw, turbine_local_ti, sorted_turbine_index,
                            turbine_velocities, turbine_ct, model_set.local_ti_model; turbine_id=downwind_turbine_id, tol=1E-6)

    end

    if velocity_only
        return turbine_velocities
    else
        return turbine_velocities, turbine_ct, turbine_ai, turbine_local_ti
    end
end

function turbine_velocities_one_direction(x, turbine_z, rotor_diameter, hub_height, turbine_yaw,
    sorted_turbine_index, ct_model, rotor_sample_points_y, rotor_sample_points_z, wind_resource,
    model_set::AbstractModelSet; wind_farm_state_id=1, velocity_only=true)

    n_turbines = Int(length(x)/2)
    println(typeof(x), n_turbines)
    turbine_x = x[1:n_turbines]
    turbine_y = x[n_turbines+1:end]
    # println(turbine_x)
    println("turbine_x type ", typeof(turbine_x))
    println("type of x ", typeof(x))

    # get number of turbines and rotor sample point
    # n_turbines = length(turbine_x)
    n_rotor_sample_points = length(rotor_sample_points_y)

    arr_type = promote_type(typeof(turbine_x[1]),typeof(turbine_y[1]),typeof(turbine_z[1]),typeof(rotor_diameter[1]),
                typeof(hub_height[1]),typeof(turbine_yaw[1]))
    turbine_velocities = zeros(arr_type, n_turbines)
    turbine_ct = zeros(arr_type, n_turbines)
    turbine_ai = zeros(arr_type, n_turbines)
    turbine_local_ti = zeros(arr_type, n_turbines)

    for d=1:n_turbines

        # get index of downstream turbine
        downwind_turbine_id = Int(sorted_turbine_index[d])

        # initialize downstream wind turbine velocity to zero
        # println("start array: ", turbine_velocities[downwind_turbine_id])
        # wind_turbine_velocity = typeof(turbine_velocities[downwind_turbine_id])(0.0)
        wind_turbine_velocity = 0.0
        # turbine_velocities[downwind_turbine_id] = 0.0

        for p=1:n_rotor_sample_points

            # scale rotor sample point coordinate by rotor diameter (in rotor hub ref. frame)
            local_rotor_sample_point_y = rotor_sample_points_y[p]*0.5*rotor_diameter[downwind_turbine_id]
            local_rotor_sample_point_z = rotor_sample_points_z[p]*0.5*rotor_diameter[downwind_turbine_id]

            locx = turbine_x[downwind_turbine_id] .+ local_rotor_sample_point_y*sin(turbine_yaw[downwind_turbine_id])
            locy = turbine_y[downwind_turbine_id] .+ local_rotor_sample_point_y*cos(turbine_yaw[downwind_turbine_id])
            locz = turbine_z[downwind_turbine_id] .+ hub_height[downwind_turbine_id] + local_rotor_sample_point_z

            # calculate the velocity at given point
            point_velocity_with_shear = point_velocity(locx, locy, locz, turbine_x, turbine_y, turbine_z, turbine_yaw, turbine_ct, turbine_ai,
                                rotor_diameter, hub_height, turbine_local_ti, sorted_turbine_index, turbine_velocities,
                                wind_resource, model_set,
                                wind_farm_state_id=wind_farm_state_id, downwind_turbine_id=downwind_turbine_id)

            # add sample point velocity to turbine velocity to be averaged later
            wind_turbine_velocity += point_velocity_with_shear

        end

        # final velocity calculation for downstream turbine (average equally across all points)
        wind_turbine_velocity /= n_rotor_sample_points

        turbine_velocities[downwind_turbine_id] = deepcopy(wind_turbine_velocity)

        # update thrust coefficient for downstream turbine
        turbine_ct[downwind_turbine_id] = calculate_ct(turbine_velocities[downwind_turbine_id], ct_model[downwind_turbine_id])

        # update axial induction for downstream turbine
        turbine_ai[downwind_turbine_id] = _ct_to_axial_ind_func(turbine_ct[downwind_turbine_id])

        # update local turbulence intensity for downstream turbine
        ambient_ti = wind_resource.ambient_tis[wind_farm_state_id]
        turbine_local_ti[downwind_turbine_id] = calculate_local_ti(turbine_x, turbine_y, ambient_ti, rotor_diameter, hub_height, turbine_yaw, turbine_local_ti, sorted_turbine_index,
                    turbine_velocities, turbine_ct, model_set.local_ti_model; turbine_id=downwind_turbine_id, tol=1E-6)

    end

    if velocity_only
        return turbine_velocities
    else
        return turbine_velocities, turbine_ct, turbine_ai, turbine_local_ti
    end

end

# turbine_velocities_one_direction!(model_set::AbstractModelSet, problem_description::AbstractWindFarmProblem; wind_farm_state_id=1) = turbine_velocities_one_direction!([0.0], [0.0],
# model_set::AbstractModelSet, problem_description::AbstractWindFarmProblem; wind_farm_state_id=1)

"""
calculate_flow_field(xrange, yrange, zrange, model_set::AbstractModelSet, turbine_x,
    turbine_y, turbine_z, turbine_yaw, turbine_ct, turbine_ai, rotor_diameter, hub_height,
    turbine_local_ti, sorted_turbine_index, wtvelocities, wind_resource; wind_farm_state_id=1)

Generates a flow field for a given state and cross section

# Arguments
- `xrange::Range`: range defining east-west locations to sample in global reference frame
- `yrange::Range`: range defining north-west locations to sample in global reference frame
- `zrange::Range`: range defining vertical locations to sample in global reference frame
- `model_set::AbstractModelSet`: defines wake-realated models to be used in analysis
- `turbine_x::Array{TF,nTurbines}`: turbine east-west locations in the global
    reference frame
- `turbine_y::Array{TF,nTurbines}`: turbine north-south locations in the global
    reference frame
- `turbine_z::Array{TF,nTurbines}`: turbine base height in the global reference frame
- `turbine_yaw::Array{TF,nTurbines}`: turbine yaw for the given wind direction in
    radians
- `turbine_ct::Array{TF,nTurbines}`: thrust coefficient of each turbine for the given state
- `turbine_ai::Array{TF,nTurbines}`: turbine axial induction for the given state
- `rotor_diameter::Array{TF,nTurbines}`: turbine rotor diameters
- `hub_height::Array{TF,nTurbines}`: turbine hub heights
- `turbine_local_ti::Array{TF,nTurbines}`: turbine local turbulence intensity for
    the given state
- `sorted_turbine_index::Array{TF,nTurbines}`: turbine north-south locations in the
    global reference frame
- `wtvelocities::Array{TF,nTurbines}`: effective inflow wind speed for given state
- `wind_resource::DiscretizedWindResource`: wind resource discreption (directions, speeds,
    frequencies, etc)
- `wind_farm_state_id::Int`: index to correct state to use from wind resource provided.
    Defaults to 1
"""
function calculate_flow_field(xrange, yrange, zrange,
    model_set::AbstractModelSet, turbine_x, turbine_y, turbine_z, turbine_yaw, turbine_ct, turbine_ai,
    rotor_diameter, hub_height, turbine_local_ti, sorted_turbine_index, wtvelocities,
    wind_resource;
    wind_farm_state_id=1)

    xlen = length(xrange)
    ylen = length(yrange)
    zlen = length(zrange)
    npoints = xlen*ylen*zlen
    point_velocities = zeros(npoints)
    point_velocities = reshape(point_velocities, (zlen, ylen, xlen))

    for zi in 1:zlen
        for yi in 1:ylen
            for xi in 1:xlen
                locx = xrange[xi]
                locy = yrange[yi]
                locz = zrange[zi]
                locx, locy = rotate_to_wind_direction(locx, locy, wind_resource.wind_directions[wind_farm_state_id])

                point_velocities[zi, yi, xi] = point_velocity(locx, locy, locz, turbine_x, turbine_y, turbine_z, turbine_yaw, turbine_ct, turbine_ai,
                    rotor_diameter, hub_height, turbine_local_ti, sorted_turbine_index, wtvelocities,
                    wind_resource, model_set,
                    wind_farm_state_id=wind_farm_state_id, downwind_turbine_id=0)

            end
        end
    end

    if zlen == 1
        return point_velocities[1,1:ylen,1:xlen]
    elseif ylen == 1
        return point_velocities[1:zlen,1,1:xlen]
    elseif xlen == 1
        return point_velocities[1:zlen,1:ylen,1]
    end

end


text:
"""
    physical_eigenvalues(x,q,[args])

Computes the physical eigenvalues of the system given the selection criteria.
Empties the tables during the run.
"""
function physical_eigenvalues(x::Float64=0.0, q::Float64=1.0,
    tabOmega=tabOmega_serial, tabKappaSqOverTwoOmega=tabKappaSqOverTwoOmega_serial,
        tabAln=tabAln_serial, tabBln=tabBln_serial, tabCln=tabCln_serial,
        tabDln=tabDln_serial, tabFln=tabFln_serial, tabGln=tabGln_serial,
        tabHln=tabHln_serial, tabTruncMln=tabTruncMln_serial,
        tabEigValsMln=tabEigValsMln_serial)
    table_function_fill!(x,q,tabOmega,tabKappaSqOverTwoOmega)
    matrix_fill!(x,q,tabAln,tabBln,tabCln,tabDln,tabFln,tabGln,tabHln,tabOmega,tabKappaSqOverTwoOmega)
    tabTruncMln!(tabTruncMln,tabAln,tabBln,tabCln,tabDln,tabFln,tabGln,tabHln)
    tabEigValsMln!(tabTruncMln,tabEigValsMln)
    physical_eig = getPhysicalEigvals(tabEigValsMln)

    # clear the temp table
    tabEigValsMln_clear!(tabEigValsMln)
    tabTruncMln_clear!(tabTruncMln)
    matrix_clear!(tabAln,tabBln,tabCln,tabDln,tabFln,tabGln,tabHln)
    table_function_clear!(tabOmega,tabKappaSqOverTwoOmega)

    # return maximum growth rate
    return physical_eig
end

ext: jl

max_stars_repo_name: KerwannTEP/LiRA

max_issues_repo_name: KerwannTEP/LiRA

max_forks_repo_name: KerwannTEP/LiRA


text:
using Revise
using Random, DataFrames, PlotlyJS, SpatialIndexing
using Printf: @sprintf

const SI = SpatialIndexing

Random.seed!(1)
pts = Vector{SI.Point{Float64, 3}}()
for i in 1:100000
    r = 1 - 1/(1+abs(randn()))
    a = rand()*pi/2
    b = rand()*pi/2
    x = (sin(a) * cos(b))^3 * r
    y = (cos(a) * cos(b))^3 * r
    z = (sin(b))^3 * r
    push!(pts, SI.Point((x, y, z)))
end
seq_tree = RTree{Float64, 3}(Int, String, leaf_capacity = 20, branch_capacity = 20)
for (i, pt) in enumerate(pts)
    insert!(seq_tree, pt, i, string(i))
end
SI.check(seq_tree)

#=
rect = SI.Rect((-1.0, -2.0, -3.0), (2.0, -0.5, 1.0))
n_in_rect = sum(br -> in(br, rect), mbrs)
SI.subtract!(seq_tree, rect);
SI.check(seq_tree)
@show seq_tree.nelems seq_tree.nelem_deletions seq_tree.nnode_reinsertions seq_tree.nnode_splits
=#

# insert a point `pt` into pareto frontier
# (i.e. only if `tree` does not contain points that dominate it)
# remove all points that are dominated by `pt`
# return true if `pt` was inserted
function pareto_insert!(tree::RTree{T,N}, pt::SI.Point{T,N}, key, val) where {T,N}
    betterbr = SI.Rect(pt.coord, ntuple(_ -> typemax(T), N))
    # check if there are points that dominate pt
    isempty(tree, betterbr) || return false
    # remove the dominated points
    worsebr = SI.Rect(ntuple(_ -> typemin(T), N), pt.coord)
    SI.subtract!(tree, worsebr)
    # insert pt
    insert!(tree, pt, key, val)
end

pareto_tree = RTree{Float64, 3}(Int, String, leaf_capacity = 8, branch_capacity = 8)
for (i, pt) in enumerate(pts)
    pareto_insert!(pareto_tree, pt, i, string(i))
end
SI.check(pareto_tree)

include(joinpath(@__DIR__, "plot_utils.jl"))

pareto_tree_plot = plot(pareto_tree);
open(joinpath(@__DIR__, "pareto3d_rtree_seq.html"), "w") do io
    PlotlyJS.savehtml(io, pareto_tree_plot, :embed)
end

bulk_pareto_tree = RTree{Float64, 3}(Int, String, leaf_capacity = 8, branch_capacity = 8)
SI.load!(bulk_pareto_tree, pareto_tree)
SI.check(bulk_pareto_tree)

bulk_pareto_tree_plot = plot(bulk_pareto_tree);
open(joinpath(@__DIR__, "pareto3d_rtree_bulk.html"), "w") do io
    PlotlyJS.savehtml(io, bulk_pareto_tree_plot, :embed)
end

using ORCA
savefig(bulk_pareto_tree_plot, joinpath(@__DIR__, "pareto3d_rtree_bulk.png"), width=900, height=1000)


text:
using LoggingExtras, Logging
import Logging: AbstractLogger, handle_message

struct AzureSinkLogger{F} <: Logging.AbstractLogger
    action_chooser::F
    io::IO
end

AzureSinkLogger(f=_->"task.logissue") = AzureSinkLogger(f, stderr)

function Logging.handle_message(logger::AzureSinkLogger, args...; kwargs...)
    #Make it a named tuple for easier working
    log = LoggingExtras.handle_message_args(args...)
    task::String = logger.action_chooser(log)

    properties=String[]
    # AZP calls it `type` not level`, and `warning` not `warn`:
    levelmap = Dict(Logging.Error => "error", Logging.Warn => "warning")
    push!(properties, "type=$(levelmap[log[:level]])")
    for key in (:_module, :group, :id, :file, :line)
        push!(properties, "$key=$(log[key])")
    end
    for key in keys(log.kwargs)
        push!(properties, "$key=$(log.kwargs[key])")
    end
    props=join(properties, ";")
    # format: ##vso[area.action property1=value;property2=value;...]message
    println(logger.io, "##vso[$task $props]$(log.message)")
end
Logging.shouldlog(::AzureSinkLogger, arg...) = true
Logging.min_enabled_level(::AzureSinkLogger) = Logging.Warn
Logging.catch_exceptions(::AzureSinkLogger) = true

function enable_azure_logging()
    # Tee-in AzureSinkLogger so that `@warn` and `@error` are printed out nicely
    global_logger(TeeLogger(
        global_logger(),
        AzureSinkLogger(),
    ))
end


text:
using SpecialFunctions: loggamma, gamma
using Printf: @printf

T = Float64

α = 0.6
β = 0.8
x = 20.0
tol = 1e-16
x_vals = T[ n*5 for n = 5:10 ]

function asymp_terms(x, α, β, tol)
    Nopt = ceil(Integer, x^(1/α))
    @printf("α = %0.4f, β = %0.4f\n", α, β)
    @printf("x = %0.4f, Nopt = %0d\n", x, Nopt)
    sgn = one(T)
    s = zero(T)
    @printf("\n%5s  %12s  %12s  %15s\n\n",
            "n", "B_n/x^n", "term", "partial sum")
    for n = 1:Nopt
        A = sinpi(n*α-β)
        B = loggamma(n*α-β+1)
        Bxn = exp(B-n*log(x))
        term = sgn * A * Bxn / π
        s += term
        sgn = -sgn
        @printf("%5d  %12.4e  %12.4e  %20.15f\n", n, Bxn, term, s)
        if  Bxn < tol
            break
        end
    end
end

function minimal_x_N(α, β, tol, minx=10, maxN=100)
    first_n = 1
    huge = convert(T, typemax(1))
    for k = minx:maxN
        x = convert(T, k)
        if log(x) < floor(α * log(huge))
            Nopt = min(maxN, ceil(Integer, x^(1/α)))
        else
            Nopt = maxN
        end
        for n = first_n:Nopt
            bound = loggamma(n*α-β+1) - n*log(x)
            @printf("%5d  %8.3f  %12.4e\n", n, x, exp(bound))
            if bound < log(tol)
                @goto done
            end
        end
        first_n = Nopt
    end
    @label done
end

asymp_terms(x, α, β, tol)


text:
module modelutils

using Flux   # must be Flux@0.90
using BSON
using ArgCheck, Dates

export save, load!, chan3cat, unsqueeze, MultiDense, mlp, BRNNenc, randn_repar

e_k(T, n, k) = begin; out = zeros(T, n); out[k] = 1; out; end
e_k(n, k) = e_k(Float32, n, k)
e1(T, n) = e_k(T, n, 1)
e1(n) = e_k(Float32, n, 1)

unsqueeze(xs, dim) = reshape(xs, (size(xs)[1:dim-1]..., 1, size(xs)[dim:end]...));

randn_repar(σ::AbstractArray, n, d, stochastic=true) = !stochastic ? zeros(Float32, n, d) : σ .* randn(Float32, n,d)

if Flux.has_cuarrays()
    using CuArrays
    # note I'm using cu(zeros(...)), cu(randn(...)), CuArrays.randn(...) in particular occasionally
    # (apparently non-deterministically) kept bugging out on the GPU side :/.
    randn_repar(σ::CuArray, n, d, stochastic=true) = !stochastic ? cu(zeros(Float32, n, d)) : σ .* cu(randn(Float32, n,d))
    randn_repar(σ::TrackedArray{Float32, N, CuArray{Float32,N}}, n, d, stochastic=true) where N =
        !stochastic ? cu(zeros(Float32, n, d)) : σ .* cu(randn(Float32, n,d))
end

# Add a third channel to a 2-channel image Tensor
chan3cat(x::AbstractArray{T,4}) where T = cat(x, zero(x)[:,:,1:1,:], dims=3)
chan3cat(x::AbstractArray{T,3}) where T = cat(x, zero(x)[:,:,1:1], dims=3)

"""
    get_strtype_wo_params(x)
Useful for models defined here. The print (or typeof) for these structs result
in the overall type, plus a load of parameters, or details in parentheses. This
function just returns the first part: the overall type, without all the other
junk. This uses the Base.show methods, but captures the output via IOBuffer.
"""
function get_strtype_wo_params(x)
    _io = IOBuffer()
    print(_io, x)    # (does not print to stdout)
    tstr = String(take!(_io))
    close(_io)
    ixparen = findfirst("(", tstr)
    ixparen === nothing && return tstr
    return tstr[1:first(ixparen)-1]
end

"""
    MultiDense(Dense(in_1, out_1), Dense(in_2, out_2))

2× Dense layer for convenience for variational (μ/σ) like application. Allows
one to keep both layers together, and can be called like:

    m(x)  # = (m.Dense1(x), m.Dense2(x))

for an instance `m`; the result is a Tuple.
"""
struct MultiDense{F1, F2, S, T}
    Dense1::Dense{F1,S,T}
    Dense2::Dense{F2,S,T}
end
Flux.@treelike MultiDense    # new Flux calls this '@functor'. May need to replace.
(a::MultiDense)(x::AbstractArray) = (a.Dense1(x), a.Dense2(x))

function Base.show(io::IO, l::MultiDense)
  fs = [d.σ == identity ? "" : ", " * string(d.σ) for d in [l.Dense1, l.Dense2]]
  print(io, "MultiDense((", size(l.Dense1.W, 2), ", ", size(l.Dense1.W, 1), fs[1], ")")
  print(io, ", (", size(l.Dense2.W, 2), ", ", size(l.Dense2.W, 1), fs[2], "))")
end

"""
    mlp(units::Tuple; activation::Function, final::Function)

A convenience wrapper for defining a fully connected feedforward (MLP)
component. `units` is a tuple which specifies the number of units in each layer,
with the first and last components taken to be the inputs / outputs of the MLP;
all intermediate units are hidden layers.

Example:
    mlp((20,128,128,40)) === Chain(Dense(20,128,relu), Dense(128,128,relu), Dense(128, 40))

The activation functions can be customized via `activation=...` (for hidden layers)
and `final=...` for the final layer.
"""
function mlp(units::Tuple; activation::Function=relu, final::Function=identity)
    last = length(units)-1
    layers = [Dense(i,j, l==last ? final : activation) for (l,(i,j)) in
            enumerate(zip(units[1:end-1], units[2:end]))]
    return last == 1 ? layers[1] : Chain(layers...)
end
mlp(x...; activation::Function=relu, final::Function=identity) = mlp(x, activation=activation, final=final)

"""
    BRNNenc(d_in, d_state)

Bidirectional RNN using an LSTM in each direction with `d_in` dimensional inputs
and `d_state` dimension for _each_. This is an encoder variant which only uses
the final state (of each LSTM), concatenated. Therefore, its usage is as
follows for `xs`, a Vector where each element corresponds to the input at time
``t``:

    m(xs)::Vector

for an instance `m`, where the output vector is the final state concatenation.
"""
struct BRNNenc{V}
    forward::V
    backward::V
end

BRNNenc(in::Integer, out::Integer) = BRNNenc(LSTM(in, out), LSTM(in, out))

function (m::BRNNenc)(xs::AbstractVector)
    m.forward.(xs)
    m.backward.(reverse(xs))
    return vcat(m.forward.state[1], m.backward.state[1])
end

Flux.@treelike BRNNenc

function Base.show(io::IO, l::BRNNenc)
  print(io, "Bidirectional RNN. (Forward / Backward):")
  show(l.forward); print(", "); show(l.backward)
end

################################################################################
##  Save and Load utils for parameters of arbitrary models.                   ##
##                                                                            ##
##  These utils can be used for *any* models as they don't save the model def ##
##  itself, but uses the fact that the parameter 'vector' can be extracted    ##
##  easily using Flux's treelike structure and similarly loaded easily.       ##
################################################################################

"""
    save(fname::String, ps::Flux.Params, opt::Flux.ADAM; timestamp::Bool=false, force::Bool=false)
    save(fname::String, ps::Flux.Params; timestamp::Bool=false, force::Bool=false)

Save a parameter vector extracted via `Flux.params(models...)` to `fname`. in
BSON format. Since the state of the optimizer is critical for continuing
optimization (should it be required), the optimizer can optionally be given too.
At the moment this is limited to Adam, but this is easily extended.

Two options:
   (i) `timestamp = true` will insert a timestamp just before the extension.
   (ii) `force=true` will overwrite an existing file of the same name. The
        default behaviour is to throw an error instead.
"""
function save(fname::String, ps::Flux.Params, opt::Flux.ADAM; timestamp::Bool=false, force::Bool=false)
    fname = _get_and_check_fname(fname, timestamp, force)
    BSON.bson(fname, ps=[cpu(Tracker.data(p)) for p in ps], etabeta=(opt.eta, opt.beta),
                opt_state=[mapleaves(cpu, opt.state[p]) for p in ps])
end

function save(fname::String, ps::Flux.Params; timestamp::Bool=false, force::Bool=false)
    fname = _get_and_check_fname(fname, timestamp, force)
    BSON.bson(fname, ps=[cpu(Tracker.data(p)) for p in ps])
end

function _get_and_check_fname(fname::String, timestamp::Bool, force::Bool)
    fn = splitext(fname)
    @argcheck fn[2] in [".bson", ""]
    timestamp && (fname = fn[1] * Dates.format(Dates.now(), "yyyymmddTHH:MM:SS") * fn[2])
    !force && (@assert !isfile(fname) "fname already exists. To overwrite, use `force=true`.")
    return fname
end

"""
    load!(ps::Flux.Params, fname::String)
    load!(ps::Flux.Params, opt::Flux.ADAM, fname::String)

The inverse operation to save; loads a parameter vector from a BSON file created
per the `save` function. This operation will overwrite the parameter values in
the `ps` parameter 'vector' and optionally also the optimizer state (and η, β)
if specified as an argument.

The function takes an _existing_ parameter vector in order to link it to some
pre-specified models; this allows the file format to be agnostic to the original
model.
"""
function load!(ps::Flux.Params, fname::String)
    tf = Flux.has_cuarrays() && Tracker.data(first(ps)) isa CuArray ? gpu : identity
    _load_pars!(ps, BSON.load(fname)[:ps], tf)
end

function load!(ps::Flux.Params, opt::Flux.ADAM, fname::String)
    tf = Flux.has_cuarrays() && Tracker.data(first(ps)) isa CuArray ? cu : identity
    f=BSON.load(fname)
    _load_pars!(ps, f[:ps], tf)
    opt.eta, opt.beta = f[:etabeta]
    for (p, p_saved) in zip(ps, f[:opt_state])
        opt.state[p] = tf(p_saved)
    end
end

function _load_pars!(ps_to::Flux.Params, ps_from::Vector, tf::Function)
    try
        for (p, p_saved) in zip(ps_to, ps_from)
            p.data .= tf(p_saved)
        end
    catch e
        # Debugging dimension mismatch:
        if e isa DimensionMismatch
            @warn "Specified file has parameters of different dimensions:"; flush(stderr)
            for (p, p_saved) in zip(ps_to, ps_from)
                txt = ("(cur/new): ", size(p.data), " <-- ", size(p_saved) , "\n")
                if size(p.data) == size(p_saved)
                    print(txt...)
                else
                    printstyled(IOContext(stdout, :color => true), txt...; bold=true)
                end
            end
        end
        rethrow(e)
    end
end

end

max_stars_repo_name: ornithos/mtds-dblpend

max_issues_repo_name: ornithos/mtds-dblpend

max_forks_repo_name: ornithos/mtds-dblpend


text:
"""
# Euclid's Algorithm for GCD

greatest common divisor (GCD) is the greatest positive integer that able to divide all provided number evenly without remainder

# Reference
- Greatest Common Divisor -- https://brilliant.org/wiki/greatest-common-divisor/

# Example
```julia
gcd_euclid_algorithm([30, 12, 6]) # returns 6
gcd_euclid_algorithm([13, 48]) # returns 1
gcd_euclid_algorithm([2, 4, 6, 8, 16]) # returns 2
```

"""
function gcd_euclid_algorithm(arr::Vector{N})where N<:Integer

    result = arr[1]

    for i in 2:length(arr)
        result = find_gcd(result,arr[i])
        if result == 1
            break
         end
    end

    return result
end

function find_gcd(number1::Int,number2::Int)::Int
    while number2!=0
        number1, number2 = number2, number1 % number2
    end

    return number1
end

max_stars_repo_name: gohpeijin/Julia

max_issues_repo_name: gohpeijin/Julia

max_forks_repo_name: gohpeijin/Julia


text:
## Metsa-Oy Forest Products Supply Chain
## Edward J. Xu <edxu96@outlook.com>
## Jan 24th, 2020

function cal_dis_profit_2(h_mt, y_mi, y_mj, y_paper_m, z_mik, z_mjk, z_paper_mk,
    alpha_t, beta_t, a_i, b_i, e_i, c_i, r_saw, r_plywood, a_j,
    c_j, r_j, a_paper_t, b_paper_j, c_paper, r_paper, gamma_ik, delta_ik,
    gamma_jk, delta_jk, gamma_paper_k, delta_paper_k, p_fuel, sigma, omega_i,
    omega_j, omega_paper, nu_saw, nu_plywood, nu_j, nu_paper, o_saw, o_plywood,
    o_j, o_paper, i_t1, i_t2, j_t2, x_saw_m, x_plywood_m, x_mj, x_paper_m; m)

  dis_profit =
    sigma^(m-1) * (
      -sum(alpha_t[t] * h_mt[m, t] + beta_t[t] * h_mt[m, t]^2
        for t = 1:6) +
      -sum(c_i[i] * y_mi[m, i] for i=1:5) +
      -sum(c_j[j] * y_mj[m, j] for j=1:2) -
      c_paper * y_paper_m[m] +
      sum(alpha_t[t] * (h_mt[m, t] - sum(a_i[i] * y_mi[m, i] for i in
        i_t1[t])) for t=1:3) +
      sum(alpha_t[t] * (h_mt[m, t] + sum(b_i[i] * y_mi[m, i] for i in
        i_t2[t-3]) - sum(a_j[j] * y_mj[m, j] for j in j_t2[t-3]) -
        a_paper_t[t-3] * y_paper_m[m]) for t in [4, 6]) +
      ## when t = 5
      alpha_t[5] * (h_mt[m, 5] + sum(b_i[i] * y_mi[m, i] for i in
        i_t2[2]) - a_paper_t[2] * y_paper_m[m]) +
      p_fuel * sum(e_i[i] * y_mi[m, i] for i=1:5) +
      sum(omega_i[i]^(m-1) * gamma_ik[i, k] * z_mik[m, i, k] -
        omega_i[i]^(m-1) * delta_ik[i, k] * z_mik[m, i, k]^2
        for i=1:5, k=1:4) +
      sum(omega_j[j]^(m-1) * gamma_jk[j, k] * z_mjk[m, j, k] -
        omega_j[j]^(m-1) * delta_jk[j, k] * z_mjk[m, j, k]^2
        for j=1:2, k=1:4) +
      sum(omega_paper^(m-1) * gamma_paper_k[k] * z_paper_mk[m, k] -
        omega_paper^(m-1) * delta_paper_k[k] * z_paper_mk[m, k]^2
        for k=1:4)
      )

  if m <= 2
    m += 1
    dis_profit += - sigma^(m-2) * (o_saw * (x_saw_m[m] - x_saw_m[m - 1]) +
      o_plywood * (x_plywood_m[m] - x_plywood_m[m-1]) +
      sum(o_j[j] * (x_mj[m, j] - x_mj[m-1, j]) for j=1:2) +
      o_paper * (x_paper_m[m] - x_paper_m[m-1]))
  end

  println("dis_profit = $(round(dis_profit; digits=4)) ;")
end

function print_result_2(result_h_mt, result_y_mi, result_y_mj, result_y_paper_m,
    result_z_mik, result_z_mjk, result_z_paper_mk, alpha_t, beta_t, a_i, b_i,
    e_i, c_i, r_saw, r_plywood, a_j, c_j, r_j, a_paper_t, b_paper_j, c_paper,
    r_paper, gamma_ik, delta_ik, gamma_jk, delta_jk, gamma_paper_k,
    delta_paper_k, p_fuel, sigma, omega_i, omega_j, omega_paper, nu_saw,
    nu_plywood, nu_j, nu_paper, o_saw, o_plywood, o_j, o_paper, i_t1, i_t2,
    j_t2, result_x_saw_m, result_x_plywood_m, result_x_mj, result_x_paper_m)

  h_mt = result_h_mt
  y_mi = result_y_mi
  y_mj = result_y_mj
  y_paper_m = result_y_paper_m
  z_mik = result_z_mik
  z_mjk = result_z_mjk
  z_paper_mk = result_z_paper_mk
  x_saw_m = result_x_saw_m
  x_plywood_m = result_x_plywood_m
  x_mj = result_x_mj
  x_paper_m = result_x_paper_m

  ## Sales distribution
  for m=1:3
    println("m = $m ;")

    for i=1:5
      println(round.(result_z_mik[m, i, :] / sum(result_z_mik[m, i, :]) * 100;
        digits=2))
    end

    for j=1:2
      println(round.(result_z_mjk[m, j, :] / sum(result_z_mjk[m, j, :]) * 100;
        digits=2))
    end

    println(round.(result_z_paper_mk[m, :] / sum(result_z_paper_mk[m, :]) * 100;
      digits=2))

    ## discounted profits
    cal_dis_profit_2(h_mt, y_mi, y_mj, y_paper_m, z_mik, z_mjk, z_paper_mk,
      alpha_t, beta_t, a_i, b_i, e_i, c_i, r_saw, r_plywood, a_j, c_j, r_j,
      a_paper_t, b_paper_j, c_paper, r_paper, gamma_ik, delta_ik, gamma_jk,
      delta_jk, gamma_paper_k, delta_paper_k, p_fuel, sigma, omega_i, omega_j,
      omega_paper, nu_saw, nu_plywood, nu_j, nu_paper, o_saw, o_plywood, o_j,
      o_paper, i_t1, i_t2, j_t2, x_saw_m, x_plywood_m, x_mj, x_paper_m, m=m)

    ## capacities
    println(round.(x_saw_m[m]; digits=2))
    println(round.(x_plywood_m[m]; digits=2))
    println(round.(x_mj[m, :]; digits=2))
    println(round.(x_paper_m[m]; digits=2))
  end
end


text:
#=
This file uses the Kaplan Yorke dimension to provide high accuracy estimates for
the fractal dimension for the Lorenz-96 system.
=#
using DrWatson
@quickactivate :FractalDimension # uses DynamicalSystems, PyPlot
using DynamicalSystems, Statistics, OrdinaryDiffEq

Dmax = 13
ics = 10 # how many random initial conditions to obtain
diffeq = (alg = Vern9(), reltol = 1e-9, abstol = 1e-9, maxiters = typemax(Int))
Ds = 5:Dmax
K = [Float64[] for _ ∈ 1:length(Ds)]
Λ = [Vector{Float64}[] for _ ∈ 1:length(Ds)]
T = 100000.0
F = 8.0

for (i, D) ∈ enumerate(Ds)
    @show D
    for j ∈ 1:ics
        ds = Systems.lorenz96(D; F = F)
        @time λs = lyapunovs(ds, T; Ttr = 100.0, diffeq...)
        push!(Λ[i], λs)
        push!(K[i], kaplanyorke_dim(λs))
    end
end

name = savename("lorenz96", @dict(Dmax, ics, T, F), "jld2")
result = @strdict Ds ics K Λ
@tagsave(datadir(name), result)

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
function printnodes(tp::Union{NodeTreeProblem,TreeProblem};
    outfile::String = "")
    output = ""
    for u in getnodes(tp.bt, tp.bt.depth), a in 1:tp.pd.npop
        level = round(sum(JuMP.getvalue(tp.assign[a,u,:])))/tp.nlevels
        if level > 0
            output = string(output, string(tp.pd.pops[a], "\t", u, "\t", level, "\n"))
        end
    end
    if outfile == ""
        print(output)
    else
        open(outfile, "w") do f
            write(f, output)
        end
    end
end

# warning: lots of magic constants here
function printtree(tp::Union{NodeTreeProblem,TreeProblem})

    pd = tp.pd
    bt = tp.bt
    depth = bt.depth

    spaces = Dict(zip(0:depth, [3 for i in 0:depth]))
    for d=(depth-1):-1:0
        spaces[d] = 2*spaces[d+1]+1
    end
    for d=0:depth
        nodes = getnodes(bt,d)
        for n in nodes
            print(" "^spaces[d])
            if d == depth
                nodeassign = round.(JuMP.getvalue(tp.assign[:,n,1]))
            else
                nodeassign = 0
            end
            if sum(nodeassign) < 1
                str = "()"
            else
                str = pd.pops[findfirst(nodeassign .> 0)][1:2]
            end
            Printf.@printf("%2s", str)
            print(" "^spaces[d])
        end
        println()
        if d < depth
            for n in nodes
                print(" "^(spaces[d]-2), " /"," "^2,"\\ "," "^(spaces[d]-2))
            end
            println()
            for n in nodes
                weight1 = JuMP.getvalue(tp.weight[(n,getchildren(bt,n)[1])])
                weight2 = JuMP.getvalue(tp.weight[(n,getchildren(bt,n)[2])])
                print(" "^(spaces[d]-3))
                if round(weight1) > 0
                    Printf.@printf("%3.0f", weight1)
                else
                    Printf.@printf("%3s", " ")
                end
                print(" "^2)
                if round(weight2) > 0
                    Printf.@printf("%3.0f", weight2)
                else
                    Printf.@printf("%3s", " ")
                end
                print(" "^(spaces[d]-3))
            end
            println()
            for n in nodes
                print(" "^(spaces[d]-2),"/ "," "^2," \\"," "^(spaces[d]-2))
            end
            println()
        end
    end

end

function Base.show(io::IO, tp::Union{NodeTreeProblem,TreeProblem}; offset::String="")
    println(io, offset, string(typeof(tp)))
end


text:
## CATEGORICAL ARRAY DECODER UTILITY

# """
#     reconstruct(A)

# For reconstructing categorical arrays from their elements alone. Here
# `A` is of type `AbstractArray{T}` where `T` is a subtype of
# `CategoricalString` or `CategoricalValue`. The function `reconstruct` has
# the property that `reconstruct(broadcast(identity, C)) == C`, whenever `C`
# is a `CategoricalArray`. In other words, `reconstruct` is a left-inverse
# for the function `C -> broadcast(identity, C)` that strips a
# CategoricalArray of its "categorical wrapper".

# Does not handle missing values.

# """
# function reconstruct(A::AbstractArray{<:CategoricalValue{T},N}) where {T,N}
#     firstnonmissing = findfirst(x->!ismissing(x), A)
#     isnothing(firstnonmissing) && error("No non-missing values encountered. ")
#     pool = A[firstnonmissing].pool
#     refs = broadcast(x -> x.level, A)
#     return CategoricalArray{T,N}(refs, pool)
# end
# function reconstruct(A::AbstractArray{<:CategoricalString,N}) where {T,N}
#     firstnonmissing = findfirst(x->!ismissing(x), A)
#     isnothing(firstnonmissing) && error("No non-missing values encountered. ")
#     pool = A[firstnonmissing].pool
#     refs = broadcast(x -> x.level, A)
#     return CategoricalArray{String,N}(refs, pool)
# end

CategoricalElement = Union{CategoricalValue,CategoricalString}

"""
    classes(x)

All the categorical values in the same pool as `x` (including `x`),
returned as a list, with an ordering consistent with the pool. Here
`x` has `CategoricalValue` or `CategoricalString` type, and
`classes(x)` is a vector of the same eltype.

Not to be confused with the levels of `x.pool` which have a
different type. In particular, while `x in classes(x)` is always
true, `x in x.pool.levels` is not true.

    julia> v = categorical([:c, :b, :c, :a])
    julia> levels(v)
    3-element Array{Symbol,1}:
     :a
     :b
     :c
    julia> classes(v[4])
    3-element Array{CategoricalValue{Symbol,UInt32},1}:
     :a
     :b
     :c

"""
function classes(x::CategoricalElement)
    p = x.pool
    return [p.valindex[p.invindex[v]] for v in p.levels]
end

raw(x::CategoricalElement) = x.pool.index[x.level] # a method just for testing

"""
   int(x)

The positional integer of the `CategoricalString` or `CategoricalValue`
`x`, in the ordering defined by the pool of `x`. The type of `int(x)`
is the refrence type of `x`.

Not to be confused with `x.ref`, which is unchanged by reordering of
the pool of `x`, but has the same type.

    int(X::CategoricalArray)
    int(W::Array{<:CategoricalString})
    int(W::Array{<:CategoricalValue})

Broadcasted versions of `int`.

    julia> v = categorical([:c, :b, :c, :a])
    julia> levels(v)
    3-element Array{Symbol,1}:
     :a
     :b
     :c
    julia> int(v)
    4-element Array{UInt32,1}:
     0x00000003
     0x00000002
     0x00000003
     0x00000001

See also: decoder
"""
int(x::CategoricalElement) = x.pool.order[x.pool.invindex[x]]
int(X::CategoricalArray) = broadcast(r -> X.pool.order[r], X.refs)
int(V::Array{<:CategoricalElement}) = broadcast(int, V)

struct CategoricalDecoder{T,R} # <: MLJType
    pool::CategoricalPool{T,R}
    invorder::Vector{Int}
end

"""
    d = decoder(x)

A callable object for decoding the integer representation of a
`CategoricalString` or `CategoricalValue` sharing the same pool as
`x`. (Here `x` is of one of these two types.) Specifically, one has
`d(int(y)) == y` for all `y in classes(x)`. One can also call `d` on
integer arrays, in which case `d` is broadcast over all elements.

    julia> v = categorical([:c, :b, :c, :a])
    julia> int(v)
    4-element Array{UInt32,1}:
     0x00000003
     0x00000002
     0x00000003
     0x00000001
    julia> d = decoder(v[3])
    julia> d(int(v)) == v
    true

See also: int, classes

"""
decoder(element::CategoricalElement) =
    CategoricalDecoder(element.pool, sortperm(element.pool.order))
## in the next lot need to skip the missing one
# decoder(X::CategoricalArray) = CategoricalDecoder(X.pool)
# function decoder(V::Array{<:CategoricalElement})
#     isempty(V) && error("Unable to extract decoder from empty array. ")
#     return X[1]
# end

(decoder::CategoricalDecoder{T,R})(i::Integer) where {T,R} =
    CategoricalValue{T,R}(decoder.invorder[i], decoder.pool)
(decoder::CategoricalDecoder{String,R})(i::Integer) where R =
    CategoricalString{R}(decoder.invorder[i], decoder.pool)
(decoder::CategoricalDecoder)(I::AbstractArray{<:Integer}) = broadcast(decoder, I)

## TABULAR DATA

const istable = Tables.istable

# hack for detecting JuliaDB.NDSparse tables without loading as dependency:
isndsparse(X) = isdefined(X, :data_buffer)

"""
    container_type(X)

Return `:table`, `:sparse`, or `:other`, according to whether `X` is a
supported table format, a supported sparse table format, or something
else.

The first two formats, together abstract vectors, support the
`MLJBase` accessor methods `selectrows`, `selectcols`, `select`,
`nrows`, `schema`, and `union_scitypes`.

"""
function container_type(X)
    if istable(X)
        return :table
    elseif isndsparse(X)
        return :sparse
    else
        return :other
    end
end

## UTILITY FOR CONVERTING BETWEEN TABULAR DATA AND MATRICES

"""
    MLJBase.matrix(X)

Convert a table source `X` into an `Matrix`; or, if `X` is
a `AbstractMatrix`, return `X`. Optimized for column-based sources.

If instead X is a sparse table, then a `SparseMatrixCSC` object is
returned. The integer relabelling of column names follows the
lexicographic ordering (as indicated by `schema(X).names`).

"""
matrix(X) = matrix(Val(container_type(X)), X)
matrix(::Val{:other}, X) = throw(ArgumentError)
matrix(::Val{:other}, X::AbstractMatrix) = X

function matrix(::Val{:table}, X)
    cols = Tables.columns(X) # property-accessible object
    mat = reduce(hcat, [getproperty(cols, ftr) for ftr in propertynames(cols)])
    # tightest eltype:
    return broadcast(identity, mat)
end

function matrix(::Val{:sparse}, X)
    K = keys(X)
    features = schema(X).names
    index_given_feature = Dict{Symbol,Int}()
    for j in eachindex(features)
        index_given_feature[features[j]] = j
    end
    I = [k[1] for k in K]
    J = [index_given_feature[k[2]] for k in K]
    V = [v[1] for v in values(X)]
    return sparse(I, J, V)
end

"""
    MLJBase.table(cols; prototype=cols)

Convert a named tuple of vectors `cols`, into a table. The table
type returned is the "preferred sink type" for `prototype` (see the
Tables.jl documentation).

    MLJBase.table(X::AbstractMatrix; names=nothing, prototype=nothing)

Convert an abstract matrix `X` into a table with `names` (a tuple of
symbols) as column names, or with labels `(:x1, :x2, ..., :xn)` where
`n=size(X, 2)`, if `names` is not specified.  If prototype=nothing,
then a named tuple of vectors is returned.

Equivalent to `table(cols, prototype=prototype)` where `cols` is the
named tuple of columns of `X`, with `keys(cols) = names`.

"""
function table(cols::NamedTuple; prototype=cols)
    istable(prototype) || error("prototype is not tabular.")
    return Tables.materializer(prototype)(cols)
end
function table(X::AbstractMatrix; names=nothing, prototype=nothing)
    if names == nothing
        _names = tuple([Symbol(:x, j) for j in 1:size(X, 2)]...)
    else
        _names = names
    end
    cols = NamedTuple{_names}(tuple([X[:,j] for j in 1:size(X, 2)]...))
    _prototype = (prototype == nothing ? cols : prototype)
    return table(cols; prototype=_prototype)
end

## UNIFIED API FOR ACCESSING TABLES, MATRICES AND VECTORS

"""
    selectrows(X, r)

Select single or multiple rows from any table, sparse table, or
abstract vector `X`.  If `X` is tabular, the object returned is a
table of the preferred sink type of `typeof(X)`, even a single row is
selected.

"""
selectrows(X, r) = selectrows(Val(container_type(X)), X, r)
selectrows(::Val{:other}, X, r) = throw(ArgumentError)

"""
    selectcols(X, c)

Select single or multiple columns from any table or sparse table
`X`. If `c` is an abstract vector of integers or symbols, then the
object returned is a table of the preferred sink type of
`typeof(X)`. If `c` is a *single* integer or column, then a `Vector`
or `CategoricalVector` is returned.

"""
selectcols(X, c) = selectcols(Val(container_type(X)), X, c)
selectcols(::Val{:other}, X, c) = throw(ArgumentError)

"""
    select(X, r, c)

Select element of a table or sparse table at row `r` and column
`c`. In the case of sparse data where the key `(r, c)`, zero or
`missing` is returned, depending on the value type.

See also: selectrows, selectcols

"""
select(X, r, c) = select(Val(container_type(X)), X, r, c)
select(::Val{:other}, X, r, c) = throw(ArgumentError)

"""
    schema(X)

Returns a struct with properties `names`, `types`
with the obvious meanings. Here `X` is any table or sparse table.

"""
schema(X) = schema(Val(container_type(X)), X)
schema(::Val{:other}, X) = throw(ArgumentError)

"""
    nrows(X)

Return the number of rows in a table, sparse table, or abstract vector.

"""
nrows(X) = nrows(Val(container_type(X)), X)
nrows(::Val{:other}, X) = throw(ArgumentError)

# project named tuple onto a tuple with only specified `labels` or indices:
project(t::NamedTuple, labels::AbstractArray{Symbol}) = NamedTuple{tuple(labels...)}(t)
project(t::NamedTuple, label::Colon) = t
project(t::NamedTuple, label::Symbol) = project(t, [label,])
project(t::NamedTuple, indices::AbstractArray{<:Integer}) =
    NamedTuple{tuple(keys(t)[indices]...)}(tuple([t[i] for i in indices]...))
project(t::NamedTuple, i::Integer) = project(t, [i,])

# multiple columns:
function selectcols(::Val{:table}, X, c::Union{Colon, AbstractArray{I}}) where I<:Union{Symbol,Integer}
    cols = Tables.columntable(X) # named tuple of vectors
    newcols = project(cols, c)
    return Tables.materializer(X)(newcols)
end

# single column:
function selectcols(::Val{:table}, X, c::I) where I<:Union{Symbol,Integer}
    cols = Tables.columntable(X) # named tuple of vectors
    return cols[c]
end

# multiple rows (using columntable):
function selectrows(::Val{:table}, X, r::Union{Colon,AbstractVector{<:Integer}})
    cols = Tables.columntable(X)
    new_cols = NamedTuple{keys(cols)}(tuple([c[r] for c in values(cols)]...))
    return Tables.materializer(X)(new_cols)
end

# single row (using columntable):
function selectrows(::Val{:table}, X, r::Integer)
    cols = Tables.columntable(X)
    new_cols = NamedTuple{keys(cols)}(tuple([c[r:r] for c in values(cols)]...))
    return Tables.materializer(X)(new_cols)
end

## ALTERNATIVE CODE FOR PREVIOUS TWO FUNCTIONS.
## ROWTABLE SELECTION OF ROWS INSTEAD OF COLUMNTABLE SELECTION
# # multiple rows:
# function selectrows(::Val{:table}, X, r::Union{Colon,AbstractVector{<:Integer}})
#     rows = Tables.rowtable(X) # vector of named tuples
#     return Tables.materializer(X)(rows[r])
# end

# # single row:
# function selectrows(::Val{:table}, X, r::Integer)
#     rows = Tables.rowtable(X) # vector of named tuples
#     return Tables.materializer(X)([rows[r]])
# end

select(::Val{:table}, X, r::Integer, c::Symbol) = selectcols(selectrows(X, r), c)[1]
select(::Val{:table}, X, r::Integer, c) = selectcols(selectrows(X, r), c)
select(::Val{:table}, X, r, c::Symbol) = selectcols(X, c)[r]
select(::Val{:table}, X, r, c) = selectcols(selectrows(X, r), c)

function schema(::Val{:table}, X)
    istable(X) || throw(ArgumentError)
    if !Tables.columnaccess(X)
        return Tables.schema(Tables.rows(X))
    else
        return Tables.schema(Tables.columns(X))
    end
end

function nrows(::Val{:table}, X)
    if !Tables.columnaccess(X)
        return length(collect(X))
    else
        cols = Tables.columntable(X)
        !isempty(cols) || return 0
        return length(cols[1])
    end
end

## ACCESSORS FOR ABSTRACT VECTORS

selectrows(::Val{:other}, v::AbstractVector, r) = v[r]
nrows(::Val{:other}, v::AbstractVector) = length(v)
selectrows(::Val{:other}, v::CategoricalVector, r) = @inbounds v[r]

## ACCESSORS FOR JULIA NDSPARSE ARRAYS (N=2)

nrows(::Val{:sparse}, X) = maximum([r[1] for r in keys(X)])
function select(::Val{:sparse}, X, r::Integer, c::Symbol)
    try
        X[r,c][1]
    catch exception
        exception isa KeyError || throw(exception)
        missing
    end
end
select(::Val{:sparse}, X, r::AbstractVector{<:Integer}, c::Symbol) = [select(X, s, c) for s in r]
select(::Val{:sparse}, X, ::Colon, c::Symbol) = [select(X, s, c) for s in 1:nrows(X)]
selectrows(::Val{:sparse}, X, r::Integer) = X[r:r,:]
selectrows(::Val{:sparse}, X, r) = X[r,:]
selectcols(::Val{:sparse}, X, c::Symbol) = select(X, :, c)
selectcols(::Val{:sparse}, X, c::AbstractVector{Symbol}) = X[:,sort(c)]
selectcols(::Val{:sparse}, X, ::Colon) = X
select(::Val{:sparse}, X, r::Integer, c::AbstractVector{Symbol}) = X[r,sort(c)]
select(::Val{:sparse}, X, r::Integer, ::Colon) = X[r,:]
select(::Val{:sparse}, X, r, c) = X[r,sort(c)]

function schema(::Val{:sparse}, X)
    names = sort(unique([r[2] for r in keys(X)]))
    types = [eltype(selectcols(X, name)) for name in names]
    return Tables.Schema(names, types)
end


text:
countries = readcsv("countries.csv")
names = filter(x -> length(x) > 0, replace.(lowercase.(countries[2:end,2]), " ", "_"))

if !isdefined(:loaded) || !loaded
    jldopen("90.jld") do file
    # jldopen("85.jld") do file
        global X, vocab, vocab2idx, idx2vocab
        X = read(file, "X")
        vocab = read(file, "vocab")
        vocab2idx = read(file, "vocab2idx")
        idx2vocab = read(file, "idx2vocab")
        global loaded = true
    end

    # @show length(vocab)
end

# note: transpoed
X_c = zeros(size(X, 2), length(names))
total_missing = 0
for (i, name) in enumerate(names)
    try
        idx = vocab2idx[name]
        X_c[:,i] = X[idx,:]
    catch e
        if isa(e, KeyError)
            warn("key $name not found in vocab")
            total_missing += 1
        else
            rethrow(e)
        end
    end
end
if total_missing > 0
    warn("Missing $total_missing")
end

size(X_c)

using HDF5
h5open("96.jl.h5", "w") do f
    write(f, "X", X_c)
end


text:
function fix_text!(ax, default=:small)
    sc = ax.scene[OldAxis]
    if default == :big
        sc[:names, :textsize] = (25, 25, 25)
        sc[:ticks, :textsize][] = (15, 20, 20)
        sc[:frame, :linewidth][] = (180, 180, 180)
    else
        sc[:names, :textsize] = (10, 10, 10)
        sc[:ticks, :textsize][] = (10, 10, 10)
        sc[:frame, :linewidth][] = (180, 180, 180)
    end
end

function default3d_plot!(fig, f::ScalarField, cmap; kwargs...)
    if !haskey(kwargs, :levels)
        kwargs = Iterators.flatten((kwargs, (:levels=>6,)))
    end

    contour(fig[1,1], f,
        alpha = 0.3, transparency = true, colormap=cmap; kwargs...)
end

function default3d_plot!(fig, f::ScalarVariable, cmap; kwargs...)
    size = get(kwargs, :size3d, 35)

    scattervariable(fig[1,1], f, size=size)
end

function section_plot!(ax, f::ScalarField, slice_dir, idx_value; kwargs...)
    dim = dir_to_idx(slice_dir)
    levels = get(kwargs, :levels, 6)
    f_section = @lift(slice(f, dim, $idx_value))
    labels = string.(filter(i->i≠slice_dir, (:x,:y,:z)))

    plt = contour!(ax, f_section, levels=levels)

    ax.xlabel = labels[1]
    ax.ylabel = labels[2]
    ax.aspect = DataAspect()

    return plt
end

function section_plot!(ax, f::ScalarVariable, slice_dir, idx_value; kwargs...)
    dim = dir_to_idx(slice_dir)
    ϵ = get(kwargs, :ϵ, 0.2)
    grid = getdomain(f)
    sort!(f, dim)

    f_section = @lift slice(f, dim, $idx_value, ϵ)
    labels = string.(filter(i->i≠slice_dir, (:x,:y,:z)))

    plt = scattervariable!(ax, f_section, size=2)

    ax.xlabel = labels[1]
    ax.ylabel = labels[2]
    ax.aspect = DataAspect()

    return plt
end

function add_legend!(fig, section_plt, ::ScalarField; kwargs...)
    levels = get(kwargs, :levels, 6)
    legend = Colorbar(fig, section_plt,
        label = "Field values",
        width = Relative(3/4),
        height = 10,
        vertical = false,
        tellheight = true,
        colormap = cgrad(:viridis, levels, categorical = true))

    fig[2,1:2] = legend
end

function add_legend!(fig, section_plt, ::ScalarVariable; kwargs...)
    # legend = Colorbar(fig, section_plt,
    #     label = "Variable values",
    #     width = Relative(3/4),
    #     height = 10,
    #     vertical = false,
    #     tellheight = true,
    #     colormap = cgrad(:viridis))

    # fig[2,1:2] = legend
end

function slider(f::ScalarField, slice_dir; kwargs...)
    dim = dir_to_idx(slice_dir)
    grid = getdomain(f)
    slider_domain = grid[dim]
    sl = Slider(axes(slider_domain, 1))
    sl_label = @lift round(slider_domain[$(sl.value)], digits=2)

    return sl, sl_label
end

function slider(f::ScalarVariable, slice_dir; kwargs...)
    dim = dir_to_idx(slice_dir)
    grid = getdomain(f)
    slider_domain = grid[dim]
    sl = Slider(1:2*10^4:length(slider_domain))
    @debug sl.value
    sl_label = @lift round(grid[dim][$(sl.value)], digits=2)

    return sl, sl_label
end

function section_widget(f, slice_dir=:x; kwargs...)
    fig = Figure(resolution=(800, 600))
    f_approx = approximate_field(f)
    grid = getdomain(f_approx)
    dim = dir_to_idx(slice_dir)
    cmap = RGBAf0.(to_colormap(:viridis, 50), 1.0)
    cmap[1:2] .= RGBAf0(0,0,0,0)

    lvl_slider = Slider(1:20)
    # Contour plot
    ax, plt = default3d_plot!(fig, f_approx, cmap; kwargs...)
    fix_text!(ax)

    # Slider
    sl, sl_label = slider(f_approx, slice_dir; kwargs...)

    section_ax = Axis(fig[1,2])
    # Section through the field
    section_plt = section_plot!(section_ax, f_approx, slice_dir, sl.value; kwargs...)

    add_legend!(fig, section_plt, f; kwargs...)

    # Transparent plane

    plane_origin = [minimum(grid)...]
    ws = maximum(grid) .- plane_origin
    ws[dim] = 0

    p = lift(sl.value) do idx
        autolimits!(section_ax)
        update!(fig.scene)

        new_origin = grid[dim][idx]
        plane_origin[dim] = new_origin
        origin = Point3f0(plane_origin...)
        widths = Point3f0(ws...)
        @debug "Plane with origin at $origin"
        FRect3D(origin, widths)
    end

    a = RGBAf0(0,0,0,0)
    c = RGBAf0(0.2, 0.2, 1.0, 1.0)
    img = AbstractPlotting.ImagePattern([c a; a c])
    mesh!(fig[1, 1], p; color = RGBAf0(0.2, 0.2, 1.0, 0.3), transparency=true)

    # Build app
    app = JSServe.App() do
        dom = md"""
        $(fig.scene)
        $(DOM.div("Slice at: ", sl_label, sl))
        """

        DOM.div(JSServe.TailwindCSS, JSServe.Styling, dom)
    end

end


text:
files = ["gll_library.f90","sem_mesh_mod.f90"]
FF = "gfortran"

for file in files
    fname = join(split(file, ".")[1:end - 1], ".")
    command = `$(FF) -fPIC -shared ./src/$(fname).f90 -o ./lib/$(fname).so`
    run(command)
end

max_stars_repo_name: ziyixi/SeisScripts

max_issues_repo_name: ziyixi/SeisScripts

max_forks_repo_name: ziyixi/SeisScripts


text:
cd(@__DIR__)
load_data("dataset.csv", dlm = ",")

nn1 = load_model("model.jld2")
nn2 = load_model("model.h5")

save_model(nn1)
save_model(nn2)


text:
#
#
#
export gradTV, gradTV2, gradIntrICTV12, forwardLogs
export gradDistance
@doc doc"""
    gradDistance(M,y,x[, p=2])

compute the (sub)gradient of the distance (squared)

```math
f(x) = \frac{1}{2} d^p_{\mathcal M}(x,y)
```

to a fixed [`MPoint`](@ref)` y` on the [`Manifold`](@ref) `M` and `p` is an
integer. The gradient reads

```math
  \nabla f(x) = -d_{\mathcal M}^{p-2}(x,y)\log_xy
```

for $p\neq 1$ or $x\neq  y$. Note that for the remaining case $p=1$,
$x=y$ the function is not differentiable. This function returns then the
[`zeroTVector`](@ref)`(M,x)`, since this is an element of the subdifferential.

# Optional

* `p` – (`2`) the exponent of the distance,  i.e. the default is the squared
  distance
"""
gradDistance(M,y,x,p::Int=2) = (p==2) ? -log(M,x,y) : -distance(M,x,y)^(p-2)*log(M,x,y)

@doc doc"""
    ∇u,∇v = gradIntrICTV12(M,f,u,v,α,β)

compute (sub)gradient of the intrinsic infimal convolution model using the mid point
model of second order differences, see [`costTV2`](@ref), i.e. for some $f\in\mathcal M$
on a [`Power`](@ref) manifold $\mathcal M$ this function computes the (sub)gradient of

```math
E(u,v) =
\frac{1}{2}\sum_{i\in\mathcal G} d_{\mathcal M}(g(\frac{1}{2},v_i,w_i),f_i)
+ \alpha
\bigl(
\beta\mathrm{TV}(v) + (1-\beta)\mathrm{TV}_2(w)
\bigr),
```
where both total variations refer to the intrinsic ones, [`gradTV`](@ref) and
[`gradTV2`](@ref), respectively.
"""
function gradIntrICTV12(M::mT,f::P,u::P,v::P,α,β) where {mT <: Manifold, P <: MPoint}
  c = midPoint(M,u,v,f)
  iL = log(M,c,f)
  return AdjDxGeo(M,u,v,1/2,iL) + α*β*gradTV(M,u), AdjDyGeo(M,u,v,1/2,iL) + α * (1-β) * gradTV2(M,v)
end
@doc doc"""
    gradTV(M,(x,y),[p=1])

compute the (sub) gradient of $\frac{1}{p}d^p_{\mathcal M}(x,y)$ with respect
to both $x$ and $y$.
"""
function gradTV(M::mT where {mT <: Manifold}, xT::Tuple{P,P} where {P <: MPoint}, p::Number=1)
  x = xT[1];
  y = xT[2];
  if p==2
      return (-log(M,x,y), -log(M,y,x))
  else
    d = distance(M,x,y);
    if d==0 # subdifferential containing zero
      return (zeroTVector(M,x),zeroTVector(M,y))
    else
      return (-log(M,x,y)/(d^(2-p)), -log(M,y,x)/(d^(2-p)))
    end
  end
end
@doc doc"""
    ξ = gradTV(M,λ,x,[p])
Compute the (sub)gradient $\partial F$ of all forward differences orrucirng,
in the power manifold array, i.e. of the function

$F(x) = \sum_{i}\sum_{j\in\mathcal I_i} d^p(x_i,x_j)$

where $i$ runs over all indices of the [`Power`](@ref) manifold `M`
and $\mathcal I_i$ denotes the forward neighbors of $i$.

# Input
* `M` – a [`Power`](@ref) manifold
* `x` – a [`PowPoint`](@ref).

# Ouput
* ξ – resulting tangent vector in $T_x\mathcal M$.
"""
function gradTV(M::Power,x::PowPoint,p::Int=1)::PowTVector
  R = CartesianIndices(M.powerSize)
  d = length(M.powerSize)
  maxInd = last(R)
  ξ = zeroTVector(M,x)
  c = costTV(M,x,p,0)
  for i in R # iterate over all pixel
    di = 0.
    for k in 1:d # for all direction combinations
      ek = CartesianIndex(ntuple(i  ->  (i==k) ? 1 : 0, d) ) #k th unit vector
      j = i+ek # compute neighbor
      if all( map(<=, j.I, maxInd.I)) # is this neighbor in range?
        if p != 1
          g = (c[i]==0 ? 1 : 1/c[i]) .* gradTV(M.manifold,(x[i],x[j]),p) # Compute TV on these
        else
          g = gradTV(M.manifold,(x[i],x[j]),p) # Compute TV on these
        end
        ξ[i] += g[1]
        ξ[j] += g[2]
      end
    end # directions
  end # i in R
  return ξ
end

@doc doc"""
    ξ = forwardLogs(M,x)

compute the forward logs $F$ (generalizing forward differences) orrucirng,
in the power manifold array, the function

```math
$F_i(x) = \sum_{j\in\mathcal I_i} \log_{x_i} x_j,\quad i \in \mathcal G,
```

where $\mathcal G$ is the set of indices of the [`Power`](@ref) manifold `M`
and $\mathcal I_i$ denotes the forward neighbors of $i$.

# Input
* `M` – a [`Power`](@ref) manifold
* `x` – a [`PowPoint`](@ref).

# Ouput
* `ξ` – resulting tangent vector in $T_x\mathcal M$ representing the logs, where
  $\mathcal N$ is thw power manifold with the number of dimensions added to `size(x)`.
"""
function forwardLogs(M::Power, x::PowPoint{P,Nt}) where {P <: MPoint, Nt}
  sX = size(x)
  R = CartesianIndices(sX)
  d = length(sX)
  maxInd = [last(R).I...] # maxInd as Array
  if d > 1
    d2 = fill(1,d+1)
    d2[d+1] = d
  else
    d2 = 1
  end
  N = Power(M.manifold, (prod(sX)*d,) )
  xT = PowPoint(repeat(getValue(x),inner=d2))
  ξ = zeroTVector(N,xT)
  for i in R # iterate over all pixel
    for k in 1:d # for all direction combinations
      I = [i.I...] # array of index
      J = I .+ 1 .* (1:d .== k) #i + e_k is j
      if all( J .<= maxInd ) # is this neighbor in range?
        j = CartesianIndex{d}(J...) # neigbbor index as Cartesian Index
        ξ[i,k] = log(M.manifold,x[i],x[j]) # Compute log and store in kth entry
      end
    end # directions
  end # i in R
  return ξ
end

@doc doc"""
    gradTV2(M,(x,y,z),p)

computes the (sub) gradient of $\frac{1}{p}d_2^p(x,y,z)$ with respect
to $x$, $y$, and $z$, where $d_2$ denotes the second order absolute difference
using the mid point model, i.e. let
```math
  \mathcal C = \bigl\{ c\in \mathcal M \ |\ g(\tfrac{1}{2};x_1,x_3) \text{ for some geodesic }g\bigr\}
```
denote the mid points between $x$ and $z$ on the manifold $\mathcal M$.
Then the absolute second order difference is defined as

```math
d_2(x,y,z) = \min_{c\in\mathcal C_{x,z}} d(c,y).
```

While the (sub)gradient with respect to $y$ is easy, the other two require
the evaluation of an [`adjointJacobiField`](@ref).
See [Illustration of the Gradient of a Second Order Difference](@ref secondOrderDifferenceGrad)
for its derivation.
"""
function gradTV2(M::mT where {mT <: Manifold}, xT::Tuple{P,P,P} where {P <: MPoint}, p::Number=1)
  x = xT[1];
  y = xT[2];
  z = xT[3];
  c = midPoint(M,x,z,y) # nearest mid point of x and z to y
  d = distance(M,y,c)
  innerLog = -log(M,c,y)
  if p==2
      return ( AdjDxGeo(M,x,z,1/2,innerLog), -log(M,y,c), AdjDyGeo(M,x,z,1/2,innerLog))
  else
    if d==0 # subdifferential containing zero
      return (zeroTVector(M,x),zeroTVector(M,y),zeroTVector(M,z))
    else
      return ( AdjDxGeo(M,x,z,1/2,innerLog/(d^(2-p))), -log(M,y,c)/(d^(2-p)), AdjDyGeo(M,x,z,1/2,innerLog/(d^(2-p))) )
    end
  end
end
@doc doc"""
    gradTV2(M,x [,p=1])

computes the (sub) gradient of $\frac{1}{p}d_2^p(x_1,x_2,x_3)$
with respect to all $x_1,x_2,x_3$ occuring along any array dimension in the
[`PowPoint`](@ref) `x`, where `M` is the corresponding [`Power`](@ref) manifold.
"""
function gradTV2(M::Power,x::PowPoint,p::Int=1)::PowTVector
  R = CartesianIndices(M.powerSize)
  d = length(M.powerSize)
  minInd, maxInd = first(R), last(R)
  ξ = zeroTVector(M,x)
  c = costTV2(M,x,p,false)
  for i in R # iterate over all pixel
    di = 0.
    for k in 1:d # for all direction combinations (TODO)
      ek = CartesianIndex(ntuple(i  ->  (i==k) ? 1 : 0, d) ) #k th unit vector
      jF = i+ek # compute forward neighbor
      jB = i-ek # compute backward neighbor
      if all( map(<=, jF.I, maxInd.I) ) && all( map(>=, jB.I, minInd.I)) # are neighbors in range?
        if p != 1
          g = (c[i] == 0 ? 1 : 1/c[i]) .* gradTV2(M.manifold,(x[jB],x[i],x[jF]),p) # Compute TV2 on these
        else
          g = gradTV2(M.manifold,(x[jB],x[i],x[jF]),p) # Compute TV2 on these
        end
        ξ[jB] += g[1]
        ξ[i] += g[2]
        ξ[jF] += g[3]
      end
    end # directions
  end # i in R
  return ξ
end


text:
include("../src/auxiliary/data_simulation_fns.jl")
include("../src/bounded_diffusion_domain.jl")
include("../src/radial_ornstein_uhlenbeck.jl")
include("../src/euler_maruyama_dom_restr.jl")
using Random
using Plots

θ = [0.05, √2.0]
Pˣ = RadialOU(θ...)

dt = 1/5000
T = 1.0
tt = 0.0:dt:T

x0 = ℝ{1}(0.5)

N = 100
samples = Vector{Any}(undef, N)
Random.seed!(4)
for i in 1:N
    XX, _ = simulate_segment(0.0, x0, Pˣ, tt)
    samples[i] = XX
end

num_obs = 1000
skip = div(length(tt), num_obs)
Time = collect(tt)[1:skip:end]
path = [x[1] for x in samples[1].yy[1:skip:end]]

p = plot(Time, path, alpha=0.3, color="steelblue", label="")
for i in 2:N
    path = [x[1] for x in samples[i].yy[1:skip:end]]
    plot!(Time, path, alpha=0.3, color="steelblue", label="")
end
p


text:
struct MLMPCSolver <: Solver
    solver
end

function solve(sol::MLMPCSolver, p::NoCrashProblem)
    mdp = NoCrashMDP{typeof(p.rmodel), typeof(p.dmodel.behaviors)}(p.dmodel, p.rmodel, p.discount, p.throw) # make sure an MDP
    return MLMPCPolicy(solve(sol.solver, mdp))
end

struct MLMPCPolicy{P<:Policy} <: Policy
    planner::P
end

srand(p::MLMPCPolicy, s) = srand(p.planner, s)
action_info(p::MLMPCPolicy, b) = action_info(p.planner, most_likely_state(b))
action(p::MLMPCPolicy, b) = first(action_info(p, b))

struct MeanMPCSolver <: Solver
    solver
end

function solve(sol::MeanMPCSolver, p::NoCrashProblem)
    mdp = NoCrashMDP{typeof(p.rmodel), typeof(p.dmodel.behaviors)}(p.dmodel, p.rmodel, p.discount, p.throw) # make sure an MDP
    return MeanMPCPlanner(solve(sol.solver, mdp))
end

struct MeanMPCPlanner{P<:Policy} <: Policy
    planner::P
end

srand(p::MeanMPCPlanner, s) = srand(p.planner, s)
action_info(p::MeanMPCPlanner, b) = action_info(p.planner, mean(b))
action(p::MeanMPCPlanner, b) = first(action_info(p, b))


text:
using Test

using OpenTelemetryAPI
using OpenTelemetrySDK
using OpenTelemetryExporterPrometheus
using HTTP

@testset "PrometheusExporter" begin
    p = MeterProvider()

    m = Meter("my_metrics"; provider = p)
    c = Counter{Int}("fruit_counter", m)

    c(; name = "apple", color = "red")
    c(2; name = "lemon", color = "yellow")
    c(1; name = "lemon", color = "yellow")
    c(2; name = "apple", color = "green")
    c(5; name = "apple", color = "red")
    c(4; name = "lemon", color = "yellow")

    e = PrometheusExporter(; host = "0.0.0.0", port = 9966)
    r = MetricReader(p, e)

    resp = HTTP.request("GET", "http://localhost:9966")
    body = split(String(resp.body), "\n"; keepempty = false)
    expected = split(
        """
        # HELP fruit_counter
        # TYPE fruit_counter counter
        fruit_counter{color="red",name="apple"} 6
        fruit_counter{color="green",name="apple"} 2
        fruit_counter{color="yellow",name="lemon"} 7
        """,
        "\n";
        keepempty = false,
    )

    for (expect, real) in zip(sort(expected), sort(body))
        @test startswith(real, expect)
    end
end