text:
# import Pkg;

# Pkg.add("JLD")
# Pkg.add("PlotlyJS")
# Pkg.add("JSON")
# Pkg.add("JuMP")
# Pkg.add("YAML")
# Pkg.add("CSV")
# Pkg.add("DataFrames")
# Pkg.add("Query")
# Pkg.add("OrderedCollections")
# Pkg.add("Missings")

import YAML

# using JLD
using OrderedCollections
using CSV
using DataFrames
using Dates
using Query
using Statistics
using JuMP
using Missings

yamlFile = YAML.load_file("../mintsDefinitionsV2.yaml"; dicttype=OrderedDict{String,Any})      # loading yaml file and arrange in order
dirctory = yamlFile["dataFolder"]                            # taking the directory of data stored
sensorNames = yamlFile["liveStack"]                         # all the sensor names
#sensorNames = ["PPD42NSDuo"]

# only 11th node data is saved
nodeIDs_array = yamlFile["nodeIDs"]                         # node's IDs array
node11 = nodeIDs_array[11]                                  # Taking 11th node
nodeId = node11["nodeID"]                                   # taking 11th nodeID

path = dirctory*"/"*nodeId

############# anaylzing and averaging data
function ss(dd, mm, yy, sensr)
    CSV_fileName_arry = readdir(path*"/"*yy*"/"*mm*"/"*dd*"/")       # taking all the file names from the folder
    num_files = length(CSV_fileName_arry)                            # number of files

    CSV_file = "MINTS_"*nodeId*"_"*sensr*"_"*yy*"_"*mm*"_"*dd*".csv"            # CSV file name
    if (CSV_file in CSV_fileName_arry)                                          # check CSV fime exist or not
        df = CSV.read(path*"/"*yy*"/"*mm*"/"*dd*"/"*CSV_file, DataFrame)        # If the CSV file exist then take the data into a dataframe
        #df.dateTime =  SubString.(string.(df.dateTime), 1, 19)                  # Removing last three decimal numbers in dateTime column

        if sensr == "GPSGPGGA2"
            df = select!(df, Not(:timestamp))
            df = select!(df, Not(:latitudeDirection))
            df = select!(df, Not(:longitudeDirection))
            df = select!(df, Not(:altitudeUnits))
            df = select!(df, Not(:undulationUnits))
            df = select!(df, Not(:age))
            df = select!(df, Not(:stationID))
        end

        if sensr == "GPSGPRMC2"
            df = select!(df, Not(:timestamp))
            df = select!(df, Not(:status))
            df = select!(df, Not(:latitudeDirection))
            df = select!(df, Not(:longitudeDirection))
            df = select!(df, Not(:dateStamp))
            df = select!(df, Not(:magVariation))
            df = select!(df, Not(:magVariationDirection))
            df = select!(df, Not(:trueCourse))
        end

        df1 = df
        df2 = df[1:2,:]

        try
            #df1 = df
            df1.dateTime =  SubString.(string.(df1.dateTime), 1, 19)
            df1.dateTime = DateTime.(df1.dateTime,"yyyy-mm-dd HH:MM:SS")

            df1.dateTime = map((x) -> round(x, Dates.Second(30)), df1.dateTime)       # Rounding dateTime for 30 seconds
            gdf1 = groupby(df1, :dateTime)                                            # making groups by same dateTime
            cgdf1 = combine(gdf1, valuecols(gdf1) .=> mean)                            # Calculate the mean of each group and then combine the groups

            colFullName = names(cgdf1)                                               # Take the column names of dataFrame into a array
            for x in 2:length(colFullName)                                          # Starting from 2nd column and go through all the columns
                num_let = length(colFullName[x])                                    # taking number of characters in one column name
                colName = SubString.(string.(colFullName[x]), 1, num_let-5)         # removing last five characters of each column name (that is _mean)
                rename!(cgdf1,colFullName[x] => sensr*"_"*colName)                   # adding sensor name to the column name
            end

            return cgdf1                             # return the dataFrame to for loop

        catch e
            println("Error on "*CSV_file)
            #df2 = df[1:10,:]

            if sensr == "GPSGPGGA2"
                clName = names(df2)
                if "Column18" in clName
                    df2 = select!(df2, Not(:Column18))
                end
                if "Column19" in clName
                    df2 = select!(df2, Not(:Column19))
                end

                try
                    df2.latitudeCoordinate = parse.(Float64, df2.latitudeCoordinate)
                    df2.longitudeCoordinate = parse.(Float64, df2.longitudeCoordinate)
                    df2.numberOfSatellites = parse.(Float64, df2.numberOfSatellites)
                catch e
                end
            end

            if sensr == "GPSGPRMC2"
                clName = names(df2)
                if "Column15" in clName
                    df2 = select!(df2, Not(:Column15))
                end
                if "Column16" in clName
                    df2 = select!(df2, Not(:Column16))
                end
                if "Column17" in clName
                    df2 = select!(df2, Not(:Column17))
                end
                if "Column18" in clName
                    df2 = select!(df2, Not(:Column18))
                end
                if "Column19" in clName
                    df2 = select!(df2, Not(:Column19))
                end
                if "Column20" in clName
                    df2 = select!(df2, Not(:Column20))
                end
                if "Column21" in clName
                    df2 = select!(df2, Not(:Column21))
                end

                try
                    #df2.latitudeCoordinate = parse.(Float64, df2.latitudeCoordinate)
                    #df2.longitudeCoordinate = parse.(Float64, df2.longitudeCoordinate)
                    df2.speedOverGround = parse.(Float64, df2.speedOverGround)
                    df2.longitude = parse.(Float64, df2.longitude)
                catch e
                end
            end

            if sensr == "PPD42NSDuo"
                clName = names(df2)
                if "Column10" in clName
                    df2 = select!(df2, Not(:Column10))
                end
                if "Column11" in clName
                    df2 = select!(df2, Not(:Column11))
                end

                try
                    df2.LPOPmMid = parse.(Float64, df2.LPOPmMid)
                catch e
                end
            end

            println(df2)
            df2.dateTime =  SubString.(string.(df2.dateTime), 1, 19)
            df2.dateTime = DateTime.(df2.dateTime,"yyyy-mm-dd HH:MM:SS")

            df2.dateTime = map((x) -> round(x, Dates.Second(30)), df2.dateTime)       # Rounding dateTime for 30 seconds
            gdf2 = groupby(df2, :dateTime)                                            # making groups by same dateTime
            cgdf2 = combine(gdf2, valuecols(gdf2) .=> mean)                            # Calculate the mean of each group and then combine the groups

            colFullName = names(cgdf2)                                               # Take the column names of dataFrame into a array
            for x in 2:length(colFullName)                                          # Starting from 2nd column and go through all the columns
                num_let = length(colFullName[x])                                    # taking number of characters in one column name
                colName = SubString.(string.(colFullName[x]), 1, num_let-5)         # removing last five characters of each column name (that is _mean)
                rename!(cgdf2,colFullName[x] => sensr*"_"*colName)                   # adding sensor name to the column name
            end

            return cgdf2
        end
    end
end

#global sensr_exsist = false
n = 0

for sensr in sensorNames                                                # Go through all the sensor names
    k = 0
    years = readdir(path)                                               # taking all the year folder names

    for y in years                                                      #
        months = readdir(path*"/"*y)

        for m in months
            days = readdir(path*"/"*y*"/"*m)

            for d in days
                filess = readdir(path*"/"*y*"/"*m*"/"*d)
                CSV_file = "MINTS_"*nodeId*"_"*sensr*"_"*y*"_"*m*"_"*d*".csv"

                if (CSV_file in filess)
                    println("Date: "*y*"_"*m*"_"*d)
                    println(sensr)
                    averaged_df = ss(d, m, y, sensr)
                    df2 = averaged_df
                    k = k + 1

                    println("done")
                    global sensr_exsist = true

                    if (k > 1)

                        #try
                        global new_df = outerjoin(new_df, df2, on = intersect(names(new_df), names(df2)))
                        #catch e
                            #println("Error of adding: "*"MINTS_"*nodeId*"_"*sensr*"_"*y*"_"*m*"_"*d*".csv")
                        #end

                    else
                        new_df = averaged_df
                    end

                else
                    println("No data of "*sensr)
                    sensr_exsist = false
                    #continue
                end
            end
        end
    end

    global n = n + 1

    if sensr_exsist

        if (n > 1)
            println("******************************************************")
            global final_df = outerjoin(final_df, new_df, on = :dateTime, makeunique=true)
        else
            final_df = new_df
        end

    end

end

final_df2 = sort!(final_df)
CSV.write("/mnt/3fbf694d-d8e0-46c0-903d-69d994241206/mintsData/avJlData/"*nodeId*"_df.csv", final_df2)


text:
"""
    abstract type AbstractCodeGenAction <: AbstractFrontendAction
Supertype for `CodeGenAction`s.
"""
abstract type AbstractCodeGenAction <: AbstractFrontendAction end

Base.unsafe_convert(::Type{CXCodeGenAction}, x::T) where {T<:AbstractCodeGenAction} = x.ptr
Base.cconvert(::Type{CXCodeGenAction}, x::T) where {T<:AbstractCodeGenAction} = x

struct LLVMOnlyAction <: AbstractCodeGenAction
    ptr::CXCodeGenAction
end

ext: jl


text:
"""
    add_current_update_v_pp!( j_dofs, p, position_old, position_new,
                              marker_charge)

Add current for one particle and update v
(according to `H_p1` part in Hamiltonian splitting)
"""
function add_current_update_v_pp!(
    j_dofs::AbstractArray,
    p::ParticleMeshCoupling1D,
    position_old::Float64,
    position_new::Float64,
    marker_charge::Float64,
)

    # Read out particle position and velocity
    # Compute index_old, the index of the last DoF on the grid the particle
    # contributes to, and r_old, its position (normalized to cell size one).

    xi = (position_old[1] - p.xmin) / p.delta_x[1]
    index_old = trunc(Int, xi)
    r_old = xi - index_old

    # Compute the new box index index_new and normalized position r_old.
    xi = (position_new[1] - p.xmin) / p.delta_x[1]
    index_new = trunc(Int, xi)
    r_new = xi - index_new

    if index_old == index_new
        update_jv_pp!(j_dofs, p, r_old, r_new, index_old, marker_charge)

    elseif index_old < index_new
        update_jv_pp!(j_dofs, p, r_old, 1.0, index_old, marker_charge)

        update_jv_pp!(j_dofs, p, 0.0, r_new, index_new, marker_charge)

        for ind in (index_old + 1):(index_new - 1)
            update_jv_pp!(j_dofs, p, 0.0, 1.0, ind, marker_charge)
        end

    else
        update_jv_pp!(j_dofs, p, 1.0, r_new, index_new, marker_charge)
        update_jv_pp!(j_dofs, p, r_old, 0.0, index_old, marker_charge)

        for ind in (index_new + 1):(index_old - 1)
            update_jv_pp!(j_dofs, p, 1.0, 0.0, ind, marker_charge)
        end
    end
end

"""
    update_jv_pp!( j_dofs, p, lower, upper, index, marker_charge)

Helper function for `add_current_update_v`.
"""
function update_jv_pp!(
    j_dofs::AbstractArray,
    p::ParticleMeshCoupling1D,
    lower::Float64,
    upper::Float64,
    index::Int64,
    marker_charge::Float64,
)
    n_cells = p.n_grid[1]

    # Evaluation of the primitive integral at the lower
    # and upper bound of the gridcell

    horner_primitive_1d(p.spline_val, p.spline_degree, p.spline_pp.poly_coeffs_fp, lower)
    horner_primitive_1d(
        p.spline_val_more, p.spline_degree, p.spline_pp.poly_coeffs_fp, upper
    )

    p.spline_val .= (p.spline_val_more .- p.spline_val) .* (p.delta_x[1])

    ind = 1
    @inbounds for i_grid in (index - p.spline_degree):index
        i_mod = mod(i_grid, n_cells) + 1
        j_dofs[i_mod] = j_dofs[i_mod] + (marker_charge * p.spline_val[ind] * p.scaling)
        ind = ind + 1
    end
end

"""
    add_current_update_v!( j_dofs, p, position_old, position_new, marker_charge)

Add current for one particle and update v (according to ``H_{p1}``
part in Hamiltonian splitting)

- Read out particle position and velocity
- Compute index_old, the index of the last DoF on the grid the
particle contributes to, and `r_old`, its position (normalized to cell size one).

"""
function add_current_update_v!(
    j_dofs::AbstractArray,
    p::ParticleMeshCoupling1D,
    position_old::Float64,
    position_new::Float64,
    marker_charge::Float64,
)
    xi = (position_old - p.xmin) / p.delta_x[1]
    index_old = trunc(Int, xi)
    r_old = xi - index_old

    # Compute the new box index index_new and normalized position r_old.

    xi = (position_new - p.xmin) / p.delta_x[1]
    index_new = trunc(Int, xi)
    r_new = xi - index_new

    if index_old == index_new
        if r_old < r_new
            update_jv!(j_dofs, p, r_old, r_new, index_old, marker_charge, 1.0)
        else
            update_jv!(j_dofs, p, r_new, r_old, index_old, marker_charge, -1.0)
        end

    elseif index_old < index_new
        update_jv!(j_dofs, p, r_old, 1.0, index_old, marker_charge, 1.0)

        update_jv!(j_dofs, p, 0.0, r_new, index_new, marker_charge, 1.0)

        for ind in (index_old + 1):(index_new - 1)
            update_jv!(j_dofs, p, 0.0, 1.0, ind, marker_charge, 1.0)
        end

    else
        update_jv!(j_dofs, p, r_new, 1.0, index_new, marker_charge, -1.0)
        update_jv!(j_dofs, p, 0.0, r_old, index_old, marker_charge, -1.0)

        for ind in (index_new + 1):(index_old - 1)
            update_jv!(j_dofs, p, 0.0, 1.0, ind, marker_charge, -1.0)
        end
    end
end

"""
    update_jv!(j_dofs, p, lower, upper, index, marker_charge, qoverm, sign)

Helper function for `add_current_update_v`.
"""
function update_jv!(
    j_dofs::AbstractArray,
    p::ParticleMeshCoupling1D,
    lower::Float64,
    upper::Float64,
    index::Int64,
    marker_charge::Float64,
    sign::Float64,
)
    n_cells = p.n_grid[1]

    c1 = 0.5 * (upper - lower)
    c2 = 0.5 * (upper + lower)

    uniform_bsplines_eval_basis!(p.spline_val, p.spline_degree, c1 * p.quad_x[1] + c2)

    p.spline_val .*= p.quad_w[1] * c1

    @inbounds for j in 2:(p.n_quad_points)
        uniform_bsplines_eval_basis!(
            p.spline_val_more, p.spline_degree, c1 * p.quad_x[j] + c2
        )

        p.spline_val .+= p.spline_val_more .* p.quad_w[j] .* c1
    end

    p.spline_val .*= sign * p.delta_x[1]

    ind = 1
    for i_grid in (index - p.spline_degree):index
        i_mod = mod(i_grid, n_cells) + 1
        j_dofs[i_mod] += marker_charge * p.spline_val[ind] * p.scaling
        ind = ind + 1
    end
end


text:
immutable UVWVertex
    position::Vec3f0
    uvw::Vec3f0
end

immutable Uniforms{SV, BF, CF, AF}
    model::Mat4f0
    modelinv::Mat4f0
    isovalue::Float32
    startvalue::SV
    breakloop::BF
    to_color::CF
    accumulation::AF
end

function vert_volume(vertex, canvas, uniforms, intensities)
    world_vert = uniforms.model * Visualize.to_vec4(vertex.position)
    pos = canvas.projectionview * world_vert
    return (
        pos,
        UVWVertex(world_vert[Vec(1, 2, 3)], vertex.uvw)
    )
end

function is_outside(position::Vec3f0)
    position[1] > 1f0 || position[2] > 1f0 ||
    position[3] > 1f0 || position[1] < 0f0 ||
    position[2] < 0f0 || position[3] < 0f0
end

function mip(accum, pos, stepdir, intensities, uniforms)
    max(accum, intensities[pos][1])
end

function raycast_loop(
        front::Vec3f0, dir::Vec3f0, num_samples, stepsize::Float32,
        uniforms, intensities
    )
    stepsize_dir = dir * stepsize
    pos = front
    pos += stepsize_dir; # apply first, to padd
    accumulator = uniforms.startvalue
    for i = 1:num_samples
        if uniforms.breakloop(accumulator) || is_outside(pos)
            break
        end
        accumulator = uniforms.accumulation(
            accumulator, pos, stepsize_dir, intensities, uniforms
        )
        pos += stepsize_dir
    end
    return uniforms.to_color(accumulator, uniforms)
end

function frag_volume(vertex_out, canvas, uniforms, intensities)
    max_distance = 1.73f0
    num_samples = 128
    step_size = max_distance / Float32(num_samples)
    dir = normalize(vertex_out.position - canvas.eyeposition)
    dir = (uniforms.modelinv * Vec4f0(dir[1], dir[2], dir[3], 0f0))[Vec(1, 2, 3)]
    color = raycast_loop(vertex_out.uvw, dir, num_samples, step_size, uniforms, intensities)
    (color, )
end


text:
module NucleosomeInference

# Write your package code here.
using Turing
using Random
using DataFrames
using StatsFuns
using Plots, StatsPlots, Measures

export randclasses, simpletwoclass, plotconfig, simplemulticlass, sample_mh_jump

include("sample.jl")
include("plotting.jl")
include("model.jl")
include("sampling.jl")

end


text:
module JuFormatting

using Base.Ryu

# Write your package code here.

export
    format,
    @f_str

include("main.jl")
include("numberformatparser.jl")
include("floatingpointformating.jl")
include("integerformating.jl")
include("fmacro.jl")

end


text:
branin(x) = branin(x[1], x[2])
branin(x1, x2; a = 1, b = 5.1/(4π^2), c = 5/π, r = 6, s = 10, t = 1/(8π),
       noiselevel = 0.) =
    a * (x2 - b*x1^2 + c*x1 - r)^2 + s*(1 - t)*cos(x1) + s + noiselevel * randn()

minima(::typeof(branin)) = [([-π, 12.275], 0.397887), ([π, 2.275], 0.397887), ([9.42478, 2.475], 0.397887)]

euclidean(x, y) = sqrt(sum((x .- y).^2))
function regret(optimizer, func)
    mins = [m[1] for m in minima(func)]
    (observed_regret = minimum(map(m -> euclidean(m, optimizer), mins)))
end

config = ConfigParameters()
config.verbose_level = 0

lb = [-5., 0.]; ub = [10., 15.]
optimizer, optimum = bayes_optimization(branin, lb, ub, config)
@test regret(optimizer, branin) < .05


text:
### A Pluto.jl notebook ###
# v0.11.12

using Markdown
using InteractiveUtils

# ╔═╡ 96bc8c74-f10a-11ea-2f86-9f75ccaa5f89
using Primes

# ╔═╡ 986d2a4a-f10c-11ea-0608-c3567591b9f4
using Statistics

# ╔═╡ bfa53152-f10c-11ea-3e60-c70874f29c08
using Gadfly  # the Julia version of the ggplot2 graphics package

# ╔═╡ e166f74e-f10c-11ea-01ec-8332273a87a1
using Colors  # for color keys and color maps in graphics

# ╔═╡ eb919846-f10c-11ea-1eaa-c59b08e5ea3e
using Cairo   # for graphics conversions, eg SVG --> PNG

# ╔═╡ fa380c72-f10c-11ea-3589-f3f3b35f6ced
using DataFrames # 2D arrays with named columns; hard to use Gadfly without them

# ╔═╡ b315f4d6-f10e-11ea-0fdc-9390980b46f1
md"""## Prime After Prime (2020 version for Pluto.jl)

This notebook explores correlations among the residues modulo $m$ of consecutive prime numbers. The notebook is an accompaniment to the article "[Prime After Prime](http://bit-player.org/2016/prime-after-prime)" published at bit-player.org on May 31, 2016. These correlations are *not* my discovery; they were first fully described and analyzed by Robert J. Lemke Oliver and Kannan Soundararajan in March 2016. For the mathematical background, see their arXiv preprint: [Unexpected Biases in the Distribution of Consecutive Primes](http://arxiv.org/abs/1603.03720).

The code is written in the [Julia programming language](http://julialang.org/). The original version was written in Julia version 0.4, which is now very much out of date. In August of 2020, inspired by an inquiry from a high school student, I set about updating the code for use with Julia 1.0 or later. This was easier than I expected, partly because I discovered I was able to throw away much of the code. Functions I had written are available in Julia packages. Also, I have only included the essentials for exploring the mathematics and plotting results. I have eliminated various false starts and dead ends that serve mainly to distract.

The 2016 version of the notebook, primeafterprime.ipynb, remains available in this repository.

All my development work was done in a Jupyter notebook, but this version has been imported into the new notebook format of [Pluto.jl](https://github.com/fonsp/Pluto.jl).

"""

# ╔═╡ b4672322-f111-11ea-149d-8be5ae03d1d9
md"""## Basic operating instructions

- Install Pluto.js via the Julia Pkg manager
- cd to the directory where this file is stored
- Do ```import Pluto``` and ```Pluto.run()```
- In your default browser, navigate to ```localhost:1234```
- Open this file
- Evaluate all the cells. (This should happen automatically, I think, but I'm just learning Pluto.js. It could take a few minutes.)

To test that the system is working, go to the last three cells in the file, uncomment the lines there, and evaluate them. The first run may take a while, but you should eventually see a pink-and-blue heatmap like the ones in the bit-player article.

For your own experiments, the protocol is as follows:

Generate a list of consecutive primes. For example, get the first 10,000 primes larger than 500,000:

```Julia
p50000 = generate_primes(500000, 10000)
```

Now examine those primes modulo m and analyze correlations between successive residues mod m. For m = 7:

```Julia
p50000mod7 = analyze_primes_mod_m(p50000, 7)
```

The matrices p50000mod7[:nz\_class\_counts] and p50000mod7[:nz\_pair\_counts] are likely to be the main focus of interest.

To produce a graphic heatmap showing the variations in the pair correlation function, do:

```Julia
hm = heatmap(p50000mod7)
```

(Or do something else with the data that you find interesting!)
"""

# ╔═╡ c6901830-f10a-11ea-0e71-9b095034a67a
# Generate a list of 'n' primes starting with the first prime
# greater than or equal to 'start'. The result is returned as an
# array of integers of the same type as the type of 'start'.

# CAUTION: If 'start' is just below a boundary between integer types,
# this procedure may throw an OverflowError or, worse, return a
# list of primes that wraps around to those at the origin of the
# number line. The boundaries in question are:
# typemax(Int64) = 9223372036854775807
# typemax(Int128) = 170141183460469231731687303715884105727
# If you are working in that neighborhood, it's best to explicitly
# coerce 'start' to the larger type, i.e. Int128(9223372036854775807) or
# BigInt(170141183460469231731687303715884105727)

function generate_primes(start, n)
    primelist = Array{typeof(start)}(undef, n)
    p = nextprime(start)
    for i = 1:n
        primelist[i] = p
        p = nextprime(p + 1)
    end
    return primelist
end

# ╔═╡ d8b4f6a2-f10a-11ea-0263-bb54f81fc3cd
# Auxiliary to 'analyze_primes_mod_m'

# Given a modulus m, return a vector of all possible
# congruence classes for p mod m, where p is a prime > m.
# When m is prime, the vector will consist of the range
# 1:m-1, but for composite m there may be missing values.
# For example, congruence_classes(10) = [1, 3, 7, 9]

function collect_congruence_classes(m)
    classes = Array{Int64}(undef, 0)
    for i = 1:m
        if gcd(i, m) == 1
            push!(classes, i)
        end
    end
    return classes
end

# ╔═╡ 7587d2c8-f10c-11ea-02fd-ad6dce0ead8a
# Auxiliary to 'analyze_primes_mod_m'

function nonzero_class_counts(class_counts, classes)
    n = length(classes)
    nz_counts = zeros(Int64, n)
    for i = 1:n
        nz_counts[i] = class_counts[classes[i] + 1]
    end
    return nz_counts
end

# ╔═╡ 82e434c0-f10c-11ea-2d53-f796bcc82a68
# Auxiliary to 'analyze_primes_mod_m'

function nonzero_pair_counts(pair_counts, classes)
    n = length(classes)
    nz_counts = zeros(Int64, n, n)
    for i = 1:n
        for j = 1:n
            nz_counts[i, j] = pair_counts[classes[i] + 1, classes[j] + 1]
        end
    end
    return nz_counts
end

# ╔═╡ 8da414dc-f10c-11ea-2511-119c66bf5293
# Auxiliary to 'analyze_primes_mod_m'

# calculate the mean value of all off-diagonal elements and the
# mean of the on-diagonal elements; return the ratio of off/on

function diagonal_ratio(count_matrix)
    n = size(count_matrix, 1)
    diagonal_sum = 0
    for i = 1:n
        diagonal_sum += count_matrix[i, i]
    end
    off_diagonal_sum = sum(count_matrix) - diagonal_sum
    return (off_diagonal_sum / (n^2 - n)) / (diagonal_sum / n)
end

# ╔═╡ b0dc1758-f10c-11ea-000f-a7b0d08cd063
# given a list of primes and an integer modulus m,
# returns a statistical analysis of the primes mod m, organized
# as a data dictionary

function analyze_primes_mod_m(primelist, m)

    # associative structure to hold all data on primes
    d = Dict()

    # basic info about the set of primes and the modulus
    d[:first_prime] = primelist[1]
    d[:last_prime] = primelist[end]
    d[:median_prime] = primelist[div(end, 2)]
    d[:prime_count] = length(primelist)
    d[:m] = m

    # set of integers that are possible residues for a prime p mod m
    d[:congruence_classes] = collect_congruence_classes(m)

    # log10 of the median will become x coordinate in certain graphs
    d[:magnitude] = Float64(log10(d[:median_prime]))
    d[:digits] = round(Int, (ceil(d[:magnitude])))

    # arrays of zeros to hold counts of residue classes and consecutive pairs
    class_counts = zeros(Int, m)
    pair_counts = zeros(Int, (m, m))

    # loop through the primes mod m, recording the number in
    # each congruence class and the number of each p_{i}, p_{i+1}
    # pair of successive primes; adjust indices to translate
    # from 0:m-1 range of modular values to 1:m range of Julia arrays.
    prev = primelist[1] % m
    class_counts[prev + 1] += 1
    for i = 2:length(primelist)
        next = primelist[i] % m
        class_counts[next + 1] += 1
        pair_counts[prev + 1, next + 1] += 1
        prev = next
    end

    # store the counts in the dictionary
    d[:class_counts] = class_counts
    d[:pair_counts] = pair_counts

    # results above include zero entries for impossible congruence
    # classes; here we separately record a vector and a matrix without
    # the zero entries
    d[:nz_class_counts] = nonzero_class_counts(class_counts, d[:congruence_classes])
    d[:nz_pair_counts] = nonzero_pair_counts(pair_counts, d[:congruence_classes])

    # normalize the results to allow comparisons of samples of different sizes
    mc = mean(d[:nz_class_counts])
    mp = mean(d[:nz_pair_counts])
    d[:norm_class_counts] = map(x -> (x - mc) / mc, d[:class_counts])
    d[:norm_pair_counts] = map(x -> (x - mp) / mp, d[:pair_counts])
    d[:norm_nz_class_counts] = map(x -> (x - mc) / mc, d[:nz_class_counts])
    d[:norm_nz_pair_counts] = map(x -> (x - mp) / mp, d[:nz_pair_counts])

    # elementary stats on the normalized results
    d[:norm_class_counts_mean] = mean(d[:norm_nz_class_counts])
    d[:norm_class_counts_std] = std(d[:norm_nz_class_counts])
    d[:norm_pair_counts_mean] = mean(d[:norm_nz_pair_counts])
    d[:norm_pair_counts_std] = std(d[:norm_nz_pair_counts])
    d[:diag_ratio] = diagonal_ratio(d[:nz_pair_counts])
    return d
end

# ╔═╡ 0218ae7e-f10d-11ea-3916-c3f4ee27d9b6
Gadfly.set_default_plot_size(400pt, 400pt)

# ╔═╡ 11310d90-f10d-11ea-35a1-4b37f1f49b81
function colorfn(x)
    if x == 0.0
        return RGB(0.5, 0.5, 0.5)
    else
        index = round(Int, min(200, 1 + 200x))
        return Colors.diverging_palette(245.0, 10.0, 200, wcolor=RGB(1,0,0), c=0.7)[index]
    end
end

# ╔═╡ 1f3ee64e-f10d-11ea-3aec-bd0e36046247
function clamp_colors(mat)
    max_value = maximum(mat)
    if max_value > 1
        adjusted_mat = similar(mat)
        for i in eachindex(mat)
            if mat[i] > 0
                adjusted_mat[i] = mat[i] / max_value
            else
                adjusted_mat[i] = mat[i]
            end
        end
        return adjusted_mat
    end
    return mat
end

# ╔═╡ 25a9d4e4-f10d-11ea-1c22-ef9ac3b8fcd4
function heatmap(d::Dict)
    mat = clamp_colors(d[:norm_nz_pair_counts])
    labels = d[:congruence_classes]
    spy(mat,
        Scale.x_discrete(labels = x -> string(labels[x])),
        Scale.y_discrete(labels = y -> string(labels[y])),
        Guide.xlabel("j = second prime mod " * string(d[:m])),
        Guide.ylabel("i = first prime mod " * string(d[:m]), orientation=:vertical),
        Guide.title("normalized counts of consecutive " * string(d[:digits]) * "-digit primes mod " * string(d[:m])),
        Guide.colorkey(""),
        Scale.color_continuous(colormap = colorfn, minvalue = -1.0, maxvalue = 1.0))
end

# ╔═╡ 2fba1106-f10d-11ea-3214-83caf8fe272c
function heatmap_z(d::Dict)
    mat = clamp_colors(d[:norm_pair_counts])
    labels = collect(0:d[:m] - 1)
    spy(mat,
        Scale.x_discrete(labels = x -> string(labels[x])),
        Scale.y_discrete(labels = y -> string(labels[y])),
        Guide.xlabel("j = second prime mod " * string(d[:m])),
        Guide.ylabel("i = first prime mod " * string(d[:m]), orientation=:vertical),
        Guide.title("normalized counts of consecutive " * string(d[:digits]) * "-digit primes mod " * string(d[:m])),
        Guide.colorkey(""),
        Scale.color_continuous(colormap = colorfn, minvalue = -1.0, maxvalue = 1.0))
end

# ╔═╡ 38afa348-f10d-11ea-1619-59dc225fa5a6
# testmod7 = generate_primes(1000000, 100000)

# ╔═╡ 49183112-f10d-11ea-19eb-c3fb0a4852d8
# test7data = analyze_primes_mod_m(testmod7, 7)

# ╔═╡ 5035f328-f10d-11ea-0029-39abf9f73636
# hm7 = heatmap(test7data)

# ╔═╡ Cell order:
# ╟─b315f4d6-f10e-11ea-0fdc-9390980b46f1
# ╟─b4672322-f111-11ea-149d-8be5ae03d1d9
# ╠═96bc8c74-f10a-11ea-2f86-9f75ccaa5f89
# ╠═c6901830-f10a-11ea-0e71-9b095034a67a
# ╠═d8b4f6a2-f10a-11ea-0263-bb54f81fc3cd
# ╠═7587d2c8-f10c-11ea-02fd-ad6dce0ead8a
# ╠═82e434c0-f10c-11ea-2d53-f796bcc82a68
# ╠═8da414dc-f10c-11ea-2511-119c66bf5293
# ╠═986d2a4a-f10c-11ea-0608-c3567591b9f4
# ╠═b0dc1758-f10c-11ea-000f-a7b0d08cd063
# ╠═bfa53152-f10c-11ea-3e60-c70874f29c08
# ╠═e166f74e-f10c-11ea-01ec-8332273a87a1
# ╠═eb919846-f10c-11ea-1eaa-c59b08e5ea3e
# ╠═fa380c72-f10c-11ea-3589-f3f3b35f6ced
# ╠═0218ae7e-f10d-11ea-3916-c3f4ee27d9b6
# ╠═11310d90-f10d-11ea-35a1-4b37f1f49b81
# ╠═1f3ee64e-f10d-11ea-3aec-bd0e36046247
# ╠═25a9d4e4-f10d-11ea-1c22-ef9ac3b8fcd4
# ╠═2fba1106-f10d-11ea-3214-83caf8fe272c
# ╠═38afa348-f10d-11ea-1619-59dc225fa5a6
# ╠═49183112-f10d-11ea-19eb-c3fb0a4852d8
# ╠═5035f328-f10d-11ea-0029-39abf9f73636

max_stars_repo_name: bit-player/primeafterprime

max_issues_repo_name: bit-player/primeafterprime

max_forks_repo_name: bit-player/primeafterprime


text:
module RocketBuilder

const datadir = joinpath(@__DIR__, "..", "data")

include("db.jl")
include("tank-editor.jl")
include("engine-editor.jl")

end


text:
using JLD
using DataFrames
using Dates
using LinearAlgebra
using DifferentialEquations
using Distributions
using Random
using CSV
using Printf
using RCall
using CPUTime
using StatsPlots
include("utils.jl")
include("sir_ode.jl")
include("abc.jl")
include("abcMCMC.jl")
include("abcSMC.jl")
include("BayesianCalibration.jl")

df_pre = DataFrame(CSV.File("../data/covid19_bc_pre_vaccine.csv"))
df_post = DataFrame(CSV.File("../data/covid19_bc_post_vaccine.csv"))

""" Pre-vaccine analysis """
tot_days = size(df_pre)[1]
time_interval = 7

y_pre = hcat([df_pre[:S],df_pre[:I],df_pre[:R]]...)[1:7:tot_days,:]
y_pre = y_pre./1e4
y_dim = size(y_pre)[1]
# test ABC
function data_generator_day(p)
    initial_state = y_pre[1,:]; time_window=(1,Float64(tot_days))
    solve_ode(initial_state,time_window,p,Float64(7))
end

# test ABC
function data_generator_week(p)
    initial_state = y_pre[1,:]; time_window=(1,Float64(tot_days)/time_interval)
    solve_ode(initial_state,time_window,p)
end

algo_parameters_abc = (prior = (Truncated(Normal(0.1,0.5),0,Inf),Truncated(Normal(0.2,0.2),0,Inf)),epsilon = 0.1,
eta= identity_mapping, d= distanceFunction)
output_abc=ABC(y_pre', data_generator_week, algo_parameters_abc, 500)

algo_parameter_mcmc = (prior = (Truncated(Normal(0.1,0.5),0,Inf),Truncated(Normal(0.2,0.2),0,Inf)),epsilon = 1, eta = identity_mapping,
                        d= compute_norm, proposal = proposal_Normal,
                        proposalRatio = proposalRatio_Normal,sd = (0.25,0.25),
                        thinning = 10000, burn_in=100, verbose=true)

output_mcmc = ABC_MCMC(y_pre', data_generator_week, algo_parameter_mcmc, 500)

algo_parameters_smc = (prior = (Truncated(Normal(0.1,0.5),0,Inf),Truncated(Normal(0.2,0.2),0,Inf)), time_final=5, eps_list = [10, 5, 1, 0.5, 0.1],
eta= identity_mapping, d= distanceFunction, kernel=proposal_Normal,
kernel_density=proposal_Normal_density,
sd=(0.5,0.5), resample_method=systematic_resample, verbose=true)
output_smc = ABC_SMC(y_pre', data_generator_week, algo_parameters_smc, 500)

abc_solutions=zeros(3,y_dim,500)
abc_mcmc_solutions=zeros(3,y_dim,500)
abc_smc_solutions=zeros(3,y_dim,500)
for i=1:500
    abc_solutions[:,:,i] = data_generator_week(output_abc[1][i,:])[1:3,:]
    abc_mcmc_solutions[:,:,i] = data_generator_week(output_mcmc[1][i,:])[1:3,:]
    abc_smc_solutions[:,:,i] = data_generator_week(output_smc[1][i,:])[1:3,:]
end

# save("../results/abc_solution_500.jld","abc_solutions",abc_solutions)
# save("../results/abc_mcmc_solution_500.jld","abc_mcmc_solutions",abc_mcmc_solutions)
# save("../results/abc_smc_solution_500.jld","abc_smc_solutions",abc_smc_solutions)

# alpha_level=1
# beta_plot = density(output_abc[1][:,1],alpha=alpha_level, label="ABC",  title="Predicted β Posterior", xlabel="β", ylabel="Density")
# density!(output_mcmc[1][:,1],alpha=alpha_level, label="ABC-MCMC")
# density!(output_smc[1][:,1],alpha=alpha_level, label="ABC-SMC")
# savefig(beta_plot,"../figs/beta_plot.png")

# gamma_plot = density(output_abc[1][:,2],alpha=alpha_level, label="ABC",  title="Predicted γ Posterior", xlabel="γ", ylabel="Density")
# density!(output_mcmc[1][:,2],alpha=alpha_level, label="ABC-MCMC")
# density!(output_smc[1][:,2],alpha=alpha_level, label="ABC-SMC")
# savefig(gamma_plot,"../figs/gamma_plot.png")

 # ess
 ess_output = zeros(3)
 ess_output_time = copy(ess_output)
 ess_output[1]  = try rcopy(R"mean(mcmcse::ess($(output_abc[1])))") catch; 0 end
 ess_output_time[1]  = ess_output[1] / output_abc[2]
 ess_output[2]  = try rcopy(R"mean(mcmcse::ess($(output_mcmc[1])))") catch; 0 end
 ess_output_time[2]  = ess_output[2] / output_mcmc[2]
 ess_output[3]  = try rcopy(R"mean(mcmcse::ess($(output_smc[1])))") catch; 0 end
 ess_output_time[3]  = ess_output[3] / output_smc[2]

 save("../results/covid/ess_pre.jld","ess_output",ess_output)
 save("../results/covid/ess_cputime_pre.jld","ess_output_time",ess_output_time)

# solution_day = data_generator_day((0.042,0.07))
# solution_week = data_generator_week((0.12,0.2))

# plot(solution_day[2:end], yaxis=:log, label = ["s_model" "i_model" "r_model"], title="Daily:β=0.042, γ=0.07");
# plot!(1:7:250, y_pre, yaxis=:log, label=["s" "i" "r"], seriestype=:scatter)

# plot(1:7:250,solution_week', yaxis=:log, label = ["s_model" "i_model" "r_model"], title="Weekly:β=0.2, γ=0.3");
# plot!(1:7:250, y_pre, yaxis=:log, label=["s" "i" "r"], seriestype=:scatter)

# plot(hcat([df_pre[:S],df_pre[:I],df_pre[:R]]...), yaxis=:log, label=["s" "i" "r"], seriestype=:scatter)


text:
a, b = parse.(Int, split(readline()))
if a < b
    println("a < b")
elseif a > b
    println("a > b")
else
    println("a == b")
end


text:
using ImageCore
using IterTools

# AbstractGray and Color3 tests should be generated seperately
function generate_test_types(number_types::AbstractArray{<:DataType}, color_types::AbstractArray{<:UnionAll})
    test_types = map(Iterators.product(number_types, color_types)) do T
        try
            T[2]{T[1]}
        catch err
            !isa(err, TypeError) && rethrow(err)
        end
    end
    test_types = filter(x->x != false, test_types)
    if isempty(filter(x->x<:Color3, test_types))
        test_types = [number_types..., test_types...]
    end
    test_types
end

_base_colorant_type(::Type{<:Number}) = Gray
_base_colorant_type(::Type{T}) where T<:Colorant = base_colorant_type(T)
"""
    test_numeric(dist, a, b, T; filename=nothing)

simply test that `dist` works for 2d image as expected, more tests go to `Distances.jl`
"""
function test_numeric(dist, a, b, T; filename=nothing)
    size(a) == size(b) || error("a and b should be the same size")
    if filename == nothing
        filename = "references/$(typeof(dist))_$(ndims(a))d"

        if eltype(a) <: Color3
            filename = filename * "_$(_base_colorant_type(T))"
        elseif eltype(a) <: Union{Number, AbstractGray}
            filename = filename * "_$(_base_colorant_type(T))"
        end
    end
    # @test_reference "$(filename)_$(eltype(a))_$(eltype(b)).txt" assess(dist, a, b)
    @test_reference "$(filename).txt" Float64(assess(dist, a, b))
end

"""
    test_cross_type(dist, a, b, type_list)

simply test if operations between `N0f8`, `Bool` and `Float32` types works as expected.
`a` and `b` should be simple enough to get rid of InexactError.
"""
function test_cross_type(dist, a, b, type_list)
    size(a) == size(b) || error("a and b should be the same size")
    rsts = [[assess(dist, Ta.(a), Tb.(b)),
             assess(dist, Tb.(a), Ta.(b))] for (Ta, Tb) in subsets(type_list, 2)]
    rsts = hcat(rsts...)
    @test all(isapprox.(rsts, rsts[1]; rtol=1e-5))
end

function test_ndarray(d, sz, T)
    x = rand(T, sz)
    y = rand(T, sz)
    @test_nowarn assess(d, x, y)

    T <: AbstractGray || return nothing
    for _ in 1:10
        x = rand(T, sz)
        y = rand(T, sz)

        @test assess(d, x, y) ≈ assess(d, channelview(x), channelview(y)) atol=1e-4
    end
end


text:
using Losers
using Test
using TestSetExtensions


text:
using DataStructures

module JMLJulia
end

function auto_to_julia(ex)
    @match ex begin
        :(!$a)          => :(to_julia($a))
        Expr(hd, tl...) => Expr(hd, map(auto_to_julia, tl)...)
        a               => a
    end
end

macro auto_to_julia(ex)
    auto_to_julia(ex) |> esc
end

@auto_to_julia to_julia(dexp :: DExp) =
    @match dexp begin
        DAssign(reg, val) =>
            let sym = Symbol(reg)
                :($sym = $(!val))
            end
        DClosure(freevars, argnames, body) =>
            # julia already has closures, so no need to create a new impl
            let argnames = map(Symbol, argnames)
                Expr(:function,
                    Expr(:tuple, argnames...),
                    !body
                )
            end
        DIf(cond, trueClause, falseClause) =>
            :($(!cond) ? $(!trueClause) : $(!falseClause))
        DConst(v :: AbstractString) => Meta.parse(v)
        DConst(v) => v
        DVar(s) => Symbol(s)
        DBlock(elts) => Expr(:block, [!x for x in elts]...)
        DCall(f, args) => Expr(:call, !f, [!arg for arg in args]...)

        DAttr(sub, attr) => Expr(:., !sub, QuoteNode(Symbol(attr)))
        DList(elts) =>
            let elts = [!elt for elt in elts]
                foldr(elts, init = :($nil())) do elt, prev
                    :($cons($elt, $prev))
                end
            end
        # after generation, we will transform the LineNumberNode to
        # set correct file
        DLoc(l, c, exp) => Expr(:block, LineNumberNode(l, c), !exp)
        DStaged(v) => v
        DImport(paths, name, actual) =>
            let actual = Symbol(actual),
                syms = map(Symbol, [paths..., name]),
                mod = foldl((a, b) -> :($a.$b), syms),
                _ = JMLJulia.eval(:(using $(syms...))),
                mod = JMLJulia.eval(mod)

                :($actual = $mod)
            end
        DModule(modname, exports, stmts) =>
            let stmts = [!stmt for stmt in stmts]
                # mod = Expr(
                #     :tuple,
                #     [   let k = Symbol(k),
                #             v = Symbol(v)
                #             :($k = $v)
                #         end
                #         for (k, v) in exports
                #     ]...
                # )
                Expr(:block, stmts...)
            end
    end

to_julia(
    source_code :: AbstractString
    ;
    scope = global_scope(),
    mode :: Symbol = :module,
    fname :: AbstractString = "<repl>",
    modules = OrderedDict("main" => new_module_spec("main", fname))
) =
    let rexp = runparser(source_code, Val(:rexp), Val(mode)),
        lexp = to_lexp(rexp)
        to_julia(scoping_analysis(scope, lexp, modules))
    end

max_forks_repo_name: thautwarm/Rupy


text:
struct Signal

    freq::Union{Vector,Number}
    amp::Union{Vector,Number}
    func::Union{Vector,Function}
    offset::Union{Vector,Number}
    info::Dict

    function Signal(fr, a, f = sin, o=0, info = Dict(); name::String = "")

        if !isempty(name)
            info["name"] = name
        end

        new(fr, a, f, o, info)

    end

end

function Signal(;freq, amp, f = sin, offset=0, info = Dict(), name::String = "")

    if !isempty(name)
        info["name"] = name
    end

    Signal(freq, amp, f, offset, info)

end

function latexify(s::Signal)

    name = "signal"

    if haskey(s.info, "name")
        name = s.info["name"]
    end

    if s.func |> typeof <: Function
        s = L"%$(name)(t)=%$(s.func)(2π %$(s.freq) (t + %$(s.offset) s)) · %$(s.amp)"
    else

        s1 = L"%$(name)(t)="
        s2 = [
            L" %$(s.func[i])(2π %$(s.freq[i]) (t + %$(s.offset[i]) s)) · %$(s.amp[i]) +"
            for _ = 1:length(s.func)
        ]
        s3 = *(s2...)

        s = s1 * s3[1:end-2]

    end
    return s

end

function (s::Signal)(t)

    if typeof(s.func) <: Function

        d = s.func(s.freq * (t.+ s.offset) * 2 * pi ) * s.amp

    else
        @assert length(s.freq) == length(s.amp) == length(s.func) # check if e.g. same # functions and freqs are given

        d = sum(
            s.func[i](s.freq[i] * t * 2 * pi + s.offset[i]) * s.amp[i] for
            i = 1:length(s.func) )
    end

    return d

end

function +(s1::Signal, s2::Signal)

    frequencies = vcat([s1.freq, s2.freq]...)
    amps = vcat([s1.amp, s2.amp]...)
    offsets = vcat([s1.offset, s2.offset]...)
    functions = vcat([s1.func, s2.func]...)
    dicts = merge(s1.info,s2.info)
    return Signal(freq=frequencies, amp = amps, offset=offsets, f = functions, info=dicts)

end

function *(n::Number, s::Signal)

    amp = s.amp*n
    return Signal(s.freq,amp,s.func,s.offset,s.info)

end


text:
## convert DSM to ReachableMatrix
function toReachableMatrix(DSM)

    I = zeros(Int64,size(DSM)) + Diagonal(ones(Int64,size(DSM))); #identity matrix
    R1 = zeros(Int64,size(DSM))
    calcmax = 1000

    for i  = 1 : calcmax

        R1 = (DSM+I)^i
        (Rindex) =findall(x-> 1<x , R1);
        R1[Rindex] = ones(size(Rindex));

        R2 = (DSM+I)^(i+1)
        (Rindex) =findall(x-> 1<x , R2);
        R2[Rindex] = ones(size(Rindex));

        if R1 == R2
            break
        end

    end

    return R1
end

##order Reachable matrix for seaqencing
function OrderReachable(Reachable)

    original = copy(Reachable);
    (DSMsize,) = size(Reachable);
    valuesDSM = collect(1:DSMsize);
    count = 0;
    level = zeros(Int,DSMsize);

    while count < DSMsize

        R = Reachable;
        A = Reachable';
        RA = R .* A;
        (RAsize,) = size(RA);
        elements = collect(1:RAsize);

        for i = 1 : RAsize

            if RA[i,:] == R[i,:];
                    count += 1;
                    level[count] = valuesDSM[i];
                    elements[i] = 0;
            end
        end

        elements = filter(x-> x > 0, elements);
        Reachable = Reachable[elements , elements];
        valuesDSM = valuesDSM[elements];

    end

    orderedReachable = original[level,level]

    return orderedReachable , level

end

## Sequencing DSM
function Sequencing(DSM,label)
    cDSM = copy(DSM);
    clabel = copy(label);

    Reachable = toReachableMatrix(cDSM);

    (OrderedReachable,level) = OrderReachable(Reachable);

    SequencedDSM = cDSM[level,level];
    Sequencedlabel = clabel[level];

    return SequencedDSM, Sequencedlabel

end


text:
# Tests for Atomistic.jl

using Atomistic
using AtomsBase
using DFTK
using InteratomicPotentials
using PeriodicTable
using StaticArrays
using Test
using Unitful
using UnitfulAtomic

@testset "Atomistic.jl" begin
    @testset "Unit Tests" begin
        include("unit/unit_convention.jl")
        @testset "Atomistic API Unit Tests" begin
            include("unit/api/exceptions.jl")
            include("unit/api/md_result.jl")
            include("unit/api/md_simulator.jl")
        end
        @testset "Integrations Unit Tests" begin
            include("unit/integrations/atomsbase_integration.jl")
            include("unit/integrations/dftk_integration.jl")
        end
        @testset "Implementations Unit Tests" begin
            @testset "NBodySimulator Implementation" begin
                include("unit/implementations/nbodysimulator/nbodysimulator_base.jl")
                include("unit/implementations/nbodysimulator/nbodysimulator_simulator.jl")
                include("unit/implementations/nbodysimulator/nbodysimulator_result.jl")
            end
        end
        @testset "Analysis Unit Tests" begin
            include("unit/analysis/plotting.jl")
            include("unit/analysis/visualization.jl")
        end
    end
    @testset "Integration Tests" begin
        include("integration/atomsbase_nbodysimulator.jl")
        include("integration/atomsbase_nbodysimulator_interatomicpotentials.jl")
        include("integration/atomsbase_nbodysimulator_dftk.jl")
    end
end


text:
# Routines related to fragmenting CAD entities in gmsh, while preserving physical groups

function process_material_hierarchy!(
        new_physical_groups::Dict{String, Vector{Tuple{Int32,Int32}}},
        material_hierarchy::Vector{String})
    # Get the material groups and the entities in each group
    groups = collect(keys(new_physical_groups))
    material_indices = findall(x->occursin("MATERIAL", x), groups)
    material_groups = groups[material_indices]
    # Ensure each material group is present in the hierarchy and warn now if it's not.
    # Otherwise the error that occurs later is not as easy to decipher
    for material in material_groups
        @assert material ∈ material_hierarchy "material_hierarchy does not contain: '$material'"
    end
    material_dict = Dict{String, Vector{Tuple{Int32,Int32}}}()
    all_material_entities = Tuple{Int32,Int32}[]
    for material in material_groups
        # Note that this is assignment by reference. Changes to material_dict are reflected
        # in new_physical_groups
        material_dict[material] = new_physical_groups[material]
        append!(all_material_entities, new_physical_groups[material])
    end
    # Remove duplicates
    all_material_entities = collect(Set(all_material_entities))
    # For each entity with a material, ensure that it only exists in one of material groups.
    # If it exists in more than one material group, apply the material hierarchy so that the
    # entity only has one material.
    numerical_material_hierarchy = Dict{String,Int32}()
    i = 1
    for material in material_hierarchy
        numerical_material_hierarchy[material] = i
        i += 1
    end
    for ent in all_material_entities
        materials = String[]
        for material in material_groups
            if ent ∈  material_dict[material]
                push!(materials, material)
            end
        end
        if 1 < length(materials)
            # Get the highest priority material
            mat_num = minimum([numerical_material_hierarchy[mat] for mat in materials])
            priority_mat = material_hierarchy[mat_num]
            # Pop ent from all other materials in dict
            deleteat!(materials, materials .== priority_mat)
            for material in materials
                deleteat!(material_dict[material], findfirst(x-> x == ent, material_dict[material]))
            end
        end
    end
end

function gmsh_group_preserving_fragment(object_dim_tags::Vector{Tuple{Signed,Int32}},
                                        tool_dim_tags::Vector{Tuple{Signed,Int32}};
                                        material_hierarchy::Vector{String} = String[])
    # Get all the physical groups
    old_physical_groups = Dict{String,Array{Tuple{Int32,Int32},1}}()
    groups = gmsh.model.get_physical_groups()
    names = [gmsh.model.get_physical_name(grp[1], grp[2]) for grp in groups]
    for (i, name) in enumerate(names)
        ents = gmsh.model.get_entities_for_physical_group(groups[i][1], groups[i][2])
        dim = groups[i][1]
        old_physical_groups[name] = [(dim, ent) for ent in ents]
    end

    # Fragment
    nents = length(object_dim_tags) + length(tool_dim_tags)
    @info "Fragmenting $nents entities"
    out_dim_tags, out_dim_tags_map = gmsh.model.occ.fragment(object_dim_tags, tool_dim_tags)

    # Create a dictionary of new physical groups using the parent child relationship
    # between input_dim_tags and out_dim_tags_map. The parent at index i of input_dim_tags
    # has children out_dim_tags_map[i]
    new_physical_groups = Dict{String, Vector{Tuple{Int32,Int32}}}()
    input_dim_tags = vcat(object_dim_tags, tool_dim_tags)
    # For each physical group
    for name in names
        new_physical_groups[name] = Tuple{Int32, Int32}[]
        # For each of the dim tags in the physical group
        for dim_tag in old_physical_groups[name]
            # If the dim_tag was one of the entities in the fragment
            if dim_tag ∈  input_dim_tags
                # Get its children
                index = findfirst(x->x == dim_tag, input_dim_tags)
                children = out_dim_tags_map[index]
                # Add the children to the new physical group
                for child in children
                    if child ∉  new_physical_groups[name]
                        push!(new_physical_groups[name], child)
                    end
                end
            else
                # If it wasn't in the fragment, no changes necessary.
                push!(new_physical_groups[name], dim_tag)
            end
        end
    end

    # Remove old groups and synchronize
    for name in names
        gmsh.model.remove_physical_name(name)
    end
    @debug "Synchronizing model"
    gmsh.model.occ.synchronize()

    # Process the material hierarchy if it exists so that each entity has one
    # or less material physical groups
    if 0 < length(material_hierarchy)
        process_material_hierarchy!(new_physical_groups, material_hierarchy)
    end

    # Create new physical groups
    for (i, name) in enumerate(names)
        dim = groups[i][1]
        tags = [dim_tag[2] for dim_tag in new_physical_groups[name]]
        ptag = gmsh.model.add_physical_group(dim, tags)
        gmsh.model.set_physical_name(dim, ptag, name)
    end

    return out_dim_tags
end

function gmsh_group_preserving_fragment(object_dim_tags::Vector{Tuple{Int32,Int32}},
                                        tool_dim_tags::Vector{Tuple{Int32,Int32}};
                                        material_hierarchy::Vector{String} = String[])
    # Get all the physical groups
    old_physical_groups = Dict{String, Vector{Tuple{Int32,Int32}}}()
    groups = gmsh.model.get_physical_groups()
    names = [gmsh.model.get_physical_name(grp[1], grp[2]) for grp in groups]
    for (i, name) in enumerate(names)
        ents = gmsh.model.get_entities_for_physical_group(groups[i][1], groups[i][2])
        dim = groups[i][1]
        old_physical_groups[name] = [(dim, ent) for ent in ents]
    end

    # Fragment
    nents = length(object_dim_tags) + length(tool_dim_tags)
    @info "Fragmenting $nents entities"
    out_dim_tags, out_dim_tags_map = gmsh.model.occ.fragment(object_dim_tags, tool_dim_tags)

    # Create a dictionary of new physical groups using the parent child relationship
    # between input_dim_tags and out_dim_tags_map. The parent at index i of input_dim_tags
    # has children out_dim_tags_map[i]
    new_physical_groups = Dict{String, Vector{Tuple{Int32,Int32}}}()
    input_dim_tags = vcat(object_dim_tags, tool_dim_tags)
    # For each physical group
    for name in names
        new_physical_groups[name] = Tuple{Int32,Int32}[]
        # For each of the dim tags in the physical group
        for dim_tag in old_physical_groups[name]
            # If the dim_tag was one of the entities in the fragment
            if dim_tag ∈  input_dim_tags
                # Get its children
                index = findfirst(x->x == dim_tag, input_dim_tags)
                children = out_dim_tags_map[index]
                # Add the children to the new physical group
                for child in children
                    if child ∉  new_physical_groups[name]
                        push!(new_physical_groups[name], child)
                    end
                end
            else
                # If it wasn't in the fragment, no changes necessary.
                push!(new_physical_groups[name], dim_tag)
            end
        end
    end

    # Remove old groups and synchronize
    for name in names
        gmsh.model.remove_physical_name(name)
    end
    @debug "Synchronizing model"
    gmsh.model.occ.synchronize()

    # Process the material hierarchy if it exists so that each entity has one
    # or less material physical groups
    if 0 < length(material_hierarchy)
        process_material_hierarchy!(new_physical_groups, material_hierarchy)
    end

    # Create new physical groups
    for (i, name) in enumerate(names)
        dim = groups[i][1]
        tags = [dim_tag[2] for dim_tag in new_physical_groups[name]]
        ptag = gmsh.model.add_physical_group(dim, tags)
        gmsh.model.set_physical_name(dim, ptag, name)
    end

    return out_dim_tags
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""ClientIPConfig represents the configurations of Client IP based session affinity.

    IoK8sApiCoreV1ClientIPConfig(;
        timeoutSeconds=nothing,
    )

    - timeoutSeconds::Int32 : timeoutSeconds specifies the seconds of ClientIP type session sticky time. The value must be >0 && <=86400(for 1 day) if ServiceAffinity == \"ClientIP\". Default value is 10800(for 3 hours).
"""
mutable struct IoK8sApiCoreV1ClientIPConfig <: SwaggerModel
    timeoutSeconds::Any # spec type: Union{ Nothing, Int32 } # spec name: timeoutSeconds

    function IoK8sApiCoreV1ClientIPConfig(;timeoutSeconds=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ClientIPConfig, Symbol("timeoutSeconds"), timeoutSeconds)
        setfield!(o, Symbol("timeoutSeconds"), timeoutSeconds)
        o
    end
end # type IoK8sApiCoreV1ClientIPConfig

const _property_map_IoK8sApiCoreV1ClientIPConfig = Dict{Symbol,Symbol}(Symbol("timeoutSeconds")=>Symbol("timeoutSeconds"))
const _property_types_IoK8sApiCoreV1ClientIPConfig = Dict{Symbol,String}(Symbol("timeoutSeconds")=>"Int32")
Base.propertynames(::Type{ IoK8sApiCoreV1ClientIPConfig }) = collect(keys(_property_map_IoK8sApiCoreV1ClientIPConfig))
Swagger.property_type(::Type{ IoK8sApiCoreV1ClientIPConfig }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1ClientIPConfig[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ClientIPConfig }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ClientIPConfig[property_name]

function check_required(o::IoK8sApiCoreV1ClientIPConfig)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ClientIPConfig }, name::Symbol, val)
end


text:
using Test, ChineseTokenizers, Pkg

@testset "Display of LCCRF is correct." begin
    matrix = Matrix{Any}(missing, 2, 2)
    matrix[2, 2] = "*"
    matrix[1, 1] = [1, 2]
    model = ChineseTokenizers.LinearChainConditionalRandomField{Float64}(["A", "B"], [1.0, 2.0], matrix)
    display(model)
end

@testset "Result of split is meaningful and IO is correct." begin
    tk = ChineseTokenizers.chmm()
    display(tk)
    segs = split("这是樊胜美的妈吗", tk)
    println()
    @show segs
    open(joinpath(dirname(pathof(ChineseTokenizers)), "..", "chmm"), "w") do io
        write(io, tk)
    end
    tk2 = open(joinpath(dirname(pathof(ChineseTokenizers)), "..", "chmm"), "r") do io
        read(io, HiddenMarkovModel)
    end
    @test tk == tk2
end

@testset "Added userdict should makes sense." begin
    tk1 = HiddenMarkovModel()
    xs  = "这是樊胜美的妈吗"
    segs1 = split(xs, tk1)
    @show segs1
    tk2 = HiddenMarkovModel{Float32, Int32}(; userdict = joinpath(dirname(pathof(ChineseTokenizers)), "..", "userdict"))
    segs2 = split(xs, tk2)
    @show segs2
    @test length(segs1) > length(segs2)
end

@testset "Added poswords may makes sense." begin
    tk1 = HiddenMarkovModel()
    xs  = "我的话你是不是不听了？"
    segs1 = split(xs, tk1)
    @show segs1
    tk2 = HiddenMarkovModel{Float32, Int32}(; poswords = joinpath(dirname(pathof(ChineseTokenizers)), "..", "poswords"))
    segs2 = split(xs, tk2)
    @show segs2
    @test length(segs1) >= length(segs2)
end

# @testset "Read and write is correct." begin
#     xs =  """
#     ( (IP (NP-TMP (NT 10月)
#     (NT 2号) )
#     (PU ，)
#     (NP-SBJ (-NONE- *pro*))
#     (VP (VE 有)
#     (IP-OBJ (NP-SBJ (NP (NN 患者) )
#     (CC 及)
#     (NP (NP (PN 其) )
#     (NP (NN 家属) ) ) )
#     (VP (VP (VP (VV 到)
#     (NP-OBJ-PN (NR 桐城市)
#     (NN 人民)
#     (NN 医院) ) )
#     (VP (VV 去)
#     (VP (VV 看病) ) ) )
#     (PU ，)
#     (VP (PP-PRP (P 因)
#     (IP (NP-SBJ (NN 护士)
#     (NN 嗓门) )
#     (VP (VA 大)
#     (ADVP (AD 些) ) ) ) )
#     (PU ，)
#     (ADVP (AD 便) )
#     (VP (VV 大打出手) ) )
#     (PU ，)
#     (VP (BA 将)
#     (IP-OBJ (NP-SBJ-1 (NP (DNP (NP (NN 医院) )
#     (DEG 的) )
#     (NP (NN 急诊科)
#     (NN 医生) ) )
#     (NP-PN (NR 杨辉) ) )
#     (VP (VP (ADVP (AD 狠) )
#     (VP (VV 揣)
#     (NP-OBJ (-NONE- *-1))) )
#     (PU ，)
#     (VP (VV 打成)
#     (NP-OBJ (NN 骨折) ) ) ) ) ) ) ) )
#     (PU 。) ) )
#     """
#     obj = ChineseTree(xs)
#     display(obj)
#     ctb = parsectb()
#     model = CykModel(ctb)
#     p = 0; t = nothing
#     @time p, t = cyk(posword(obj), model)
#     @show p
#     display(decnf(t))
# end


text:
#!/usr/bin/env julia

using DipoleChain

using Printf
using SparseArrays

using ArgParse
using Arpack

s = ArgParseSettings()
s.autofix_names = true
@add_arg_table s begin
    "-R"
        metavar = "R"
        help = "separation distance"
        arg_type = Float64
        required = true
    "-N"
        metavar = "N"
        help = "number of rotors"
        arg_type = Int
        required = true
    "--l-max"
        metavar = "L"
        help = "local basis truncation"
        arg_type = Int
        required = true
    "--l-total_max"
        metavar = "L"
        help = "many-body basis truncation (default: no truncation)"
        arg_type = Int
    "--sym-lp"
        metavar = "LP"
        help = "total l parity of symmetry block"
        arg_type = Int
        default = 0
    "--sym-m"
        metavar = "M"
        help = "total m value of symmetry block"
        arg_type = Int
        default = 0
    "--A-start"
        metavar = "jA"
        help = "first rotor of subsystem A"
        arg_type = Int
        default = 1
    "--A-size"
        metavar = "NA"
        help = "length of subsystem A"
        arg_type = Int
        default = 1
end
c = parse_args(ARGS, s, as_symbols=true)

R = c[:R]
N = c[:N]
l_max = c[:l_max]
l_total_max = c[:l_total_max]
if isnothing(l_total_max)
    l_total_max = l_max * N
end
sym_lp = c[:sym_lp]
sym_m = c[:sym_m]
A_start = c[:A_start]
A_size = c[:A_size]

println("[ ] Constructing basis.")
@time basis = SingleBlockBasis(N, l_max, l_total_max, sym_lp, sym_m)
println("[+] Constructed basis: $(basis.size).")

println("[ ] Constructing matrix elements.")
@time K = sparse_rot(basis)
@time V = sparse_pot(basis)
@time H = K + V/R^3
# Number of non-zero elements.
println("[+] Constructed matrix elements: $(nnz(H)).")

println("[ ] Diagonalizing.")
@time d = eigs(H, which=:SR, nev=1)
E0 = d[1][1]
wf = d[2][:,1]
# Number of iterations.
println("[+] Done diagonalizing: $(d[4]).")

# Ground state energy.
println_result("E0 = $(E0)")

println("[ ] Calculating correlations.")
@time for (op1_name, op1_f, op1_pow, op2_name, op2_f, op2_pow) in
        [("x", dense_x, 1, "x", dense_x, 1),
         ("z", dense_z, 1, "z", dense_z, 1),
         ("z", dense_z, 2, "z", dense_z, 2)]
    op1_fullname = "$(op1_name)_i"
    if op1_pow != 1
        op1_fullname = "$(op1_fullname)^$(op1_pow)"
    end
    op2_fullname = "$(op2_name)_j"
    if op2_pow != 1
        op2_fullname = "$(op2_fullname)^$(op2_pow)"
    end
    println_result("<$(op1_fullname) $(op2_fullname)>")
    op1 = op1_f(basis)^op1_pow
    op2 = op2_f(basis)^op2_pow
    for i in 1:N
        for j in 1:(i-1)
            @printf("%18s ", "")
        end
        @printf(" % .15f", expval1(basis, op1*op2, i, wf))
        for j in (i+1):N
            @printf(" % .15f", expval2(basis, op1, i, op2, j, wf))
        end
        println()
    end
end
println("[+] Calculated correlations.")

println("[ ] Constructing subsystem bases.")
@time basis_A = MultiBlockBasis(A_size, l_max, l_total_max)
@time basis_B = MultiBlockBasis(N-A_size, l_max, l_total_max)
println("[+] Constructed subsystem bases: $(basis_A.size), $(basis_B.size).")

println("[ ] Calculating reduced eigenvalues.")
@time eigvals = reduced_eigenvalues(basis, basis_A, basis_B, A_start, wf)
println("[+] Calculated reduced eigenvalues.")

# Von Neumann entropy.
println_result("SvN = $(S_vn(eigvals))")
# Order-2 Rényi entropy.
println_result("S2 = $(S_renyi(eigvals, 2))")


text:
using Documenter, OptimTensorBoard

makedocs(;
    modules=[OptimTensorBoard],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/tkf/OptimTensorBoard.jl/blob/{commit}{path}#L{line}",
    sitename="OptimTensorBoard.jl",
    authors="Takafumi Arakaki <aka.tkf@gmail.com>",
    assets=String[],
)

deploydocs(;
    repo="github.com/tkf/OptimTensorBoard.jl",
)


text:
# colouring the points

"""
    getcolour(s::RANSAC.FittedShape) = get(ColorSchemes.rainbow, rand())

Don't fail for any shape, just return a random value.
"""
getcolour(s::RANSAC.FittedShape) = get(colorschemes[:rainbow], rand())

# plane orange
# cylinder red
# sphere green
# cone blue

getcolour(s::RANSAC.FittedPlane) = get(colorschemes[:seaborn_bright], 0.1)
getcolour(s::RANSAC.FittedCylinder) = get(colorschemes[:seaborn_bright], 0.33)
getcolour(s::RANSAC.FittedSphere) = get(colorschemes[:seaborn_bright], 0.2)
getcolour(s::RANSAC.FittedCone) = get(colorschemes[:seaborn_bright], 0.0)


text:
"""
This module provides some linear fitting methods.
"""
module LinearFitting

using LinearAlgebra: dot
using Polynomials: polyder, polyfit, degree, coeffs, Poly
using MLStyle: @match

using ..FiniteStrains

export energy_strain_expansion,
       energy_strain_derivative,
       strain_volume_derivative,
       energy_volume_expansion,
       energy_volume_derivatives,
       energy_volume_derivative_at_order

energy_strain_expansion(f::Vector{<:Real}, e::Vector{<:Real}, n::Int)::Poly =
    polyfit(f, e, n)

energy_strain_derivative(p::Poly, m::Int)::Poly = polyder(p, m)

function strain_volume_derivative(s::EulerianStrain, v0::Real, v::Real, m::Int)
    m == 1 && return -1 / (3v) * cbrt(v0 / v)^2
    -(3m + 2) / (3v) * strain_volume_derivative(s, v0, v, m - 1)
end
function strain_volume_derivative(s::LagrangianStrain, v0::Real, v::Real, m::Int)
    m == 1 && return -1 / (3v) * cbrt(v / v0)^2
    -(3m - 2) / (3v) * strain_volume_derivative(s, v0, v, m - 1)
end
function strain_volume_derivative(s::NaturalStrain, v0::Real, v::Real, m::Int)
    m == 1 && return 1 / (3v)
    -m / v * strain_volume_derivative(s, v0, v, m - 1)
end
function strain_volume_derivative(s::InfinitesimalStrain, v0::Real, v::Real, m::Int)
    m == 1 && return (1 - getstrain(s, v0, v))^4 / 3 / v0
    -(3m + 1) / (3v) * strain_volume_derivative(s, v0, v, m - 1)
end

function energy_volume_expansion(
    s::FiniteStrain,
    v0::Real,
    v::Real,
    p::Poly,
    highest_order::Int = degree(p),
)
    # The zeroth order value plus values from the first to the ``highest_order`.
    p(v) + dot(
        energy_volume_derivatives(s, v0, v, p, highest_order),
        getstrain(s, v0, v) .^ collect(1:highest_order),
    )
end

function energy_volume_derivatives(
    s::FiniteStrain,
    v0::Real,
    v::Real,
    p::Poly,
    highest_order::Int,
)
    0 ≤ highest_order ≤ degree(p) ? (x = 1:highest_order) :
    throw(DomainError("The `highest_order` must be within 0 to $(degree(p))!"))
    strain_derivatives = map(m -> strain_volume_derivative(s, v0, v, m), x)
    energy_derivatives = map(f -> f(v), map(m -> energy_strain_derivative(p, m), x))
    map(
        m -> energy_volume_derivative_at_order(m)(strain_derivatives, energy_derivatives),
        x,
    )
end

function energy_volume_derivative_at_order(m::Int)::Function
    function (f::Vector{<:Real}, e::Vector{<:Real})
        @match m begin
            1 => e[1] * f[1]
            2 => e[2] * f[1]^2 + e[1] * f[1]
            3 => e[3] * f[1]^3 + 3 * f[1] * f[2] * e[2] + e[1] * f[3]
            4 => e[4] * f[1]^4 + 6 * f[1]^2 * f[2] * e[3] +
                 (4 * f[1] * f[3] + 3 * f[3]^2) * e[2] + e[1] * f[3]
            _ => error("Expansion is not defined at order = $(m)!")
        end
    end
end

end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 67c8fac0-974e-11eb-2898-d14b79387954
begin
	init_notebook = true
	using OpticSim, OpticSim.Geometry, OpticSim.Emitters
	import OpticSim.NotebooksUtils as NB

	NB.InitNotebook()
end

# ╔═╡ a891a210-87bc-11eb-311d-4d020986fe19
begin
	using DataFrames

	# defining the optical system
	sys = AxisymmetricOpticalSystem(
		DataFrame(SurfaceType = ["Object", "Standard", "Standard", "Standard", "Stop", "Standard", "Standard", "Image"],
				  Radius = [Inf, 26.777, 66.604, -35.571, 35.571, 35.571, -26.777, Inf],
				  Thickness = [Inf, 4.0, 2.0, 4.0, 2.0, 4.0, 44.748, missing],
				  Material = [OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_SK16, OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_SF2, OpticSim.GlassCat.Air, OpticSim.GlassCat.SCHOTT.N_SK16, OpticSim.GlassCat.Air, missing],
				  SemiDiameter = [Inf, 8.580, 7.513, 7.054, 6.033, 7.003, 7.506, 15.0]))
	@show sys
end

# ╔═╡ d6c83680-879e-11eb-31d4-7dbda7e93e48
begin
	init_notebook

	import PlutoUI

	using Makie
	using Makie.AbstractPlotting
	using Makie.AbstractPlotting.MakieLayout

	defs = OpticSim.NotebooksUtils.Defs("ran")

	NB.DefsClearHTML(defs)

	# this function is needed to allow the visualization scene to be displayed inside a pluto notebook
	function Makie.display(obj)
		# @info "RG: $obj"
		return obj
	end

	#text of this cell to appear in the notebook
	NB.DefsAddHTML(defs, NB.HTMLFromObj(md"Basic Initialization (code)"))

	PlutoUI.Show(MIME"text/html"(), NB.DefsHTML(defs))
end

# ╔═╡ 6ce01760-879e-11eb-2b13-4d3d07c4b4ce
md"# Introductions To Emitters"

# ╔═╡ 87954830-8854-11eb-2baf-d3b79b804e07
md"## Basic Emitters"

# ╔═╡ 9608ed90-8854-11eb-080d-4d1f60d596f5
md"""
Emitters are defined by Pixels and Spatial Layouts. An emitter has a spectrum, and an optical power distribution over the hemisphere. These are intrinsic physical properties of the emitter.

The **basic emitter** is constructed as a combination of 4 basic elements and a 3D transform. The basic elements include:
- Spectrum
- Angular Power Distribution
- Rays Origins Distribution
- Rays Directions Distribution

The OpticSim package comes with various implementations of each of these basic elements:
- Spectrum - the **generate** interface returns a tuple (power, wavelength)
  * **Uniform** - A flat spectrum bounded (default: from 450nm to 680nm). the range is sampled uniformly.
  * **DeltaFunction** - Constant wave length.
  * **Measured** - measured spectrum to compute emitter power and wavelength (created by reading CSV files – more details will follow).
- Angular Power Distribution - the interface **apply** returns an OpticalRay with modified power
  * **Lambertian**
  * **Cosine**
  * **Gaussian**
- Rays Origins Distribution - the interface **length** returns the number of samples, and **generate** returns the n'th sample.
  * **Point** - a single point
  * **RectUniform** - a uniformly sampled rectangle with user defined number of samples
  * **RectGrid** - a rectangle sampled in a grid fashion
  * **Hexapolar** - a circle (or an ellipse) sampled in an hexapolar fasion (rings)
- Rays Directions Distribution - the interface **length** returns the number of samples, and **generate** returns the n'th sample.
  * **Constant**
  * **RectGrid**
  * **UniformCone**
  * **HexapolarCone**

"""

# ╔═╡ d970a5b0-8858-11eb-1425-cdfb68e222ba
md"### Examples of basic emitters"

# ╔═╡ f21787a0-8858-11eb-19b1-914e835922eb
md"#### Point Origin Samples"

# ╔═╡ 5efca830-885b-11eb-3f76-3d9ac33369d8
md"#### Rectangle and Ellipse Origins samples"

# ╔═╡ 4d265a60-885c-11eb-2b6d-7b244ccfd172
md"#### Angular Power Distribution"

# ╔═╡ 4553fb50-88d7-11eb-3fc3-6ddfeb3f1910
md"In the following drawings, the length of the sampled rays is relative to its power."

# ╔═╡ fd3287f0-885f-11eb-279d-990c4222ba83
md"## Composite Emitters"

# ╔═╡ 03fc33b0-8860-11eb-2a20-0f8c2b1f79a3
md"""
**Composite Emitters** are emitters that can hold a list of basic emitters or composite ones. Composite emitters allow you to combine multiple basic emitters, such as an R, G and B emitters in to a Pixel emitter, which then can be combined in to a "Display" composite emitters.
The Transform data structure allow the user to define the geometric arrangment of the composited elements.

"""

# ╔═╡ 94bffc60-8860-11eb-37ca-b3203b46693a
md"### Display example"

# ╔═╡ b39c3952-8860-11eb-13a4-0b5004282a09
begin
	# construct the emitter's basic components
	local S = Spectrum.Uniform()
	local P = AngularPower.Lambertian()
	local O = Origins.RectGrid(1.0, 1.0, 3, 3)
	local D = Directions.HexapolarCone(deg2rad(5.0), 3)

	# construct the source. in this example a "pixel" source will contain only one source as we are simulating a "b/w" display.
	# for RGB displays we can combine 3 sources to simulate "a pixel".
	local Tr = Transform(Geometry.Vec3(0.5, 0.5, 0.0))
	local source1 = Sources.Source(Tr, S, O, D, P)

	# create a list of pixels - each one is a composite source
	local pixels = Vector{Sources.CompositeSource{Float64}}(undef, 0)
	for y in 1:10 # image_height
		for x in 1:10 # image_width
			# pixel position relative to the display's origin
			local pixel_position = Geometry.Vec3((x-1) * 1.1, (y-1) * 1.5, 0.0)
			local Tr = Transform(pixel_position)

			# constructing the "pixel"
			pixel = Sources.CompositeSource(Tr, [source1])

			push!(pixels, pixel)
		end
	end

	local Tr = Transform(Geometry.Vec3(0.0, 0.0, 0.0))
	my_display = Sources.CompositeSource(Tr, pixels)
	md"###### CONSTRUCTION CODE: Display is composed $(length(pixels)) pixels, esch one generating $(length(pixels[1])) rays, resulting in $(length(my_display)) total."
end

# ╔═╡ 04d97450-8860-11eb-2bda-036792395ee4
md"----"

# ╔═╡ 0b6553c2-8860-11eb-0c5e-89632a8dd612
md"## Cooke Triplet Example"

# ╔═╡ 5f77cfd0-8854-11eb-377c-ef975f9abf63
md"## Initialization Stuff - Can be ignored"

# ╔═╡ 68e9b210-87ad-11eb-0f3a-5bb2dbf7d86c
begin

	NB.DefsClearHTML(defs)

	default_resolution_info = NB.GetVarInfo(
		@bind default_resolution PlutoUI.Select(["Small", "Medium", "Large"], default="Medium")
	)
	makie_backend_info = NB.GetVarInfo(
		@bind makie_backend PlutoUI.Select(["Static", "Web"], default="Static")
	)
	show_toc_info = NB.GetVarInfo(
		@bind show_toc PlutoUI.CheckBox(default=true)
	)

	fb3_items = [
		"Options",
		"@  Drawings Size $(default_resolution_info.html)",
		"Document Options",
		"@  $(show_toc_info.html) Show Table Of Content",
		"@  Makie Backend $(makie_backend_info.html) (Leave as Static for now)",
		]

	NB.DefsAddHTML(defs,
		NB.HTMLFloatingBox(
			fb3_items,
			name="plutoui-docmenu",
			header="Document Menu",
			width="25%",
			top="45%",
		)
	)

	# wide document layout
	NB.DefsAddHTML(defs, NB.HTMLNewDocLayout())

	# add TOC
	NB.DefsAddHTML(defs,
		NB.HTMLFromObj(
			PlutoUI.TableOfContents(title = "Document Table of Content", depth = 4)
		)*
		NB.HTMLFixTOC()
	)

	#text of this cell to appear in the notebook
	NB.DefsAddHTML(defs, NB.HTMLFromObj(md"Define Floating Menu, TOC and New Layout (code)"))

	PlutoUI.Show(MIME"text/html"(), NB.DefsHTML(defs))
end

# ╔═╡ 3a5d3ba0-87ae-11eb-1717-93be0b802cab
begin
	drawing = 1
	OpticSim.NotebooksUtils.SetBackend(defs, makie_backend)
	md"Set Makie backend"
end

# ╔═╡ 96e423a0-885a-11eb-02a3-8704e8dbdab6
begin
	function Vis.scene(resolution = (1000, 1000))
		# @info "RG: Vis.Scene Replacement"
		scene, layout = MakieLayout.layoutscene(resolution = resolution)
		Vis.set_current_main_scene(scene)
		lscene = layout[1, 1] = MakieLayout.LScene(scene, scenekw = (camera = Makie.cam3d_cad!, axis_type = Makie.axis3d!, raw = false))
		Vis.set_current_3d_scene(lscene)
		return scene, lscene
	end

	function resolution()
		if (default_resolution == "Small")
			return (300, 300)
		elseif (default_resolution == "Medium")
			return (500, 500)
		else
			return (1000, 1000)
		end
	end
end

# ╔═╡ e92e18de-8847-11eb-2fdf-b90425c7327b
begin
	drawing	# notebook only - create dependency on drawing backend - comment if running in REPL

	local s = Sources.Source(origins=Origins.Point(), directions=Directions.RectGrid(π/4, π/4, 15, 15))
	Vis.draw(s, resolution=resolution(), debug=true)
end

# ╔═╡ 544d5870-8848-11eb-385e-afcce49d099d
begin
	drawing	# notebook only - create dependency on drawing backend - comment if running in REPL

	local s = Sources.Source(origins=Origins.Point(), directions=Directions.UniformCone(π/6, 1000))
	Vis.draw(s, resolution=resolution(), debug=true)
end

# ╔═╡ 78875510-8848-11eb-18d6-eba32897a99e
begin
	drawing	# notebook only - create dependency on drawing backend - comment if running in REPL

	local s = Sources.Source(origins=Origins.Point(), directions=Directions.HexapolarCone(π/6, 10))
	Vis.draw(s, resolution=resolution(), debug=true)
end

# ╔═╡ 69f94f40-885b-11eb-08b0-e9faf7969dda
begin
	drawing	# notebook only - create dependency on drawing backend - comment if running in REPL

	local s = Sources.Source(origins=Origins.RectGrid(1.0, 1.0, 10, 10), directions=Directions.Constant())
	Vis.draw(s, resolution=resolution(), debug=true)
end

# ╔═╡ ac0ee1b0-885b-11eb-0395-a19b37c5a472
begin
	drawing	# notebook only - create dependency on drawing backend - comment if running in REPL

	local s = Sources.Source(origins=Origins.RectGrid(1.0, 1.0, 3, 3), directions=Directions.HexapolarCone(π/6, 10))
	Vis.draw(s, resolution=resolution(), debug=true)
end

# ╔═╡ e14e85c0-8846-11eb-1f4a-2768e655321c
begin
	drawing	# notebook only - create dependency on drawing backend - comment if running in REPL

	local s = Sources.Source(origins=Origins.Hexapolar(2, 8.0, 8.0), directions=Directions.RectGrid(π/8, π/8, 5, 5))
	Vis.draw(s, resolution=resolution(), debug=true)
end

# ╔═╡ 69d61a60-885c-11eb-3f97-75e98121ec8c
begin
	drawing	# notebook only - create dependency on drawing backend - comment if running in REPL

	local s = Sources.Source(
		origins=Origins.Hexapolar(1, 8.0, 8.0),
		directions=Directions.RectGrid(π/6, π/6, 15, 15),
		# try to comment and un-comment the following lines to see the diffrence
		power=AngularPower.Cosine(10.0)
	)
	Vis.draw(s, resolution=resolution(), debug=true)
end

# ╔═╡ 0ae26f80-885d-11eb-2c2a-ad23b34daac1
begin
	drawing	# notebook only - create dependency on drawing backend - comment if running in REPL

	local s = Sources.Source(
		origins=Origins.RectGrid(1.0, 1.0, 3, 3),
		directions=Directions.HexapolarCone(π/6, 10),
		# try to comment and un-comment the following lines to see the diffrence
		power=AngularPower.Gaussian(2.0, 2.0)
	)
	Vis.draw(s, resolution=resolution(), debug=true)
end

# ╔═╡ 941bb0a0-8861-11eb-09fd-87e2022e068e
begin
	drawing	# notebook only - create dependency on drawing backend - comment if running in REPL

	Vis.draw(my_display;  resolution=resolution(), color=:red, debug=false)
	local rays = AbstractArray{OpticalRay{Float64, 3}}(collect(my_display))
	Vis.draw!(rays; color=:red, debug=false)
end

# ╔═╡ def239f0-87bc-11eb-2edb-2f859ac41bee
begin
	drawing	# notebook only - create dependency on drawing backend - comment if running in REPL

	local s1 = Sources.Source(sourcenum = 1, origins=Origins.Hexapolar(5, 8.0, 8.0), directions=Directions.Constant())

	# create a second source and rotate it by -5 degs
	local s2 = Sources.Source(sourcenum = 2, origins=Origins.Hexapolar(5, 8.0, 8.0), directions=Directions.Constant(), transform=Transform(zeros(Geometry.Vec3), rotationY(deg2rad(-8)) * unitZ3()))

	# create the "ray generator"
	local combined_sources = Sources.CompositeSource(Transform(Geometry.Vec3(0.0, 0.0, 10.0), unitZ3() * -1), [s1 s2])

	# and draw the system + the generated rays
	Vis.drawtracerays(sys, raygenerator = combined_sources,  resolution=resolution(), test = true, trackallrays = true, colorbysourcenum = true, drawgen = false)
end

# ╔═╡ Cell order:
# ╟─6ce01760-879e-11eb-2b13-4d3d07c4b4ce
# ╟─87954830-8854-11eb-2baf-d3b79b804e07
# ╟─9608ed90-8854-11eb-080d-4d1f60d596f5
# ╟─d970a5b0-8858-11eb-1425-cdfb68e222ba
# ╟─f21787a0-8858-11eb-19b1-914e835922eb
# ╠═e92e18de-8847-11eb-2fdf-b90425c7327b
# ╠═544d5870-8848-11eb-385e-afcce49d099d
# ╠═78875510-8848-11eb-18d6-eba32897a99e
# ╠═5efca830-885b-11eb-3f76-3d9ac33369d8
# ╠═69f94f40-885b-11eb-08b0-e9faf7969dda
# ╠═ac0ee1b0-885b-11eb-0395-a19b37c5a472
# ╠═e14e85c0-8846-11eb-1f4a-2768e655321c
# ╟─4d265a60-885c-11eb-2b6d-7b244ccfd172
# ╟─4553fb50-88d7-11eb-3fc3-6ddfeb3f1910
# ╠═69d61a60-885c-11eb-3f97-75e98121ec8c
# ╠═0ae26f80-885d-11eb-2c2a-ad23b34daac1
# ╟─fd3287f0-885f-11eb-279d-990c4222ba83
# ╟─03fc33b0-8860-11eb-2a20-0f8c2b1f79a3
# ╟─94bffc60-8860-11eb-37ca-b3203b46693a
# ╠═b39c3952-8860-11eb-13a4-0b5004282a09
# ╠═941bb0a0-8861-11eb-09fd-87e2022e068e
# ╟─04d97450-8860-11eb-2bda-036792395ee4
# ╟─0b6553c2-8860-11eb-0c5e-89632a8dd612
# ╟─a891a210-87bc-11eb-311d-4d020986fe19
# ╠═def239f0-87bc-11eb-2edb-2f859ac41bee
# ╟─5f77cfd0-8854-11eb-377c-ef975f9abf63
# ╟─67c8fac0-974e-11eb-2898-d14b79387954
# ╟─d6c83680-879e-11eb-31d4-7dbda7e93e48
# ╟─68e9b210-87ad-11eb-0f3a-5bb2dbf7d86c
# ╟─3a5d3ba0-87ae-11eb-1717-93be0b802cab
# ╟─96e423a0-885a-11eb-02a3-8704e8dbdab6


text:
using Test
using MIPVerify

@testset "sequential.jl" begin

    @testset "Sequential" begin
        nnparams = Sequential(
            [
                Conv2d(ones(1, 4, 4, 16), ones(16), 2),
                ReLU(),
                Flatten([3, 4, 1, 2]),
                Linear(ones(1000, 20), ones(20)),
                MaskedReLU([-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0]),
                Linear(ones(20, 10), ones(10)),
            ],
            "testnet",
        )
        io = IOBuffer()
        Base.show(io, nnparams)
        @test String(take!(io)) == """
        sequential net testnet
          (1) Conv2d(4, 16, kernel_size=(1, 4), stride=(2, 2), padding=same)
          (2) ReLU()
          (3) Flatten(): flattens 4 dimensional input, with dimensions permuted according to the order [3, 4, 1, 2]
          (4) Linear(1000 -> 20)
          (5) MaskedReLU with expected input size (20,). (10 zeroed, 7 as-is, 3 rectified).
          (6) Linear(20 -> 10)
        """
    end

end


text:
using ParallelTreeSampling
using Test

@testset "ParallelTreeSampling.jl" begin
    # Write your tests here.
end


text:
using Test
using LogicCircuits

include("../helper/validate_sdd.jl")

@testset "CNF compiler tests" begin

   cnf = zoo_cnf("easy/C17_mince.cnf")
   vtree = zoo_vtree("easy/C17_mince.min.vtree");
   mgr = SddMgr(vtree)

   @test_throws Exception compile(mgr.left, cnf)
   @test_throws Exception compile(right_most_descendent(mgr), cnf)

   r = mgr(cnf)
   @test compile(mgr, cnf) === r

   @test respects_vtree(r)
   @test respects_vtree(r, vtree)
   @test !respects_vtree(r, Vtree(num_variables(vtree), :balanced))

   cnfs = [
            ("easy","C17_mince",32,92,45)
            ("easy","majority_mince",32,132,61)
            ("easy","b1_mince",8,169,84)
            ("easy","cm152a_mince",2048,127,62)
            ("iscas89","s208.1.scan",262144,1942,927)
          ]

   for (suite, name, count, size, nodes) in cnfs

      cnf = zoo_cnf("$suite/$name.cnf")
      vtree = zoo_vtree("$suite/$name.min.vtree");

      mgr = SddMgr(vtree)
      # cnfΔ = @time compile_cnf(mgr, cnf)
      cnfΔ = compile(mgr, cnf)

      validate(cnfΔ)

      @test model_count(cnfΔ) == count
      @test sdd_size(cnfΔ) == size
      @test sdd_num_nodes(cnfΔ) == nodes

   end

end


text:
# Extract the data
flex_data = getflexibilitydata(m)
flex_data2 = getflexibilitydata(m2)
@test_throws ErrorException getflexibilitydata(Model())

# Check data type
@test isa(flex_data, FlexibilityAnalysis.FlexibilityData)
@test isa(flex_data2, FlexibilityAnalysis.FlexibilityData)

# Check F index extractor function
@test getflexibilityindex(m) == nothing

# Check prelimnary errors
@test_throws ErrorException getconfidencelevel(m)

# checking the model printing
print_to_string(f::Function, args...) = begin
    io = IOBuffer()
    f(io,args...)
    String(take!(io))
end
@test print_to_string(JuMP.show, m.ext[:FlexData].flexibility_constraints[1]) == "0.67*Qc + -2*T[2] + -x + 100 >= 0"
@test print_to_string(JuMP.show, m.ext[:FlexData].flexibility_constraints[4]) == "1*Qc + -1.5*T[1] + -2*T[2] + -1*T[3] + 2044 <= 0"


text:
export iauFw2m
"""
Form rotation matrix given the Fukushima-Williams angles.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
    gamb     double         F-W angle gamma_bar (radians)
    phib     double         F-W angle phi_bar (radians)
    psi      double         F-W angle psi (radians)
    eps      double         F-W angle epsilon (radians)

Returned:
    r        double[3][3]   rotation matrix

Notes:

    1. Naming the following points:

        e = J2000.0 ecliptic pole,
        p = GCRS pole,
        E = ecliptic pole of date,
    and   P = CIP,

    the four Fukushima-Williams angles are as follows:

    gamb = gamma = epE
    phib = phi = pE
    psi = psi = pEP
    eps = epsilon = EP

    2. The matrix representing the combined effects of frame bias,
    precession and nutation is:

    NxPxB = R_1(-eps).R_3(-psi).R_1(phib).R_3(gamb)

    3. Three different matrices can be constructed, depending on the
    supplied angles:

    .  To obtain the nutation x precession x frame bias matrix,
    generate the four precession angles, generate the nutation
    components and add them to the psi_bar and epsilon_A angles,
    and call the present function.

    .  To obtain the precession x frame bias matrix, generate the
    four precession angles and call the present function.

    .  To obtain the frame bias matrix, generate the four precession
    angles for date J2000.0 and call the present function.

    The nutation-only and precession-only matrices can if necessary
    be obtained by combining these three appropriately.

Called:
    iauIr        initialize r-matrix to identity
    iauRz        rotate around Z-axis
    iauRx        rotate around X-axis

Reference:

    Hilton, J. et al., 2006, Celest.Mech.Dyn.Astron. 94, 351

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauFw2m(double gamb, double phib, double psi, double eps,
#              double r[3][3])

function iauFw2m(gamb::Real, phib::Real, psi::Real, eps::Real)
    r = zeros(Float64, 3, 3)

    ccall((:iauFw2m, libsofa_c), Cvoid,
          (Cdouble, Cdouble, Cdouble, Cdouble,
          Ptr{Cdouble}),
           convert(Float64, gamb), convert(Float64, phib),
           convert(Float64, psi), convert(Float64, eps),
           r)

    return SMatrix{3,3}(r')
 end


text:
g_distance = basic_osm_graph_stub()
g_time = basic_osm_graph_stub(:time)

# Basic use
edge = iterate(keys(g_distance.edge_to_highway))[1]
node1_id = edge[1]
node2_id = edge[2]
path = shortest_path(g_distance, node1_id, node2_id)
@test path[1] == node1_id
@test path[2] == node2_id

# Test using nodes rather than node IDs get same results
path_from_nodes = shortest_path(g_distance, g_distance.nodes[node1_id], g_distance.nodes[node2_id])
@test path_from_nodes == path

# Pass in weights directly
path_with_weights = shortest_path(g_distance, g_distance.nodes[node1_id], g_distance.nodes[node2_id]; weights=g_distance.weights)
@test path_with_weights == path

# Also test astar doesn't error
path_astar = shortest_path(g_distance, node1_id, node2_id, algorithm=:astar)
@test path_astar==path

# Test edge weight sum equals the weight in g_distance.weights
@test total_path_weight(g_distance, path) == g_distance.weights[g_distance.node_to_index[node1_id],g_distance.node_to_index[node2_id]]
@test total_path_weight(g_distance, path) == sum(weights_from_path(g_distance, path))
n_nodes = length(g_distance.nodes)
ones_weights = ones(n_nodes, n_nodes)
@test total_path_weight(g_distance, path, weights=ones_weights) == 1 * (length(path) - 1)
@test all(weights_from_path(g_distance, path, weights=ones_weights) .== 1)

# Test time weights
path_time_weights = shortest_path(g_time, node1_id, node2_id)
path_time_weights_astar = shortest_path(g_time, node1_id, node2_id, algorithm=:astar)
@test path_time_weights[1] == node1_id
@test path_time_weights[2] == node2_id
@test path_time_weights == path_time_weights_astar

edge_speed = g_distance.highways[g_distance.edge_to_highway[edge]].tags["maxspeed"]
@test isapprox(total_path_weight(g_distance, path) / total_path_weight(g_time, path), edge_speed)

# Test paths we know the result of from the stub graph
path = shortest_path(g_time, 1001, 1004)
@test path == [1001, 1006, 1007, 1004] # this highway is twice the speed so should be quicker
path = shortest_path(g_distance, 1001, 1004)
@test path == [1001, 1002, 1003, 1004]

# Test restriction (and bug fixed in PR #42). Restriction in stub stops 1007 -> 1004 -> 1003 right turn
path = shortest_path(g_distance, 1007, 1003, restrictions=nothing)
@test path == [1007, 1004, 1003]
path = shortest_path(g_distance, 1007, 1003) # using g.indexed_restrictions
@test path == [1007, 1006, 1001, 1002, 1003]

# Test bug fixed in PR #42
g_temp = deepcopy(g_distance)
g_temp.weights[g_temp.node_to_index[1004], g_temp.node_to_index[1003]] = 100
path = shortest_path(g_temp, 1007, 1003, restrictions=nothing)
@test path == [1007, 1006, 1001, 1002, 1003]


text:
# This file defines the const varibles

const ρ = 1.225 #空气密度 量纲kg*m^-3
const ν = 1.46e-4 #运动黏性系数 无量纲
const g = 9.8 #重力加速度  量纲kg*m*s^-2
const v_sound = 343.2 #音速 量纲m/s
const NR = 2  #桨叶片数
const Nbe = 10 #旋翼分段数量
const R = 1.05  #旋翼半径 量纲m
const ecut = 0.1 #桨叶根切比例 无量纲
const m_ = 0.35/R #桨叶质量密度 量纲kg/m
const Ω = 100  #旋翼转速 量纲rad/s
const αs = -2*π/180.0  #旋翼轴倾角  量纲rad
const Kβ = 0.0  #桨叶根部挥舞弹簧刚度 量纲？？？
const Iβ = m_/3*(1-ecut^3)*R^3  #桨叶挥舞惯量 量纲kg*m^2
# const R0 = 0.382  #lenght of root cut 量纲m
# const R1 = 0.136  #length of tip loss 量纲m
const vair = 50 # 来流速度 量纲 m/s
const v_air = [vair*cos(αs),0.0,vair*sin(αs)] #forward wind speed  量纲m/s
const μ_air = v_air[1]/(Ω*R)  #来流在桨盘edgewise方向分量  无量纲
const λ_air = v_air[3]/(Ω*R)  #来流在桨盘轴向分量  无量纲
const A = π*R^2 #参考面积 量纲为m^2
const fnonc = ρ*A*Ω*R^2*Ω #力的无量纲化参数 量纲kg*m*s^-2
const mnonc = ρ*A*Ω^2*R^3 #力矩的无量纲化参数 量纲kg*m^2/s^2
const T = 250 #飞行器重量 (量纲为kg*m*s^-2)
const CT = T/fnonc  #无量纲重量系数
const dpsideg = 10  # 方位角不仅长度（量纲为deg）
const dψ = dpsideg*π/180 #方位角步进步长 (量纲为rad)
const npsi = 360/dpsideg # 周向分割步数
const βang0 = 0.0 # 挥舞角初值
const dβ0 = 0.0 # 挥舞角速度初值
const θ7 = 2.0/180*π # 旋翼总距
const θtw = -5./180*π # 桨叶扭转角
const thelat = -0/180*π # 横向周期变距
const thelon = 0/180*π # 纵向周期变距

ch = Array{Float64}(NR,Nbe) #桨叶弦长
rb = Array{Float64}(NR,Nbe)  #桨叶叶素点分段径向长度值
dr = Array{Float64}(NR,Nbe)  #桨叶叶素分段长度值

for k in 1:NR
  for i in 1:Nbe
    ch[k,i] = 0.35
    dr[k,i] = R/Nbe
    rb[k,i] = R/Nbe*(i-1)+dr[k,i]/2
  end
end

max_stars_repo_name: dingzhw/RotorTrim

max_issues_repo_name: helicopterliu/RotorTrim

max_forks_repo_name: helicopterliu/RotorTrim


text:
using Base:@pure

export AnalyticTransform, ExpTransform, PolyTransform,
       IdTransform, CoulombTransform, MorseTransform
import Base: convert, Dict, ==

# TODO: rename this AnalyticTransform
struct AnalyticTransform{FT, FDT, VT} <: SpaceTransform
   id::String
   f::FT
   f_d::FDT
   valid::VT
end

@pure transform(t::AnalyticTransform, r::Number) = t.f(r)
@pure transform_d(t::AnalyticTransform, r::Number) = t.f_d(r)

AnalyticTransform(id, f, f_d) = AnalyticTransform(id, f, f_d, Val(Symbol(id)))

==(T1::AnalyticTransform, T2::AnalyticTransform) = (T1.id == T2.id)

function AnalyticTransform(strans::String; fwrap = true)
   strans0 = strans
   # if @analytic is a substring then we don't do anything
   if !occursin(r"@analytic", strans)
      # but if not, then we next look for ->
      if !occursin(r"->", strans)
         # if -> is not a substring then we assume that strans is of the form
         # "(r0/r)^4" or similar i.e. explicitly uses r as the variable.
         strans = "@analytic r -> " * strans
      else
         # @analytic is not a substring but -> is a substring. e.g.
         # r -> (r0/r)^3 or s -> (r0/s)^3. we add @analytic to compile it
         strans = "@analytic " * strans
      end
   end
   ftrans = eval(Meta.parse(strans))
   if fwrap
      return AnalyticTransform(strans0, F64fun(ftrans.f), F64fun(ftrans.f_d))
   end
   return AnalyticTransform(strans0, ftrans.f, ftrans.f_d)
end

Dict(t::AnalyticTransform) = Dict( "__id__" => "NBodyIPs_AnalyticTransform",
                                "defn" => t.id )
AnalyticTransform(D::Dict) = AnalyticTransform(D["defn"])
convert(::Val{:NBodyIPs_AnalyticTransform}, D::Dict) = AnalyticTransform(D)
hash(::BASIS, t::AnalyticTransform) = hash((AnalyticTransform, t.id))

"""
Implements the space transform
```
r -> exp( - A * (r/r0 - 1) )
```

Constructor:
```
ExpTransform(A, r0)
```
"""
struct ExpTransform{TA, T} <: SpaceTransform
   A::TA
   r0::T
end

@pure transform(t::ExpTransform, r::Number) = exp( - t.A * (r/t.r0 - 1))
@pure transform_d(t::ExpTransform, r::Number) = (-t.A/t.r0) * exp( - t.A * (r/t.r0 - 1) )

# e^{-A (r/r0 - 1)) = x
# r/r0 - 1 = - (log x) / A
# r = r0 * (1 - (log x) / A)
# TODO: write a test that checks this inverse
inv_transform(t::ExpTransform, x::Number) = t.r0 * (1 - log(x) / t.A)

Dict(t::ExpTransform) =
   Dict("__id__" => "NBodyIPs_ExpTransform", "A" => t.A, "r0" => t.r0)
ExpTransform(D::Dict) = ExpTransform(D["A"], D["r0"])
convert(::Val{:NBodyIPs_ExpTransform}, D::Dict) = ExpTransform(D)
hash(::BASIS, t::ExpTransform) = hash(t)

"""
Implements the space transform
```
r -> (r0/r)^p
```

Constructor:
```
PolyTransform(p, r0)
```
"""
struct PolyTransform{TP, T} <: SpaceTransform
   p::TP
   r0::T
end

@pure transform(t::PolyTransform, r::Number) = @fastmath((t.r0/r)^t.p)
@pure transform_d(t::PolyTransform, r::Number) = @fastmath((-t.p/t.r0) * (t.r0/r)^(t.p+1))

# x = (r0/r)^p
# r x^{1/p} = r0
inv_transform(t::PolyTransform, x::Number) = t.r0 / x^(1/t.p)

Dict(t::PolyTransform) =
   Dict("__id__" => "NBodyIPs_PolyTransform",
        "p" => t.p,
        "r0" => t.r0)
PolyTransform(D::Dict) = PolyTransform(D["p"], D["r0"])
convert(::Val{:NBodyIPs_PolyTransform}, D::Dict) = PolyTransform(D)
hash(::BASIS, t::PolyTransform) = hash(t)

"""
Implements the space transform
```
r -> 1/r
```
Constructor: `CoulombTransform()`
"""
struct CoulombTransform <: SpaceTransform
end
@pure transform(t::CoulombTransform, r::Number) = @fastmath(1/r)
@pure transform_d(t::CoulombTransform, r::Number) = @fastmath(-1/(r*r))
inv_transform(t::CoulombTransform, x::Number) = 1/x

Dict(t::CoulombTransform) = Dict("__id__" => "NBodyIPs_CoulombTransform")
CoulombTransform(D::Dict) = CoulombTransform()
convert(::Val{:NBodyIPs_CoulombTransform}, D::Dict) = CoulombTransform(D)
hash(::BASIS, t::CoulombTransform) = hash(t)

"""
Implements the space transform
```
r -> exp(-r)
```
Constructor: `MorseTransform()`
"""
struct MorseTransform <: SpaceTransform
end
@pure transform(t::MorseTransform, r::Number) = @fastmath(exp(-r))
@pure transform_d(t::MorseTransform, r::Number) = @fastmath(-exp(-r))
inv_transform(t::MorseTransform, x::Number) = - log(x)

Dict(t::MorseTransform) = Dict("__id__" => "NBodyIPs_MorseTransform")
MorseTransform(D::Dict) = MorseTransform()
convert(::Val{:NBodyIPs_MorseTransform}, D::Dict) = MorseTransform(D)
hash(::BASIS, t::MorseTransform) = hash(t)

"""
Implements the space transform
```
r -> r
```
Constructor: `IdTransform()`
"""
struct IdTransform <: SpaceTransform
end

@pure transform(t::IdTransform, r::Number) = r
@pure transform_d(t::IdTransform, r::Number) = one(typeof(r))
inv_transform(t::IdTransform, x::Number) = x

Dict(t::IdTransform) = Dict("__id__" => "NBodyIPs_IdTransform")
IdTransform(D::Dict) = IdTransform()
convert(::Val{:NBodyIPs_IdTransform}, D::Dict) = IdTransform(D)
hash(::BASIS, t::IdTransform) = hash(t)

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
#######################################################################
# Utils
#######################################################################

include("equality.jl")
include("status-info.jl")

max_stars_repo_name: julbinb/julia-pkgsrc-downloader

max_issues_repo_name: julbinb/julia-pkgsrc-downloader

max_forks_repo_name: julbinb/julia-pkgsrc-downloader