text:
# This file contains utilities to facilities the use of
# DataFrames.jl for generating benchmark settings
# The idea is to benchmark different "features" and obtain
# "measurements". Dependent features are supported.

using FileIO, DataFrames;

@doc """
    generate_all_settings( features_dict, dependent_features =  nothing )

    Return a DataFrame, the columns of which correspond to every possible
    combination of feature values.

    ## Arguments

    * `features_dict` is a Dict or NamedTuple. The keys must be Strings or
      Symbols and provide a feature name. This is also the column name in the
      returned DataFrame. The values must be Vectors of possible feature values.
      E.g. if feature `"x"` can take values `[1,2,3]` provide an entry
      `"x" => [1,2,3]`.
    * `dependent_features` (optional). A Vector of Dicts or NamedTuples, one
      for each **dependent** feature. A dependent feature depends on some features
      from `features_dict`; we call these **arguments**.
      Each `dependent_feature` is defined by an entry in `dependent_features` with
      keys `"name"`, `"depends_on"` and `"values"`.
      The value of `"values"` must either be a function that takes as keyword arguments
      the arguments and returns a vector of possible values.\n
      It can also be a Dict/NamedTuple with keys being Tuples corresponding to all
      possible value combinations of features in `depends_on` (in the order given there).

    ## Example

    ```jldoctest
    features = Dict(
        :X => [1,],
        :Y => [1,2,]
    );
    dep_features = [
        Dict(
            "name" => :Z,
            "depends_on" => [:X, :Y],
            "values" => function(;X,Y) fill( X+Y, X+Y ) end
        ),
    ];
    s = generate_all_settings(features, dep_features)

    # output

    5×3 DataFrame
     Row │ Y      X      Z
         │ Int64  Int64  Int64
    ─────┼─────────────────────
       1 │     1      1      2
       2 │     1      1      2
       3 │     2      1      3
       4 │     2      1      3
       5 │     2      1      3
    ```
"""
function generate_all_settings( features_dict :: Union{Dict, NamedTuple},
        dependent_features :: Union{Nothing,Vector{D}} where D<:Union{Dict, NamedTuple} = nothing
    )

    # turn entries of `features_dict` into Dict with key=feature name and
    # value = single column DataFrame of values
    simple_dfs = Dict{Union{String},DataFrame}()
    for (feature, feature_values) ∈ pairs(features_dict)
        simple_dfs[string(feature)] = DataFrame( feature => feature_values )
    end

    # combine to first settings DataFrame
    # rows = all possible combinations of feature values from `features_dict`
    settings_df = crossjoin( (simple_dfs[string(k)] for k in keys(features_dict))... )

    if !isnothing( dependent_features )
        for dependent_feature ∈ dependent_features
            # for each possible value combination of feature values of
            # features in `dependent_feature["depends_on"]` put a dataframe
            # into `dependent_dfs`; these are later stacked vertically
            dependent_dfs = DataFrame[];
            if length( dependent_feature["depends_on"] ) > 1
                input_args = crossjoin(
                    (simple_dfs[ string(dependency) ]
                        for dependency ∈ dependent_feature["depends_on"] )...
                );
            else
                input_args = simple_dfs[ string(dependent_feature["depends_on"][1]) ]
            end

            for arg_row = eachrow( input_args )
                # a 1-col DF with colname=dependent_feature["name"] and values=f(;arg_row...)
                row_df = DataFrame(
                    dependent_feature["name"] => let dp_vals = dependent_feature["values"];
                        if isa(dp_vals, Function)
                            dp_vals(; arg_row... )
                        elseif isa( dp_vals, Union{Dict,NamedTuple} )
                            dp_vals[ Tuple(arg_row) ];
                        else
                            error(
                                "Field `values` of dependent feature $(dependent_features["name"])
                                must be a function with kwargs or a Dict/NamedTuple with keys corresponding
                                to all possible value combinations of $(dependent_feature["depends_on"])."
                            );
                        end#if
                    end#let
                );
                # add redundant columns containing input_args information
                # and colnames for joining with settings_df
                push!( dependent_dfs, crossjoin( DataFrame(arg_row), row_df ) )
            end

            settings_df = innerjoin( settings_df, vcat( dependent_dfs... );
                on = dependent_feature["depends_on"]
            );
        end
    end
    return settings_df;
end

#=
features = Dict(
    :X => [1,],
    :Y => [1,2,]
);
dep_features = [
    Dict(
        "name" => :Z,
        "depends_on" => [:X, :Y],
        "values" => function(;X,Y) fill( X+Y, X+Y ) end
    ),
];
s = generate_all_settings(features, dep_features)
=#

function feature_names( features_dict :: Union{Dict, NamedTuple},
        dependent_features :: Union{Nothing,Vector{D}} where D<:Union{Dict, NamedTuple} = nothing
    )
    return [
        string.(collect( keys( features_dict ) ));
        [ string(dep_feature["name"]) for dep_feature ∈ dependent_features ]
    ]
end

#%%
function unpack_feature( feature_name :: String, some_dict )
    feature_symbol = Symbol( feature_name );
    feature_symbol_symbol = :(Symbol( $feature_name ));
    @eval $feature_symbol = getindex($some_dict, $feature_symbol_symbol);
end

#%%
@doc """
    add_observation_columns!(df, observations ::Union{Dict, NamedTuple} )

Add one or several new empty column(s) to DataFrame `df` or overwrite if exists.
Columns are specified by `observations` where each key gives a column name and
each value the corresponding data type.
Each column will have then have this data type in union with `Nothing`.
"""
function add_observation_columns!(df :: DataFrame, observations ::Union{Dict, NamedTuple} )
    n_rows = size(df, 1);
    for (obs_name, obs_type) ∈ pairs(observations)
        df[!, string(obs_name)] = Vector{Union{obs_type,Missings.Missing}}(Missings.missing, n_rows)
    end
    nothing
end

#=
observations = (;
    :ω => Float64,
)
add_observation_columns!( s, observations );
=#

#%%
function load_previous_results( filename; result_key = "results" )
    file_data = load( filename );

    if !haskey( file_data, result_key )
        error("Cannot retrieve previous data.");
    end

    return file_data[result_key];
end

function scan_feature_values( df :: DataFrame, feature_name :: Union{Symbol, String})
    if feature_name in names(df)
        return unique( df[ !, feature_name ] );
    else
        @warn "No column $(feature_name) found in DataFrame."
        return []
    end
end

#%%
function scan_dependent_feature_values( df :: DataFrame, feature_name :: Union{Symbol, String},
        depends_on :: Union{S, Vector{S}} where S<:Union{Symbol, String}
    )

    feature_name = string(feature_name);

    if feature_name in names(df)
        dep_args = isa( depends_on, Union{Symbol, String} ) ? [ string( depends_on ), ] : string.(depends_on);
        new_dep_dict = Dict()
        for sub_df in groupby( df, dep_args )
            new_dep_dict[ Tuple( sub_df[1, dep_args]) ] = unique( sub_df[:, feature_name ] );
        end
        return new_dep_dict
    else
        @warn "No column $(feature_name) found in DataFrame."
        return []
    end
end

function scan_dependent_feature_values(df :: DataFrame, feature :: Union{Dict, NamedTuple})
    return scan_dependent_feature_values( feature["name"], feature["depends_on"] )
end

#=
z_vals = scan_dependent_feature_values( s, "Z", ["X","Y"])
=#

#%%
function find_observation_rows( df :: DataFrame,
        observations :: Vector{<:Symbol} = Symbol[]
    )
    return findall(
        .!(vec(
            prod( hcat( ( ismissing.( df[:,obs] ) for obs in observations)... ), dims = 2 )
        ))
    );
end

find_observation_rows(df :: DataFrame, obs :: Union{Dict, NamedTuple}) = find_observation_rows( df, collect(keys(obs)) )

function fill_from_partial_results!(target :: DataFrame, source :: DataFrame,
        feature_names :: Vector{ S } where S<:Union{Symbol, String},
        observations :: Union{Nothing, Union{Dict, NamedTuple} } = nothing,
        required_observations :: Union{Nothing, Union{Dict, NamedTuple} } = nothing;
    )

    if isnothing( observations )
        observations = intersect( setdiff( names(target), string.(feature_names) ), names( source ) )
    else
        observations = collect(keys(observations));
    end

    if !isnothing( observations )
        if isnothing( required_observations )
            required_observations = observations;
        end
        investigated_row_indices = Int[];

        non_missing_indices = find_observation_rows( source, required_observations );

        for non_missing_index in non_missing_indices
            src_row = source[ non_missing_index, : ];
            for target_row ∈ eachrow( target )
                if Tuple(target_row[feature_names]) == Tuple(src_row[feature_names])
                    target_row[observations] = src_row[observations];
                    push!( investigated_row_indices, rownumber(target_row) )
                    # break here if we expect/assume unique features
                end
            end
        end
        to_do_row_indices = setdiff( 1 : size(target,1), investigated_row_indices )
        println("Filled in data from $(length(investigated_row_indices)) rows.")
        return investigated_row_indices, to_do_row_indices
    end

    return nothing
end

function fill_from_partial_results!( target :: DataFrame, source :: DataFrame,
        features_dict :: Union{Dict, NamedTuple},
        dependent_features :: Union{Nothing,Vector{D}} where D<:Union{Dict, NamedTuple} = nothing,
        observations :: Union{Nothing, Union{Dict, NamedTuple} } = nothing
    )
    feature_names = collect(keys( features_dict ));
    if !isnothing( dependent_features )
        push!( feature_names, [dep_feat["name"] for dep_feat ∈ dependent_features]... )
    end
    return fill_from_partial_results!( target, source, feature_names, observations)
end
#=
new_df = DataFrame( X=[1,1,2,2], Y=[1,2,1,2] )
add_observation_columns!( new_df, (; :Z => Int,))
src = DataFrame( X = [1,2], Y= [2,2], Z = [10,11])
fill_from_partial_results!( new_df, src, ["X", "Y"] )
=#


text:
using Ipopt
using DelimitedFiles
using PyPlot
import ForwardDiff

using CSV
using DataFrames
# using YAML

# set up boundary constraint wrapper function
function boundary_wrapper(x, params)
    # include relevant globals
    params.boundary_vertices_a
    params.boundary_normals_a
    params.boundary_vertices_b
    params.boundary_normals_b
    params.boundary_vertices_c
    params.boundary_normals_c
    params.boundary_vertices_d
    params.boundary_normals_d
    params.boundary_vertices_e
    params.boundary_normals_e

    # get number of turbines
    nturbines = Int(length(x)/2)

    # extract x and y locations of turbines from design variables vector
    turbine_x = x[1:nturbines]
    turbine_y = x[nturbines+1:end]

    boundcon_a = ff.ray_trace_boundary(boundary_vertices_a, boundary_normals_a, turbine_x[1:31], turbine_y[1:31])
    boundcon_b = ff.ray_trace_boundary(boundary_vertices_b, boundary_normals_b, turbine_x[32:42], turbine_y[32:42])
    boundcon_c = ff.ray_trace_boundary(boundary_vertices_c, boundary_normals_c, turbine_x[43:58], turbine_y[43:58])
    boundcon_d = ff.ray_trace_boundary(boundary_vertices_d, boundary_normals_d, turbine_x[59:72], turbine_y[59:72])
    boundcon_e = ff.ray_trace_boundary(boundary_vertices_e, boundary_normals_e, turbine_x[73:81], turbine_y[73:81])

    # get and return boundary distances
    return [boundcon_a; boundcon_b; boundcon_c; boundcon_d; boundcon_e]
end

# set up spacing constraint wrapper function
function spacing_wrapper(x, params)
    # include relevant globals
    params.rotor_diameter

    # get number of turbines
    nturbines = Int(length(x)/2)

    # extract x and y locations of turbines from design variables vector
    turbine_x = x[1:nturbines]
    turbine_y = x[nturbines+1:end]

    spacecon_a = 2.0*rotor_diameter[1] .- ff.turbine_spacing(turbine_x[1:31], turbine_y[1:31])
    spacecon_b = 2.0*rotor_diameter[1] .- ff.turbine_spacing(turbine_x[32:42], turbine_y[32:42])
    spacecon_c = 2.0*rotor_diameter[1] .- ff.turbine_spacing(turbine_x[43:58], turbine_y[43:58])
    spacecon_d = 2.0*rotor_diameter[1] .- ff.turbine_spacing(turbine_x[59:72], turbine_y[59:72])
    spacecon_e = 2.0*rotor_diameter[1] .- ff.turbine_spacing(turbine_x[73:81], turbine_y[73:81])

    # get and return spacing distances
    return [spacecon_a; spacecon_b; spacecon_c; spacecon_d; spacecon_e]
end

# set up objective wrapper function
function aep_wrapper(x, params)
    # include relevant globals
    params.turbine_z
    params.rotor_diameter
    params.hub_height
    params.turbine_yaw
    params.ct_model
    params.generator_efficiency
    params.cut_in_speed
    params.cut_out_speed
    params.rated_speed
    params.rated_power
    params.windresource
    params.power_models
    params.model_set
    params.rotor_points_y
    params.rotor_points_z
    params.obj_scale

    # get number of turbines
    nturbines = Int(length(x)/2)

    # extract x and y locations of turbines from design variables vector
    turbine_x = x[1:nturbines]
    turbine_y = x[nturbines+1:end]

    # calculate AEP
    AEP = obj_scale*ff.calculate_aep(turbine_x, turbine_y, turbine_z, rotor_diameter,
                hub_height, turbine_yaw, ct_model, generator_efficiency, cut_in_speed,
                cut_out_speed, rated_speed, rated_power, windresource, power_models, model_set,
                rotor_sample_points_y=rotor_points_y,rotor_sample_points_z=rotor_points_z)

    # return the objective as an array
    return [AEP]
end

# objective function
function obj(x)
    global funcalls_AEP
    AEP = -aep_wrapper(x)[1]
    append!(funcalls_AEP, -AEP)
    return AEP
end

# constraint function
function con(x, g)
    # calculate spacing constraint value
    spacing_con = spacing_wrapper(x)

    # calculate boundary constraint
    boundary_con = boundary_wrapper(x)

    # combine spacing and boundary constraint values
    g[:] = [spacing_con; boundary_con]
end

# objective gradient function
function obj_grad(x, grad_f)
    grad_f[:] = -ForwardDiff.jacobian(aep_wrapper,x)
end

# constraint gradients function
function con_grad(x, mode, rows, cols, values)
    if mode == :Structure
        # report the sparcity structure of the jacobian
        for i = 1:prob.m
            rows[(i-1)*prob.n+1:i*prob.n] = i*ones(Int,prob.n)
            cols[(i-1)*prob.n+1:i*prob.n] = 1:prob.n
        end
        return rows, cols
    else
        # calculate spacing constraint jacobian
        ds_dx = ForwardDiff.jacobian(spacing_wrapper, x)

        # calculate boundary constraint jacobian
        db_dx = ForwardDiff.jacobian(boundary_wrapper, x)

        # combine constraint jacobians into an overall array of constraint derivatives
        for i = 1:prob.m
            for j = 1:prob.n
                values[(i-1)*prob.n+j] = [ds_dx; db_dx][i, j]
            end
        end
    end
    return values
end

# this function will be called at each iteration of the optimization
function intermediate(
    alg_mod::Int,
    iter_count::Int,
    obj_value::Float64,
    inf_pr::Float64, inf_du::Float64,
    mu::Float64, d_norm::Float64,
    regularization_size::Float64,
    alpha_du::Float64, alpha_pr::Float64,
    ls_trials::Int)
    global iter_AEP
    append!(iter_AEP, -obj_value)

    return true
end

# set globals for iteration history
iter_AEP = zeros(Float64, 0)
funcalls_AEP = zeros(Float64, 0)
global iter_AEP
global funcalls_AEP

# import model set with wind farm and related details
include("./model_sets/model_set_7_ieacs4_reduced_wind_rose.jl")

# scale objective to be between 0 and 1
obj_scale = 1E-12

# set wind farm boundary parameters
boundary_vertices_a = [10363.8 6490.3; 9449.7 1602.2; 9387.0 1056.6; 9365.1 625.5; 9360.8 360.2; 9361.5 126.9; 9361.3 137.1; 7997.6 1457.9; 6098.3 3297.5;
    8450.3 6455.3; 8505.4 6422.3; 9133.0 6127.4; 9332.8 6072.6; 9544.2 6087.1; 9739.0 6171.2; 9894.9 6316.9; 10071.8 6552.5; 10106.9 6611.1]
boundary_normals_a = [0.9829601758936983 -0.1838186405319916; 0.9934614633172962 -0.11416795042154541; 0.9987121579438882 -0.050734855622757584;
    0.9998686751666075 -0.01620593781838486; 0.9999954987444023 0.0030004151269687495; -0.9998078216567232 -0.019604074934516894; -0.6957179389375846 -0.718315076718037;
    -0.6957275377423737 -0.7183057797532565; -0.8019887481131871 0.5973391397688945; 0.5138086803485797 0.8579047965820281; 0.4252760929807897 0.905063668886888;
    0.2645057513093967 0.9643841078762402; -0.0684295708121141 0.9976559496331737; -0.39636379138742883 0.9180935381958544; -0.6828023205475376 0.7306031693435896;
    -0.7996740386176392 0.6004343694034798; -0.8578802011411015 0.5138497450520954; 0.42552559023380465 0.9049463918134445]
boundary_vertices_b = [5588.4 3791.3; 4670.7 4680.2; 7274.9 7940.8; 7369.9 7896.2; 7455.1 7784.3; 7606.5 7713.0; 7638.9 7708.4; 8297.1 7398.9]
boundary_normals_b = [-0.6957460043611584 -0.7182878931288504; -0.7813688797257963 0.6240694462926818; 0.4249708760634733 0.9052070230051488; 0.7956275395848184 0.6057861159305391;
    0.4260560153872896 0.9046967844268629; 0.14056568619461773 0.9900713549359138; 0.4255255464063141 0.9049464124220882; 0.7996806883794807 -0.6004255129763556]
boundary_vertices_c = [3267.1 10100.6; 4164.1 9586.6; 5749.8 9068.6; 6054.7 8925.3; 1468.5 7781.7; 107.4 9100.0]
boundary_normals_c = [0.49718026396417986 0.8676472699919642; 0.31052117525343714 0.9505664625470563; 0.42535384615162936 0.9050271297392228; 0.24194817066179167 -0.9702891747893577;
    -0.6957228969594285 -0.7183102746351193; -0.30189947425802094 0.9533397649540959]
boundary_vertices_d = [6764.9 8399.7; 4176.8 5158.6; 2047.8 7220.7]
boundary_normals_d = [0.7814306689309158 -0.6239920749930895; -0.6957310325444781 -0.7183023947855072; -0.24248239299288069 0.9701558066045093]
boundary_vertices_e = [8953.7 11901.5; 7048.3 9531.5; 6127.7 9962.7; 4578.1 10464.9; 4524.1 10498.7]
boundary_normals_e = [0.7793586677376737 -0.6265780613955122; -0.4241667101838764 -0.9055841219742026; -0.30829751674447764 -0.9512899879475178; -0.5305632140423848 -0.847645371546978; -0.3019099610801309 0.9533364439695956]

# set globals for use in wrapper functions
struct params_struct{}
    model_set
    rotor_points_y
    rotor_points_z
    turbine_z
    ambient_ti
    rotor_diameter
    boundary_vertices_a
    boundary_normals_a
    boundary_vertices_b
    boundary_normals_b
    boundary_vertices_c
    boundary_normals_c
    boundary_vertices_d
    boundary_normals_d
    boundary_vertices_e
    boundary_normals_e
    obj_scale
    hub_height
    turbine_yaw
    ct_model
    generator_efficiency
    cut_in_speed
    cut_out_speed
    rated_speed
    rated_power
    windresource
    power_models
end

params = params_struct(model_set, rotor_points_y, rotor_points_z, turbine_z, ambient_ti,
    rotor_diameter, boundary_vertices_a, boundary_normals_a, boundary_vertices_b, boundary_normals_b,
    boundary_vertices_c, boundary_normals_c, boundary_vertices_d, boundary_normals_d, boundary_vertices_e,
    boundary_normals_e, obj_scale, hub_height, turbine_yaw, ct_model, generator_efficiency, cut_in_speed,
    cut_out_speed, rated_speed, rated_power, windresource, power_models)

# initialize design variable array
x_initial = [copy(turbine_x);copy(turbine_y)]

# add initial turbine location to plot
plt.clf()
for i = 1:length(turbine_x)
    plt.gcf().gca().add_artist(plt.Circle((turbine_x[i],turbine_y[i]), rotor_diameter[1]/2.0, fill=false,color="C0"))
end

# get number of design variables
n_designvariables = length(x_initial)

# get number of constraints
n_spacingconstraints = length(spacing_wrapper(x_initial, params))
n_boundaryconstraints = length(boundary_wrapper(x_initial, params))
n_constraints = n_spacingconstraints + n_boundaryconstraints

# set general lower and upper bounds for design variables
lb = ones(n_designvariables) * -Inf
ub = ones(n_designvariables) * Inf

# set lower and upper bounds for constraints
lb_g = ones(n_constraints) * -Inf
ub_g = zeros(n_constraints)

# create the problem
prob = createProblem(n_designvariables, lb, ub, n_constraints, lb_g, ub_g, n_designvariables*n_constraints, 0,
    obj, con, obj_grad, con_grad)
addOption(prob, "hessian_approximation", "limited-memory")
prob.x = deepcopy(x_initial)

# set intermediate function to be ran at every iteration
setIntermediateCallback(prob, intermediate)

# generate wrapper function surrogates
spacing_wrapper(x) = spacing_wrapper(x, params)
aep_wrapper(x) = aep_wrapper(x, params)
boundary_wrapper(x) = boundary_wrapper(x, params)

# print initial objective value
println("nturbines: ", nturbines)
println("rotor diameter: ", rotor_diameter[1])
println("starting AEP value (MWh): ", aep_wrapper(x_initial)[1]*1E5)

# run and time optimization
t1 = time()
status = solveProblem(prob)
t2 = time()
clkt = t2-t1
xopt = prob.x
fopt = prob.obj_val
info = Ipopt.ApplicationReturnStatus[status]

# print optimization results
println("Finished in : ", clkt, " (s)")
println("info: ", info)
println("end AEP value (MWh): ", aep_wrapper(xopt)[1]*1E5)

# extract final turbine locations
turbine_x = copy(xopt[1:nturbines])
turbine_y = copy(xopt[nturbines+1:end])

# add final turbine locations to plot
for i = 1:length(turbine_x)
    plt.gcf().gca().add_artist(plt.Circle((turbine_x[i],turbine_y[i]), rotor_diameter[1]/2.0, fill=false,color="C1", linestyle="--"))
end

# add wind farm boundary to plot
plt.gcf().gca().plot([boundary_vertices_a[:,1];boundary_vertices_a[1,1]],[boundary_vertices_a[:,2];boundary_vertices_a[1,2]], color="C2")
plt.gcf().gca().plot([boundary_vertices_b[:,1];boundary_vertices_b[1,1]],[boundary_vertices_b[:,2];boundary_vertices_b[1,2]], color="C2")
plt.gcf().gca().plot([boundary_vertices_c[:,1];boundary_vertices_c[1,1]],[boundary_vertices_c[:,2];boundary_vertices_c[1,2]], color="C2")
plt.gcf().gca().plot([boundary_vertices_d[:,1];boundary_vertices_d[1,1]],[boundary_vertices_d[:,2];boundary_vertices_d[1,2]], color="C2")
plt.gcf().gca().plot([boundary_vertices_e[:,1];boundary_vertices_e[1,1]],[boundary_vertices_e[:,2];boundary_vertices_e[1,2]], color="C2")

# set up and show turbine layout plot
axis("square")
xlim(0, 11000)
ylim(-500, 13000)
plt.show()
savefig("ipoptresults_turbinelayout.png")

# create obj value history plots
plt.clf()
plt.gcf().gca().plot(1:length(funcalls_AEP), funcalls_AEP)
title("Objective Value History")
xlabel("Function Call")
ylabel("AEP (TWh)")
savefig("ipoptresults_plot_AEP_history_by_funcall.png")

plt.clf()
plt.gcf().gca().plot(1:length(iter_AEP), iter_AEP)
title("Objective Value at Each Iteration")
xlabel("Iteration")
ylabel("AEP (TWh)")
savefig("ipoptresults_plot_AEP_history_by_iter.png")

# write text file with optimization summary info
io = open("ipoptresults_summary.txt", "a")
write(io, "number of turbines: ")
writedlm(io, nturbines)
write(io, "rotor diameter: ")
writedlm(io, rotor_diameter[1])
write(io, "starting AEP value (MWh): ")
writedlm(io, aep_wrapper(x_initial)[1]*1E5)
write(io, "Finished in : ")
writedlm(io, clkt)
write(io, "end objective value: ")
writedlm(io, aep_wrapper(xopt)[1]*1E5)
write(io, "\noptimal turbine coordinates:\n")
writedlm(io, [turbine_x turbine_y])
write(io, "\neach function call value\n")
writedlm(io, funcalls_AEP)
close(io)

# write results to csv files
dataforcsv_xopt = DataFrame(xopt = xopt)
dataforcsv_funceval = DataFrame(function_value = funcalls_AEP)
dataforcsv_funceval_at_iteration = DataFrame(function_value = iter_AEP)
CSV.write("cs4_xopt.csv", dataforcsv_xopt)
CSV.write("cs4_functionvalue_log.csv", dataforcsv_funceval)
CSV.write("cs4_functionvalue_at_iteration_log.csv", dataforcsv_funceval_at_iteration)

# # write results to a YAML file (not complete)
# dataforyaml = [
#     "optimal_turbine_coordinates" => [turbine_x turbine_y]
#     "function_evaluation_log" => funcalls_AEP
#     "function_evaluation_at_iterations_log" => iter_AEP
#     ]
# YAML.write_file("cs4_results.yml", dataforyaml)

ext: jl


text:
struct CrankNicholson <: Propagator
    dt::Float64
    H_mid::Hamiltonian
end

function (U::CrankNicholson)(ψ::Ket)
    Δ = U.H_mid * (im/2.0) * U.dt
    lhs = I .+ Δ
    rhs = (I .- Δ) * ψ
    lhs / rhs # solves the linear equation lhs*X = rhs
end


text:
using StatsBase
using RDatasets
using Distributions
d = Dirichlet([1.0, 3.0, 5.0])
arr=rand(d)
entropy(arr,2)
entropy(arr)

max_stars_repo_path: Chapter 3/6.jl

max_stars_repo_name: PacktPublishing/Julia-for-data-science

max_issues_repo_path: Chapter 3/6.jl

max_issues_repo_name: PacktPublishing/Julia-for-data-science

max_forks_repo_path: Chapter 3/6.jl

max_forks_repo_name: PacktPublishing/Julia-for-data-science


text:
@testset "Inversion of second and fourth order tensors" begin
    for T in (Float32, Float64), dim in 1:4
        Random.seed!(1234)
        # second order
        x = rand(SecondOrderTensor{dim, T})
        y = rand(SymmetricSecondOrderTensor{dim, T})
        @test (@inferred inv(x))::typeof(x) ⋅ x ≈ one(x)
        @test (@inferred inv(y))::typeof(y) ⋅ y ≈ one(y)
        # fourth order
        dim > 3 && continue
        x = rand(FourthOrderTensor{dim, T})
        y = rand(SymmetricFourthOrderTensor{dim, T})
        @test (@inferred inv(x))::typeof(x) ⊡ x ≈ one(x)
        @test (@inferred inv(y))::typeof(y) ⊡ y ≈ one(y)
    end
    @testset "Fast version using block matrix algorithm" begin
        for T in (Float32, Float64), dim in 1:11
            Random.seed!(1234)
            x = rand(SecondOrderTensor{dim, T})
            y = rand(SymmetricSecondOrderTensor{dim, T})
            if T == Float64
                @test (@inferred Tensorial.fastinv(x))::typeof(x) ⋅ x ≈ one(x)
                @test (@inferred Tensorial.fastinv(y))::typeof(y) ⋅ y ≈ one(y)
            else
                @test_throws Exception Tensorial.fastinv(x)
                @test_throws Exception Tensorial.fastinv(y)
            end
        end
    end
end

@testset "Solving linear system" begin
    for T in (Float32, Float64), dim in 1:4
        Random.seed!(1234)
        A = rand(SecondOrderTensor{dim, T})
        S = rand(SymmetricSecondOrderTensor{dim, T})
        b = rand(Vec{dim, T})
        @test (@inferred A \ b)::Vec{dim, T} ≈ Array(A) \ Array(b)
        @test (@inferred S \ b)::Vec{dim, T} ≈ Array(S) \ Array(b)
    end
    @testset "Fast version using block matrix algorithm" begin
        for T in (Float32, Float64), dim in 1:11
            Random.seed!(1234)
            A = rand(SecondOrderTensor{dim, T})
            S = rand(SymmetricSecondOrderTensor{dim, T})
            b = rand(Vec{dim, T})
            if T == Float64
                @test (@inferred Tensorial.fastsolve(A, b))::Vec{dim, T} ≈ Array(A) \ Array(b)
                @test (@inferred Tensorial.fastsolve(S, b))::Vec{dim, T} ≈ Array(S) \ Array(b)
            else
                @test_throws Exception Tensorial.fastsolve(A, b)
                @test_throws Exception Tensorial.fastsolve(S, b)
            end
        end
    end
end


text:
using DataDeps

register(DataDep(
  "computer-hardware",
  "https://archive.ics.uci.edu/ml/datasets/Computer+Hardware",
  "https://archive.ics.uci.edu/ml/machine-learning-databases/cpu-performance/machine.data",
  "c95f0e52fb8381a14a24aa7e9f4b50f72a08639ecfd8cc619ab779919dbd448b",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=10,
        feature_indices=[1; 3:9],
        id_indices=2,
        categoric_indices=1,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
using Zygote
using RayTracer: get_params, set_params!

function loss_fn(θ, color)
    rendered_color = raytrace(origin, direction, θ, light, eye_pos, 0)
    loss = sum(abs2.(rendered_color.x .- color.x) .+
               abs2.(rendered_color.y .- color.y) .+
               abs2.(rendered_color.z .- color.z))
    return loss
end

# Define the Scene Parameters

screen_size = (w = 400, h = 300)

light = PointLight(
    Vec3(1.0),
    1000.0,
    Vec3(0.15, 0.5, -110.5)
)

eye_pos = Vec3(0.0, 0.0, -5.0)

cam = Camera(eye_pos, Vec3(0.0), Vec3(0.0, 1.0, 0.0), 45.0, 1.0, screen_size...)

origin, direction = get_primary_rays(cam);

@testset "Triangle" begin
    scene = [
        Triangle(Vec3(-1.7, 1.0, 0.0), Vec3(1.0, 1.0, 0.0), Vec3(1.0, -1.0, 0.0),
                 Material(color_ambient = Vec3(1.0),
                          color_diffuse = Vec3(1.0),
                          color_specular = Vec3(1.0),
                          specular_exponent = 50.0,
                          reflection = 1.0))
    ]

    scene_new = [
        Triangle(Vec3(-1.9, 1.3, 0.1), Vec3(1.2, 1.1, 0.3), Vec3(0.8, -1.2, -0.15),
                 Material(color_ambient = Vec3(1.0),
                          color_diffuse = Vec3(1.0),
                          color_specular = Vec3(1.0),
                          specular_exponent = 50.0,
                          reflection = 1.0))
    ]

    color = raytrace(origin, direction, scene, light, eye_pos, 0)

    zygote_grads = get_params(gradient(x -> loss_fn(x, color), scene_new)[1][1])

    numerical_grads = get_params(numderiv(x -> loss_fn([x], color), scene_new[1]))

    @test isapprox(numerical_grads, zygote_grads, rtol = 1.0e-3)
end

@testset "Sphere" begin
    scene = [
        Sphere(Vec3(-1.7, 1.0, 0.0), [0.6],
               Material(color_ambient = Vec3(1.0),
                        color_diffuse = Vec3(1.0),
                        color_specular = Vec3(1.0),
                        specular_exponent = 50.0,
                        reflection = 1.0))
    ]

    scene_new = [
        Sphere(Vec3(-1.9, 1.3, 0.1), [0.9],
	       Material(color_ambient = Vec3(1.0),
                        color_diffuse = Vec3(1.0),
                        color_specular = Vec3(1.0),
                        specular_exponent = 50.0,
                        reflection = 1.0))
    ]

    color = raytrace(origin, direction, scene, light, eye_pos, 0)

    zygote_grads = get_params(gradient(x -> loss_fn(x, color), scene_new)[1][1])

    numerical_grads = get_params(numderiv(x -> loss_fn([x], color), scene_new[1]))

    @test isapprox(numerical_grads[1:4], zygote_grads[1:4], rtol = 1.0e-3)
end


text:
################################################################################
# This part is from Calus et al, JDS, 2016, 99:6403-6419
################################################################################

"""
    function DRP()
---
# Some terms
- DYD: Daughter Yield Deviation, VanRaden and Wiggans, 1991
- DRP: DeRegressed proofs
- EDC: Effective Daughter contributions, Fikse and Banos, 2001.

This DRP is kind of matrix deregression.
"""
function DRP(W, A, ng, nu, y, λ; conv=1.)
    # A is the inverse of thee A matrix.
    # R is the weight of DRP
    # l = 1:ng
    # r = ng+1:ng+nu
    #
    # lhs = [ng          ones(ng)'W   zeros(nu)'
    #        R*ones(ng)  W+A[l, l].*λ A[l, r].*λ
    #        zeros(r, r) A[r, l].*λ   A[r, r].*λ]
    # rhs = [ ones(ng)'W*y
    #         W*y
    #         zeros(nu) ]
    # b = lhs \ rhs
    μ = 0
end

function ERC(λ, Rel)
    λ * Rel / (1 - Rel)
end

function EDC(h², Rel)
    λₛ = (4 - h²) / h²
    λₛ * Rel / (1 - Rel)
end


text:
## JDF

@info "JDF.jl is available - including functionality to read / write JDF files"

using .JDF
using DataFrames

function read_table(::TableIOInterface.JDFFormat, filename:: AbstractString; kwargs...)
    return JDF.load(filename; kwargs...)
end

function write_table!(::TableIOInterface.JDFFormat, filename:: AbstractString, table:: DataFrame; kwargs...)
    JDF.save(filename, table; kwargs...)
    nothing
end

# JDF supports only DataFrames, not arbitrary Tables.jl inputs. For export, the table is converted to a DataFrame first.
write_table!(::TableIOInterface.JDFFormat, filename:: AbstractString, table; kwargs...) = write_table!(TableIOInterface.JDFFormat(), filename, DataFrame(table); kwargs...)


text:
using ..Continuations: ContinuationOptions
import ..Continuations: init

init(prob::Codim2Problem; kwargs...) =
    BifurcationSolver(prob, ContinuationOptions(; kwargs...))


text:
function save_parameters_to_file( p::Param , fid::HDF5.File  )

    #Time variables
    # -> saved at end of simulation

    #Spatial variables
    HDF5.attributes(fid)["rmin"]      = p.rmin
    HDF5.attributes(fid)["rmax"]      = p.rmax
    HDF5.attributes(fid)["rnodes"]    = p.rnodes

    HDF5.attributes(fid)["thnodes"]    = p.thnodes

    #Initial condition
    HDF5.attributes(fid)["A0"]      = p.A0
    HDF5.attributes(fid)["sigma"]   = p.σ
    HDF5.attributes(fid)["r0"]      = p.r0
    HDF5.attributes(fid)["omega"]   = p.ω
    HDF5.attributes(fid)["m_mode"]  = p.m

     #Cylinder data
     HDF5.attributes(fid)["epsilon"]    = p.ϵ
     HDF5.attributes(fid)["Omega"]      = p.Ω
     HDF5.attributes(fid)["Z"]          = p.Z0

end


text:
module CTIO

Base.Experimental.@optlevel 0

export load_data, sinog_extract, sinog_load
export load_yaml, write_yaml
export write_ct_image, read_ct_image
export load_image, write_image
export load_sinogram, write_sinogram
export load_tomogram, write_tomogram

import YAML, Mmap
using IterTools: imap
using ..Monads
using ..CTImages
using ..Utils: yaml_repr
using ..TypeDict: standardize_type

struct DataFileHelper
    header::Int
    size::Int
    width::Int
    height::Int
    tot_len::Int
    offset::Int
    dtype::DataType
    stype::DataType
end

function DataFileHelper(kwargs...)
    arg_dict = Dict(kwargs)
    header = get(arg_dict, :header, 0)
    width = get(arg_dict, :width, 0)
    tot_len = get(arg_dict, :length, 0)
    @assert (width != 0 || tot_len != 0) "Specify at least either width or length"
    width != 0 && tot_len == 0 && (tot_len = width)
    height = get(arg_dict, :height, 0)
    offset = get(arg_dict, :offset, 0)
    dtype = standardize_type(get(arg_dict, :dtype, UInt16))
    stype = standardize_type(get(arg_dict, :stype, Float32))
    size = tot_len * height
    DataFileHelper(header, size, width, height, tot_len, offset, dtype, stype)
end

function load_data(
    io::IO;
    header::Int = 0,
    dtype::Union{String,DataType} = UInt16,
)
    if header != 0
        seek(io, header)
    end
    Mmap.mmap(io, Vector{standardize_type(dtype)})
end

function load_data(f::AbstractString; kwargs...)
    open(f) do s
        load_data(s; kwargs...)
    end
end

function sinog_extract(data::Vector{<:Real}, dhp::DataFileHelper)
    if dhp.size == 0
        eff_size = length(data)
    else
        eff_size = dhp.size
    end
    width = dhp.width
    height = dhp.height
    if height == 0
        height = eff_size ÷ dhp.tot_len
        eff_size = width * height
    end
    conv_data = convert(Vector{dhp.stype}, data[1:eff_size])
    reshape(conv_data, height, width)
end

sinog_extract(data; kwargs...) = sinog_extract(data, DataFileHelper(kwargs...))

function sinog_load(s::IO; kwargs...)
    dhp = DataFileHelper(kwargs...)
    raw_data = load_data(s, header = dhp.header, dtype = dhp.dtype)
    sinog_extract(raw_data, dhp)
end

function sinog_load(f::AbstractString; kwargs...)
    open(f) do s
        sinog_load(s; kwargs...)
    end
end

function load_yaml(s::IO)
    yaml_doc = YAML.load(s; dicttype = Dict{Symbol,Any})
    (; pairs(yaml_doc)...)
end

function load_yaml(f::AbstractString)
    open(f) do s
        load_yaml(s)
    end
end

function write_yaml(io::IO, obj)
    YAML.write(io, yaml_repr(obj))
end

function write_yaml(io::IO; kwargs...)
    write_yaml(io, Dict(imap(yaml_repr, kwargs)))
end

function write_yaml(f::String, obj)
    open(f, "w") do s
        write_yaml(s, obj)
    end
end

function write_yaml(f::String; kwargs...)
    open(f, "w") do s
        write_yaml(s; kwargs...)
    end
end

"""
    read_ct_image(io::IO[, T=Float32]; kwargs...) where {T<:Real}

Read CT image from the stream `io`.

# Keyword arguments

- `header`: read rows and columns from header in the file. [default: false]
- `rows`: the image rows.
- `cols`: the image columns.
- `header_type`: type for the two values of the header. [default: `Int64`]
- `row_major`: assume the image is stored in row major ordering.

See also: [`write_ct_image`](@ref)
"""
function read_ct_image(
    io::IO,
    ::Type{T} = Float32;
    header = false,
    rows = 0,
    cols = 0,
    header_type = Int64,
    row_major = false,
) where {T<:Real}
    if header
        rows = read(io, header_type)
        cols = read(io, header_type)
    elseif cols == 0
        cols = rows
    end
    if row_major
        data = Matrix{T}(undef, cols, rows)
        return read!(io, data)' |> collect
    end
    data = Matrix{T}(undef, rows, cols)
    read!(io, data)
end

"""
    read_ct_image(f::AbstractString[, T=Float32]; <keyword arguments>) where {T<:Real}

Read CT image from file `f`.
"""
function read_ct_image(
    f::AbstractString,
    ::Type{T} = Float32;
    kwargs...,
) where {T<:Real}
    open(f) do s
        read_ct_image(s, T; kwargs...)
    end
end

"""
    write_ct_image(io::IO, image::AbstractMatrix{T}; kwargs...) where {T<:Real}

Write CT image to the stream `io`.

# Keyword arguments

- `header`: if `true`, write header information with rows and columns. [default: false]
- `header_type`: type for the two values of the header. [default: `Int64`]
- `row_major`: write image in row major ordering.

See also: [`read_ct_image`](@ref)
"""
function write_ct_image(
    io::IO,
    image::AbstractMatrix{T};
    header = false,
    header_type = Int64,
    row_major = false,
) where {T<:Real}
    if header
        write(io, header_type(rows))
        write(io, header_type(cols))
    end
    if row_major
        write(io, image')
        return
    end
    write(io, image)
end

"""
    write_ct_image(f::AbstractString, image::AbstractMatrix{T}; <keyword arguments>) where {T<:Real}

Write CT image to file `f`.
"""
function write_ct_image(
    f::AbstractString,
    image::AbstractMatrix{T};
    kwargs...,
) where {T<:Real}
    open(f, "w") do s
        write_ct_image(s, image; kwargs...)
    end
end

"""
    load_image(io::IO[, T=Float32]; <keyword arguments>) where {T<:Real}

Load image from stream `io`.

# Arguments
- `io`: the stream where to load the image from.
- `header=false`: specify if the file contains the header.
- `nϕ=0`: number of projections (rows in the file for row major ordering).
- `nd=0`: number of detectors (columns in the file for row major ordering).
- `row_major=false`: tells if the data is stored in row major ordering or not.
- `header_type=Int64`: the type of the header data.

See Also: [`write_image`](@ref), [`read_ct_image`](@ref), [`write_ct_image`](@ref)
"""
function load_image(
    io::IO,
    ::Type{T} = Float32;
    header = false,
    rows = 0,
    cols = 0,
    row_major = false,
    header_type = Int64,
) where {T<:Real}
    read_ct_image(
        io,
        T;
        header = header,
        rows = rows,
        cols = cols,
        row_major = false,
        header_type = header_type,
    ) |> CTImage
end

"""
    load_image(f::AbstractString[, T=Float32]; <keyword arguments>) where {T<:Real}

Load image from file `f`.
"""
function load_image(
    f::AbstractString,
    ::Type{T} = Float32;
    kwargs...,
) where {T<:Real}
    open(f) do s
        load_image(s, T; kwargs...)
    end
end

"""
    write_image(io::IO, image::AbstractMatrix{T}; <keyword arguments>) where {T<:Real}

Write input image to stream `io`.

# Arguments
- `io`: the stream where the image should be written to.
- `image`: the image data matrix.
- `header`: specify if the image header should be written.
- `row_major`: specify if the image should be stored in row major order.
- `header_type=Int64`: the type of the header data.

See Also: [`read_image`](@ref), [`write_ct_image`](@ref), [`read_ct_image`](@ref).
"""
function write_image(
    io::IO,
    image::AbstractMatrix{T};
    header = false,
    row_major = false,
    header_type = Int64,
) where {T<:Real}
    write_ct_image(
        io,
        image;
        header = header,
        row_major = row_major,
        header_type = header_type,
    )
end

"""
    write_image(f::AbstractString, image::AbstractMatrix{T}; <keyword arguments>) where {T<:Real}

Write input image to file `f`.
"""
function write_image(
    f::AbstractString,
    image::AbstractMatrix{T};
    kwargs...,
) where {T<:Real}
    open(f, "w") do s
        write_image(s, image; kwargs...)
    end
end

"""
    load_sinogram(io::IO[, ::Type{T}=Float32]; <keyword arguments>) where {T<:Real}

Load sinogram from stream `io`.

This function is specifically designed to load a sinogram. Do not
expect consistent result with other type of images. Use more general
function instead, refer to the `See Also` section.

# Arguments
- `io`: the stream where to load the image from.
- `header=false`: specify if the file contains the header.
- `nϕ=0`: number of projections (rows in the file for row major ordering).
- `nd=0`: number of detectors (columns in the file for row major ordering).
- `row_major=false`: tells if the data is stored in row major ordering or not.
- `header_type=Int64`: the type of the header data.

See Also: [`write_sinogram`](@ref), [`read_ct_image`](@ref), [`write_ct_image`](@ref)
"""
function load_sinogram(
    io::IO,
    ::Type{T} = Float32;
    header = false,
    nϕ = 0,
    nphi = nothing,
    nd = 0,
    row_major = false,
    header_type = Int64,
) where {T<:Real}
    nϕ = maybe(nϕ, nphi)
    read_ct_image(
        io,
        T;
        header = header,
        rows = nd,
        cols = nϕ,
        row_major = false,
        header_type = header_type,
    ) |> CTSinogram
end

"""
    load_sinogram(f::AbstractString[, T=Float32]; <keyword arguments>) where {T<:Real}

Load sinogram from file `f`.
"""
function load_sinogram(
    f::AbstractString,
    ::Type{T} = Float32;
    kwargs...,
) where {T<:Real}
    open(f) do s
        load_sinogram(s, T; kwargs...)
    end
end

"""
    write_sinogram(io::IO, sinog::AbstractMatrix{T}; <keyword arguments>) where {T<:Real}

Write sinogram `sinog` to stream `io`.

This function is specifically designed to deal with a sinogram image.

# Arguments
- `io`: the stream where to load the image from.
- `header=false`: specify if the file contains the header.
- `row_major=false`: tells if the data is stored in row major ordering or not.
- `header_type=Int64`: the type of the header data.

See Also: [`load_sinogram`](@ref), [`read_ct_image`](@ref), [`write_ct_image`](@ref).
"""
function write_sinogram(
    io::IO,
    sinog::AbstractMatrix{T};
    header = false,
    row_major = false,
    header_type = Int64,
) where {T<:Real}
    write_ct_image(
        io,
        sinog;
        header = header,
        row_major = row_major,
        header_type = header_type,
    )
end

"""
    write_sinogram(f::AbstractString, sinog::AbstractMatrix{T}) where {T<:Real}

Write sinogram to file `f`.
"""
function write_sinogram(
    f::AbstractString,
    sinog::AbstractMatrix{T};
    kwargs...,
) where {T<:Real}
    open(f, "w") do s
        write_sinogram(s, sinog; kwargs...)
    end
end

"""
    load_tomogram(io::IO[, T=Float32]; <keyword arguments>) where {T<:Real}

Load reconstructed image from stream `io`.

This function is specifically designed to load a tomogram. Do not
expect consistent result with other type of images. Use more general
function instead, refer to the `See Also` section.

# Arguments
- `io`: the stream where to load the image from.
- `header=false`: specify if the file contains the header.
- `nϕ=0`: number of projections (rows in the file for row major ordering).
- `nd=0`: number of detectors (columns in the file for row major ordering).
- `row_major=false`: tells if the data is stored in row major ordering or not.
- `header_type=Int64`: the type of the header data.

See Also: [`write_tomogram`](@ref), [`read_ct_image`](@ref), [`write_ct_image`](@ref)
"""
function load_tomogram(
    io::IO,
    ::Type{T} = Float32;
    header = false,
    rows = 0,
    cols = 0,
    row_major = false,
    header_type = Int64,
) where {T<:Real}
    load_image(
        io,
        T;
        header = header,
        rows = rows,
        cols = cols,
        row_major = row_major,
        header_type = header_type,
    ) |> CTTomogram
end

"""
    load_tomogram(f::AbstractString[, T=Float32]; <keyword arguments>) where {T<:Real}

Load reconstructed image from file `f`.
"""
function load_tomogram(
    f::AbstractString,
    ::Type{T} = Float32;
    kwargs...,
) where {T<:Real}
    open(f) do s
        load_tomogram(s, T; kwargs...)
    end
end

"""
    write_tomogram(io::IO, tomog::AbstractMatrix{T}; <keyword arguments>) where {T<:Real})

Write reconstructed image to stream `io`.

# Arguments
- `io`: the stream where to load the image from.
- `header=false`: specify if the file contains the header.
- `row_major=false`: tells if the data is stored in row major ordering or not.
- `header_type=Int64`: the type of the header data.
"""
function write_tomogram(
    io::IO,
    tomog::AbstractMatrix{T};
    header = false,
    row_major = false,
    header_type = Int64,
) where {T<:Real}
    write_image(
        io,
        tomog;
        header = header,
        row_major = row_major,
        header_type = header_type,
    )
end

"""
    write_tomogram(f::AbstractString, tomog::AbstractMatrix{T}; <keyword arguments>) where {T<:Real}

Write reconstructed image to file `f`.
"""
function write_tomogram(
    f::AbstractString,
    tomog::AbstractMatrix{T};
    kwargs...,
) where {T<:Real}
    open(f, "w") do s
        write_tomogram(s, tomog; kwargs...)
    end
end

const _io_write_funcs = (:write_image, :write_sinogram, :write_tomogram)

for nm ∈ _io_write_funcs
    @eval begin
        $nm(io::Union{IO,AbstractString}; kwargs...) =
            x -> $nm(io, x; kwargs...)
        $nm(io::Union{IO,AbstractString}, img::AbstractCTImage; kwargs...) =
            img ↣ $nm(io; kwargs...)
    end
end

end  # module


text:
include("./superres.jl")

using Plots
pyplot();

let ntests=100,ℓ=2,R=3,clusterFun=Ζ,amplitudeFun = opposite1s
    let pertFun=𝔉ᵪ(ℓ)
        let ϵR = 10.0.^LinRange(-14,-3,150), ΩR = NrangeLog(1,3,150), ΔR = 10.0.^LinRange(-4,-2,150)
            let res1 = [runSingle_Acc(ℓ,R,clusterFun, pertFun,amplitudeFun, ϵR,ΩR,ΔR) for _ in 1:ntests]
                let res=mapreduce(transpose,vcat,res1)
                    let ΔX=reduce(vcat,res[:,2]), ΔC = reduce(vcat,res[:,3]),Succ=reduce(vcat,res[:,4])
                        p1=scatter(log10.(res[:,1]),log10.(ΔX./Succ),
                                    shape=[:circle :x :square],color=[:red :black :blue],
                                    markeralpha=[0.6 1.0 0.6],markerstrokewidth=[0.1 1 0.1],
                                    msc=[:white :black :white],ms=[8 5 5],
                                    xlabel="log(SRF)",ylabel="log(\$ \\cal K \$)",
                                    lab=cat(["MP \$ {\\cal K}_{x,$(i)} \$" for i=1:R]...,dims=2));
                        plot!([t->(2*ℓ-2)*t+1],w=1,l=:dash,lc=[:black],lab="\$ SRF^{2p-2} \$ ");
                        p2=scatter(log10.(res[:,1]),log10.(ΔC./Succ),
                                shape=[:circle :x :square],color=[:red :black :blue],
                                markeralpha=[0.6 1.0 0.6],markerstrokewidth=[0.1 1 0.1],
                                msc=[:white :black :white],ms=[8 5 5],
                                xlabel="log(SRF)",ylabel="log(\$ \\cal K \$)",
                                lab=cat(["MP \$ {\\cal K}_{a,$(i)} \$" for i=1:R]...,dims=2));
                        plot!([t->(2*ℓ-1)*t+1],w=1,l=:dash,lc=[:black],lab="\$ SRF^{2p-1} \$ ");
                        titleStr = "p=$(ℓ),d=$(R), \$ {\\bf S_2} \$"
                        p=plot(p1,p2,layout=(1,2),size=(800,300),title=titleStr);
                        savefig(p,"./figures/fig4.pdf")
                    end
                end
            end
        end
    end
end

max_stars_repo_name: batenkov-team/iai-code

max_issues_repo_name: batenkov-team/iai-code

max_forks_repo_name: batenkov-team/iai-code


text:
module NeuriViz

"""

        fdt_parser(path::String, dims::Array{Int64, 1})

`fdt` files are a type of binary file format that stores "float data".
This function parses them to a Julia readable format provided the path to a `fdt` file and the dimensions of the output array.

# Arguments

- `path::String`: path to an `fdt` file
- `dims::Array{Int64, 1}`: array representing [rows, columns] of the desired output

# Returns

- `Array{Float32}`: Julia readable array containing parsed `fdt` data

"""
function fdt_parser(path, dims)
    parsed_data = open(path) do io
        read!(io, Array{Float32}(undef, (dims[1], dims[2])))
    end

    return parsed_data

end

export fdt_parser

end

max_stars_repo_name: TheCedarPrince/NeuriViz

max_issues_repo_name: TheCedarPrince/NeuriViz

max_forks_repo_name: TheCedarPrince/NeuriViz


text:
# A form is something that ends up as geometry in the graphic.

abstract type FormPrimitive end

const empty_tag = Symbol("")

struct Form{P <: FormPrimitive} <: ComposeNode
    primitives::Vector{P}
    tag::Symbol

    Form{P}(prim, tag::Symbol=empty_tag) where P = new{P}(prim, tag)
end

Form(primitives::Vector{P}, tag::Symbol=empty_tag) where P <: FormPrimitive =
        Form{P}(primitives, tag)

isempty(f::Form) = isempty(f.primitives)

isscalar(f::Form) =
        length(f.primitives) == 1

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, form::Form) =
        Form([resolve(box, units, t, primitive) for primitive in form.primitives])

Base.similar(f::Form{T}) where T = Form{T}(T[])

form_string(::Form) = "FORM"  # fallback definition

# Polygon
# -------

struct SimplePolygonPrimitive{P <: Vec} <: FormPrimitive
    points::Vector{P}
end

const SimplePolygon{P<:SimplePolygonPrimitive} = Form{P}

const Polygon = SimplePolygon
const PolygonPrimitive = SimplePolygonPrimitive

polygon() = Form([PolygonPrimitive(Vec[])])

"""
    polygon(points)

Define a polygon. `points` is an array of `(x,y)` tuples
that specify the corners of the polygon.
"""
function polygon(points::AbstractArray{T}, tag=empty_tag) where T <: XYTupleOrVec
    XM, YM = narrow_polygon_point_types(Vector[points])
    if XM == Any
        XM = Measure
    end
    if YM == Any
        YM = Measure
    end
    VecType = Tuple{XM, YM}

    return Form([PolygonPrimitive(VecType[(x_measure(point[1]), y_measure(point[2]))
                    for point in points])], tag)
end

function polygon(point_arrays::AbstractArray, tag=empty_tag)
    XM, YM = narrow_polygon_point_types(point_arrays)
    VecType = XM == YM == Any ? Vec : Tuple{XM, YM}
    PrimType = XM == YM == Any ? PolygonPrimitive : PolygonPrimitive{VecType}

    polyprims = Array{PrimType}(length(point_arrays))
    for (i, point_array) in enumerate(point_arrays)
        polyprims[i] = PrimType(VecType[(x_measure(point[1]), y_measure(point[2]))
                                        for point in point_array])
    end
    return Form{PrimType}(polyprims, tag)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::PolygonPrimitive) =
        PolygonPrimitive{AbsoluteVec2}( AbsoluteVec2[
            resolve(box, units, t, point) for point in p.points])

function boundingbox(form::PolygonPrimitive, linewidth::Measure,
                     font::AbstractString, fontsize::Measure)
    x0 = minimum([p[1] for p in form.points])
    x1 = maximum([p[1] for p in form.points])
    y0 = minimum([p[2] for p in form.points])
    y1 = maximum([p[2] for p in form.points])
    return BoundingBox(x0 - linewidth,
                       y0 - linewidth,
                       x1 - x0 + linewidth,
                       y1 - y0 + linewidth)
end

form_string(::SimplePolygon) = "SP"

struct ComplexPolygonPrimitive{P <: Vec} <: FormPrimitive
    rings::Vector{Vector{P}}
end

const ComplexPolygon{P<:ComplexPolygonPrimitive} = Form{P}

complexpolygon() = ComplexPolygon([ComplexPolygonPrimitive(Vec[])])

function complexpolygon(rings::Vector{Vector}, tag=empty_tag)
    XM, YM = narrow_polygon_point_types(rings)
    if XM == Any
        XM = Length{:cx, Float64}
    end
    if YM == Any
        YM = Length{:cy, Float64}
    end
    VecType = Tuple{XM, YM}

    return ComplexPolygon([
        ComplexPolygonPrimitive([VecType[(x_measure(point[1]), y_measure(point[2]))
                                         for point in points]])], tag)
end

function complexpolygon(ring_arrays::Vector{Vector{Vector}}, tag=empty_tag)
    XM, YM = narrow_polygon_point_types(coords)
    VecType = XM == YM == Any ? Vec : Tuple{XM, YM}
    PrimType = XM == YM == Any ? ComplexPolygonPrimitive : ComplexPolygonPrimitive{VecType}

    ComplexPolygon([PrimType[[(x_measure(x), y_measure(y)) for (x, y) in ring]
                    for ring in ring_array] for ring_array in ring_arrays], tag)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::ComplexPolygonPrimitive) =
        ComplexPolygonPrimitive( [AbsoluteVec2[
            resolve(box, units, t, point, t) for point in ring] for ring in p.rings])

form_string(::ComplexPolygon) = "CP"

# Rectangle
# ---------

struct RectanglePrimitive{P <: Vec, M1 <: Measure, M2 <: Measure} <: FormPrimitive
    corner::P
    width::M1
    height::M2
end

const Rectangle{P<:RectanglePrimitive} = Form{P}

"""
    rectangle()

Define a rectangle that fills the current context completely.
"""
function rectangle()
    prim = RectanglePrimitive((0.0w, 0.0h), 1.0w, 1.0h)
    return Rectangle{typeof(prim)}([prim])
end

"""
    rectangle(x0, y0, width, height)

Define a rectangle of size `width`x`height` with its top left corner at the point (`x`, `y`).
"""
function rectangle(x0, y0, width, height, tag=empty_tag)
    corner = (x_measure(x0), y_measure(y0))
    width = x_measure(width)
    height = y_measure(height)
    prim = RectanglePrimitive(corner, width, height)
    return Rectangle{typeof(prim)}([prim], tag)
end

"""
    rectangle(x0s, y0s, widths, heights)

Arguments can be passed in arrays in order to perform multiple drawing operations at once.
"""
rectangle(x0s::AbstractArray, y0s::AbstractArray,
                   widths::AbstractArray, heights::AbstractArray, tag=empty_tag) =
        @makeform (x0 in x0s, y0 in y0s, width in widths, height in heights),
            RectanglePrimitive{Vec2, Measure, Measure}((x_measure(x0), y_measure(y0)),
                                                    x_measure(width), y_measure(height)) tag

function resolve(box::AbsoluteBox, units::UnitBox, t::Transform,
                 p::RectanglePrimitive)
    corner = resolve(box, units, t, p.corner)
    width = resolve(box, units, t, p.width)
    height = resolve(box, units, t, p.height)

    if isxflipped(units) && hasunits(Length{:cx}, p.corner[1])
        # if coordinates are flipped we end up with the corner on the other end
        # of the rectangle, which is fix here
        x = corner[1] - width
    else
        x = corner[1]
    end

    if isyflipped(units) && hasunits(Length{:cy}, p.corner[2])
        y = corner[2] - height
    else
        y = corner[2]
    end

    return RectanglePrimitive{AbsoluteVec2, AbsoluteLength, AbsoluteLength}(
        (x, y), width, height)
end

boundingbox(form::RectanglePrimitive, linewidth::Measure,
                     font::AbstractString, fontsize::Measure) =
        BoundingBox(form.corner.x - linewidth,
                    form.corner.y - linewidth,
                    form.width + 2*linewidth,
                    form.height + 2*linewidth)

form_string(::Rectangle) = "R"

# Circle
# ------

struct CirclePrimitive{P <: Vec, M <: Measure} <: FormPrimitive
    center::P
    radius::M
end

CirclePrimitive(center::P, radius::M) where {P, M} = CirclePrimitive{P, M}(center, radius)
CirclePrimitive(x, y, r) = CirclePrimitive((x_measure(x), y_measure(y)), x_measure(r))

const Circle{P<:CirclePrimitive} = Form{P}

"""
    circle()

Define a circle in the center of the current context with a diameter equal to the width of the context.
"""
function circle()
    prim = CirclePrimitive((0.5w, 0.5h), 0.5w)
    return Circle{typeof(prim)}([prim])
end

"""
    circle(x, y, r)

Define a circle with its center at (`x`,`y`) and a radius of `r`.
"""
function circle(x, y, r, tag=empty_tag)
    prim = CirclePrimitive(x, y, r)
    return Circle{typeof(prim)}([prim], tag)
end

"""
    circle(xs, ys, rs)

Arguments can be passed in arrays in order to perform multiple drawing operations.
"""
function circle(xs::AbstractArray, ys::AbstractArray, rs::AbstractArray, tag=empty_tag)
    if isempty(xs) || isempty(ys) || isempty(rs)
        prima = CirclePrimitive[]
        return Circle{eltype(prima)}(prima, tag)
    end

    return @makeform (x in xs, y in ys, r in rs),
        CirclePrimitive{Vec2, Measure}((x_measure(x), y_measure(y)), x_measure(r)) tag
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::CirclePrimitive) =
        CirclePrimitive{AbsoluteVec2, AbsoluteLength}(
            resolve(box, units, t, p.center),
            resolve(box, units, t, p.radius))

boundingbox(form::CirclePrimitive, linewidth::Measure, font::AbstractString, fontsize::Measure) =
        BoundingBox(form.center[1] - form.radius - linewidth,
            form.center[2] - form.radius - linewidth,
            2 * (form.radius + linewidth),
            2 * (form.radius + linewidth))

form_string(::Circle) = "C"

# Ellipse
# -------

struct EllipsePrimitive{P1<:Vec, P2<:Vec, P3<:Vec} <: FormPrimitive
    center::P1
    x_point::P2
    y_point::P3
end

const Ellipse{P<:EllipsePrimitive} = Form{P}

function ellipse()
    prim = EllipsePrimitive((0.5w, 0.5h),
                            (1.0w, 0.5h),
                            (0.5w, 1.0h))
    return Ellipse{typeof(prim)}([prim])
end

function ellipse(x, y, x_radius, y_radius, tag=empty_tag)
    xm = x_measure(x)
    ym = y_measure(y)
    prim = EllipsePrimitive((xm, ym),
                            (xm + x_measure(x_radius), ym),
                            (xm, ym + y_measure(y_radius)))
    return Ellipse{typeof(prim)}([prim], tag)
end

function ellipse(xs::AbstractArray, ys::AbstractArray,
                 x_radiuses::AbstractArray, y_radiuses::AbstractArray, tag=empty_tag)
        return @makeform (x in xs, y in ys, x_radius in x_radiuses, y_radius in y_radiuses),
            EllipsePrimitive((x_measure(x), y_measure(y)),
                     (x_measure(x) + x_measure(x_radius), y_measure(y)),
                     (x_measure(x), y_measure(y) + y_measure(y_radius))) tag
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::EllipsePrimitive) =
        EllipsePrimitive{AbsoluteVec2, AbsoluteVec2, AbsoluteVec2}(
            resolve(box, units, t, p.center),
            resolve(box, units, t, p.x_point),
            resolve(box, units, t, p.y_point))

function boundingbox(form::EllipsePrimitive, linewidth::Measure,
                     font::AbstractString, fontsize::Measure)
    x0 = min(form.x_point.x, form.y_point.x)
    x1 = max(form.x_point.x, form.y_point.x)
    y0 = min(form.x_point.y, form.y_point.y)
    y1 = max(form.x_point.y, form.y_point.y)
    xr = x1 - x0
    yr = y1 - y0
    return BoundingBox(x0 - linewidth - xr,
                       y0 - linewidth - yr,
                       2 * (xr + linewidth),
                       2 * (yr + linewidth))
end

form_string(::Ellipse) = "E"

# Text
# ----

abstract type HAlignment end
struct HLeft   <: HAlignment end
struct HCenter <: HAlignment end
struct HRight  <: HAlignment end

const hleft   = HLeft()
const hcenter = HCenter()
const hright  = HRight()

abstract type VAlignment end
struct VTop    <: VAlignment end
struct VCenter <: VAlignment end
struct VBottom <: VAlignment end

const vtop    = VTop()
const vcenter = VCenter()
const vbottom = VBottom()

struct TextPrimitive{P<:Vec, R<:Rotation, O<:Vec} <: FormPrimitive
    position::P
    value::AbstractString
    halign::HAlignment
    valign::VAlignment

    # Text forms need their own rotation field unfortunately, since there is no
    # way to give orientation with just a position point.
    rot::R

    offset::O
end

const Text{P<:TextPrimitive} = Form{P}

"""
    text(x, y, value [,halign::HAlignment [,valign::VAlignment [,rot::Rotation]]])

Draw the text `value` at the position (`x`,`y`) relative to the current context.

The default alignment of the text is `hleft` `vbottom`. The vertical and horizontal
alignment is specified by passing `hleft`, `hcenter` or `hright` and `vtop`,
`vcenter` or `vbottom` as values for `halgin` and `valgin` respectively.
"""
function text(x, y, value::AbstractString,
              halign::HAlignment=hleft, valign::VAlignment=vbottom,
              rot=Rotation(), offset::Vec2=(0mm,0mm);
              tag::Symbol=empty_tag)
    moffset = (x_measure(offset[1]), y_measure(offset[2]))
    prim = TextPrimitive((x_measure(x), y_measure(y)), value, halign, valign, rot, moffset)
    Text{typeof(prim)}([prim], tag)
end

function text(x, y, value,
              halign::HAlignment=hleft, valign::VAlignment=vbottom,
              rot=Rotation(), offset::Vec2=(0mm,0mm);
              tag::Symbol=empty_tag)
    moffset = (x_measure(offset[1]), y_measure(offset[2]))
    prim = TextPrimitive((x_measure(x), y_measure(y)), string(value), halign, valign, rot, moffset)
    Text{typeof(prim)}([prim], tag)
end

"""
    text(xs, ys, values [,haligns::HAlignment [,valigns::VAlignment [,rots::Rotation]]])

Arguments can be passed in arrays in order to perform multiple drawing operations at once.
"""
text(xs::AbstractArray, ys::AbstractArray, values::AbstractArray{AbstractString},
              haligns::AbstractArray=[hleft], valigns::AbstractArray=[vbottom],
              rots::AbstractArray=[Rotation()], offsets::AbstractArray=[(0mm,0mm)];
              tag::Symbol=empty_tag) =
    @makeform (x in xs, y in ys, value in values, halign in haligns, valign in valigns, rot in rots, offset in offsets),
        TextPrimitive((x_measure(x), y_measure(y)), value, halign, valign, rot, (x_measure(offset[1]), y_measure(offset[2]))) tag

text(xs::AbstractArray, ys::AbstractArray, values::AbstractArray,
              haligns::AbstractArray=[hleft], valigns::AbstractArray=[vbottom],
              rots::AbstractArray=[Rotation()], offsets::AbstractArray=[(0mm,0mm)];
              tag::Symbol=empty_tag) =
    @makeform (x in xs, y in ys, value in values, halign in haligns, valign in valigns, rot in rots, offset in offsets),
        TextPrimitive((x_measure(x), y_measure(y)), value, halign, valign, rot, (x_measure(offset[1]), y_measure(offset[2]))) tag

function resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::TextPrimitive{P,R,O}) where {P,R,O}
    rot = resolve(box, units, t, p.rot)
    return TextPrimitive{AbsoluteVec2, typeof(rot), O}(
                resolve(box, units, t, p.position),
                p.value, p.halign, p.valign, rot, p.offset)
end

function boundingbox(form::TextPrimitive, linewidth::Measure,
                     font::AbstractString, fontsize::Measure)

    width, height = text_extents(font, fontsize, form.value)[1]

    if form.halign == hleft
        x0 = form.position.x
    elseif form.halign == hcenter
        x0 = form.position.x - width/2
    elseif form.halign == hright
        x0 = form.position.x - width
    end

    if form.valign == vbottom
        y0 = form.position.y - height
    elseif form.valign == vcenter
        y0 = form.position.y - height/2
    elseif form.valign == vtop
        y0 = form.position.y
    end

    return BoundingBox(x0 - linewidth + form.offset[1],
                       y0 - linewidth + form.offset[2],
                       width + linewidth,
                       height + linewidth)
end

form_string(::Text) = "T"

# Line
# ----

struct LinePrimitive{P<:Vec} <: FormPrimitive
    points::Vector{P}
end

const Line{P<:LinePrimitive} = Form{P}

function line()
    prim = LinePrimitive(Vec[])
    return Line{typeof(prim)}([prim])
end

function line(points::AbstractArray{T}, tag=empty_tag) where T <: XYTupleOrVec
    XM, YM = narrow_polygon_point_types(Vector[points])
    VecType = XM == YM == Any ? Vec2 : Tuple{XM, YM}
    prim = LinePrimitive(VecType[(x_measure(point[1]), y_measure(point[2])) for point in points])
    return Line{typeof(prim)}([prim], tag)
end

function line(point_arrays::AbstractArray, tag=empty_tag)
    XM, YM = narrow_polygon_point_types(point_arrays)
    VecType = XM == YM == Any ? Vec2 : Tuple{XM, YM}
    PrimType = XM == YM == Any ? LinePrimitive : LinePrimitive{VecType}

    lineprims = Array{PrimType}(length(point_arrays))
    for (i, point_array) in enumerate(point_arrays)
        p = PrimType(VecType[(x_measure(point[1]), y_measure(point[2]))
                             for point in point_array])
        lineprims[i] = p
    end
    return Form{PrimType}(lineprims, tag)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::LinePrimitive) =
        LinePrimitive{AbsoluteVec2}(
            AbsoluteVec2[resolve(box, units, t, point) for point in p.points])

function boundingbox(form::LinePrimitive, linewidth::Measure,
                     font::AbstractString, fontsize::Measure)
    x0 = minimum([p.x for p in form.points])
    x1 = maximum([p.x for p in form.points])
    y0 = minimum([p.y for p in form.points])
    y1 = maximum([p.y for p in form.points])
    return BoundingBox(x0 - linewidth,
                       y0 - linewidth,
                       x1 - x0 + linewidth,
                       y1 - y0 + linewidth)
end

form_string(::Line) = "L"

# Curve
# -----

struct CurvePrimitive{P1<:Vec, P2<:Vec, P3<:Vec, P4<:Vec} <: FormPrimitive
    anchor0::P1
    ctrl0::P2
    ctrl1::P3
    anchor1::P4
end

const Curve{P<:CurvePrimitive} = Form{P}

function curve(anchor0::XYTupleOrVec, ctrl0::XYTupleOrVec,
               ctrl1::XYTupleOrVec, anchor1::XYTupleOrVec, tag=empty_tag)
    prim = CurvePrimitive((x_measure(anchor0[1]), y_measure(anchor0[2])),
                                 (x_measure(ctrl0[1]), y_measure(ctrl0[2])),
                                 (x_measure(ctrl1[1]), y_measure(ctrl1[2])),
                                 (x_measure(anchor1[1]), y_measure(anchor1[2])))
    return Curve{typeof(prim)}([prim], tag)
end

curve(anchor0s::AbstractArray, ctrl0s::AbstractArray,
               ctrl1s::AbstractArray, anchor1s::AbstractArray, tag=empty_tag) =
        @makeform (anchor0 in anchor0s, ctrl0 in ctrl0s, ctrl1 in ctrl1s, anchor1 in anchor1s),
            CurvePrimitive((x_measure(anchor0[1]), y_measure(anchor0[2])),
                (x_measure(ctrl0[1]), y_measure(ctrl0[2])),
                (x_measure(ctrl1[1]), y_measure(ctrl1[2])),
                (x_measure(anchor1[1]), y_measure(anchor1[2]))) tag

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::CurvePrimitive) =
        CurvePrimitive{AbsoluteVec2, AbsoluteVec2, AbsoluteVec2, AbsoluteVec2}(
            resolve(box, units, t, p.anchor0),
            resolve(box, units, t, p.ctrl0),
            resolve(box, units, t, p.ctrl1),
            resolve(box, units, t, p.anchor1))

form_string(::Curve) = "CV"

# Bitmap
# ------

struct BitmapPrimitive{P <: Vec, XM <: Measure, YM <: Measure} <: FormPrimitive
    mime::AbstractString
    data::Vector{UInt8}
    corner::P
    width::XM
    height::YM
end

const Bitmap{P<:BitmapPrimitive} = Form{P}

function bitmap(mime::AbstractString, data::Vector{UInt8}, x0, y0, width, height, tag=empty_tag)
    corner = (x_measure(x0), y_measure(y0))
    width = x_measure(width)
    height = y_measure(height)
    prim = BitmapPrimitive(mime, data, corner, width, height)
    return Bitmap{typeof(prim)}([prim], tag)
end

bitmap(mimes::AbstractArray, datas::AbstractArray,
                x0s::AbstractArray, y0s::AbstractArray,
                widths::AbstractArray, heights::AbstractArray, tag=empty_tag) =
        @makeform (mime in mimes, data in datas, x0 in x0s, y0 in y0s, width in widths, height in heigths),
            BitmapPrimitive(mime, data, x0, y0, x_measure(width), y_measure(height)) tag

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::BitmapPrimitive) =
        BitmapPrimitive{AbsoluteVec2, AbsoluteLength, AbsoluteLength}(
            p.mime, p.data,
            resolve(box, units, t, p.corner),
            resolve(box, units, t, p.width),
            resolve(box, units, t, p.height))

boundingbox(form::BitmapPrimitive, linewidth::Measure, font::AbstractString, fontsize::Measure) =
        BoundingBox(form.corner.x, form.corner.y, form.width, form.height)

form_string(::Bitmap) = "B"

# Path
# ----

# An implementation of the SVG path mini-language.

abstract type PathOp end

struct MoveAbsPathOp <: PathOp
    to::Vec
end

function assert_pathop_tokens_len(op_type, tokens, i, needed)
    provided = length(tokens) - i + 1
    provided < needed &&
            error("In path $(op_type) requires $(needed) argumens but only $(provided) provided.")
end

function parsepathop(::Type{MoveAbsPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(MoveAbsPathOp, tokens, i, 2)
    op = MoveAbsPathOp((x_measure(tokens[i]), y_measure(tokens[i + 1])))
    return (op, i + 2)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::MoveAbsPathOp) =
        MoveAbsPathOp(resolve(box, units, t, p.to))

struct MoveRelPathOp <: PathOp
    to::Vec
end

function parsepathop(::Type{MoveRelPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(MoveRelPathOp, tokens, i, 2)
    op = MoveRelPathOp((tokens[i], tokens[i + 1]))
    return (op, i + 2)
end

function resolve_offset(box::AbsoluteBox, units::UnitBox, t::Transform, p::Vec)
    absp = resolve(box, units, t, p)
    zer0 = resolve(box, units, t, (0w, 0h))
    return (absp[1] - zer0[1], absp[2] - zer0[2])
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::MoveRelPathOp) =
        MoveRelPathOp(resolve_offset(box, units, t, p.to))

struct ClosePathOp <: PathOp
end

parsepathop(::Type{ClosePathOp}, tokens::AbstractArray, i) = (ClosePathOp(), i)

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::ClosePathOp) = p

struct LineAbsPathOp <: PathOp
    to::Vec
end

function parsepathop(::Type{LineAbsPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(LineAbsPathOp, tokens, i, 2)
    op = LineAbsPathOp((x_measure(tokens[i]), y_measure(tokens[i + 1])))
    return (op, i + 2)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::LineAbsPathOp) =
        LineAbsPathOp(resolve(box, units, t, p.to))

struct LineRelPathOp <: PathOp
    to::Vec
end

function parsepathop(::Type{LineRelPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(LineRelPathOp, tokens, i, 2)
    op = LineRelPathOp((x_measure(tokens[i]), y_measure(tokens[i + 1])))
    return (op, i + 2)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::LineRelPathOp) =
        LineRelPathOp(resolve(box, units, t, p.to))

struct HorLineAbsPathOp <: PathOp
    x::Measure
end

function parsepathop(::Type{HorLineAbsPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(HorLineAbsPathOp, tokens, i, 1)
    op = HorLineAbsPathOp(x_measure(tokens[i]))
    return (op, i + 1)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::HorLineAbsPathOp) =
        HorLineAbsPathOp(resolve(box, units, t, (p.x, 0mm))[1])

struct HorLineRelPathOp <: PathOp
    Δx::Measure
end

function parsepathop(::Type{HorLineRelPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(HorLineRelPathOp, tokens, i, 1)
    op = HorLineRelPathOp(x_measure(tokens[i]))
    return (op, i + 1)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::HorLineRelPathOp) =
        HorLineRelPathOp(resolve(box, units, t, p.Δx))

struct VertLineAbsPathOp <: PathOp
    y::Measure
end

function parsepathop(::Type{VertLineAbsPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(VertLineAbsPathOp, tokens, i, 1)
    op = VertLineAbsPathOp(y_measure(tokens[i]))
    return (op, i + 1)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::VertLineAbsPathOp) =
        VertLineAbsPathOp(resolve(box, units, t, (0mm, p.y))[2])

struct VertLineRelPathOp <: PathOp
    Δy::Measure
end

function parsepathop(::Type{VertLineRelPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(VertLineRelPathOp, tokens, i, 1)
    op = VertLineRelPathOp(y_measure(tokens[i]))
    return (op, i + 1)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::VertLineRelPathOp) =
        VertLineAbsPathOp(resolve(box, units, t, (0mmm, p.Δy))[2])

struct CubicCurveAbsPathOp <: PathOp
    ctrl1::Vec
    ctrl2::Vec
    to::Vec
end

function parsepathop(::Type{CubicCurveAbsPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(CubicCurveAbsPathOp, tokens, i, 6)
    op = CubicCurveAbsPathOp((tokens[i],     tokens[i + 1]),
                             (tokens[i + 2], tokens[i + 3]),
                             (tokens[i + 4], tokens[i + 5]))
    return (op, i + 6)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::CubicCurveAbsPathOp) =
        CubicCurveAbsPathOp(
            resolve(box, units, t, p.ctrl1),
            resolve(box, units, t, p.ctrl2),
            resolve(box, units, t, p.to))

struct CubicCurveRelPathOp <: PathOp
    ctrl1::Vec
    ctrl2::Vec
    to::Vec
end

function parsepathop(::Type{CubicCurveRelPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(CubicCurveRelPathOp, tokens, i, 6)
    op = CubicCurveRelPathOp((tokens[i],     tokens[i + 1]),
                             (tokens[i + 2], tokens[i + 3]),
                             (tokens[i + 4], tokens[i + 5]))
    return (op, i + 6)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::CubicCurveRelPathOp) =
        CubicCurveRelPathOp(
            resolve(box, units, t, p.ctrl1),
            resolve(box, units, t, p.ctrl2),
            resolve(box, units, t, p.to))

struct CubicCurveShortAbsPathOp <: PathOp
    ctrl2::Vec
    to::Vec
end

function parsepathop(::Type{CubicCurveShortAbsPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(CubicCurveShortAbsPathOp, tokens, i, 4)
    op = CubicCurveShortAbsPathOp((x_measure(tokens[i]),     y_measure(tokens[i + 1])),
                                  (x_measure(tokens[i + 2]), y_measure(tokens[i + 3])))
    return (op, i + 4)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::CubicCurveShortAbsPathOp) =
        CubicCurveShortAbsPathOp(
            resolve_offset(box, units, t, p.ctrl2),
            resolve_offset(box, units, t, p.to))

struct CubicCurveShortRelPathOp <: PathOp
    ctrl2::Vec
    to::Vec
end

function parsepathop(::Type{CubicCurveShortRelPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(CubicCurveShortRelPathOp, tokens, i, 4)
    op = CubicCurveShortRelPathOp((x_measure(tokens[i]),     y_measure(tokens[i + 1])),
                                  (x_measure(tokens[i + 2]), y_measure(tokens[i + 3])))
    return (op, i + 4)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::CubicCurveShortRelPathOp) =
        CubicCurveShortRelPathOp(
            resolve(box, units, t, p.ctrl2),
            resolve(box, units, t, p.to))

struct QuadCurveAbsPathOp <: PathOp
    ctrl1::Vec
    to::Vec
end

function parsepathop(::Type{QuadCurveAbsPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(QuadCurveAbsPathOp, tokens, i, 4)
    op = QuadCurveAbsPathOp((tokens[i],     tokens[i + 1]),
                            (tokens[i + 2], tokens[i + 3]))
    return (op, i + 4)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::QuadCurveAbsPathOp) =
        QuadCurveAbsPathOp(
            resolve(box, units, t, p.ctrl1),
            resolve(box, units, t, p.to))

struct QuadCurveRelPathOp <: PathOp
    ctrl1::Vec
    to::Vec
end

function parsepathop(::Type{QuadCurveRelPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(QuadCurveRelPathOp, tokens, i, 4)
    op = QuadCurveRelPathOp((tokens[i],     tokens[i + 1]),
                            (tokens[i + 2], tokens[i + 3]))
    return (op, i + 4)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::QuadCurveRelPathOp) =
        QuadCurveRelPathOp(
            (resolve(box, units, t, p.ctrl1[1]),
             resolve(box, units, t, p.ctrl1[2])),
            (resolve(box, units, t, p.to[1]),
             resolve(box, units, t, p.to[2])))

struct QuadCurveShortAbsPathOp <: PathOp
    to::Vec
end

function parsepathop(::Type{QuadCurveShortAbsPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(QuadCurveShortAbsPathOp, tokens, i, 2)
    op = QuadCurveShortAbsPathOp((tokens[i], tokens[i + 1]))
    return (op, i + 2)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::QuadCurveShortAbsPathOp) =
        QuadCurveShortAbsPathOp(resolve(box, units, t, p.to))

struct QuadCurveShortRelPathOp <: PathOp
    to::Vec
end

function parsepathop(::Type{QuadCurveShortRelPathOp}, tokens::AbstractArray, i)
    assert_pathop_tokens_len(QuadCurveShortRelPathOp, tokens, i, 2)
    op = QuadCurveShortRelPathOp((tokens[i], tokens[i + 1]))
    return (op, i + 2)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::QuadCurveShortRelPathOp) =
        QuadCurveShortRelPathOp(
            (resolve(box, units, t, p.to[1]),
             resolve(box, units, t, p.to[2])))

struct ArcAbsPathOp <: PathOp
    rx::Measure
    ry::Measure
    rotation::Float64
    largearc::Bool
    sweep::Bool
    to::Vec
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::ArcAbsPathOp) =
        ArcAbsPathOp(
            resolve(box, units, t, p.rx),
            resolve(box, units, t, p.ry),
            p.rotation,
            p.largearc,
            p.sweep,
            resolve(box, units, t, p.to))

struct ArcRelPathOp <: PathOp
    rx::Measure
    ry::Measure
    rotation::Float64
    largearc::Bool
    sweep::Bool
    to::Vec
end

function parsepathop(::Type{T}, tokens::AbstractArray, i) where T <: Union{ArcAbsPathOp, ArcRelPathOp}
    assert_pathop_tokens_len(T, tokens, i, 7)

    if isa(tokens[i + 3], Bool)
        largearc = tokens[i + 3]
    elseif tokens[i + 3] == 0
        largearc = false
    elseif tokens[i + 3] == 1
        largearc = true
    else
        error("largearc argument to the arc path operation must be boolean")
    end

    if isa(tokens[i + 4], Bool)
        sweep = tokens[i + 4]
    elseif tokens[i + 4] == 0
        sweep = false
    elseif tokens[i + 4] == 1
        sweep = true
    else
        error("sweep argument to the arc path operation must be boolean")
    end

    isa(tokens[i + 2], Number) || error("path arc operation requires a numerical rotation")

    op = T(x_measure(tokens[i]),
           y_measure(tokens[i + 1]),
           convert(Float64, tokens[i + 2]),
           largearc, sweep,
           (x_measure(tokens[i + 5]), y_measure(tokens[i + 6])))

    return (op, i + 7)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::ArcRelPathOp) =
        ArcRelPathOp(
            resolve(box, units, t, p.rx),
            resolve(box, units, t, p.ry),
            p.rotation,
            p.largearc,
            p.sweep,
            (resolve(box, units, t, p.to[1]),
             resolve(box, units, t, p.to[2])))

const path_ops = Dict(
     :M => MoveAbsPathOp,
     :m => MoveRelPathOp,
     :Z => ClosePathOp,
     :z => ClosePathOp,
     :L => LineAbsPathOp,
     :l => LineRelPathOp,
     :H => HorLineAbsPathOp,
     :h => HorLineRelPathOp,
     :V => VertLineAbsPathOp,
     :v => VertLineRelPathOp,
     :C => CubicCurveAbsPathOp,
     :c => CubicCurveRelPathOp,
     :S => CubicCurveShortAbsPathOp,
     :s => CubicCurveShortRelPathOp,
     :Q => QuadCurveAbsPathOp,
     :q => QuadCurveRelPathOp,
     :T => QuadCurveShortAbsPathOp,
     :t => QuadCurveShortRelPathOp,
     :A => ArcAbsPathOp,
     :a => ArcRelPathOp
)

# A path is an array of symbols, numbers, and measures following SVGs path
# mini-language.
function parsepath(tokens::AbstractArray)
    ops = PathOp[]
    last_op_type = nothing
    i = 1
    while i <= length(tokens)
        tok = tokens[i]
        strt = i
        if isa(tok, Symbol)
            if !haskey(path_ops, tok)
                error("$(tok) is not a valid path operation")
            else
                op_type = path_ops[tok]
                i += 1
                op, i = parsepathop(op_type, tokens, i)
                push!(ops, op)
                last_op_type = op_type
            end
        else
            op, i = parsepathop(last_op_type, tokens, i)
            push!(ops, op)
        end
    end

    return ops
end

struct PathPrimitive <: FormPrimitive
    ops::Vector{PathOp}
end

const Path = Form{PathPrimitive}

path(tokens::AbstractArray, tag=empty_tag) = Path([PathPrimitive(parsepath(tokens))], tag)

path(tokens::AbstractArray{T}, tag=empty_tag) where T <: AbstractArray =
        Path([PathPrimitive(parsepath(ts)) for ts in tokens], tag)

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::PathPrimitive) =
        PathPrimitive([resolve(box, units, t, op) for op in p.ops])

# TODO: boundingbox


text:
  module pkgs

  using Reexport

  import Pkg
Pkg.status()

import Pkg
Pkg.add("Luxor")
Pkg.add("Clipper")
Pkg.add("ProgressMeter")
Pkg.add("JSON")
Pkg.add("FFTW")
Pkg.add("HTTP")
Pkg.add("ForwardDiff")

Pkg.add("CUDA")

import CUDA
# This will download CUDA
@info CUDA.version()
@info CUDA.has_cuda()

Pkg.add("Flux")
Pkg.add("Zygote")

  end

max_stars_repo_name: bhdl/bhdl

max_issues_repo_name: bhdl/bhdl

max_forks_repo_name: bhdl/bhdl


text:
include("pkg.jl")
include("resolve.jl")

max_stars_repo_name: nreigl/julia

max_issues_repo_name: BenChung/Julia-BTI

max_forks_repo_name: BenChung/Julia-BTI