text:
module MeasurementsWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const AbstractMeasurement = @jimport org.orekit.estimation.measurements.AbstractMeasurement
const AngularAzEl = @jimport org.orekit.estimation.measurements.AngularAzEl
const AngularRaDec = @jimport org.orekit.estimation.measurements.AngularRaDec
const Class = @jimport java.lang.Class
const ComparableMeasurement = @jimport org.orekit.estimation.measurements.ComparableMeasurement
const EOPHistory = @jimport org.orekit.frames.EOPHistory
const EstimatedEarthFrameProvider = @jimport org.orekit.estimation.measurements.EstimatedEarthFrameProvider
const EstimatedMeasurement = @jimport org.orekit.estimation.measurements.EstimatedMeasurement
const EstimatedMeasurement_Status = @jimport org.orekit.estimation.measurements.EstimatedMeasurement$Status
const EstimationModifier = @jimport org.orekit.estimation.measurements.EstimationModifier
const EstimationsProvider = @jimport org.orekit.estimation.measurements.EstimationsProvider
const FieldAbsoluteDate = @jimport org.orekit.time.FieldAbsoluteDate
const FieldTransform = @jimport org.orekit.frames.FieldTransform
const FieldVector3D = @jimport org.hipparchus.geometry.euclidean.threed.FieldVector3D
const Frame = @jimport org.orekit.frames.Frame
const GeodeticPoint = @jimport org.orekit.bodies.GeodeticPoint
const GroundStation = @jimport org.orekit.estimation.measurements.GroundStation
const InterSatellitesRange = @jimport org.orekit.estimation.measurements.InterSatellitesRange
const JString = @jimport java.lang.String
const List = @jimport java.util.List
const Map = @jimport java.util.Map
const MultiplexedMeasurement = @jimport org.orekit.estimation.measurements.MultiplexedMeasurement
const Object = @jimport java.lang.Object
const ObservableSatellite = @jimport org.orekit.estimation.measurements.ObservableSatellite
const ObservedMeasurement = @jimport org.orekit.estimation.measurements.ObservedMeasurement
const PV = @jimport org.orekit.estimation.measurements.PV
const ParameterDriver = @jimport org.orekit.utils.ParameterDriver
const Position = @jimport org.orekit.estimation.measurements.Position
const Range = @jimport org.orekit.estimation.measurements.Range
const RangeRate = @jimport org.orekit.estimation.measurements.RangeRate
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const SpacecraftState = @jimport org.orekit.propagation.SpacecraftState
const StationDisplacement = @jimport org.orekit.models.earth.displacement.StationDisplacement
const Stream = @jimport java.util.stream.Stream
const TimeStamped = @jimport org.orekit.time.TimeStamped
const TimeStampedFieldPVCoordinates = @jimport org.orekit.utils.TimeStampedFieldPVCoordinates
const TimeStampedPVCoordinates = @jimport org.orekit.utils.TimeStampedPVCoordinates
const TopocentricFrame = @jimport org.orekit.frames.TopocentricFrame
const Transform = @jimport org.orekit.frames.Transform
const TurnAroundRange = @jimport org.orekit.estimation.measurements.TurnAroundRange
const UT1Scale = @jimport org.orekit.time.UT1Scale
const Vector3D = @jimport org.hipparchus.geometry.euclidean.threed.Vector3D

export AbstractMeasurement
export AngularAzEl
export AngularRaDec
export ComparableMeasurement
export EstimatedEarthFrameProvider
export EstimatedMeasurement
export EstimationModifier
export EstimationsProvider
export GroundStation
export InterSatellitesRange
export MultiplexedMeasurement
export ObservableSatellite
export ObservedMeasurement
export PV
export Position
export Range
export RangeRate
export TurnAroundRange
export add_modifier
export compare_to
export equals
export estimate
export get_base_frame
export get_base_weight
export get_class
export get_clock_offset_driver
export get_coordinates
export get_correlation_coefficients_matrix
export get_count
export get_covariance_matrix
export get_date
export get_derivatives_drivers
export get_dimension
export get_displacements
export get_east_offset_driver
export get_estimated_earth_frame
export get_estimated_measurement
export get_estimated_measurements
export get_estimated_ut1
export get_estimated_value
export get_iteration
export get_master_station
export get_measurements
export get_modifiers
export get_north_offset_driver
export get_number
export get_observed_measurement
export get_observed_value
export get_offset_geodetic_point
export get_offset_to_inertial
export get_parameter_derivatives
export get_parameters_drivers
export get_participants
export get_polar_drift_x_driver
export get_polar_drift_y_driver
export get_polar_offset_x_driver
export get_polar_offset_y_driver
export get_position
export get_prime_meridian_drift_driver
export get_prime_meridian_offset_driver
export get_propagator_index
export get_reference_frame
export get_satellites
export get_slave_station
export get_state_derivatives
export get_state_size
export get_states
export get_station
export get_status
export get_theoretical_standard_deviation
export get_time_offset
export get_transform
export get_velocity
export get_zenith_offset_driver
export hash_code
export is_enabled
export is_two_way
export modify
export notify
export notify_all
export set_enabled
export set_estimated_value
export set_parameter_derivatives
export set_state_derivatives
export set_status
export signal_time_of_flight
export to_string
export wait

include("FilteringWrapper/FilteringWrapper.jl")
include("GenerationWrapper/GenerationWrapper.jl")
include("GnssWrapper/GnssWrapper.jl")
include("ModifiersWrapper/ModifiersWrapper.jl")
include("abstract_measurement.jl")
include("angular_az_el.jl")
include("angular_ra_dec.jl")
include("comparable_measurement.jl")
include("estimated_earth_frame_provider.jl")
include("estimated_measurement.jl")
include("estimation_modifier.jl")
include("estimations_provider.jl")
include("ground_station.jl")
include("inter_satellites_range.jl")
include("multiplexed_measurement.jl")
include("observable_satellite.jl")
include("observed_measurement.jl")
include("position.jl")
include("pv")
include("range.jl")
include("range_rate.jl")
include("turn_around_range.jl")

end


text:
using JSON
using DataFrames
using Gadfly

include("utils.jl")

mfile = open("motifs.json")
hashes = JSON.parse(mfile)
close(mfile)

function getIWDB(url::ASCIIString)
  matrix = download(url)
  am = open(matrix)
  A = int(readdlm(am))
  close(am)
  return A
end

function cleanmatrix(A::Array{Int64, 2})
  cumulative_degree = vec(sum(A, 1)) .+ sum(A, 2)
  have_degree = filter((i) -> cumulative_degree[i] > 0, 1:length(cumulative_degree))
  A = A[have_degree, have_degree]
end

function countmotifs(A::Array{Int64, 2}, msize) # TODO motif dict
  nodes_id = combinations(1:size(A, 1), msize)
  motifs = ASCIIString[]
  for nbunch in nodes_id
    sub_A = A[collect(nbunch),collect(nbunch)]
    if nonempty(sub_A)
      V = size(sub_A)[1]
      E = sum(sub_A)
      if(E >= (V-1))
         m_sub_A = mhash(sub_A)
         m_id = string(V)*"_"*string(E)*"_"*string(hashes[string(V)][string(E)][m_sub_A])
         push!(motifs, m_id)
      end
    end
  end
  mcounts = {x => sum(motifs .== x) for x in unique(motifs)}
  return mcounts
end

function buildIWDBname(name::ASCIIString)
  return "https://www.nceas.ucsb.edu/interactionweb/data/predator_prey/text_matrices/" * name * "txt.txt"
end

function wrapIWDB(x::ASCIIString)
  ok_mat = x |> buildIWDBname |> getIWDB |> cleanmatrix
  return countmotifs(ok_mat, 3)
end

data_names = [
   "AkatoreA",
   "AkatoreB",
   "Catlins",
   "Berwick",
   "Narrowdale",
   "Blackrock",
   "Broad",
   "DempstersSu",
   "DempstersAu",
   "DempstersSp",
   "German",
   "Healy",
   "LilKyeburn",
   "Stony"
]

IWDBmotifs = {x => wrapIWDB(x) for x in data_names}

nrows = sum([length(v) for (k,v) in IWDBmotifs])
iwdb = DataFrame([Symbol, Symbol, Symbol, Int64], [:fw, :motif, :edges, :count], nrows)

cr = 1
for (fw, mc) in IWDBmotifs
  for (mn, cn) in mc
    iwdb[:fw][cr] = symbol(fw)
    iwdb[:motif][cr] = symbol(mn)
    iwdb[:edges][cr] = symbol(split(mn,"_")[2])
    iwdb[:count][cr] = cn
    cr += 1
  end
end

# Functions to deal with proportions
function prop(x::DataArray{Symbol, 1})
  return x
end

function prop(x::DataArray{Int64, 1})
  return x ./ sum(x)
end

# Agregate the dataframe
a_iwdb = aggregate(iwdb, :fw, prop)

draw(
  PDF("motif_iwdb.pdf", 30cm, 10cm),
  plot(a_iwdb, x="fw", xgroup="edges_prop", color="motif_prop", y="count_prop", Geom.subplot_grid(Geom.bar))
)

ext: jl


text:
function event_occurred(obj::FieldODEEventHandler, arg0::FieldODEStateAndDerivative, arg1::jboolean)
    return jcall(obj, "eventOccurred", Action, (FieldODEStateAndDerivative, jboolean), arg0, arg1)
end

function g(obj::FieldODEEventHandler, arg0::FieldODEStateAndDerivative)
    return jcall(obj, "g", RealFieldElement, (FieldODEStateAndDerivative,), arg0)
end

function init(obj::FieldODEEventHandler, arg0::FieldODEStateAndDerivative, arg1::RealFieldElement)
    return jcall(obj, "init", void, (FieldODEStateAndDerivative, RealFieldElement), arg0, arg1)
end

function reset_state(obj::FieldODEEventHandler, arg0::FieldODEStateAndDerivative)
    return jcall(obj, "resetState", FieldODEState, (FieldODEStateAndDerivative,), arg0)
end


text:
using CRF

include("util.jl")
include("features.jl")

const output_succ = colorize(:green, ".")
const output_fail = colorize(:red, "!")

# These parameters were estimated using sequences 1:5 of the weather data
const weights = [
    -0.28026, -0.54288,  0.04901, -0.25603,  1.09337,  0.49870,
     0.13389, -0.23739, -0.24375, -1.08821, -0.06286,  0.89308,
    -0.37372,  0.38654,  0.38167, -3.53123,  1.76168,  3.07147,
    -1.50942,  0.51421, -4.12856, -0.64580,  4.74698, -0.31943,
     2.64340, -0.58929, -2.62894,  0.10829,  0.18280, -0.54266,
    -0.47681,  2.13713, -1.66658, -2.48922,  1.82307, -0.30672,
     1.82972,  0.71780, -2.07108, -0.51137,  0.60169, -0.13208,
    -1.85513,  0.43814, -1.48759,  0.66490,  2.30342,  0.42272,
     1.62713,  0.70396, -1.13422,  3.62780,  4.14362,  1.62987,
    -0.85902, -3.39346, -3.32714, -2.12237, -1.70612, -0.31701,
     0.28682, -0.53939, -3.27933, -4.02713, -2.32070,  1.53830,
     2.50143,  1.38448,  1.65289, -0.06920,  0.63677,  0.05326,
     1.14765, -0.65981, -0.98722, -0.32634, -0.84214,  0.73448,
     1.06528,  0.95661,  0.84304, -0.59334, -0.36241, -0.66310,
     1.45682,  6.45608,  1.91501,  4.23405,  3.87607, -0.46321,
    -1.01568, -2.34954, -1.60510, -0.80973, -2.09774, -0.23512,
    -0.96984,  0.29281, -1.69349, -1.77548, -0.18955, -1.13785,
     0.13768, -2.34054,  0.25428, -1.02766, -0.30154, -0.35736,
    -0.48482,  0.05892, -0.01486, -0.04475,  1.87549, -4.50467,
    -1.51220, -0.85258, -1.91109,  3.69303,  1.07203,  4.04240,
     1.01121,  2.16412, -0.29784,  0.23954, -0.66587,  0.89715,
    -0.31191, -0.08533, -0.15711, -1.60186, -0.02909, -0.69339,
    -0.40224, -0.50869,  0.00055, -0.99869,  0.05103, -0.60559,
     0.44725,  0.13268, -2.81193, -1.61552, -1.24948, -2.67640,
    -1.28247, -2.80358, -0.34114, -0.74070,  1.51590, -1.77853,
     1.70655, -0.54594,  1.30247, -0.46603,  2.28510,  1.73859,
    -0.34424,  2.82764, -0.23264,  3.28881, -0.02248,  2.36526,
     0.10379,  1.30489,  0.37751,  0.86419,  0.42223,  0.25484,
]

# Load weather data
X, Y = load("weather.csv")

# Remove sequences we used for parameter estimation
splice!(X, 1:5)
splice!(Y, 1:5)

# Create sequences
crfs = Sequence[ Sequence(x, features, Θ=weights, labels=labels) for x in X ]

# Calling label on the array of sequences does the trick
Y_true = Y
Y_pred = label(crfs)

# Let's compare the predicted labels with the hidden labels
error = 0
total = 0
for (y_true, y_pred) in zip(Y_true, Y_pred)
    for (lt, lp) in zip(y_true, y_pred)
        print((lt == lp) ? output_succ : output_fail)
        total += 1
        error += (lt != lp)
    end
    print(" ")
end
print("\n")

@printf "Total: %6d\n" total
@printf "Error: %6d (%.2f %%)\n" error ((error / total) * 100.0)


text:
using Test, DeviationsLH
# Requires deviation_test_setup

dlh = DeviationsLH;

function dev_vector_test()
    @testset "DevVector" begin
        d = DevVector()
        @test length(d) == 0
        dev1 = make_deviation(1);
        dlh.append!(d, dev1);
        @test length(d) == 1

        # Scalar deviation
        dev2 = make_scalar_deviation(2);
        dlh.append!(d, dev2);
        @test length(d) == 2
        show_deviations(d; sorted = false);
        show_deviations(d; sorted = true);
        show(d);

        # Retrieve
        dev22 = dlh.retrieve(d, :d2);
        @test !isempty(dev22)
        @test dev22.dataV ≈ dev2.dataV
        @test dev_exists(d, :d2)
        @test !dev_exists(d, :notThere)

        @test isempty(retrieve(d, :notThere))
        dev = retrieve(d, :d1);
        @test dev.name == :d1

        # Set model values
        modelV = dev22.dataV .+ 1.3;
        set_model_values(d, :d2, modelV);
        dev22 = dlh.retrieve(d, :d2);
        @test dev22.modelV ≈ modelV
        @test isapprox(get_model_values(d, :d2), modelV)
        @test isapprox(get_data_values(d, :d2), get_data_values(dev22))

        # Set weights
        dev3 = make_deviation(3);
        dlh.append!(d, dev3);
        wtV = dev3.dataV .+ 2.3;
        dlh.set_weights!(d, :d3, wtV);
        dev3 = dlh.retrieve(d, :d3);
        @test dev3.wtV ≈ wtV
    end
end

function iterate_test()
    @testset "Iteration" begin
        d = make_deviation_vector(5);
        j = 0;
        for dev in d
            @test isa(dev, AbstractDeviation)
            j += 1;
        end
        @test j == length(d)
    end
end

function scalar_dev_devvector_test()
    @testset "Scalar deviation" begin
        n = 5;
        d = make_deviation_vector(n);
        devV = scalar_devs(d);
        @test length(devV) == n;
        # Not sorted, so test is weak
        dev1 = d[1];
        scalarDev1, _ = scalar_dev(dev1);
        @test any(x -> isapprox(x, scalarDev1), devV)

        scalarDev = scalar_deviation(d);
        @test isa(scalarDev, Float64)
        @test scalarDev >= 0.0

        @test isapprox(sum(devV), scalarDev, rtol = 1e-4)
    end
end

function devvec_dict_test()
    @testset "Deviation Dict" begin
        n = 5;
        d = make_deviation_vector(5);
        sds = scalar_dev_dict(d);
        @test isa(sds, Dict{Symbol, dlh.DevType})
        @test length(sds) == 5

        dev1 = d[1];
        scalarDev1, _ = scalar_dev(dev1);
        @test sds[dlh.name(dev1)] ≈ scalarDev1

        devV = scalar_devs(d);
        @test sort(devV) ≈ sort(collect(values(sds)))

        scalarDev = scalar_deviation(d);
        @test isapprox(sum(devV), scalarDev, rtol = 1e-4)
    end
end

@testset "DevVector" begin
    dev_vector_test()
    iterate_test()
    scalar_dev_devvector_test()
    devvec_dict_test()
end

# ----------------


text:
# --- PUBLIC METHODS ----------------------------------------------------------- #
function get_reactions_for_ec_number(ec_number::String)::Union{Array{Reaction,1}, Nothing}

    # Check: does this contain a ec: prefix? Yes, chop it off -
    local_ec_number = ec_number
    if (occursin("ec:",ec_number) == true)
        local_ec_number = ec_number[4:end]
    end

    # initalize -
    record_array = Reaction[]

    # formulate url -
    linkage_url_string = "http://rest.kegg.jp/link/rn/$(local_ec_number)"

    # call to get linkage array -
    ec_rn_linkage_array = http_get_call_with_url(linkage_url_string)

    # if the server returned something, it should be in the form ec# rn#
    # we need the rn# to pull down the reaction string -
    if (length(ec_rn_linkage_array)>1)

        # ec_rn_linkage can have more than one association -
        # split around the \n, and then process each item
        P = reverse(split(chomp(ec_rn_linkage_array),'\n'))
        while (!isempty(P))

            # pop -
            local_record = pop!(P)

            # The second fragment is the rn# -
            # split around the \t
            rn_number = string(split(local_record,'\t')[2])

            # ok, run a second query to pull the reaction string down -
            reaction_string = chomp(http_get_call_with_url("http://rest.kegg.jp/list/$(rn_number)"))

            # if we have something in reaction string, create a string record -
            if (length(reaction_string)!=0)

                # create a new Reaction wrapper -
                reaction_wrapper = Reaction()

                # populate w/easy stuff -
                reaction_wrapper.ec_number = ec_number
                reaction_wrapper.kegg_reaction_number = rn_number

                # parse the body string -
                reaction_wrapper.kegg_enzyme_name = string(split(reaction_string,";")[1])
                reaction_wrapper.kegg_reaction_markup = chomp(string(split(reaction_string,";")[2]))

                # cache -
                push!(record_array,reaction_wrapper)
            end
        end
    end

    # check - zero entries?
    if (length(record_array) == 0)
        return nothing
    end

    # return list of records -
    return record_array
end

# get_ec_number_for_gene
function get_ec_number_for_gene(gene_location::String)::Union{Array{String,1},Nothing}

    # TODO: check gene location string -

    # initialize -
    ec_number_array = String[]

    # formulate the url -
    url_string = "http://rest.kegg.jp/link/ec/$(gene_location)"

    # get the sequence -
    ecdata = http_get_call_with_url(url_string)

    # check: do we have only the new line? Yes, then return nothing
    if (ecdata == "\n")
        return nothing
    end

    # ok, parse the response which is of the form: {gene_location}\t{ec_number}\n*
    # we can have more than 1 record -
    record_array = split(ecdata,"\n")

    # process the records -
    for record in record_array

        # split along the \t -
        if (record!="")
            ec_number = string(split(record,"\t")[2])
            push!(ec_number_array,ec_number)
        end
    end

    # return this record -
    return ec_number_array
end
# ------------------------------------------------------------------------------ #

function write_ec_numbers_to_disk(path_to_input_gene_file::String, path_to_output_file::String, kegg_organism_code::String)

    # initalize -
    ec_number_buffer = String[]

    # get -
    gene_list = readdlm(path_to_input_gene_file,',')

    # how many genes do we have?
    (number_of_genes,number_of_cols) = size(gene_list)

    # create a progress meter -
    p = Progress(number_of_genes,color=:yellow)

    # start -
    @info "Downloading gene = ecnumber pairs from KEGG ..."

    # iterate through the list of tags -
    for gene_index = 1:number_of_genes

        # get the gene tag -
        gene_tag = String(gene_list[gene_index,2])

        # get the buffer -
        line_buffer = download_ec_number_from_kegg(gene_tag,kegg_organism_code)

        # check - does this have multiple records?
        P = split(line_buffer,'\n')
        number_of_fragments = length(P)
        for fragment_index = 1:number_of_fragments

            fragment_string = string(P[fragment_index])
            if (length(fragment_string)>0)

                # add line to overall buffer -
                push!(ec_number_buffer,fragment_string)

                # message -
                msg = "$(gene_tag) => $(fragment_string) ($(gene_index) of $(number_of_genes))"

                # update the progress bar -
                ProgressMeter.next!(p; showvalues = [(:status,msg)])
            end
        end
    end

    # dump to disk -
    open("$(path_to_output_file)/ec_numbers.dat", "w") do f

        for line_item in ec_number_buffer
            write(f,"$(line_item)\n")
        end
    end

    @info "Completed ...\r"

end

function write_protein_sequences_to_disk(path_to_input_gene_file::String, path_to_output_file::String, kegg_organism_code::String)

    # get -
    gene_list = readdlm(path_to_input_gene_file,',')

    # how many genes do we have?
    (number_of_genes, number_of_cols) = size(gene_list)

    # start -
    @info "Downloading protein sequences from KEGG ..."

    # create a progress meter -
    p = Progress(number_of_genes,color=:yellow)

    # iterate through the list of tags -
    for gene_index = 1:number_of_genes

        # get the gene tag -
        gene_tag = String(gene_list[gene_index,2])

        # get the buffer -
        buffer = download_protein_sequence_from_kegg(gene_tag, kegg_organism_code)

        # dump to disk -
        file_name = "$(path_to_output_file)/P_$(gene_tag).seq"
        writedlm(file_name,buffer)

        # update the progress bar -
        ProgressMeter.next!(p; showvalues = [(:gene_index,("$(gene_index) of $(number_of_genes)")), (:gene,gene_tag)])
    end

    @info "Completed ...\r"
end

function write_gene_sequences_to_disk(path_to_input_gene_file::String,path_to_output_file::String,kegg_organism_code::String)

    # get -
    gene_list = readdlm(path_to_input_gene_file,',')

    # how many genes do we have?
    (number_of_genes, number_of_cols) = size(gene_list)

    # start -
    @info "Downloading gene sequences from KEGG ..."

    # create a progress meter -
    p = Progress(number_of_genes,color=:yellow)

    # iterate through the list of tags -
    for gene_index in 1:number_of_genes

        # get the gene tag -
        gene_tag = String(gene_list[gene_index,2])

        # get the buffer -
        buffer = download_gene_sequence_from_kegg(gene_tag,kegg_organism_code)

        # dump to disk -
        file_name = "$(path_to_output_file)/$(gene_tag).seq"
        writedlm(file_name,buffer)

        # update the progress bar -
        ProgressMeter.next!(p; showvalues = [(:gene_index,("$(gene_index) of $(number_of_genes)")), (:gene,gene_tag)])

        # message -
        # msg = "$(gene_tag) => downloaded ($(gene_index) of $(number_of_genes))"
        # println(msg)
    end

    @info "Completed ...\r"
end

function write_gene_symbols_to_disk(cobra_dictionary, path_to_output_gene_file::String, organism_id::String)

    # get list of genes from the cobra_dictionary -
    gene_list = cobra_dictionary["genes"]

    # initialize buffer -
    gene_index = 1
    gene_name_buffer = String[]
    for gene_name in gene_list
        line = "$(gene_index),$(organism_id):$(gene_name)"
        push!(gene_name_buffer,line)
        gene_index = gene_index + 1
    end

    # dump to disk -
    open("$(path_to_output_gene_file)", "w") do f

        for line_item in gene_name_buffer
            write(f,"$(line_item)\n")
        end
    end

    # start -
    @info "Wrote gene symbols to $(path_to_output_gene_file)"
end


text:
"""
    make_bed(dir, lmap, sid, nf, target)
---
# Description
This function will merge the genotypes files of GSGT ver-2 format in `dir` into one
`vcf` file according linkage map `lmap`.  The sample ID in `dir` will be replaced
with animal ID according to file `sid`.  The animal ID is in the 2nd column of `sid`.
The sample ID is in `nf`th field, and not a `NA`.

# dir
In the future, put all the files from a same platform in the same folder.  This
folder should contains the genotype files of GSGT ver-2 format.  Nothing else.

# lmap
This file contains 3 columns:
1. SNP name
2. Chromsome number.  Here they 1-26 autosomes only.
3. Base pair position.

The SNP are ordered on their chromosome numbers and bp positions.

# sid
ID info file is partial flexible.  But the first six columns are fixed, namely:
1. Herdbook_number
2. AnimalID
3. BirthYear
4. BreedGroup
5. Breed
6. Gender

In the future, when new platform arrives, just append another column, indicating
which ID are genotyped on this platform.  The current platforms are:
7. SampleID_LD
8. SampleID_HD
9. SampleID_17Kbeta
10. SampleID_17K
11. SampleID_17Kgamma

# Example
    NSG.make_bed("dat/raw/a17k", "dat/maps/a17k.map", "dat/idinfo/id.info", 10, "dat/plink/a17k")

Above will merge all the files in `dat/raw/a17k` into `dat/plink/a17k.{bed,bim,fam}`.
"""
function make_bed(dir, lmap, sid, nf, target)
    title("Merge files in $dir into $target.bed")
    item("Check parameters")
    if !isdir(dir)
        warning("$dir doesn't exist")
        return
    end

    dir = abspath(dir)
    files = readdir(dir, join=true)
    if length(files) == 0
        warning("No file in $dir")
    end

    if !isfile(lmap)
        warning("Map $lmap doesn't exist")
        return
    end
    lmap = abspath(lmap)

    if !isfile(sid)
        warning("ID description $sid doesn't exist")
        return
    end

    target =  begin
        pt, target = splitdir(target)
        if length(pt) == 0
            warning("You are saving the bed files in the current dir")
            pt = pwd()
        else
            pt = abspath(pt)
            isdir(pt) || mkpath(pt)
        end
        joinpath(pt, target)
    end

    done("OK")

    item("Create an ID reference")
    tmp = mktempdir(".")
    open(joinpath(tmp, "id.txt"), "w") do io
        buffer = readlines(sid)[2:end] # read in all ID descriptions and skip header
        for line in buffer
            id, nm, yr, brd = split(line)[[2, nf, 3, 5]]
            yr = parse(Int, yr)
            if nm ≠ "NA" && brd=="10" && yr > 1999
                println(io, "$nm $id")
            end
        end
    end
    done()

    item("Create beagle files")
    cd(tmp)
    run(`$bin_dir/mrg2bgl id.txt $lmap $files`)
    done()

    item("Create VCF files")
    for chr in 1:26
        run(pipeline(`java -jar $bin_dir/beagle2vcf.jar $chr $chr.mrk $chr.bgl -`,
                     "$chr.vcf"))
    end

    mv("1.vcf", "ori.vcf")
    open("ori.vcf", "a") do io
        for chr in 2:26
            for line in eachline("$chr.vcf")
                line[1] ≠ '#' && println(io, line)
            end
        end
    end
    done()

    item("Make $target.bed")
    _ = read(run(`$plink --sheep --vcf ori.vcf --const-fid --out $target`), String)
    cd("..")
    rm(tmp, recursive=true, force=true)
    done()
end

"""
    id_desc(raw, tgt)
---
This funciton is to simplify ID descriptions in `raw` to `tgt`.

# Descriptions

## Originally, the genotyped ID were described as below:
  1. Herdbook_number
  2. AnimalID
  3. SampleID_LD
  4. SampleID_HD
  5. SampleID_17Kbeta
  6. SampleID_17K
  7. BirthYear
  8. BreedGroup
  9. Breed
  10. Gender
  11. SampleID_17Kgamma

## New descriptions will be:

  1. `HBN` => Herdbook_number
  2. `ID` => AnimalID
  3. `SID` => SampleID
  4. `PF` => Platform, e.g, LD, a17k, b17k, c17k and so on
  5. `YoB` => BirthYear
  6. `Group` => BreedGroup
  7. `Breed` => Breed
  8. `Sex` => Gender

In the future, descriptions of new genotyped ID will just simply append to the previous one.
If an ID were genotyped on more than one platforms, then the PF field can be "LDHD".

**OBS!: dropped, as an ID may have more than one `SID`, as well as more than one `PF`.**
"""
function id_desc(raw, tgt)
    # dic = Dict(3=>"LD", 4=>"HD", 5=>"b17k", 6=>"a17k", 11=>"c17k")
    # open(tgt, "w") do oo
    #     println(oo, "HBN ID SID PF YoB Group Breed Sex")
    #     open(raw, "r") do ii
    #         _ = readline(raw)   # skip the header
    #         for line in eachline(ii)
    #             rcd = split(line)
    #             pf = ""
    #             for i in [3, 4, 5, 6, 11]
    #                 rcd[i] ≠ "NA" && (pf*=dic[i])
    #             end
    #             println(oo, join([rcd[1, 2],
end


text:
struct Node
    id::Union{Integer, Nothing}
    alias::Union{String, Nothing}
    label::String
    properties::Dict
end

Node(alias::String, label::String, properties::Dict) = Node(nothing, alias, label, properties)
Node(label::String, properties::Dict) = Node(nothing, nothing, label, properties)
Node(alias::String, label::String) = Node(nothing, alias, label, Dict())
Node(label::String) = Node(nothing, nothing, label, Dict())
Node(id::Integer) = Node(id, nothing, "", Dict())

# TODO: should be moved to some common module
function prop_value_to_string(prop_value::Any) return "$prop_value" end
function prop_value_to_string(prop_value::String) return "\"$prop_value\"" end

function string(n::Node)
    alias = ""
    label = n.label
    props_str = ""

    if n.alias != nothing
        alias = n.alias
    end

    if length(n.properties) != 0
        props = ["$prop_name: " * prop_value_to_string(prop_value) for (prop_name, prop_value) in pairs(n.properties)]
        props_str = "{" * join(props, ",") * "}"
    end
    "($alias:$label $props_str)"
end

function isequal(x::Node, y::Node)
    if x.id != nothing && y.id != nothing && x.id == y.id
        return true
    else
        if x.label == y.label && x.properties == y.properties
            return true
        end
    end
    return false
end

Base.:(==)(x::Node, y::Node) = isequal(x, y)


text:
# Floating point type for simulations
T = Float64

# Geometrical setup
setup = PlateSetup{T}(;
    name = "plates/someplates",
    width = 50.0,
    depth = 50.0,
    heights = fill(5.0, 5),
    bend = 0.0,
    twist = 0.0,
    # refinement = 5.0,
)

ncell = length(setup.heights)

# Get compartimentalized coefficient vectors
coeffs = coefficients(
    setup;
    D = [0.002 * I(3) for _ = 1:ncell],
    T₂ = fill(Inf, ncell),
    ρ = fill(1.0, ncell),
    κ = (; interfaces = fill(1e-4, ncell - 1), boundaries = fill(0.0, ncell)),
    γ = 2.67513e-4,
)


text:
module PointCloud
using jInv.Mesh
using jInv.Utils
using jInv.InverseSolve
using ShapeReconstructionPaLS.Utils
using ShapeReconstructionPaLS.ParamLevelSet
using MAT
using SparseArrays
using Distributed
import jInv.ForwardShare.getData
import jInv.ForwardShare.getSensTMatVec
import jInv.ForwardShare.getSensMatVec
import jInv.ForwardShare.ForwardProbType

## The rotation planes are:
# first rotation (theta) is for planes x1 and x2
# second rotation (phi) is for planes x1 and x3
# Sampling is done by summing planes of x1 and x2 for all values of x3.

export PointCloudParam, getPointCloudParam
mutable struct PointCloudParam <: ForwardProbType
    Mesh      				:: RegularMesh
	P						:: Array{Array{Int64,1}}
	Normals					:: Array{Float64,2}
	npcAll					:: Int64
	margin					:: Float64
	workerSubIdxs
    theta_phi_rad			:: Array{Float64,2}
	b						:: Array{Float64,2}
	samplingBinning			:: Int64
	method					:: String
	Jacobian				:: SparseMatrixCSC{Float32,Int32}
	#S						:: SparseMatrixCSC{Float32,Int32}
end

function getPointCloudParamInternal(Mesh::RegularMesh, P:: Array{Array{Int64,1}}, Normals::Array{Float64,2},npcAll::Int64,margin::Float64,workerSubIdxs,theta_phi_rad::Array{Float64,2} ,b::Array{Float64,2},samplingBinning::Int64,method)
	return PointCloudParam(Mesh, P:: Array{Array{Int64,1}}, Normals::Array{Float64,2}, npcAll, margin, workerSubIdxs, theta_phi_rad::Array{Float64,2} ,b::Array{Float64,2}, samplingBinning, method,spzeros(Float32,Int32,0,0));
end

function getPointCloudParam(Mesh::RegularMesh,P::Array{Array{Int64,1}}, Normals::Array{Float64,2}, margin::Float64,theta_phi_rad::Array{Float64,2} ,b::Array{Float64,2},samplingBinning::Int64,numWorkers::Int64,method = MATBased)
	## This function does use the parallel mechanism of jInv (i.e., returns a RemoteChannel), even if numWorkers=1.
	if numWorkers > nworkers()
		numWorkers = nworkers();
	end
	#SourcesSubInd = Array{Array{Int64,1}}(numWorkers);
	SourcesSubInd = Array{Array{Int64}}(undef,numWorkers);
	ActualWorkers = workers();
	if numWorkers < nworkers()
		ActualWorkers = ActualWorkers[1:numWorkers];
	end
	#pFor   = Array{RemoteChannel}(numWorkers)
	pFor   = Array{RemoteChannel}(undef,numWorkers);
	i = 1; nextidx() = (idx=i; i+=1; idx)
	idx=i;
	ndips  = length(P);
	# send out jobs
	@sync begin
		for w = ActualWorkers
			@async begin
				while true
					idx = nextidx();
					if idx > numWorkers
						break
					end
					I_p = getIndicesOfKthWorker(numWorkers,idx,ndips);
					tmp = initRemoteChannel(getPointCloudParamInternal,w, Mesh,P:: Array{Array{Int64,1}}, Normals::Array{Float64,2},ndips, margin ,I_p,theta_phi_rad::Array{Float64,2} ,b::Array{Float64,2},samplingBinning,method);
					pFor[1]  = initRemoteChannel(getPointCloudParamInternal,w, Mesh,P:: Array{Array{Int64,1}}, Normals::Array{Float64,2},ndips, margin,I_p,theta_phi_rad::Array{Float64,2} ,b::Array{Float64,2},samplingBinning,method);
					wait(pFor[1]);

				end
			end
		end
	end
	return pFor # Array of Remote Refs
end

import jInv.Utils.clear!
function clear!(pFor::PointCloudParam)
	pFor.SampleMat = spzeros(0);
	clear!(pFor.Mesh);
	return pFor;
end

#include("generateDipMatrix.jl");
include("prepareSyntheticPointCloudData.jl");
include("getData.jl");
include("sensitivityFuncs.jl");
end

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
export hmm_train

import ..HMMModel

using mlpack._Internal.io

import mlpack_jll
const hmm_trainLibrary = mlpack_jll.libmlpack_julia_hmm_train

# Call the C binding of the mlpack hmm_train binding.
function hmm_train_mlpackMain()
  success = ccall((:hmm_train, hmm_trainLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module hmm_train_internal
  import ..hmm_trainLibrary

import ...HMMModel

# Get the value of a model pointer parameter of type HMMModel.
function IOGetParamHMMModel(paramName::String)::HMMModel
  HMMModel(ccall((:IO_GetParamHMMModelPtr, hmm_trainLibrary), Ptr{Nothing}, (Cstring,), paramName))
end

# Set the value of a model pointer parameter of type HMMModel.
function IOSetParamHMMModel(paramName::String, model::HMMModel)
  ccall((:IO_SetParamHMMModelPtr, hmm_trainLibrary), Nothing, (Cstring, Ptr{Nothing}), paramName, model.ptr)
end

# Serialize a model to the given stream.
function serializeHMMModel(stream::IO, model::HMMModel)
  buf_len = UInt[0]
  buf_ptr = ccall((:SerializeHMMModelPtr, hmm_trainLibrary), Ptr{UInt8}, (Ptr{Nothing}, Ptr{UInt}), model.ptr, Base.pointer(buf_len))
  buf = Base.unsafe_wrap(Vector{UInt8}, buf_ptr, buf_len[1]; own=true)
  write(stream, buf)
end
# Deserialize a model from the given stream.
function deserializeHMMModel(stream::IO)::HMMModel
  buffer = read(stream)
  HMMModel(ccall((:DeserializeHMMModelPtr, hmm_trainLibrary), Ptr{Nothing}, (Ptr{UInt8}, UInt), Base.pointer(buffer), length(buffer)))
end
end # module

"""
    hmm_train(input_file; [batch, gaussians, input_model, labels_file, seed, states, tolerance, type, verbose])

This program allows a Hidden Markov Model to be trained on labeled or unlabeled
data.  It supports four types of HMMs: Discrete HMMs, Gaussian HMMs, GMM HMMs,
or Diagonal GMM HMMs

Either one input sequence can be specified (with `input_file`), or, a file
containing files in which input sequences can be found (when
`input_file`and`batch` are used together).  In addition, labels can be provided
in the file specified by `labels_file`, and if `batch` is used, the file given
to `labels_file` should contain a list of files of labels corresponding to the
sequences in the file given to `input_file`.

The HMM is trained with the Baum-Welch algorithm if no labels are provided.  The
tolerance of the Baum-Welch algorithm can be set with the `tolerance`option.  By
default, the transition matrix is randomly initialized and the emission
distributions are initialized to fit the extent of the data.

Optionally, a pre-created HMM model can be used as a guess for the transition
matrix and emission probabilities; this is specifiable with `output_model`.

# Arguments

 - `input_file::String`: File containing input observations.
 - `batch::Bool`: If true, input_file (and if passed, labels_file) are
      expected to contain a list of files to use as input observation sequences
      (and label sequences).  Default value `false`.

 - `gaussians::Int`: Number of gaussians in each GMM (necessary when type
      is 'gmm').  Default value `0`.

 - `input_model::HMMModel`: Pre-existing HMM model to initialize training
      with.
 - `labels_file::String`: Optional file of hidden states, used for labeled
      training.  Default value ``.

 - `seed::Int`: Random seed.  If 0, 'std::time(NULL)' is used.  Default
      value `0`.

 - `states::Int`: Number of hidden states in HMM (necessary, unless
      model_file is specified).  Default value `0`.

 - `tolerance::Float64`: Tolerance of the Baum-Welch algorithm.  Default
      value `1e-05`.

 - `type_::String`: Type of HMM: discrete | gaussian | diag_gmm | gmm.
      Default value `gaussian`.

 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `output_model::HMMModel`: Output for trained HMM.

"""
function hmm_train(input_file::String;
                   batch::Union{Bool, Missing} = missing,
                   gaussians::Union{Int, Missing} = missing,
                   input_model::Union{HMMModel, Missing} = missing,
                   labels_file::Union{String, Missing} = missing,
                   seed::Union{Int, Missing} = missing,
                   states::Union{Int, Missing} = missing,
                   tolerance::Union{Float64, Missing} = missing,
                   type_::Union{String, Missing} = missing,
                   verbose::Union{Bool, Missing} = missing,
                   points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, hmm_trainLibrary), Nothing, ());

  IORestoreSettings("Hidden Markov Model (HMM) Training")

  # Process each input argument before calling mlpackMain().
  IOSetParam("input_file", input_file)
  if !ismissing(batch)
    IOSetParam("batch", convert(Bool, batch))
  end
  if !ismissing(gaussians)
    IOSetParam("gaussians", convert(Int, gaussians))
  end
  if !ismissing(input_model)
    hmm_train_internal.IOSetParamHMMModel("input_model", convert(HMMModel, input_model))
  end
  if !ismissing(labels_file)
    IOSetParam("labels_file", convert(String, labels_file))
  end
  if !ismissing(seed)
    IOSetParam("seed", convert(Int, seed))
  end
  if !ismissing(states)
    IOSetParam("states", convert(Int, states))
  end
  if !ismissing(tolerance)
    IOSetParam("tolerance", convert(Float64, tolerance))
  end
  if !ismissing(type_)
    IOSetParam("type", convert(String, type_))
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("output_model")
  # Call the program.
  hmm_train_mlpackMain()

  return hmm_train_internal.IOGetParamHMMModel("output_model")
end


text:
#=
Author: Daniel Carrera (dcarrera@gmail.com)

Compute the core radius using the table from Zeng et al (2016).
http://adsabs.harvard.edu/abs/2016ApJ...819..127Z

In general I'll use a bi-linear interpolation on the table from
Zeng et al. But if the mass is out of range, I'll use R ~ M^(1/3.7)
(also from Zeng et al.) I don't expect that to happen very often.
=#
using DataFrames
using Missings
using CSV

#####################
#
# COMPUTE PLANET RADII
#
#####################
#
# I specify the column names by hand because CSV.read keeps adding spaces
# to the column names.
#
mr_names = Symbol[]
append!(mr_names, [:Mearth])
append!(mr_names, map(p -> Symbol("$(p)%fe"), 100:-5:5))
append!(mr_names, [:rocky])
append!(mr_names, map(p -> Symbol("$(p)%h2o"), 5:5:100))

mr_table = "$(@__DIR__)/../data/Zeng_2016_MR_table.tsv"
mr_table = CSV.read(mr_table, datarow=2, header=mr_names, delim='\t')

"""
Compute the core radius for either a rock-iron or water-rock core.

Examples:

	#
	# Radius of a 3.0 M_earth core with 10% water + 90% rock:
	#
	radius = core_radius(3.0, h2o=0.1) # In Earth radii.

	#
	# Radius of a 3.0 M_earth core with 10% iron + 90% rock:
	#
	radius = core_radius(3.0, fe=0.1) # In Earth radii.

Citation:

	Zeng et al (2016)
	http://adsabs.harvard.edu/abs/2016ApJ...819..127Z
"""
function core_radius(mass; h2o=0.0, fe=0.0)::Float64
	global mr_table
	@assert(mass > 0)
	@assert(h2o >= 0)
	@assert(h2o <= 1)
	@assert(fe >= 0)
	@assert(fe <= 1)

	if (h2o > 0) && (fe > 0)
		@error("h2o and fe fractions cannot both be positive.")
	end

	#
	# Column names.
	#
	local p1, p2, c1, c2
	if (h2o == 0) && (fe == 0)
		p1 = p2 = 0
		c1 = c2 = :rocky
	else
		p   = (h2o > 0 ?  h2o  :  fe ) * 100
		mat = (h2o > 0 ? "h2o" : "fe")

		p1 = floor(Int, p / 5) * 5
		p2 = ceil( Int, p / 5) * 5

		c1 = p1 == 0 ? :rocky : Symbol("$(p1)%$(mat)")
		c2 = p2 == 0 ? :rocky : Symbol("$(p2)%$(mat)")
	end
	#=
	NOTE:	I considered collapsing the `if` below and just making sure
			that p1 != p2. But when you add the extra lines to deal with
			the :rocky column and fe vs h2o columns, I'm not sure that
			the code would look any better.
	=#

	#
	# (Single column) ?
	#
	if c1 == c2
		col = c1
		#
		# (Single colum) AND (Mass out of bounds) ?
		# 		=> Extrapolate.
		#
		if mass <= mr_table[1,:Mearth]
			return mr_table[1,col] * (mass/mr_table[1,:Mearth])^(1/3.7)
		end
		if mass >= mr_table[end,:Mearth]
			return mr_table[end,col] * (mass/mr_table[end,:Mearth])^(1/3.7)
		end
		#
		# (Single colum) AND NOT (Mass out of bounds)
		# 		=> Interpolate.
		#
		row1 = sum(mr_table[:Mearth] .< mass)
		row2 = row1 + 1

		mass1 = mr_table[row1, :Mearth]
		mass2 = mr_table[row2, :Mearth]

		R1 = mr_table[row1, col]
		R2 = mr_table[row2, col]

		return R1 * (mass - mass2)/(mass1 - mass2) + R2 * (mass - mass1)/(mass2 - mass1)
	else
		#
		# NOT (Single column)
		#
		col1 = c1
		col2 = c2
		#
		# NOT (Single column) AND (Mass out of bounds) ?
		# 		=> Interpolate + Extrapolate.
		#
		if mass <= mr_table[1,:Mearth]
			R1 = mr_table[1,col1]
			R2 = mr_table[1,col2]

			R = R1 * (p - p2)/(p1 - p2) + R2 * (p - p1)/(p2 - p1)
			return R * (mass/mr_table[1,:Mearth])^(1/3.7)
		end
		if mass >= mr_table[end,:Mearth]

			R1 = mr_table[end,col1]
			R2 = mr_table[end,col2]

			R = R1 * (p - p2)/(p1 - p2) + R2 * (p - p1)/(p2 - p1)
			return R * (mass/mr_table[end,:Mearth])^(1/3.7)
		end
		#
		# NOT (Single column) AND NOT (Mass out of bounds)
		# 		=> Bi-linear interpolation.
		#
		row1 = sum(mr_table[:Mearth] .< mass)
		row2 = row1 + 1

		mass1 = mr_table[row1, :Mearth]
		mass2 = mr_table[row2, :Mearth]

		R11 = mr_table[row1, col1] # (row,col) -- (mass, h2o or fe)
		R12 = mr_table[row1, col2] # (row,col) -- (mass, h2o or fe)
		R21 = mr_table[row2, col1] # (row,col) -- (mass, h2o or fe)
		R22 = mr_table[row2, col2] # (row,col) -- (mass, h2o or fe)

		R = 0.0
		R += R11 * (mass - mass2)/(mass1 - mass2) * (p - p2)/(p1 - p2)
		R += R12 * (mass - mass2)/(mass1 - mass2) * (p - p1)/(p2 - p1)
		R += R21 * (mass - mass1)/(mass2 - mass1) * (p - p2)/(p1 - p2)
		R += R22 * (mass - mass1)/(mass2 - mass1) * (p - p1)/(p2 - p1)
		return R
	end
end


text:
function jamsd_set_modeltype(m::JAMSDMathProgModel, idx)
    discrete = any((m.vartypes .== :Int) + (m.vartypes .== :Bin) .> 0)
    if discrete
        if m.model_type == qcp
            m.model_type = miqcp
        elseif m.model_type == nlp
             m.model_type = minlp
        end
    end
    res = ccall((:ctx_setmodeltype, libjamsd), Cint, (Ptr{context}, Cint), m.jamsd_ctx, m.model_type)
    res != 0 && error("return code $res from JAMSD")
    res = ccall((:ctx_setobjsense, libjamsd), Cint, (Ptr{context}, Cint), m.jamsd_ctx, sense_to_jamsd[m.sense])
    res != 0 && error("return code $res from JAMSD")
end

function jamsd_declare_mathprgm(mp, ctx::Ptr{context}, emp::Ptr{empinfo})
    # TODO(xhub) that is JuMP-specific
    m = mp.emp.model_ds.internalModel.inner

    jamsd_mp = emp_mp_alloc(emp, ctx)

    CONFIG[:debug] && println("DEBUG: jamsd_declare_mathprgm: tackling MP sense = $(mp.sense) objequ = $(mp.objequ) equs = $(mp.equs) vars = $(mp.vars)")

    if mp.objequ > 0 || mp.objvar > 0
        typ = 0

        emp_mp_start(jamsd_mp, typ)
        emp_mp_objdir(jamsd_mp, mp.sense)
        emp_mp_objequ(jamsd_mp, mp.objequ-1)

        if mp.objvar > 0
            emp_mp_objvar(jamsd_mp, m.v_index_map[mp.objvar])
        else
            emp_mp_objvar(jamsd_mp, mp.objvar-1)
        end

        for eidx in mp.equs
            if eidx != mp.objequ
                emp_mp_constraint(jamsd_mp, eidx-1)
            end
        end

        for vidx in mp.vars
            emp_mp_var(jamsd_mp, m.v_index_map[vidx])
        end
    else
        typ = 2

        emp_mp_start(jamsd_mp, typ)
        VIvars = keys(mp.matching)
        # TODO this should be an array of booleans ... --xhub
        equ_seen = fill(-1, length(mp.equs))
        sidx = 1
        for vidx in mp.vars
            if vidx in VIvars
                eidx = mp.matching[vidx]
                emp_mp_vipair(jamsd_mp, eidx-1, m.v_index_map[vidx])
                equ_seen[sidx] = eidx
                sidx += 1
                CONFIG[:debug] && println("DEBUG: eqn $eidx perp x[$vidx]")
            else
                # QVI case
                emp_mp_var(jamsd_mp, vidx-1)
            end
        end

        for idx in mp.equs
            if !(idx in equ_seen)
                CONFIG[:debug] && println("DEBUG: eqn $eidx is a constraint")
                emp_mp_constraint(jamsd_mp, idx-1)
            end
        end
    end

    return jamsd_mp
end


text:
module TransformManager

export Transform

using ..Management
using CoreExtended
using MatrixMath

""" TODO """
type Transform
	id::Symbol
	enabled::Bool

  position::Vec3f
  rotation::Vec3f
	scaling::Vec3f

	Transform(id=:_) = new(id,true,zeros(Vec3f),zeros(Vec3f),ones(Vec3f))
end

# preset for manager: includes code here
#presetManager(Transform)
m=Management.presetManager(Transform)

end #TransformManager

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
#=
P19 (**) Rotate a list N places to the left.
Examples:
* (rotate '(a b c d e f g h) 3)
(D E F G H A B C)

* (rotate '(a b c d e f g h) -2)
(G H A B C D E F)

Hint: Use the predefined functions length and append, as well as the result of problem P17.
=#

function rotate(coll, n)
    if n > 0
        vcat(coll[n + 1:end], coll[1:n])
    else
        vcat(coll[end + n + 1:end], coll[1:end + n])
    end
end

@assert rotate([:a, :b, :c, :d, :e, :f, :g, :h], 3) == [:d, :e, :f, :g, :h, :a, :b, :c]
@assert rotate([:a, :b, :c, :d, :e, :f, :g, :h], -2) == [:g, :h, :a, :b, :c, :d, :e, :f]
@assert rotate([:a, :b, :c], 0) == [:a, :b, :c]

println("Tests passed: JL-19.jl")


text:
function unnormedCategorical(p̃::AbstractVector)
    normalize!(p̃,1)
    return rand(Categorical(p̃))
end

function urban_transmission(u,p,t)
    @unpack T,β = p
    S_urb = u[:,1,1]
    I_urb = u[:,3,1] .+ u[:,4,1]
    I_rural = u[:,3,2] .+ u[:,4,2]
    Î = T*I_urb + I_rural
    λ_urb = β*T'*(Î ./N̂)
    return sum(S_urb.*λ_urb)
end

function affect_urb_transmission!(integrator)
    @unpack T,β = integrator.p
    u = integrator.u
    S_urb = u[:,1,1]
    I_urb = u[:,3,1] .+ u[:,4,1]
    I_rural = u[:,3,2] .+ u[:,4,2]
    Î = T*I_urb + I_rural
    λ_urb = β*T'*(Î ./N̂)
    # inf_rate = S_urb.*λ_urb
    inf_rate = [S_urb[i]*λ_urb[i] for i = 1:n]
    area = unnormedCategorical(inf_rate)
    u[area,1,1] -= 1
    u[area,2,1] += 1
end
jump_urb_trans = ConstantRateJump(urban_transmission,affect_urb_transmission!)

function rural_transmission(u,p,t)
    @unpack T,β = p
    S_rural = u[:,1,2]
    I_urb = u[:,3,1] .+ u[:,4,1]
    I_rural = u[:,3,2] .+ u[:,4,2]
    Î = T*I_urb + I_rural
    λ_rural = β*(Î./N̂)
    return β*sum(S_rural.*λ_rural)
end
function affect_rural_transmission!(integrator)
    @unpack β,T = integrator.p
    u = integrator.u
    S_rural = u[:,1,2]
    I_urb = u[:,3,1] .+ u[:,4,1]
    I_rural = u[:,3,2] .+ u[:,4,2]
    Î = (T*I_urb) .+ I_rural
    λ_rural = β*(Î./N̂)
    inf_rate = [S_rural[i]*λ_rural[i] for i = 1:n]
    area = unnormedCategorical(inf_rate)
    u[area,1,2] -= 1
    u[area,2,2] += 1
end
jump_rural_trans = ConstantRateJump(rural_transmission,affect_rural_transmission!)

function incubation(u,p,t)
    return p.σ*sum(u[:,2,:])
end

function affect_incubation!(integrator)
    @unpack δ,σ = integrator.p
    u = integrator.u
    E_urb = u[:,2,1]
    E_rural = u[:,2,2]
    incubationrates = σ*(E_urb .+ E_rural)
    choose_disease = 3 + rand(Bernoulli(δ))
    area = unnormedCategorical(incubationrates)
    urb_vs_rural = 1+ rand(Bernoulli(E_rural[area]/(E_urb[area] + E_rural[area])))
    u[area,choose_disease,urb_vs_rural] += 1
    u[area,2,urb_vs_rural] -= 1
    u[area,choose_disease+4,urb_vs_rural] += 1#Add to cumulative cases
end
jump_incubation= ConstantRateJump(incubation,affect_incubation!)

function recovery(u,p,t)
    return p.γ*sum(u[:,3:5,:])
end

function affect_recovery!(integrator)
    @unpack γ = integrator.p
    u = integrator.u
    sub_I = sum(u, dims = 3)[:,3]
    diseased_I = sum(u, dims = 3)[:,4]
    H = sum(u, dims = 3)[:,5]
    area = unnormedCategorical(γ*(sub_I .+ diseased_I .+H))
    group = 2 + unnormedCategorical(γ*[sub_I[area],diseased_I[area],H[area]])
    urb_vs_rural = 1 + rand(Bernoulli(u[area,group,2]/(u[area,group,1]+u[area,group,2])))
    u[area,group,urb_vs_rural] -= 1
    u[area,6,urb_vs_rural] += 1
end
jump_recovery = ConstantRateJump(recovery,affect_recovery!)

function hospitalisation(u,p,t)
    return p.τ*sum(u[:,4,:])
end

function affect_hospitalisation!(integrator)
    @unpack τ = integrator.p
    u = integrator.u
    diseased_I_urb = u[:,4,1]
    diseased_I_rural = u[:,4,2]
    area = unnormedCategorical(τ*(diseased_I_urb .+ diseased_I_rural))
    urb_vs_rural = 1 + rand(Bernoulli(u[area,4,2]/(u[area,4,1]+u[area,4,2])))
    u[area,4,urb_vs_rural] -= 1
    u[area,5,urb_vs_rural] += 1
end
jump_hosp = ConstantRateJump(hospitalisation,affect_hospitalisation!)

function death(u,p,t)
    return p.μ₁*sum(u[:,5,:])
end

function affect_death!(integrator)
    @unpack μ₁ = integrator.p
    u = integrator.u
    H_urb = u[:,5,1]
    H_rural = u[:,5,2]
    area = unnormedCategorical(μ₁*(H_urb .+ H_rural))
    urb_vs_rural = 1 + rand(Bernoulli(u[area,5,2]/(u[area,5,1] +u[area,5,2])))
    u[area,5,urb_vs_rural] -=1
    u[area,9,urb_vs_rural] +=1
end
jump_death = ConstantRateJump(death,affect_death!)

function nairobi_import(u,p,t)
    return sum(p*u[31,1,1])
end

function affect_nairobi_import!(integrator)
    @unpack ϵ_nai = integrator.p
    u = integrator.u
    u[31,1,1] -= 1
    u[31,2,1] += 1
end
jump_nairobi_trans = ConstantRateJump(nairobi_transmission,affect_nairobi_transmission!)

function mombassa_import(u,p,t)
    return sum(p*u[29,1,1])
end

function affect_mombassa_import!(integrator)
    @unpack ϵ_mom = integrator.p
    u = integrator.u
    u[29,1,1] -= 1
    u[29,2,1] += 1
end
jump_mombassa_trans = ConstantRateJump(mombassa_transmission,affect_mombassa_transmission!)


text:
"""
    easter_rata(y::Year) → Int

Returns Easter date as a *[Rata Die](https://en.wikipedia.org/wiki/Rata_Die)* number.

Based on *[Algo R](http://www.linuxtopia.org/online_books/programming_books/python_programming/python_ch38.html)*.
"""
function easter_rata(y::Year)

    local c::Int
    local e::Int
    local p::Int

     # Algo R only works after 1582
     if y.value < 1582
          # Are you using this? Send me a postcard!
          error("Year cannot be less than 1582. Provided: $(y.value).")
     end

    # Century
     c = div(y.value , 100) + 1

     # Shifted Epact
     e = mod(14 + 11*(mod(y.value, 19)) - div(3*c, 4) + div(5+8*c, 25), 30)

     # Adjust Epact
     if (e == 0) || ((e == 1) && ( 10 < mod(y.value, 19) ))
        e += 1
     end

     # Paschal Moon
     p = Date(y.value, 4, 19).instant.periods.value - e

     # Easter: locate the Sunday after the Paschal Moon
     return p + 7 - mod(p, 7)
end

"""
    easter_date(y::Year) → Date

Returns result of `easter_rata` as a `Date` instance.
"""
function easter_date(y::Year)
    # Compute the gregorian date for Rata Die number
     return Date(Dates.rata2datetime( easter_rata(y) ))
end

# weekday_target values:
# const Monday,Tuesday,Wednesday,Thursday,Friday,Saturday,Sunday = 1,2,3,4,5,6,7
# See query.jl on Dates module
# See also dayofweek(dt) function.
"""
    findweekday(weekday_target::Int, yy::Int, mm::Int, occurrence::Int, ascending::Bool) → Date

Given a year `yy` and month `mm`, finds a date where a choosen weekday occurs.

`weekday_target` values are declared in module `Base.Dates`:
`Monday,Tuesday,Wednesday,Thursday,Friday,Saturday,Sunday = 1,2,3,4,5,6,7`.

If `ascending` is true, searches from the beggining of the month. If false, searches from the end of the month.

If `occurrence` is `2` and `weekday_target` is `Monday`, searches the 2nd Monday of the given month, and so on.
"""
function findweekday(weekday_target::Int, yy::Int, mm::Int, occurrence::Int, ascending::Bool)

    local dt::Date = Date(yy, mm, 1)
    local dt_dayofweek::Integer
    local offset::Integer

    if occurrence <= 0
        error("occurrence must be >= 1. Provided $(occurrence).")
    end

    if ascending
        dt_dayofweek = dayofweek(dt)
        offset = rem(weekday_target + 7 - dt_dayofweek, 7) # rem = MOD function
    else
        dt = lastdayofmonth(dt)
        dt_dayofweek = dayofweek(dt)
        offset = rem(dt_dayofweek + 7 - weekday_target, 7)
    end

    if occurrence > 1
        offset += 7 * (occurrence - 1)
    end

    if ascending
        return dt + Dates.Day(offset)
    else
        return dt - Dates.Day(offset)
    end
end


text:
using WordCloud
function drawjuliacircle(sz)
    juliacirclessvg = WordCloud.Render.Drawing(sz, sz, :svg)
    WordCloud.Render.origin()
    WordCloud.Render.background(0, 0, 0, 0)
    WordCloud.Render.juliacircles(sz ÷ 4)
    WordCloud.Render.finish()
    juliacirclessvg
end

docs = (readdir(joinpath(dirname(Sys.BINDIR), "share/doc/julia/html/en", dir), join=true) for dir in ["manual", "base", "stdlib"])
docs = docs |> Iterators.flatten

words, weights = processtext(maxnum=400, maxweight=1) do
    counter = Dict{String,Int}()
    for doc in docs
        content = html2text(open(doc))
        countwords(content, counter=counter)
    end
    counter
end

wc = wordcloud(
    [words..., "∴"], # add a placeholder for julia-logo
    [weights..., weights[1]],
    density=0.65,
    mask=shape(box, 900, 300, cornerradius=0, color=0.95),
    colors=((0.796, 0.235, 0.20), (0.584, 0.345, 0.698), (0.22, 0.596, 0.149)),
    angles=(0, -45, 45),
    # font = "Verdana Bold",
)
setangles!(wc, "julia", 0)
# setangles!(wc, "function", 45)
# initword!(wc, "function")
setcolors!(wc, "julia", (0.796, 0.235, 0.20))
# setfonts!(wc, "julia", "forte")
initword!(wc, "julia")
juliacircles = drawjuliacircle(getfontsizes(wc, "∴") |> round)
setsvgimages!(wc, "∴", juliacircles) # replace image
sz1 = size(getimages(wc, "∴"))
sz2 = size(getimages(wc, "julia"))
y1, x1 = (size(wc.mask) .- (sz1[1], sz1[2] + sz2[2])) .÷ 2
y2 = (size(wc.mask, 1) - sz2[1]) ÷ 2
x1 = round(Int, x1 * 0.9)
setpositions!(wc, "∴", (x1, y1))
setpositions!(wc, "julia", (x1 + sz1[2], y2))

pin(wc, ["julia", "∴"]) do
    placewords!(wc)
    generate!(wc, 2000)
end
println("results are saved to juliadoc.svg")
# paint(wc, "juliadoc.png")
paint(wc, "juliadoc.svg")
wc
#eval# runexample(:juliadoc)
#md# ![](juliadoc.svg)


text:
using PdfFraud
using GenerativeModels
using Flux
using ValueHistories
using ArgParse
using BSON
using GMExtensions
using CuArrays
using DelimitedFiles, DataFrames, CSV
using PyPlot

# this is for running the script from REPL - default values
if length(ARGS) == 0
	push!(ARGS, "vae_runs/2")
end
# parse arguments
s = ArgParseSettings()
@add_arg_table s begin
    "modelpath"
    	default = ""
    	arg_type = String
    	help = "a master dir of where the models are saved"
    "--imodels"
    	default = [0]
    	help = "only compute models given by indices"
    	arg_type = Int
    	nargs = '+'
end
parsed_args = parse_args(ARGS, s)

# extract settings from the parsed arguments
modelpath = parsed_args["modelpath"]
modelpath = (length(modelpath) > 0 && modelpath[1] == "/") ? modelpath :
	joinpath(PdfFraud.experiment_path, modelpath, "models")
reprpath = joinpath(dirname(modelpath), "representations")
mkpath(reprpath)
imodels = parsed_args["imodels"]

# get labeled data
data, labels, sha256, page_nums = PdfFraud.load_pdf_data()
data = data |> gpu;

# save the page info
lf = joinpath(reprpath, "index.csv")
index = DataFrame()
index[!,:sha256] = vcat(map(x->repeat([x[1]], x[2]), zip(sha256, page_nums))...)
index[!,:page] = vcat(map(x->1:x, page_nums)...)
CSV.write(lf, index)

# load model(s)
model_files = joinpath.(modelpath, readdir(modelpath))
@info "Found $(length(model_files)) saved models."
if imodels[1] != 0
	model_files = model_files[imodels]
end
@info "Continuing with $(length(model_files)) models."

# now get the representations

# now lets check the encodings
zfs = joinpath.(reprpath, filter(x->x!="index.csv", readdir(reprpath)))
Zs = map(zf -> readdlm(zf, ',') ,zfs)

linds = map(x->x!=nothing,labels)
Y = labels[linds];

figure(figsize=(10,25))
for iz in 41:50#min(10,length(Zs))
	subplot(5,2,iz%10+1)
	Z = Zs[iz][linds,:];
	for y in unique(Y)
		inds = Y .== y
		scatter(Z[inds,1], Z[inds,2])
	end
	title(basename(zfs[iz]), fontsize=8)
end
tight_layout()


text:
"""
    Polynomial

one of functions mapping the set of real numbers is the algebraic polynomials, the set of functions of the form:
``p_n(x) = a_n x^n + a_{n-1} x^{n-1} + ... +a_{1}x + a_0``
- nth Larange interpolating polynomial
- Neville’s Iterated Interpolation
- Newton’s Divided-Difference Formula
- Natural Cubic Spline
- Clamped Cubic Spline
"""
module Polynomial
using SymPy
using LinearAlgebra
using Latexify

@inline function product(x::Vector)
    res = 1
    for i in 1:size(x)[1]
        res = res*x[i]
    end
    return res
end

"""
    Lagrange(X::Vector, Y::Vector; a::String="x")

use the ``(x_0, y_0), (x_1, y_1), ..., (x_n, y_n)``, to caculate a polynomial.
input X and Y, default use x to express functions.
output: function.
"""
@inline function Lagrange(X::Vector, Y::Vector; a::String="x")
    a = symbols(a)
    m = size(X)[1]
    A = ones(m) * X'
    An = A - A'
    An[An .== 0] .= 1
    up = a .- X
    upr = product(up)
    sum = 0
    for i in 1:m
        down = product(An[:,i])
        u = upr/up[i]
        sum = sum + Y[i] * (u/down)
    end
    return simplify(sum)
end

"""
    Neville(x::Vector, y::Vector, x₀::Real; tab::Bool=false)

input two vector, x and y. satisfy the function `y=f(x)`, then input x₀, which you want
to  approximate f(x₀). The default tab setting is false, just output the estimate. you can set true
to output a result table.
"""
@inline function Neville(x::Vector, y::Vector, x₀::Real; tab::Bool=false)
    n, = size(x)
    res = zeros(n, n)
    for i in 1:n
        res[i, 1] = y[i]
    end
    for i in 2:n
        for j in 2:i
            res[i,j] = ((x₀ - x[i-j+1])*res[i, j-1] - (x₀ - x[i])*res[i-1,j-1])/(x[i] - x[i-j+1])
        end
    end
    if tab
        return res
    else
        return res[n, n]
    end
end

"""
    NDDF(x::Vector, y::Vector; a::String="x", simple::Bool=true, tab::Bool=false, backward::Bool=false)

input ``x_1, x_2, x_3, ...x_n```,values ``f(x_1), f(x_2), ...,f(x_n)``, set Args `simpel = true`
output simplify ans, set tab to output a table. set `backward=true` use the Newton Backward–Difference Formula.
"""
@inline function NDDF(x::Vector, y::Vector; a::String="x", simple::Bool=true, tab=false, backward::Bool=false)
    a = symbols(a)
    n, = size(x)
    res = zeros(n, n)
    for i in 1:n
        res[i,1] = y[i]
    end
    for i in 2:n
        for j in 2:i
            res[i,j] = (res[i, j-1] - res[i-1, j-1])/(x[i] - x[i-j+1])
        end
    end
    sum = res[1,1]
    for i in 2:n
        p = 1
        for j in 1:i-1
            p = p*(a - x[j])
        end
        if backward
            sum = sum + res[n, i]*p
        else
            sum = sum + res[i, i]*p
        end
    end
    if tab
        return res
    else
        if simple
            return simplify(sum)
        else
            return sum
        end
    end
end

"""
    NCSpline( NCSpline(x::Vector, y::Vector; Latex::Bool=false, a::String="t")
input ``x_1, x_2, x_3, ...x_n```,values ``f(x_1), f(x_2), ...,f(x_n)``, set Args `Latex = true`  to
output information about the function.
"""
@inline function NCSpline(x::Vector, y::Vector; Latex::Bool=false, a::String="t")
    a = symbols(a)
    n, = size(x)
    h, α = zeros(n), zeros(n-1)

    for i in 1:n-1
        h[i] = x[i+1] - x[i]
    end
    for i in 2:n-1
        α[i] = 3/h[i] * (y[i+1] - y[i]) - 3/h[i-1] *(y[i] - y[i-1])
    end

    l = ones(n)
    μ = zeros(n)
    z = zeros(n)

    for i in 2:n-1
        l[i] = 2*(x[i+1] - x[i-1]) - h[i-1]*μ[i-1]
        μ[i] = h[i]/l[i]
        z[i] = (α[i] - h[i-1]*z[i-1])/l[i]
    end

    c, b, d = zeros(n), zeros(n), zeros(n)
    for j in reverse(1:n-1)
        c[j] = z[j] - μ[j]*c[j+1]
        b[j] = (y[j+1] - y[j])/h[j] - h[j]*(c[j+1]+2*c[j])/3
        d[j] = (c[j+1] - c[j])/(3*h[j])
    end
    A = [b[1:end-1] c[1:end-1] d[1:end-1]]
    X = []
    for i in 2:4
        for j in 1:3
            push!(X,(a - x[i-1])^(j))
        end
    end
    res = diag(A * reshape(X, (n-1,n-1))) + y[1:end-1]

    if Latex
        for i in 1:n-1
            print(res[i], ",","t ∈","[$(x[i]), $(x[i+1])]")
            print("\n")

        end
    else
        return res
    end
end

"""
    CCSpline(x::Vector, y::Vector, endpoint::Vector; Latex::Bool=false, a::String="t")

input ``x_1, x_2, x_3, ...x_n```,values ``f(x_1), f(x_2), ...,f(x_n)`` and ``f'(x_1), f'(x_n)``, set Args `Latex = true`  to
output information about the function.
"""
@inline function CCSpline(x::Vector, y::Vector, endpoint::Vector; Latex::Bool=false, a::String="t")
    n, = size(x)
    h = ones(n)
    α = ones(n)
    a = symbols(a)

    for i in 1:n-1
        h[i] = x[i+1] - x[i]
    end
    α[1] = 3(y[2] - y[1])/h[1] - 3*endpoint[1]
    α[n] = 3*endpoint[2] - 3*(y[n] - y[n-1])/h[n-1]
    for i in 2:n-1
        α[i] = 3/h[i] * (y[i+1] - y[i]) - 3/h[i-1] *(y[i] - y[i-1])
    end

    l = ones(n)
    l[1] = 2*h[1]
    μ = fill(0.5, n)
    z = fill(α[1]/l[1], n)
    for i in 2:n-1
        l[i] = 2*(x[i+1] - x[i-1]) - h[i-1]*μ[i-1]
        μ[i] = h[i]/l[i]
        z[i] = (α[i] - h[i-1]*z[i-1])/l[i]
    end
    l[n] = h[n-1]*(2 - μ[n-1])
    z[n] = (α[n] - h[n-1]*z[n-1])/l[n]
    c = fill(z[n], n)
    b, d = zeros(n), zeros(n)
    for j in reverse(1:n-1)
        c[j] = z[j] - μ[j]*c[j+1]
        b[j] = (y[j+1] - y[j])/h[j] - h[j]*(c[j+1]+2*c[j])/3
        d[j] = (c[j+1] - c[j])/(3*h[j])
    end
    A = [b[1:end-1] c[1:end-1] d[1:end-1]]
    X = []
    for i in 2:4
        for j in 1:3
            push!(X,(a - x[i-1])^(j))
        end
    end
    res = diag(A * reshape(X, (n-1,n-1))) + y[1:end-1]
    if Latex
        for i in 1:n-1
            print(res[i], ",","t ∈","[$(x[i]), $(x[i+1])]")
            print("\n")

        end
    else
        return res
    end

end

end


text:
##############################################################################################################################
# Call functions (used in map/pmap etc ...)
##############################################################################################################################

# DataCell(...)
(x::T where T<:CellData)(f::S where S<:CellData) = x |> f;
(x::T where T<:CellData)(f::S where S<:CellFun) = x |> f;
(x::T where T<:CellData)(p::S where S<:PipeGeneric) = x |> p

# FunctionCell(...)
(f::T where T<:CellFun)(x::S where S<:CellData) =  x |> f
(f1::T where T<:CellFun)(f2::S where S<:CellFun) = f1 + f2
(f::T where T<:CellFun)(p::S where S<:PipeGeneric) = f + p
(f::T where T<:CellFun)(x::S where S) =  x |> f

# PipeCell(...)
(p::T where T<:PipeGeneric)(x::S where S<:CellData) = x |> p
(p::T where T<:PipeGeneric)(f::S where S<:CellFun) = p + f
(p1::T where T<: PipeGeneric)(p2::S where S<:PipeGeneric) = p1 + p2
(p::T where T<:PipeGeneric)(x::S where S) = x |> p

# Generic calls
(f::T where T<:CellFun)(x::S where S<:AbstractCell) = begin
	try
		getx!(x) |> f
	catch
		warn("[call] Passing generic cell content into function cell failed, will make a pipe.")
		x+f
	end

end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
using Test

# -------------------------------------------------------------------------------------------------
# Test constructing and freeing an instance.
@testset "TessInst constructor" begin
    @test_nowarn TessInst("eng", datadir)

    local inst = TessInst("eng", datadir)
    @suppress @test_nowarn show(inst)
    tess_delete!(inst)
    @suppress @test_nowarn show(inst)
end

# -------------------------------------------------------------------------------------------------
# Test tess_delete!()
@testset "TessInst tess_delete!" begin
    local tess = TessInst("eng", datadir)

    @test_nowarn tess_delete!(tess)
    @test_nowarn tess_delete!(tess)
end

# -------------------------------------------------------------------------------------------------
# Test tess_recognize()
@testset "tess_recognize" begin
    local inst = TessInst("eng", datadir)
    local pix  = pix_with()

    # ---------------------------------------------------------------------------------------------
    # Normal tests.
    @suppress @test tess_recognize(inst) == false
    @test_nowarn tess_image(inst, pix)
    @test_nowarn tess_resolution(inst, 72)
    @test tess_recognize(inst) == true

    # ---------------------------------------------------------------------------------------------
    # Test being called with a freed instance.
    local err = "Instance has been freed."
    tess_delete!(inst)
    @test (@test_logs (:error, err) tess_recognize(inst)) == false
end


text:
struct DistanceSphereToSphere{M,R1,S,m} <: TaskMap{M,R1,S}
    center::SVector{m,S}
    radius1::S
    radius2::S
end
DistanceSphereToSphere{M,R1}(c::SVector{m,S}, r1::S, r2::S) where {M,R1,S,m} =
    DistanceSphereToSphere{M,R1,S,dim(M)}(c, r1, r2)

struct DistanceSphereToSphereT{M,R1,S,m} <: TaskMapT{M,R1,S}
    base_map::DistanceSphereToSphere{M,R1,S,m}
end
DistanceSphereToSphereT{M,R1}(c::SVector{m,S}, r1::S, r2::S) where {M,R1,S,m} =
    DistanceSphereToSphereT{M,R1,S,m}(DistanceSphereToSphere{M,R1}(c, r1, r2))

task_map_emb(::EmbRep, ::EmbRep, xme, task_map::DistanceSphereToSphere{ℝ{m},R1}) where m =
    SA[norm(xme - task_map.center) - task_map.radius1 - task_map.radius2]


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

# default maximum lag to be used in various methods
default_maxlag(data) = 0.1diagonal(boundingbox(data))

"""
    EmpiricalVariogram(data, var₁, var₂=var₁; [optional parameters])

Computes the empirical (a.k.a. experimental) omnidirectional
(cross-)variogram for variables `var₁` and `var₂` stored in
spatial `data`.

## Parameters

  * nlags    - number of lags (default to `20`)
  * maxlag   - maximum lag (default to half of maximum lag of data)
  * distance - custom distance function (default to `Euclidean` distance)
  * algo     - accumulation algorithm (default to `:ball`)

Available algorithms:

  * `:full` - loop over all pairs of points in the data
  * `:ball` - loop over all points inside maximum lag ball

All implemented algorithms produce the exact same result.
The `:ball` algorithm is considerably faster when the
maximum lag is much smaller than the bounding box of
the data.

The function [`values`](@ref) can be used to retrieve
the abscissa, ordinate and bin count of an empirical
variogram:

```julia
julia> x, y, n = values(γ)
```

See also: [`DirectionalVariogram`](@ref)

## References

* Chilès, JP and Delfiner, P. 2012. [Geostatistics: Modeling Spatial Uncertainty]
  (https://onlinelibrary.wiley.com/doi/book/10.1002/9781118136188)
* Hoffimann, J and Zadrozny, B. 2019. [Efficient variography with partition variograms]
  (https://www.sciencedirect.com/science/article/pii/S0098300419302936)
"""
struct EmpiricalVariogram{V,D}
  abscissa::Vector{Float64}
  ordinate::Vector{V}
  counts::Vector{Int}
  distance::D
end

function EmpiricalVariogram(data, var₁::Symbol, var₂::Symbol=var₁;
                            nlags=20, maxlag=default_maxlag(data),
                            distance=Euclidean(), algo=:ball)

  # retrieve table and domain
  𝒯 = values(data)
  𝒟 = domain(data)

  # retrieve the column names of data values
  vars = Tables.columnnames(𝒯)

  # sanity checks
  @assert nelements(data) > 1 "variogram requires at least 2 elements"
  @assert (var₁, var₂) ⊆ vars "invalid variable names"
  @assert algo ∈ (:full, :ball) "invalid accumulation algorithm"
  @assert nlags  > 0 "number of lags must be positive"
  @assert maxlag > 0 "maximum lag distance must be positive"

  # ball search with NearestNeighbors.jl requires AbstractFloat and MinkowskiMetric
  # https://github.com/KristofferC/NearestNeighbors.jl/issues/13
  isfloat     = coordtype(data) <: AbstractFloat
  isminkowski = distance isa MinkowskiMetric

  # warn users requesting :ball option with invalid parameters
  (algo == :ball && !isfloat) && @warn ":ball algorithm requires floating point coordinates"
  (algo == :ball && !isminkowski) && @warn ":ball algorithm requires Minkowski metric"

  # empirical variograms are defined on point sets
  𝒫 = PointSet([centroid(𝒟, i) for i in 1:nelements(𝒟)])
  𝒮 = georef(𝒯, 𝒫)

  # choose accumulation algorithm
  if algo == :ball && isfloat && isminkowski
    xsums, ysums, counts = ball_search_accum(𝒮, var₁, var₂, maxlag, nlags, distance)
  else
    xsums, ysums, counts = full_search_accum(𝒮, var₁, var₂, maxlag, nlags, distance)
  end

  # bin (or lag) size
  δh = maxlag / nlags
  lags = range(δh/2, stop=maxlag - δh/2, length=nlags)

  # variogram abscissa
  abscissa = @. xsums / counts
  abscissa[counts .== 0] .= lags[counts .== 0]

  # variogram ordinate
  ordinate = @. (ysums / counts) / 2
  ordinate[counts .== 0] .= zero(eltype(ordinate))

  EmpiricalVariogram(abscissa, ordinate, counts, distance)
end

"""
    values(γ)

Returns the abscissa, the ordinate, and the bin counts
of the empirical variogram `γ`.
"""
Base.values(γ::EmpiricalVariogram) = γ.abscissa, γ.ordinate, γ.counts

"""
    distance(γ)

Retun the distance used to compute the empirical variogram `γ`.
"""
distance(γ::EmpiricalVariogram) = γ.distance

"""
    merge(γα, γβ)

Merge the empirical variogram `γα` with the empirical variogram `γβ`
assuming that both variograms have the same number of lags.
"""
function merge(γα::EmpiricalVariogram{D}, γβ::EmpiricalVariogram{D}) where {D}
  xα = γα.abscissa
  xβ = γβ.abscissa
  yα = γα.ordinate
  yβ = γβ.ordinate
  nα = γα.counts
  nβ = γβ.counts

  n = nα + nβ
  x = @. (xα*nα + xβ*nβ) / n
  y = @. (yα*nα + yβ*nβ) / n

  # adjust empty bins
  x[n .== 0] .= xα[n .== 0]
  y[n .== 0] .= 0

  d = γα.distance

  EmpiricalVariogram(x, y, n, d)
end

# -----------
# IO METHODS
# -----------

function Base.show(io::IO, ::EmpiricalVariogram)
  print(io, "EmpiricalVariogram")
end

function Base.show(io::IO, ::MIME"text/plain", γ::EmpiricalVariogram)
  println(io, γ)
  println(io, "  abscissa: ", extrema(γ.abscissa))
  println(io, "  ordinate: ", extrema(γ.ordinate))
  print(  io, "  N° pairs: ", sum(γ.counts))
end

# ------------------------
# ACCUMULATION ALGORITHMS
# ------------------------

function full_search_accum(data, var₁, var₂, maxlag, nlags, distance)
  # retrieve table and point set
  𝒯 = values(data)
  𝒫 = domain(data)

  # collect vectors for variables
  cols = Tables.columns(𝒯)
  Z₁   = Tables.getcolumn(cols, var₁)
  Z₂   = Tables.getcolumn(cols, var₂)

  # lag size
  δh = maxlag / nlags

  # accumulation type
  V = typeof((Z₁[1] - Z₂[1]) ⋅ (Z₁[1] - Z₂[1]))

  # lag sums and counts
  xsums  = zeros(nlags)
  ysums  = zeros(V, nlags)
  counts = zeros(Int, nlags)

  # loop over all pairs of points
  @inbounds for j in 1:nelements(𝒫)
    pⱼ  = 𝒫[j]
    z₁ⱼ = Z₁[j]
    z₂ⱼ = Z₂[j]
    for i in j+1:nelements(𝒫)
      pᵢ  = 𝒫[i]
      z₁ᵢ = Z₁[i]
      z₂ᵢ = Z₂[i]

      # evaluate spatial lag
      h = evaluate(distance, coordinates(pᵢ), coordinates(pⱼ))
      h > maxlag && continue # early exit if out of range

      # evaluate (cross-)variance
      v = (z₁ᵢ - z₁ⱼ) ⋅ (z₂ᵢ - z₂ⱼ)

      # bin (or lag) where to accumulate result
      lag = ceil(Int, h / δh)
      lag == 0 && @warn "duplicate coordinates found"

      if 0 < lag ≤ nlags && !ismissing(v)
        xsums[lag] += h
        ysums[lag] += v
        counts[lag] += 1
      end
    end
  end

  xsums, ysums, counts
end

function ball_search_accum(data, var₁, var₂, maxlag, nlags, distance)
  # retrieve table and point set
  𝒯 = values(data)
  𝒫 = domain(data)

  # collect vectors for variables
  cols = Tables.columns(𝒯)
  Z₁   = Tables.getcolumn(cols, var₁)
  Z₂   = Tables.getcolumn(cols, var₂)

  # lag size
  δh = maxlag / nlags

  # accumulation type
  V = typeof((Z₁[1] - Z₂[1]) ⋅ (Z₁[1] - Z₂[1]))

  # lag sums and counts
  xsums  = zeros(nlags)
  ysums  = zeros(V, nlags)
  counts = zeros(Int, nlags)

  # fast ball search
  ball = MetricBall(maxlag, distance)
  searcher = BallSearch(𝒫, ball)

  # loop over points inside ball
  @inbounds for j in 1:nelements(𝒫)
    pⱼ  = 𝒫[j]
    z₁ⱼ = Z₁[j]
    z₂ⱼ = Z₂[j]
    for i in search(pⱼ, searcher)
      i ≤ j && continue # avoid double counting
      pᵢ  = 𝒫[i]
      z₁ᵢ = Z₁[i]
      z₂ᵢ = Z₂[i]

      # evaluate spatial lag
      h = evaluate(distance, coordinates(pᵢ), coordinates(pⱼ))

      # evaluate (cross-)variance
      v = (z₁ᵢ - z₁ⱼ) ⋅ (z₂ᵢ - z₂ⱼ)

      # bin (or lag) where to accumulate result
      lag = ceil(Int, h / δh)
      lag == 0 && @warn "duplicate coordinates found"

      if 0 < lag ≤ nlags && !ismissing(v)
        xsums[lag] += h
        ysums[lag] += v
        counts[lag] += 1
      end
    end
  end

  xsums, ysums, counts
end

# temporary fix for ⋅ with missing values
# https://github.com/JuliaLang/julia/issues/40743
⋅(::Missing, ::Missing) = missing


text:
using MacroTools: @q

struct Loop end

function DataFlow.toexpr(f::Loop, λ, out, is...)
  exs = MacroTools.block(λ).args
  @capture(exs[end], (args__,) -> body_) || return :($(Loop())($λ, $out, $(is...)))
  isempty(is) && (is = Iterators.repeated(Symbol("<undef>")))
  iters = Expr(:block, [:($x=$y) for (x, y) in zip(args, is)]...)
  loop = Expr(:for, iters, body)
  @q begin
    $(exs[1:end-1]...)
    $loop
    $out
  end
end

# Inline all temporary arrays
function rmtemps(v::IVertex)
  v = λopen(v)
  deps = dependents(v)
  v = DataFlow.postwalk!(v) do v
    (iscall(v, getindex) && value(v[2]) isa Loop && length(deps[v[2]]) == 1) || return v
    body = λclose(v[2,1])
    body = DataFlow.spliceinputs(body.value.body, body.inputs..., v.inputs[3:end]...)
    @assert iscall(body, setindex!)
    return body[3]
  end
  λclose(v)
end

function accum!(op, var, val)
  var[] = op(var[], val)
  return
end

@tk function reduce(f, v0, xs)
  [i] -> (out = Ref{Any}(v0)) => accum!(f, out, xs[i])
end

max_stars_repo_name: JuliaTagBot/Tokamak

max_issues_repo_name: JuliaTagBot/Tokamak

max_forks_repo_name: FluxML/Tokamak


text: const polyCoeffs240 =  [1, 1547827339295914963, 750097436599629445, 604784063937344339, 1271955066541750687, 1470694949862035313, 1243845019994552891, 1650351806820104120, 2130371604632259008, 1212189416326483249, 769753200018148741, 2296686257144156462, 859810657562839084, 333487777687686012, 1057292136184957112, 363693033707135167, 1584773882785214289, 838875284686242637, 1093679796858682774, 1316510027393789289, 231506995470207753, 237402142406808300, 1841929466371482366, 267038393335434283, 1673541231392068341, 2132729040348167035, 1323611672972245071, 2251387572802242477, 1849474270259047765, 1439565750301017841, 1081538560132929511, 90770412441711845, 42017911305177385, 723182840872907252, 917426827997895711, 1662687298091796154, 1810900253425438293, 324016555806047026, 1749745103146061464, 1054148991637098924, 350458865663042008, 1779016697074766349, 1211307133118578020, 867512037603801700, 288724544640042135, 1126969982050969258, 1839217483181604223, 518757867725905742, 1676441379238930084, 1476311745810790030, 1480861890393073092, 556637090378129049, 455342924775446017, 887497743742028692, 797007306523672598, 830670442326060780, 451286479655240195, 1759989367550666568, 86915997953418891, 426949377800588473, 2285824022782725194, 1151027105792548185, 550252000900659228, 354091764183966270, 2032795167925932002, 2199221944139777637, 780471592146373291, 2231921365164269566, 1851901169598918075, 140113195605090396, 1353628688323112462, 162537673241432467, 1592721132148202971, 1744837105826825987, 1397302899870616114, 1157875316102389319, 2176533313906351715, 120389027538781101, 351486096123753686, 224417537071118550, 2277350329548037899, 1949088665271210168, 1976516773334723536, 959309607677407011, 1916978920654628344, 2212414480778710795, 355392287161828980, 783935032974299243, 1713981722173263421, 654887489236103022, 925582323427035027, 707517829869537958, 1083562472932415398, 445019679064415640, 2192194756151073659, 1386043007147474441, 159843582405934225, 1141062712108981653, 1762343867915050320, 2251574297665438284, 2023542522369988614, 2191648684154694604, 1184036338684177649, 1285034448333455030, 290485676696016822, 213839997608644114, 2297762437359529086, 2043374560060668135, 274205219642385478, 1078271998162098290, 1549258621467896508, 158642893168718536, 475157220527607674, 700112152720402965, 1282474771641140547, 1670132011286277521, 107842274436726759, 388101247182244832, 2085380049483372198, 1935729167444287315, 2170300308000700348, 496986891144163950, 1770501265920493620, 2296965632479061796, 720357698320081269, 717859941225115125, 2010726538763802712, 2047338551148458926, 168518658098489881, 1743815482675216709, 394482365539068770, 1784202883994446021, 2221301488379738228, 2194868754897163475, 53615190408416683, 422659426207623140, 497057675777612163, 388747796955029337, 745209608479495978, 1404501305911586160, 1373617901449301043, 1730038721698058250, 369267192725121047, 324374980050137665, 515017073541639506, 1412750200631374939, 1022446327954834517, 1293575611388039917, 1496163824512106942, 1339377782028162474, 289461904370031844, 443363836550641481, 1756211325868169709, 1168360282606641827, 1670965127236174407, 285886342464666033, 320685127801992871, 1829121914815803042, 999316482512205945, 562472117383593959, 905412501732872461, 1295269365455611835, 1125026401768115222, 894594128082204542, 1583591029710780012, 665206848563181410, 46943504701469505, 632185758202463829, 1327802460807221826, 1328263496101079979, 2200012762604574954, 1539105106426148896, 997365122099163575, 2185200321605649001, 1664657111101620907, 142014510944189293, 1158114054654491853, 935788610046390263, 659394534870007553, 1075761873900949604, 984706060682777015, 1475412463424495849, 1063562063908425150, 1127555948209467738, 930319747378944022, 435366201697401635, 2128864348361181526, 284871845180688067, 666609654061752675, 309680133825086297, 1303885603752481458, 491658597388373745, 998123405889782846, 19200982783974250, 451821309786076120, 1858721772114928664, 715147319681428375, 472810999624273854, 2187117497156092536, 272571963996680260, 78239773931032697, 2097855418581091296, 1994241269917186343, 1121276429600727953, 249684125427174941, 1221890798949056842, 664778794280501631, 524711698228357758, 1060729160058730257, 273795353724894827, 1483753628583389638, 28765196066229742, 1306265322032484129, 1702136928232713186, 770557375707640426, 2105292271365658580, 2215227204923328982, 321876906940937052, 487529166448461396, 2115700404867419492, 174418936197175522, 1753016082221364638, 1415256826652670517, 931799821254113175, 1263201645408095639, 2073683068646729292, 1548978890552512492, 1975160376345701791, 774874220929733277, 300798243423740172, 247943417452409322, 981209458062432143, 1266382982589353078, 1800780586730386109, 899047367627219020, 236522531257453874, 574917585520884546, 1298755622256427695, 2296070755750404630, 2305843009213693472, 1]


text:
using Test

@testset "relu.jl" begin

@testset "ReLU" begin
    @testset "Initialize without tightening algorithm" begin
        p = ReLU()
        @test p.tightening_algorithm === nothing
    end

    @testset "Initialize with tightening algorithm" begin
        p = ReLU(MIPVerify.mip)
        @test !(p.tightening_algorithm === nothing)
        @test p.tightening_algorithm == MIPVerify.mip
    end

    @testset "Base.show" begin
        p = ReLU()
        io = IOBuffer()
        Base.show(io, p)
        @test String(take!(io)) == "ReLU()"
    end
end

end


text:
using Plots
using RecipesBase
using LinearAlgebra
using SparseArrays
using Distributed

include("C:/Users/lucas/OneDrive/Documents/GitHub/GDB.jl/src/prim/walls.jl")
include("C:/Users/lucas/OneDrive/Documents/GitHub/GDB.jl/src/prim/grids.jl")
include("C:/Users/lucas/OneDrive/Documents/GitHub/GDB.jl/src/prim/matrices.jl")
include("C:/Users/lucas/OneDrive/Documents/GitHub/GDB.jl/src/prim/splines.jl")
include("C:/Users/lucas/OneDrive/Documents/Dartmouth/Research/Code/GDB.jl/src/prim/projection.jl")
include("C:/Users/lucas/OneDrive/Documents/GitHub/GDB.jl/src/prim/multigrid.jl")
include("C:/Users/lucas/OneDrive/Documents/GitHub/GDB.jl/src/ideal_mhd/setup.jl")

# derivatives from matrices without projections
# call simple to put into multigrid_grids
# have v_cycle_general spit out residual after one time through
# norm of residual divided by norm of rhs
# return x and E_fine, if norm(E_fine)/norm(b) stop
# A - laplacian
# b - f_to_grid, first argument is function (say, (x,y) -> sin(x) ), second argument is grid
# vec to mesh spits out vector that you can apply a heatmap to using above notation
# if singular, b = transpose(A) * b, A = transpose(A) * A

grid, wall, deltas, corners = simple(60, 60, 10, 1, 20, 20, 3)
multiresgrids = multigrid_grids(grid, [wall], deltas, corners, [3,3])
A = laplacian(grid)
b = f_to_grid((x,y) -> sin(x), grid)
x = A \ b

plotvar(x, grid, :heatmap)


text:
using Makie
import DataStructures: DefaultDict

export Port

@recipe(Steps, t, y) do scene
    Theme(
        color = :black,
        fill = nothing,
    )
end

function Makie.plot!(myplot::Steps)
    t = myplot[1]
    y = myplot[2]

    line = Node(Point2f0[])
    function update!(t, y)
        empty!(line[])
        y_old = y[1]
        for (_t,_y) in zip(t,y)
            push!(line[], Point2f0(_t,y_old))
            push!(line[], Point2f0(_t,_y))
            y_old = _y
        end
    end

    Makie.Observables.onany(update!, t, y)
    update!(t[],y[])

    poly!(myplot, line, color=myplot[:fill])
    lines!(myplot, line, color=myplot[:color])
    myplot
end

@recipe(StateTrace, t, state, color, offset, height) do scene
    Theme(
    )
end

function Makie.plot!(plt::StateTrace)
    t,state,color,offset,height=plt[:t],plt[:state],plt[:color],plt.offset,plt.height
    rects = Node(Rect2D{Float64}[])
    colors = Node(Any[])

    Makie.Observables.onany(t, state, color, offset, height) do t,state,color, offset, height
        t₀,s₀ = first(t),first(state)
        for (t₁,s₁) in zip(t[2:end],state[2:end])
            push!(rects[], Rect(float(t₀), float(offset), float(t₁-t₀), float(height)))
            push!(colors[], color[s₀])
            t₀,s₀ = t₁,s₁
        end
    end

    # Force update
    t[] = t[]

    poly!(plt, rects, color=colors)

    plt
end

struct Port{W,I}
    name::I
end
Port{W}(name::I) where {W,I} = Port{W,I}(name)

function Makie.default_theme(scene::Makie.SceneLike, ::Type{<: Makie.Plot(Neuron)})
    Theme(
        color=RGBAf0(0.9,0.9,0.9,1.0),
        branch_width=0.1,
        branch_length=1.0,
        angle_between=10/180*pi,
        xticks = [],
        root_position = Point2f0(0,0),
        ports = Dict{Symbol,Vector{Symbol}}(),
        fxaa = false,
    )
end

function Makie.plot!(treeplot::Makie.Plot(Neuron))
    # get the actual object to plot
    tree = to_value(treeplot[1])
    # get position of the root = offset of all nodes
    offset = treeplot[:root_position]

    # if obj is an observable dict-type holding named values, get named value ...
    maybe_get(obj::Node, key) = if isa(obj[],Union{DefaultDict,Dict})
        @lift $obj[key]
    else
        obj
    end

    # if obj is a dict-type holding named values, get named value ...
    maybe_get(obj::Union{DefaultDict,Dict}, key) = obj[key]
    # ... otherwise return the obj itself
    maybe_get(obj, key) = obj

    function serialize_tree(tree::NeuronOrSegment{ID,T,WT,IT}, l, w, ω, given_ports) where {ID,T,WT,IT}
        sectors = Float64[]
        depths = Float64[]
        all_points = Vector{Pair{ID,Point2f0}}[]
        all_ports = Vector{Pair{ID,Point2f0}}[]
        all_points_flat = Pair{ID,Point2f0}[]
        all_ports_flat = Pair{ID,Point2f0}[]
        all_parents_flat = Pair{ID,ID}[]

        # go through all branches once to collect information
        for subtree in tree.next_upstream
            # get angle and depth of sector
            α,d,points,ports,parents = serialize_tree(subtree, l, w, ω, given_ports)
            # compute depth of new sector
            ll = maybe_get(l, subtree.id)
            c = √(ll^2+d^2-2*ll*d*cos(π-α/2))
            # compute angle of new sector
            β=2acos((ll^2+c^2-d^2)/(2ll*c))

            push!(sectors, β+ maybe_get(ω, subtree.id))
            push!(depths, c)
            push!(all_points, points)
            push!(all_ports, ports)
            append!(all_parents_flat, parents)
            push!(all_parents_flat, subtree.id=>tree.id)
        end

        # calculate total angle and depth for current (sub-)tree
        sector = sum(sectors)
        depth = isempty(depths) ? 0.0 : maximum(depths)

        # calculate angles for all branches
        branch_angles = cumsum(sectors) .- sectors./2 .- sector/2

        for (branch, branch_angle, points, ports) in zip(tree.next_upstream, branch_angles, all_points, all_ports)
            # shift by ll and rotate by `branch_angle`
            ll = maybe_get(l, branch.id)
            c=cos(branch_angle)
            s=sin(branch_angle)
            transform = ((x,y),) -> Point2f0(c*x-s*(y+ll),s*x+c*(y+ll))

            branch_end = transform(Point2f0(0.0,0.0))
            # add branch for branch
            push!(all_points_flat, branch.id=>branch_end)

            # add ports for branch
            branch_ports = get(given_ports, branch.id, ID[])
            append!(all_ports_flat, Pair.(branch_ports, LinRange(Point2f0(0,0), branch_end, length(branch_ports)+2)[2:end-1]))

            # keep branches from branches' subtree
            if !isempty(points)
                append!(all_points_flat, map(((name,pts),) -> name=>transform(pts), points))
            end

            # keep ports from branches' subtree
            if !isempty(ports)
                append!(all_ports_flat, map(((prt,pos),) -> prt=>transform(pos), ports))
            end
        end

        if isa(tree, Neuron)
            neuron_ports = get(given_ports, tree.id, ID[])

            ww = maybe_get(w, tree.id)
            append!(all_ports_flat, Pair.(neuron_ports, LinRange(Point2f0(0,-ww), Point2f0(0,ww), length(neuron_ports)+2)[2:end-1]))
        end

        return (sector=sector, depth=depth, nodes=all_points_flat, ports=all_ports_flat, parents=all_parents_flat)
    end

    # get all nodes, their ports and parents in flat format
    serialized = lift(treeplot[:branch_length],treeplot[:branch_width],treeplot[:angle_between], treeplot[:ports]) do l,w,ω,prts
        ser=serialize_tree(tree,l,w,ω,prts)
        push!(ser.nodes, tree.id => Point2f0(0.0,0.0))
        ser
    end
    parent = DefaultDict(tree.id, serialized[].parents...)
    # get all the positions of the ports
    ports = lift(serialized,offset) do ser,off
        [k=>v+off for (k,v) in ser.ports]
    end

    # plot branches
    for (name,parent_name) in serialized[].parents
        c  = maybe_get(treeplot[:color], name)
        w1 = maybe_get(treeplot[:branch_width], name)
        w2 = maybe_get(treeplot[:branch_width], parent_name)

        # dynamically recompute polygon
        branch_poly = lift(w1, w2, serialized, offset) do w1,w2,ser,off
            node_dict = Dict(ser.nodes...)
            b1 = node_dict[name]
            b2 = node_dict[parent_name]
            normal = [0 1; -1 0]*(b2-b1)
            normal /= sqrt(normal'*normal)*2
            Ref(off) .+ Point2f0[
                -normal*w1+b1,
                 normal*w1+b1,
                 normal*w2+b2,
                -normal*w2+b2
            ]
        end

        # draw polygon for the branch
        poly!(treeplot, branch_poly, color=c, strokewidth = 1, strokecolor=c)

        # draw circles to cap off branch
        poly!(treeplot, lift(w1,serialized, offset) do w1,ser,off
            node_dict = Dict(ser.nodes...)
            b2 = node_dict[name]
            decompose(Point2f0, Circle(off + b2, w1/2))
        end, color=c, strokewidth = 1, strokecolor=c)

        poly!(treeplot, lift(w1,offset) do w1,off
            decompose(Point2f0, Circle(off, w1/2))
        end, color=c, strokewidth = 1, strokecolor=c)
    end

    # draw polygon for the root
    c = maybe_get(treeplot[:color], tree.id)
    w = maybe_get(treeplot[:branch_width], tree.id)
    root_poly = lift(w,offset) do w,offset
        Point2f0[(-2w/2*√(3), -2w/2), (0,2w), (2w/2*√(3), -2w/2)].+offset
    end
    poly!(treeplot, root_poly, color=c, strokewidth = 1, strokecolor=c)

    treeplot.attributes[:ports] = ports
    treeplot
end


text:
using Unitful

"""
    Spectrum <: AbstractSpectrum

A 1-dimensional spectrum stored as vectors of real numbers. The wavelengths are assumed to be in angstrom.
"""
mutable struct Spectrum{W <: Number,F <: Number} <: AbstractSpectrum{W,F}
    wave::Vector{W}
    flux::Vector{F}
    meta::Dict{Symbol,Any}
end

Spectrum(wave, flux, meta::Dict{Symbol,Any}) = Spectrum(collect(wave), collect(flux), meta)

function Base.show(io::IO, spec::Spectrum)
    print(io, "Spectrum($(eltype(spec.wave)), $(eltype(spec.flux)))")
    for (key, val) in spec.meta
        print(io, "\n  $key: $val")
    end
end


text:
using Unitful, UnitfulAstro
using PhysicalConstants.CODATA2018: h, c_0, k_B

export blackbody

"""
    blackbody(wave::Vector{<:Quantity}, T::Quantity)
    blackbody(wave::Vector{<:Real}, T::Real)

Create a blackbody spectrum using Planck's law. The curve follows the mathematical form

``B_\\lambda(T) = \\frac{2hc^2}{\\lambda^5}\\frac{1}{e^{hc/\\lambda k_B T} - 1}``

If `wave` and `T` are not `Unitful.Quantity`, they are assumed to be in angstrom and Kelvin, and the returned flux will be in units `W m^-2 Å^-1`.

The physical constants are calculated using [PhysicalConstants.jl](https://github.com/juliaphysics/physicalconstants.jl), specifically the CODATA2018 measurement set.

# References
[Planck's Law](https://en.wikipedia.org/wiki/Planck%27s_law)

# Examples
```jldoctest
julia> using Spectra, Unitful, UnitfulAstro

julia> wave = range(1, 3, length=100)u"μm"
(1.0:0.020202020202020204:3.0) μm

julia> bb = blackbody(wave, 2000u"K")
UnitfulSpectrum (100,)
  λ (μm) f (W μm^-1 m^-2)
  T: 2000 K
  name: Blackbody

julia> blackbody(ustrip.(u"angstrom", wave), 6000)
Spectrum (100,)
  T: 6000
  name: Blackbody

julia> bb.wave[argmax(bb)]
1.4444444444444444 μm

julia> 2898u"μm*K" / bb.T # See if it matches up with Wien's law
1.449 μm
```
"""
function blackbody(wave::AbstractVector{<:Quantity}, T::Quantity)
    out_unit = u"W/m^2" / unit(eltype(wave))
    flux = _blackbody(wave, T) .|> out_unit
    return spectrum(wave, flux, name = "Blackbody", T = T)
end

function blackbody(wave::AbstractVector{<:Real}, T::Real)
    flux = ustrip.(u"W/m^2/angstrom", _blackbody(wave * u"angstrom", T * u"K"))
    return spectrum(wave, flux, name = "Blackbody", T = T)
end

_blackbody(wave::AbstractVector{<:Quantity}, T::Quantity) = blackbody(T).(wave)

"""
  blackbody(T::Quantity)

Returns a function for calculating blackbody curves.
"""
blackbody(T::Quantity) = w->2h * c_0^2 / w^5 / (exp(h * c_0 / (w * k_B * T)) - 1)


text:
module user

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
using Random, Simplices, Test
using Simplices: Vertex, Simplex, windingnumber, windingangle
using Simplices: centre, assessconvergence, position, value
using Simplices: bestvertex, issortedbyangle, hypervolume

@testset "Simplices tests" begin

  Random.seed!(0)

  @testset "Vertices" begin
    v1 = Vertex([0.0, 1.0], 2.0 - 3.0im)
    v2 = Vertex([0.0, 1.0], 2.0 - 3.0im)
    @test !(v1 == v2)
    @test isequal(v1, v2)
  end

  @testset "Simplex winding number" begin

    irrelevant = [0.0, 0.0]
    @testset "Simplex encloses zero" begin
      v1 = Vertex(irrelevant, 1.0 - im)
      v2 = Vertex(irrelevant, 0.0 + im)
      v3 = Vertex(irrelevant, -1.0 - im)
      encloseszero = Simplex([v1, v2, v3])
      @test isapprox(1, abs(windingangle(encloseszero)) / (2π))
      @test windingnumber(encloseszero) == 1
    end

    @testset "Simplex doesn't enclose zero" begin
      v1 = Vertex(irrelevant, 1.0 + im)
      v2 = Vertex(irrelevant, 2.0 + im)
      v3 = Vertex(irrelevant, 1.0 + im * 2)
      doesntenclosezero = Simplex([v1, v2, v3])
      @test abs(windingangle(doesntenclosezero)) / (2π) < 1.0e-3
      @test windingnumber(doesntenclosezero) == 0
      v1 = Vertex(irrelevant, 1.0 + im)
      v2 = Vertex(irrelevant, 0.0 + 2*im)
      v3 = Vertex(irrelevant, -1.0 + im)
      doesntenclosezero = Simplex([v1, v2, v3])
      @test abs(windingangle(doesntenclosezero)) / (2π) < 1.0e-3
      @test windingnumber(doesntenclosezero) == 0
      v1 = Vertex(irrelevant, 1.0 + 0im)
      v2 = Vertex(irrelevant, 0.0 + im)
      v3 = Vertex(irrelevant, 1.0 + 0im)
      doesntenclosezero = Simplex([v1, v2, v3])
      @test abs(windingangle(doesntenclosezero)) / (2π) < 1.0e-3
      @test windingnumber(doesntenclosezero) == 0
    end
  end

  @testset "Simplex convergence" begin
    @testset "Simplices with identical vertices is converged" begin
      dim = 2
      T = Float64
      U = ComplexF64
      pos = T[1.0, 1.0]
      val = rand(U)
      v1 = Vertex(pos, val)
      v2 = Vertex(pos, val)
      v3 = Vertex(pos, val)
      s = Simplex([v1, v2, v3])

      defaults = Simplices.convergenceconfig(dim, T)
      returncode = assessconvergence(s, defaults)
      @test returncode == :XTOL_REACHED
    end

    @testset "Simplices in a simplex eps apart are converged" begin
      dim = 2
      T = Float64
      U = ComplexF64
      v1 = Vertex([one(T), one(T)], rand(U))
      v2 = Vertex([one(T), one(T) + eps(T)], rand(U))
      v3 = Vertex([one(T) + eps(T), one(T)], rand(U))
      s = Simplex([v1, v2, v3])
      defaults = Simplices.convergenceconfig(dim, T)
      returncode = assessconvergence(s, defaults)
      @test returncode == :XTOL_REACHED
    end

    @testset "Simplices in a chain eps apart are converged" begin
      dim = 2
      T = Float64
      U = ComplexF64
      v1 = Vertex([one(T), one(T)], rand(U))
      v2 = Vertex([one(T), one(T) + eps(T)], rand(U))
      v3 = Vertex([one(T), one(T) + eps(T) + eps(T)], rand(U))
      s = Simplex([v1, v2, v3])
      defaults = Simplices.convergenceconfig(dim, T)
      returncode = assessconvergence(s, defaults)
      @test returncode == :XTOL_REACHED
      @assert all(isapprox(position(v1)[d], position(v2)[d],
                        atol=defaults[:xtol_abs][d],
                        rtol=defaults[:xtol_rel][d]) for d ∈ 1:dim)
      @assert all(isapprox(position(v2)[d], position(v3)[d],
                        atol=defaults[:xtol_abs][d],
                        rtol=defaults[:xtol_rel][d]) for d ∈ 1:dim)
      @assert !all(isapprox(position(v1)[d], position(v3)[d],
                        atol=defaults[:xtol_abs][d],
                        rtol=defaults[:xtol_rel][d]) for d ∈ 1:dim)
      @test returncode == :XTOL_REACHED
    end

    @testset "Simplices vertices are in order" begin
      for _ ∈ 1:10
        dim = 2
        T = Float64
        U = ComplexF64
        v1 = Vertex(rand(dim), rand(U))
        v2 = Vertex(rand(dim), rand(U))
        v3 = Vertex(rand(dim), rand(U))
        s = Simplex([v1, v2, v3])
        p = s.permabs
        @test abs(value(s[p[1]])) < abs(value(s[p[2]]))
        @test abs(value(s[p[2]])) < abs(value(s[p[3]]))
        @test issortedbyangle(s)
      end
    end
  end

  @testset "Shapes" begin
    @testset "centres" begin
      v1 = Vertex([0.0, 0.0], one(ComplexF64))
      v2 = Vertex([1.0, 0.0], one(ComplexF64))
      v3 = Vertex([0.0, 1.0], one(ComplexF64))
      s = Simplex([v1, v2, v3])
      @test all(centre(s) .== [1/3, 1/3])
    end

    @testset "extremas" begin
      s = Simplex(x->im, [1.0, 3.0], 1.0)
      exs = Simplices.extrema(s)
      @test exs[1] == (1.0, 2.0)
      @test exs[2] == (3.0, 4.0)
    end

    @testset "hypervolumes" begin
      x0 = rand(2)
      a, b = rand(2)
      v1 = Vertex(x0 .+ [0.0, 0.0], one(ComplexF64))
      v2 = Vertex(x0 .+ [a, 0.0], one(ComplexF64))
      v3 = Vertex(x0 .+ [0.0, b], one(ComplexF64))
      s = Simplex([v1, v2, v3])
      @test hypervolume(s) ≈ (a * b) / 2
      x0 = rand(3)
      a, b, c = rand(3)
      v1 = Vertex(x0 .+ [0, 0, 0], one(ComplexF64))
      v2 = Vertex(x0 .+ [a, 0, 0], one(ComplexF64))
      v3 = Vertex(x0 .+ [0, b, 0], one(ComplexF64))
      v4 = Vertex(x0 .+ [0, 0, c], one(ComplexF64))
      s = Simplex([v1, v2, v3, v4])
      @test hypervolume(s) ≈ (a * b * c) / 6
    end
  end

  @testset "Partition of hypercubes" begin
    @testset "unity" begin
      for dims in 1:6
        positions = Simplices.partitionunitypositions(dims, Float64)
        ss = [Simplex(x->one(ComplexF64), [p for p in ps]) for ps in positions]
        for i ∈ 1:10
          x = rand(dims)
          xins = [x ∈ s for s in ss]
          @test sum(xins) == 1 # can only exist in 1 simplex
        end
      end
    end
  @testset "random hypercube" begin
      for dims in 1:6
        a = rand(dims) .- 0.5
        b = a .+ rand(dims)
        ss = Simplices.partitionhypercube(x->1, a, b)
        for i ∈ 1:10
          x = a .+ rand(dims) .* (b .- a)
          xins = [x ∈ s for s in ss]
          @test sum(xins) == 1 # can only exist in 1 simplex
        end
      end
    end
   @testset "partition simplices don't overlap" begin
      for dims in 1:4 # anything higher is too expensive
        ss = Simplices.partitionhypercube(x->1, zeros(dims), ones(dims))
        for a in ss, b in ss
          if isequal(a, b)
            @test a ∈ b
          else
            @test !(a ∈ b)
          end
        end
      end
    end
  end

end