text:
module QJuliaRegisters

# Plain baseline registers
const half2 = NTuple{2, VecElement{Float16}}
const half4 = NTuple{4, VecElement{Float16}}
#
const float2  = NTuple{2,  VecElement{Float32}}
const float4  = NTuple{4,  VecElement{Float32}}
const float8  = NTuple{8,  VecElement{Float32}}
const float16 = NTuple{16, VecElement{Float32}}
#
const double2  = NTuple{2,  VecElement{Float64}}
const double4  = NTuple{4,  VecElement{Float64}}
const double8  = NTuple{8,  VecElement{Float64}}

# Helper functions for plain registers

function register_type(T)
  if     T == Float16 || T == ComplexF16 || T == half2 || T == half4
    return Float16
  elseif T == Float32 || T == ComplexF32 || T == float2 || T == float4 || T == float8 || T == float16
    return Float32
  elseif T == Float64 || T == ComplexF64 || T == double2 || T == double4 || T == double8
    return Float64
  elseif T == BigFloat || T == Complex{BigFloat}
    return BigFloat
  else
    error("Cannot deduce a type.")
  end

  return nothing
end

function register_size(T)

  if     T == Float16 || T == Float32 || T == Float64 || T == BigFloat
    return 1
  elseif T == ComplexF16 || T == ComplexF32 || T == ComplexF64 || T == Complex{BigFloat} || T == half2 || T == float2 || T == double2
    return 2
  elseif T == double4 || T == float4 || T == half4
    return 4
  elseif T == double8 || T == float8
    return 8
  elseif T == float16
    return 16
  else
    error("Cannot deduce the register size (type is not supported).")
  end

  return nothing
end

# More generic registers

abstract type GenericRegister{T<:Real,N} end

mutable struct FloatN{T<:AbstractFloat,N} <: GenericRegister{T,N}

  val::NTuple{N,  VecElement{T}}

  FloatN{T,N}() where {T} where {N} = new( NTuple{N,  VecElement{T}}(ntuple(i->0, N)) )
  FloatN{T,N}(src::FloatN) where {T} where {N}  = new( src.val  )
  FloatN{T,N}(reg::NTuple{N, VecElement{T}}) where {T} where {N}  = new( reg )

end

# re-implement plain complex type for compatibility

mutable struct Complex2{T<:AbstractFloat, N} <: GenericRegister{T,N}

  val::Complex{T}

  Complex2{T,N}() where {T} where {N} = new( Complex{T}(0.0, 0.0) )
  Complex2{T,N}(src::Complex2) where {T} where {N} = new( src.val  )
  Complex2{T,N}(reg::Complex{T}) where {T} where {N} = new( reg )
  Complex2{T,N}(rea::T, img::T) where {T} where {N} = new( Complex{T}(rea, img) )
  Complex2{T,N}(tpl::NTuple{2, VecElement{T}}) where {T} where {N} = new( Complex{T}(tpl[1], tpl[2]) )

end

mutable struct IntN{T<:Integer,N} <: GenericRegister{T,N}

  val::NTuple{N,  VecElement{T}}

  IntN{T,N}() where {T} where {N} = new( NTuple{N,  VecElement{T}}(ntuple(i->0, N)) )
  IntN{T,N}(src::IntN) where {T} where {N}  = new( src.val  )
  IntN{T,N}(reg::NTuple{N, VecElement{T}}) where {T} where {N}  = new( reg )

end

# Floating point registers
Half2    = FloatN{Float16, 2}
Half4    = FloatN{Float16, 4}
Half8    = FloatN{Float16, 8}
Half16   = FloatN{Float16, 16}
Single2  = FloatN{Float32, 2}
Single4  = FloatN{Float32, 4}
Single8  = FloatN{Float32, 8}
Single16 = FloatN{Float32, 16}
Double2  = FloatN{Float64, 2}
Double4  = FloatN{Float64, 4}
Double8  = FloatN{Float64, 8}
# Big float type
BigDouble2  = FloatN{BigFloat, 2}

# Additional complex types (are they compatible with Half2, Single2, Double2?)
ComplexH = Complex2{Float16,2}
ComplexS = Complex2{Float32,2}
ComplexD = Complex2{Float64,2}
ComplexB = Complex2{BigFloat,2}

# Signed: BigInt Int128 Int16 Int32 Int64 Int8
# Unsigned: UInt128 UInt16 UInt32 UInt64 UInt8

Int2     = IntN{Int32, 2}
Int4     = IntN{Int32, 4}
UInt2    = IntN{UInt32, 2}
UInt4    = IntN{UInt32, 4}
LongInt2 = IntN{Int64, 2}
LongInt4 = IntN{Int64, 4}

# Helper methods
function register_type(reg::GenericRegister{T,N}) where T where N; return T; end
function register_size(reg::GenericRegister{T,N}) where T where N; return N; end

end #QJuliaRegisters

max_stars_repo_name: alexstrel/QJulia

max_issues_repo_name: alexstrel/QJulia

max_forks_repo_name: alexstrel/QJulia


text:
module ReinforcementLearningEnvironmentDiscrete
using ReinforcementLearningBase, SparseArrays, POMDPs, POMDPModels, Random, LinearAlgebra
import StatsBase: sample, wsample
import GR: imshow
import ReinforcementLearningBase: interact!, getstate, reset!, plotenv, actionspace

const ENV_RNG = MersenneTwister(0)

include("mdp.jl")
include("randommdp.jl")
include("cliffwalking.jl")
include("maze.jl")
include("pomdps.jl")

export MDP, POMDPEnv, MDPEnv, DiscreteMaze, treeMDP, DetMDP, DetTreeMDP,
DetTreeMDPwithinrew, StochMDP, StochTreeMDP, AbsorbingDetMDP, CliffWalking,
DeterministicNextStateReward, DeterministicStateActionReward, NormalNextStateReward,
NormalStateActionReward

end # module

ext: jl


text:
using Pkg
Pkg.activate(pwd())

# # Optimal control
# ## Computing trajectories with no control

struct PMSM{T<:Real}
    ρ::T
    ω::T
    A::Matrix{T}
    invA::Matrix{T}

    function PMSM(ρ, ω)
        A = -ρ*[1 0; 0 1] -ω*[0 -1; 1 0]
        return new{eltype(A)}(ρ, ω, A, inv(A))
    end
end

#+

function expA(p::PMSM, t)
    ρ, ω = p.ρ, p.ω
    return exp(-ρ*t)*[cos(ω*t) sin(ω*t); -sin(ω*t) cos(ω*t)]
end

#+

ρ = 0.1
ω = 2
x0 = [0;-0.5]
q = [1;0]

ps = PMSM(ρ, ω)

# ### Exercise:
#
# Verify that the matrix exponential is computed correctly and that it is different from
# the elementwise exponential.
#
# **Hint**: The matrix exponential can also be computed directly by the `exp` function from
# the `LinearAlgebra` package.
#
# ---
# ### Solution:

# ---
#
# ### Exercise:
#
# Write two function `trajectory_fin_diff` and `trajectory_exact` which compute the
# trajectory. The first one should use the finite difference method to discretize the
# time, while the second one should use the closed-form formula.
#
# Plot both trajectories on time interval ``[0,10]`` with time discretization step
# ``\Delta t=0.01``. Since ``x(t)`` is a two-dimensional vector, plot each component
# on one axis.
#
# ---
# ### Solution:

# ---

ts = 0:0.0001:10

xs1 = trajectory_fin_diff(ps, x0, ts, q)
xs2 = trajectory_exact(ps, x0, ts, q)

plot(xs1[1,:], xs1[2,:], label="Finite differences");
plot!(xs2[1,:], xs2[2,:], label="True value")

# ## Solving the optimal control problem
# ### Exercise:
#
# Solve the optimal time for ``x_{\rm tar}= (0.25, -0.5)`` with the maximum voltage
# ``U_{\rm max} = 0.1``.
#
# **Hint**: To solve the equation above for ``t``, use the bisection method.

function bisection(f, a, b; tol=1e-6)
    fa = f(a)
    fb = f(b)
    fa == 0 && return a
    fb == 0 && return b
    fa*fb > 0 && error("Wrong initial values for bisection")
    while b-a > tol
        c = (a+b)/2
        fc = f(c)
        fc == 0 && return c
        if fa*fc > 0
            a = c
            fa = fc
        else
            b = c
            fb = fc
        end
    end
    return (a+b)/2
end

# ---
# ### Solution:

# ---

function trajectory_control(p::PMSM, x0, ts, q, U_max, p0)
    xs = zeros(length(x0), length(ts))

    for (i, t) in enumerate(ts)
        eAt  = expA(p, t)
        emAt = expA(p, -t)
        xs[:, i] = eAt*(x0 + p.invA * (I - emAt)*q + U_max/ρ*(exp(ρ*t) - 1)*p0)
    end
    return xs
end

#+

p0 = ps.ρ/(U_max*(exp(ps.ρ*τ)-1))*(expA(ps, -τ)*x_t - x0 - ps.invA*(I-expA(ps, -τ))*q)
p0 /= norm(p0)

#+

ts = range(0, τ; length=100)
traj = trajectory_control(ps, x0, ts, q, U_max, p0)

#+

plot(traj[1,:], traj[2,:], label="Optimal trajectory");
scatter!([x0[1]], [x0[2]], label="Starting point");
scatter!([x_t[1]], [x_t[2]], label="Target point")

#+

ts = 0:0.01:10
plt = plot();

for α = 0:π/4:2*π
    trj = trajectory_control(ps, x0, ts, q, U_max, [sin(α); cos(α)])
    plot!(plt, trj[1,:], trj[2,:], label="")
end

display(plt)


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: iotdeviceadvisor
using AWS.Compat
using AWS.UUIDs

"""
    create_suite_definition()
    create_suite_definition(params::Dict{String,<:Any})

Creates a Device Advisor test suite.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"suiteDefinitionConfiguration"`: Creates a Device Advisor test suite with suite
  definition configuration.
- `"tags"`: The tags to be attached to the suite definition.
"""
create_suite_definition(; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("POST", "/suiteDefinitions"; aws_config=aws_config)
create_suite_definition(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("POST", "/suiteDefinitions", params; aws_config=aws_config)

"""
    delete_suite_definition(suite_definition_id)
    delete_suite_definition(suite_definition_id, params::Dict{String,<:Any})

Deletes a Device Advisor test suite.

# Arguments
- `suite_definition_id`: Deletes a Device Advisor test suite with defined suite Id.

"""
delete_suite_definition(suiteDefinitionId; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("DELETE", "/suiteDefinitions/$(suiteDefinitionId)"; aws_config=aws_config)
delete_suite_definition(suiteDefinitionId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("DELETE", "/suiteDefinitions/$(suiteDefinitionId)", params; aws_config=aws_config)

"""
    get_suite_definition(suite_definition_id)
    get_suite_definition(suite_definition_id, params::Dict{String,<:Any})

Gets information about a Device Advisor test suite.

# Arguments
- `suite_definition_id`: Requests suite definition Id with GetSuiteDefinition API call.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"suiteDefinitionVersion"`: Requests the suite definition version of a test suite.
"""
get_suite_definition(suiteDefinitionId; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/suiteDefinitions/$(suiteDefinitionId)"; aws_config=aws_config)
get_suite_definition(suiteDefinitionId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/suiteDefinitions/$(suiteDefinitionId)", params; aws_config=aws_config)

"""
    get_suite_run(suite_definition_id, suite_run_id)
    get_suite_run(suite_definition_id, suite_run_id, params::Dict{String,<:Any})

Gets information about a Device Advisor test suite run.

# Arguments
- `suite_definition_id`: Requests the information about Device Advisor test suite run based
  on suite definition Id.
- `suite_run_id`: Requests the information about Device Advisor test suite run based on
  suite run Id.

"""
get_suite_run(suiteDefinitionId, suiteRunId; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/suiteDefinitions/$(suiteDefinitionId)/suiteRuns/$(suiteRunId)"; aws_config=aws_config)
get_suite_run(suiteDefinitionId, suiteRunId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/suiteDefinitions/$(suiteDefinitionId)/suiteRuns/$(suiteRunId)", params; aws_config=aws_config)

"""
    get_suite_run_report(suite_definition_id, suite_run_id)
    get_suite_run_report(suite_definition_id, suite_run_id, params::Dict{String,<:Any})

Gets a report download link for a successful Device Advisor qualifying test suite run.

# Arguments
- `suite_definition_id`: Device Advisor suite definition Id.
- `suite_run_id`: Device Advisor suite run Id.

"""
get_suite_run_report(suiteDefinitionId, suiteRunId; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/suiteDefinitions/$(suiteDefinitionId)/suiteRuns/$(suiteRunId)/report"; aws_config=aws_config)
get_suite_run_report(suiteDefinitionId, suiteRunId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/suiteDefinitions/$(suiteDefinitionId)/suiteRuns/$(suiteRunId)/report", params; aws_config=aws_config)

"""
    list_suite_definitions()
    list_suite_definitions(params::Dict{String,<:Any})

Lists the Device Advisor test suites you have created.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: Request the list of all the Device Advisor test suites.
- `"nextToken"`: Requests the Device Advisor test suites next token.
"""
list_suite_definitions(; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/suiteDefinitions"; aws_config=aws_config)
list_suite_definitions(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/suiteDefinitions", params; aws_config=aws_config)

"""
    list_suite_runs()
    list_suite_runs(params::Dict{String,<:Any})

Lists the runs of the specified Device Advisor test suite. You can list all runs of the
test suite, or the runs of a specific version of the test suite.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: MaxResults for list suite run API request.
- `"nextToken"`: Next pagination token for list suite run request.
- `"suiteDefinitionId"`: Lists the runs of the specified Device Advisor test suite based on
  suite definition Id.
- `"suiteDefinitionVersion"`: Lists the runs of the specified Device Advisor test suite
  based on suite definition version.
"""
list_suite_runs(; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/suiteRuns"; aws_config=aws_config)
list_suite_runs(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/suiteRuns", params; aws_config=aws_config)

"""
    list_tags_for_resource(resource_arn)
    list_tags_for_resource(resource_arn, params::Dict{String,<:Any})

Lists the tags attached to an IoT Device Advisor resource.

# Arguments
- `resource_arn`: The ARN of the IoT Device Advisor resource.

"""
list_tags_for_resource(resourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/tags/$(resourceArn)"; aws_config=aws_config)
list_tags_for_resource(resourceArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/tags/$(resourceArn)", params; aws_config=aws_config)

"""
    list_test_cases()
    list_test_cases(params::Dict{String,<:Any})

Lists all the test cases in the test suite.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"intendedForQualification"`: Lists all the qualification test cases in the test suite.
- `"maxResults"`: Requests the test cases max results.
- `"nextToken"`: Requests the test cases next token.
"""
list_test_cases(; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/testCases"; aws_config=aws_config)
list_test_cases(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("GET", "/testCases", params; aws_config=aws_config)

"""
    start_suite_run(suite_definition_id)
    start_suite_run(suite_definition_id, params::Dict{String,<:Any})

Starts a Device Advisor test suite run.

# Arguments
- `suite_definition_id`: Request to start suite run based on suite definition Id.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"suiteDefinitionVersion"`: Request to start suite run based on suite definition version.
- `"suiteRunConfiguration"`: Request to start suite run based on suite configuration.
- `"tags"`: The tags to be attached to the suite run.
"""
start_suite_run(suiteDefinitionId; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("POST", "/suiteDefinitions/$(suiteDefinitionId)/suiteRuns"; aws_config=aws_config)
start_suite_run(suiteDefinitionId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("POST", "/suiteDefinitions/$(suiteDefinitionId)/suiteRuns", params; aws_config=aws_config)

"""
    tag_resource(resource_arn, tags)
    tag_resource(resource_arn, tags, params::Dict{String,<:Any})

Adds to and modifies existing tags of an IoT Device Advisor resource.

# Arguments
- `resource_arn`: The resource ARN of an IoT Device Advisor resource.
- `tags`: The tags to be attached to the IoT Device Advisor resource.

"""
tag_resource(resourceArn, tags; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("POST", "/tags/$(resourceArn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resourceArn, tags, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("POST", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), params)); aws_config=aws_config)

"""
    untag_resource(resource_arn, tag_keys)
    untag_resource(resource_arn, tag_keys, params::Dict{String,<:Any})

Removes tags from an IoT Device Advisor resource.

# Arguments
- `resource_arn`: The resource ARN of an IoT Device Advisor resource.
- `tag_keys`: List of tag keys to remove from the IoT Device Advisor resource.

"""
untag_resource(resourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("DELETE", "/tags/$(resourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("DELETE", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), params)); aws_config=aws_config)

"""
    update_suite_definition(suite_definition_id)
    update_suite_definition(suite_definition_id, params::Dict{String,<:Any})

Updates a Device Advisor test suite.

# Arguments
- `suite_definition_id`: Updates a Device Advisor test suite with suite definition id.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"suiteDefinitionConfiguration"`: Updates a Device Advisor test suite with suite
  definition configuration.
"""
update_suite_definition(suiteDefinitionId; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("PATCH", "/suiteDefinitions/$(suiteDefinitionId)"; aws_config=aws_config)
update_suite_definition(suiteDefinitionId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = iotdeviceadvisor("PATCH", "/suiteDefinitions/$(suiteDefinitionId)", params; aws_config=aws_config)


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## Test case from Issue #445
##
## Added threaded implementation (2015-04-06)

using Base.Threads

#STOCKCORR - The original, unoptimised code that simulates two correlated assets
function stockcorr()
    ## Correlated asset information
    CurrentPrice = [78. 102.]     # Initial Prices of the two stocks
    Corr = [1. 0.4; 0.4 1.]       # Correlation Matrix
    T = 500                       # Number of days to simulate = 2years = 500days
    n = 1000000                   # Number of simulations
    dt = 1/250                    # Time step (1year = 250days)
    Div=[0.01 0.01]               # Dividend
    Vol=[0.2 0.3]                 # Volatility

    ## Market Information
    r = 0.03                      # Risk-free rate

    ## Define storages
    SimulPriceA = zeros(T,n)      # Simulated Price of Asset A
    SimulPriceA[1,:] = CurrentPrice[1]
    SimulPriceB = zeros(T,n)      # Simulated Price of Asset B
    SimulPriceB[1,:] = CurrentPrice[2]

    ## Generating the paths of stock prices by Geometric Brownian Motion
    UpperTriangle=chol(Corr)      # UpperTriangle Matrix by Cholesky decomposition

    for i = 1:n
       Wiener = randn(T-1,2)
       CorrWiener = Wiener*UpperTriangle
       for j = 2:T
          SimulPriceA[j,i] = SimulPriceA[j-1,i]*exp((r-Div[1]-Vol[1]^2/2)*dt+Vol[1]*sqrt(dt)*CorrWiener[j-1,1])
          SimulPriceB[j,i] = SimulPriceB[j-1,i]*exp((r-Div[2]-Vol[2]^2/2)*dt+Vol[2]*sqrt(dt)*CorrWiener[j-1,2])
       end
    end

    return (SimulPriceA, SimulPriceB)
end

# Run paths in parallel
# NOTE: this has to be in its own function due to #10718
function runpath!(n, Wiener, CorrWiener, SA, SB, T, UpperTriangle, k11, k12, k21, k22, rngs)
    @threads for i = 1:n
    #for i = 1:n
        randn!(rngs[threadid()], Wiener)
        #randn!(rngs[1], Wiener)
        A_mul_B!(CorrWiener, Wiener, UpperTriangle)
        @simd for j = 2:T
            @inbounds SA[j, i] = SA[j-1, i] * exp(k11 + k12*CorrWiener[j-1, 1])
            @inbounds SB[j, i] = SB[j-1, i] * exp(k21 + k22*CorrWiener[j-1, 2])
        end
    end
end

#precompile(runpath!, (Int64, Array{Float64,2}, Array{Float64,2}, Array{Float64,2}, Int64, Array{Float64, 2}, Float64, Float64, Float64, Float64))

# Threaded version
function pstockcorr(n)
    ## Correlated asset information
    const CurrentPrice = [78. 102.]     # Initial Prices of the two stocks
    const Corr = [1. 0.4; 0.4 1.]       # Correlation Matrix
    const T = 500                       # Number of days to simulate = 2years = 500days
    const dt = 1/250                    # Time step (1year = 250days)
    const Div=[0.01 0.01]               # Dividend
    const Vol=[0.2 0.3]                 # Volatility

    ## Market Information
    const r = 0.03                      # Risk-free rate

    ## Define storages
    SimulPriceA = zeros(T,n)            # Simulated Price of Asset A
    SimulPriceA[1,:] = CurrentPrice[1]
    SimulPriceB = zeros(T,n)            # Simulated Price of Asset B
    SimulPriceB[1,:] = CurrentPrice[2]

    ## Generating the paths of stock prices by Geometric Brownian Motion
    const UpperTriangle = full(chol(Corr))    # UpperTriangle Matrix by Cholesky decomposition

    # Optimization: pre-allocate these for performance
    # NOTE: the new GC will hopefully fix this, but currently GC time
    # kills performance if we don't do in-place computations
    Wiener = Array{Float64}(T-1, 2)
    CorrWiener = Array{Float64}(T-1, 2)

    # Runtime requirement: need per-thread RNG since it stores state
    rngs = [MersenneTwister(777+x) for x in 1:nthreads()]

    # Optimization: pre-computable factors
    # NOTE: this should be automatically hoisted out of the loop
    # NOTE: applied this to the MATLAB version as well
    k11 = (r-Div[1]-Vol[1]^2/2)*dt
    k12 = Vol[1]*sqrt(dt)
    k21 = (r-Div[2]-Vol[2]^2/2)*dt
    k22 = Vol[2]*sqrt(dt)

    # NOTE: this reduces reported GC allocation to 576 bytes (due to?)
    #runpath!(1, Wiener, CorrWiener, SimulPriceA, SimulPriceB, T, UpperTriangle, k11, k12, k21, k22, rngs)

    # run
    @time runpath!(n, Wiener, CorrWiener, SimulPriceA, SimulPriceB, T, UpperTriangle, k11, k12, k21, k22, rngs)

    return (SimulPriceA, SimulPriceB)
end

@time pstockcorr(1000000)
#ccall(:jl_threading_profile, Void, ())

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
"""
    AbstractEntry{V,T<:Real}

Supertype of all entry types, which store changes of a specified variable of type `V`
with timestamps of type `T`.
"""
abstract type AbstractEntry{V,T<:Real} end

Base.zero(::Type{E}) where {E<:AbstractEntry} = E()
Base.valtype(::Type{E}) where {V,E<:AbstractEntry{V}} = V

return_type(::Type{E}, v, t) where {E<:AbstractEntry} =
    _return_type(E, _typeof(v), _typeof(t))

_typeof(::T) where {T} = T
_typeof(::Type{T}) where {T} = T
_typeof(::AbstractArray{T}) where {T} = T
_typeof(::AbstractClock{T}) where {T} = T

_return_type(::Type{E}, ::Type{V}, ::Type{T}) where {V,T,E<:AbstractEntry} = E{V,T}
_return_type(::Type{E}, ::Type{V}, ::Type{T}) where {V,T,VE,E<:AbstractEntry{VE}} = E{T}
_return_type(
    ::Type{E},
    ::Type{V},
    ::Type{T},
) where {V,T,VE,TE<:Real,E<:AbstractEntry{VE,TE}} = E

Base.:(==)(e1::AbstractEntry, e2::AbstractEntry) =
    getts(e1) == getts(e2) && getvs(e1) == getvs(e2)

# Show methods
function Base.show(io::IO, ::MIME"text/plain", e::AbstractEntry)
    summary(io, e)
    ts = getts(e)
    isempty(ts) && return nothing
    print(io, " with timestamps:\n ")
    join(io, ts, "\n ")
    return nothing
end

# Recipes
@userplot struct TimeSeries{T}
    es::T
end
@userplot struct PhasePortrait{T}
    es::T
end

@recipe function f(args::TimeSeries)
    es = _unpack(args.es...)
    for e in es
        @series begin
            getts(e), getvs(e)
        end
    end
    return nothing
end
@recipe function f(args::PhasePortrait)
    es = _unpack(args.es...)
    length(es) < 2 && error("Must have at least two entries")
    length(es) > 3 && error("Must have at most three entries")
    ts = copy(getts(es[1]))
    for e in es[2:end]
        append!(ts, getts(e))
    end
    sort!(unique!(ts))
    return Tuple(map(ei -> gettime(ei, ts), es)) # this is not type stable for array arguments
end

_unpack(es::AbstractArray{<:AbstractEntry}) = vec(es)
_unpack(A::DOKSparseArray{<:AbstractEntry}) = collect(values(A.dok))
_unpack(es::AbstractEntry...) = es
# Tools
"""
    store!(e::AbstractEntry, v, t::Union{Real,AbstractClock})

Store a entry with value `v` at time `t`.
"""
store!(e::AbstractEntry, v, c::AbstractClock) = store!(e, v, currenttime(c))

"""
    del!(e::AbstractEntry, t::Union{Real,AbstractClock})

Marks a entry `e` deleted at time `t`.
If an AbstractClock `c` is given, `t = currenttime(c)`.
"""
del!(e::AbstractEntry, c::AbstractClock) = del!(e, currenttime(c))

"""
    getts(e::AbstractEntry{V,T}) -> Vector{T}

Get timestamps of an entry `e`.
"""
getts

"""
    getvs(e::AbstractEntry{V,T}) -> Vector{V}

Get values at each timestamp of an entry `e`.
"""
getvs

"""
    tspan(e::AbstractEntry{V,T}) -> T

Get last time of an entry `e`.
"""
tspan(e::AbstractEntry) = (ts = getts(e); ts[end] - ts[1])

# search algorithm
abstract type AbstractSearch end
struct LinearSearch <: AbstractSearch end
struct BinarySearch <: AbstractSearch end

"""
    gettime([alg::AbstractSearch], e::AbstractEntry, t::Real)
    gettime([alg::AbstractSearch], e::AbstractEntry, ts)
    gettime([alg::AbstractSearch], es::AbstractArray{<:AbstractEntry}, t::Real)
    gettime([alg::AbstractSearch], es::AbstractArray{<:AbstractEntry}, ts)

Get the value of `e` at `t::Real` or values at each time `t` in an iterate `ts`.
If `t` is not a timestamp  in `getts(e)`, return value at time `getts(e)[i]` where
`getts(e)[i] < t < getts(e)[i+1]`. The `alg` is a search algorithm that finds
the index `i` of a target time `t`, which can be `LinearSearch` or `BinarySearch`
(by default).
For `es`, an array of entries, return values at `t` as a similar array.

!!! note

    `ts` must be monotonically increasing.
"""
gettime(e::AbstractEntry, t) = gettime(BinarySearch(), e, t) # t or ts
gettime(alg::AbstractSearch, e::AbstractEntry{V}, ts) where {V} =
    @inbounds gettime!(alg, Vector{V}(undef, length(ts)), e, ts)
function gettime!(alg::AbstractSearch, dst, e::AbstractEntry{V}, ts) where {V}
    Base.@_propagate_inbounds_meta
    state = _initstate(e)
    v = zero(V)
    for (i, t) in enumerate(ts)
        v::V, state = _gettime_itr(alg, e, t, v, state)
        dst[i] = v
    end
    return dst
end
_gettime_itr(::AbstractSearch, ::AbstractEntry, ::Real, v, ::Nothing) = v, nothing

gettime(es::AbstractArray{<:AbstractEntry}, t) = gettime(BinarySearch(), es, t) # t or ts
function gettime(alg::AbstractSearch, es::AbstractArray{<:AbstractEntry}, t::Real)
    ret = similar(es, valtype(eltype(es)))
    @inbounds for i in 1:length(es)
        if isassigned(es, i)
            ret[i] = gettime(alg, es[i], t)
        else
            ret[i] = zero(eltype(ret))
        end
    end
    return ret
end
function gettime(alg::AbstractSearch, es::AbstractArray{<:AbstractEntry}, ts)
    tlen = length(ts)
    ret = Array{valtype(eltype(es))}(undef, tlen, size(es)...)
    @inbounds for i in 1:length(es)
        if isassigned(es, i)
            gettime!(alg, view(ret, (i-1)*tlen+1:i*tlen), es[i], ts)
        else
            ret[(i-1)*tlen+1:i*tlen] .= zero(eltype(ret))
        end
    end
    return ret
end

"""
    DynamicEntry{V,T} <: AbstractEntry{V,T}

Entry type to store changing history of a variable whose value changing overtime.
"""
struct DynamicEntry{V,T<:Real} <: AbstractEntry{V,T}
    vs::Vector{V}
    ts::Vector{T}
    function DynamicEntry{V,T}(v::V, t::T) where {V,T}
        return new{V,T}([v], [t])
    end
    function DynamicEntry{V,T}() where {V,T}
        return new{V,T}(V[], T[])
    end
end
DynamicEntry{V,T}(v, t) where {V,T} = DynamicEntry{V,T}(convert(V, v), convert(T, t))
DynamicEntry{V}(v, t::T) where {V,T} = DynamicEntry{V,T}(v, t)
DynamicEntry(v::V, t::T) where {V,T} = DynamicEntry{V,T}(v, t)
DynamicEntry{V,T}(v, c::C) where {V,T,C<:AbstractClock} =
    DynamicEntry{V,T}(v, currenttime(c))
DynamicEntry{V}(v, c::AbstractClock{T}) where {V,T} = DynamicEntry{V,T}(v, currenttime(c))
DynamicEntry(v::V, c::AbstractClock{T}) where {V,T} = DynamicEntry{V,T}(v, currenttime(c))

store!(e::DynamicEntry, v, t::Real) = (push!(getvs(e), v); push!(getts(e), t); e)
del!(e::DynamicEntry, ::Real) = e

getts(e::DynamicEntry) = e.ts
getvs(e::DynamicEntry) = e.vs

function gettime(::LinearSearch, e::DynamicEntry, t::Real)
    ts = getts(e)
    vs = getvs(e)
    @inbounds ts[1] > t && return zero(valtype(e))
    @inbounds for i in eachindex(ts)
        ts[i] > t && return vs[i-1]
    end
    return vs[end]
end

function gettime(::BinarySearch, e::DynamicEntry{V}, t::Real) where {V}
    match = searchsortedlast(getts(e), t)
    return match == 0 ? zero(V) : getvs(e)[match]
end

_initstate(e::DynamicEntry) = 1, length(getvs(e))

function _gettime_itr(
    ::LinearSearch,
    e::DynamicEntry{V},
    t::Real,
    ::V,
    state::Tuple{Int,Int},
) where {V}
    ts = getts(e)
    vs = getvs(e)
    lo, hi = state
    lo == 1 && ts[1] > t && return zero(V), state
    @inbounds for i in lo:hi
        ts[i] > t && return vs[i-1], (i - 1, hi)
    end
    return vs[hi], nothing
end
function _gettime_itr(
    ::BinarySearch,
    e::DynamicEntry{V},
    t::Real,
    ::V,
    state::Tuple{Int,Int},
) where {V}
    ts = getts(e)
    vs = getvs(e)
    lo, hi = state
    m = searchsortedlast(ts, t, lo, hi, Base.Order.Forward)
    if m + 1 == lo
        return zero(V), state
    elseif m == hi
        return vs[hi], nothing
    else
        return vs[m], (m, hi)
    end
end

"""
    StaticEntry{V,T} <: AbstractEntry{V,T}

Entry type to store changing history of a variable whose value not changing overtime.
"""
mutable struct StaticEntry{V,T<:Real} <: AbstractEntry{V,T}
    init::Bool         # initialized or not
    v::V               # value
    s::T               # start time
    delete::Bool       # deleted or not
    e::T               # end time
    function StaticEntry{V,T}(v::V, t::T) where {V,T}
        return new{V,T}(true, v, t, false) # e.e is not assigned
    end
    function StaticEntry{V,T}() where {V,T}
        return new{V,T}(false) # most of fields are not assigned
    end
end
StaticEntry{V,T}(v, t) where {V,T} = StaticEntry{V,T}(convert(V, v), convert(T, t))
StaticEntry{V}(v, t::T) where {V,T} = StaticEntry{V,T}(v, t)
StaticEntry(v::V, t::T) where {V,T} = StaticEntry{V,T}(v, t)
StaticEntry{V,T}(v, c::C) where {V,T,C<:AbstractClock} = StaticEntry{V,T}(v, currenttime(c))
StaticEntry{V}(v, c::AbstractClock{T}) where {V,T} = StaticEntry{V,T}(v, currenttime(c))
StaticEntry(v::V, c::AbstractClock{T}) where {V,T} = StaticEntry{V,T}(v, currenttime(c))

function store!(e::StaticEntry, v, t::Real)
    e.init && error("the StaticEntry has be initialized")
    e.init = true
    e.v = v
    e.s = t
    e.delete = false
    return e
end
function del!(e::StaticEntry, t::Real)
    e.init && e.delete && error("the StaticEntry has be deleted")
    e.delete = true
    e.e = t
    return e
end

getts(e::StaticEntry) = e.delete ? [e.s, e.e] : [e.s]
getvs(e::StaticEntry) = e.delete ? [e.v, e.v] : [e.v]

gettime(::AbstractSearch, e::StaticEntry{V}, t::Real) where {V} =
    t < e.s ? zero(V) : !e.delete ? e.v : t <= e.e ? e.v : zero(V)

_initstate(::StaticEntry) = true

function _gettime_itr(
    ::AbstractSearch,
    e::StaticEntry{V},
    t::Real,
    ::V,
    state::Bool,
) where {V}
    state && e.s > t && return zero(V), true
    e.delete || return e.v, nothing
    if e.e >= t
        return e.v, false
    else
        return zero(V), nothing
    end
end


text:
module ContinuumArrays
using IntervalSets, LinearAlgebra, LazyArrays, FillArrays, BandedMatrices, QuasiArrays, Infinities, InfiniteArrays, StaticArrays, BlockArrays, RecipesBase
import Base: @_inline_meta, @_propagate_inbounds_meta, axes, size, getindex, convert, prod, *, /, \, +, -, ==, ^,
                IndexStyle, IndexLinear, ==, OneTo, _maybetail, tail, similar, copyto!, copy, diff,
                first, last, show, isempty, findfirst, findlast, findall, Slice, union, minimum, maximum, sum, _sum,
                getproperty, isone, iszero, zero, abs, <, ≤, >, ≥, string, summary, to_indices, view
import Base.Broadcast: materialize, BroadcastStyle, broadcasted, Broadcasted
import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupport, most, combine_mul_styles, AbstractArrayApplyStyle,
                        adjointlayout, arguments, _mul_arguments, call, broadcastlayout, layout_getindex, UnknownLayout,
                        sublayout, sub_materialize, ApplyLayout, BroadcastLayout, combine_mul_styles, applylayout,
                        simplifiable, _simplify, AbstractLazyLayout, PaddedLayout
import LinearAlgebra: pinv, dot, norm2, ldiv!, mul!
import BandedMatrices: AbstractBandedLayout, _BandedMatrix
import BlockArrays: block, blockindex, unblock, blockedrange, _BlockedUnitRange, _BlockArray
import FillArrays: AbstractFill, getindex_value, SquareEye
import ArrayLayouts: mul, ZerosLayout, ScalarLayout, AbstractStridedLayout
import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclusion, SubQuasiArray,
                    QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat, QuasiArrayLayout,
                    ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle, AbstractQuasiLazyLayout,
                    LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle, _factorize,
                    AbstractQuasiFill, UnionDomain, __sum, _cumsum, __cumsum, applylayout, _equals, layout_broadcasted, PolynomialLayout
import InfiniteArrays: Infinity, InfAxes

export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, ℵ₁, Inclusion, Basis, grid, plotgrid, affine, ..

const QMul2{A,B} = Mul{<:Any, <:Any, <:A,<:B}
const QMul3{A,B,C} = Mul{<:AbstractQuasiArrayApplyStyle, <:Tuple{A,B,C}}

cardinality(::AbstractInterval) = ℵ₁

Inclusion(d::AbstractInterval{T}) where T = Inclusion{float(T)}(d)
first(S::Inclusion{<:Any,<:AbstractInterval}) = leftendpoint(S.domain)
last(S::Inclusion{<:Any,<:AbstractInterval}) = rightendpoint(S.domain)

norm2(x::Inclusion{T,<:AbstractInterval}) where T = sqrt(dot(x,x))

function dot(x::Inclusion{T,<:AbstractInterval}, y::Inclusion{V,<:AbstractInterval}) where {T,V}
    x == y || throw(DimensionMismatch("first quasivector has axis $(x) which does not match the axis of the second, $(y)."))
    TV = promote_type(T,V)
    isempty(x) && return zero(TV)
    a,b = endpoints(x.domain)
    convert(TV, b^3 - a^3)/3
end

include("maps.jl")

const QInfAxes = Union{Inclusion,AbstractAffineQuasiVector}

sub_materialize(_, V::AbstractQuasiArray, ::Tuple{QInfAxes}) = V
sub_materialize(_, V::AbstractQuasiArray, ::Tuple{QInfAxes,QInfAxes}) = V
sub_materialize(_, V::AbstractQuasiArray, ::Tuple{Any,QInfAxes}) = V
sub_materialize(_, V::AbstractQuasiArray, ::Tuple{QInfAxes,Any}) = V

# ambiguity error
sub_materialize(_, V::AbstractQuasiArray, ::Tuple{InfAxes,QInfAxes}) = V
sub_materialize(_, V::AbstractQuasiArray, ::Tuple{QInfAxes,InfAxes}) = V
sub_materialize(::ApplyLayout{typeof(hcat)}, V::AbstractQuasiArray, ::Tuple{QInfAxes,Any}) = V
sub_materialize(::ApplyLayout{typeof(hcat)}, V::AbstractQuasiArray, ::Tuple{QInfAxes,InfAxes}) = V

#
# BlockQuasiArrays

BlockArrays.blockaxes(::Inclusion) = blockaxes(Base.OneTo(1)) # just use 1 block
function BlockArrays.blockaxes(A::AbstractQuasiArray{T,N}, d) where {T,N}
    @_inline_meta
    d::Integer <= N ? blockaxes(A)[d] : Base.OneTo(1)
end

@inline to_indices(A::AbstractQuasiArray, inds, I::Tuple{Block{1}, Vararg{Any}}) =
    (unblock(A, inds, I), to_indices(A, _maybetail(inds), tail(I))...)
@inline to_indices(A::AbstractQuasiArray, inds, I::Tuple{BlockRange{1,R}, Vararg{Any}}) where R =
    (unblock(A, inds, I), to_indices(A, _maybetail(inds), tail(I))...)
@inline to_indices(A::AbstractQuasiArray, inds, I::Tuple{BlockIndex{1}, Vararg{Any}}) =
    (inds[1][I[1]], to_indices(A, _maybetail(inds), tail(I))...)
@inline to_indices(A::AbstractQuasiArray, I::Tuple{BlockRange, Vararg{Any}}) = to_indices(A, axes(A), I)

@inline to_indices(A::AbstractQuasiArray, inds, I::Tuple{AbstractArray{<:BlockIndex{1}}, Vararg{Any}}) =
    (inds[1][I[1]], to_indices(A, _maybetail(inds), tail(I))...)

checkpoints(x::Number) = x
checkpoints(d::AbstractInterval{T}) where T = width(d) .* SVector{3,float(T)}(0.823972,0.01,0.3273484) .+ leftendpoint(d)
checkpoints(d::UnionDomain) = mapreduce(checkpoints,union,d.domains)
checkpoints(x::Inclusion) = checkpoints(x.domain)
checkpoints(A::AbstractQuasiMatrix) = checkpoints(axes(A,1))

include("operators.jl")
include("bases/bases.jl")
include("basisconcat.jl")

include("plotting.jl")

end


text:
using Pkg;Pkg.activate(".");Pkg.instantiate();

using DifferentialEquations,ComponentArrays
using LinearAlgebra
using DelimitedFiles,CSV,DataFrames
using Parameters
using Distributed
using Printf

# Load the file that contains the Polyharmonic interpolator function
include("polyHarmonicInterp.jl")

"""

	plodyModel()

	This function is used in DifferentialEquations and returns the RHS of the
	differential equation needed to solve ploidy evolution of the discrete
	system.

	INPUTS
		du:		A vector of the RHS of y' = du
		u:		A vector of the solution variable
		p:		A list of parameters that we pass
		t:		The current time point
	OUTPUTS
		du:		A vector of the RHS of y' = du

"""
function ploidy_wellmixed_model(du,u,pars,t)

	# Grab the parameters
	@unpack (focal_birthRate,nChrom,minChrom,maxChrom,stepChrom,misRate,deathRate) = pars

	# Iterate over each compartment
	for i in CartesianIndices(u)

		#=
		1) calculateParents(focal, minChrom, maxChrom, stepChrom)
		2) updatePopulation (a struct?)
			sum(inflow(focal,parent,u[parent]) for parent in parentList)
			- outflow(focal,u[focal])
			a. inflow needs to calculate birth and flow rate
			b. outflow needs to calculate birth, death and flow rate
		=#

		focalCN =  collect(i[k] for k in 1:nChrom);

		parentCNList = calculateParents(focalCN, minChrom,maxChrom,stepChrom)

		inflow = 0.0

		# Inflow from the parents
		for parentCN in parentCNList

			# Get birth rate from precomputed
			birthRate = focal_birthRate[Int.(parentCN)...]

			# Get flow rate from parentCN -> focalCN
			flowRate = q(parentCN,focalCN,misRate,nChrom)

			# Will need to have this be the index in the future
			inflow += birthRate*flowRate*u[Int.(parentCN)...] # *(1 - 1/avg(CN))

		end

		# Grab the focal cells birth rate
		birthRate = focal_birthRate[i]

		# Add it to the inflow to the focal compartment
		inflow += birthRate*(1.0 - 2*misRate)*u[i]

		# Death of the focal compartment
		outflow = deathRate*u[i]

		# Update the RHS
		du[i] = inflow - outflow

	end

end

function ploidy_spatial_model(du,u,pars,t)

	# Grab the parameters
	@unpack (focal_birthRate,nChrom,chromArray,misRate,deathRate,Γ,Γₑ,ϕ,Ξ,χ,δ,Np,dp,k,E_vessel,
	domain_Dict,boundary_Dict,max_cell_cycle_duration,max_birthRate,debugging) = pars

	# FIX ME:: Assume that all chroms have same lower, upper and step size
	minChrom,stepsize,maxChrom = Int.([chromArray[1] |> f for f in (first,step,length)])

	# Convert chrom array to integer
	intChromArray = [Int(first(x)):Int(step(x)):Int(last(x)) for x in chromArray]
	s, E = u.s, u.E

	coordsList = [1:num_points for num_points in Np]
	# dp = (Np .- 1).^(-1)
	sum_dp_invsq = sum(dp.^(-2))

	# cflow = 0.0
	# dflow = 0.0

	# Iterate over each coordinate
	for coord in Iterators.product((1:length(coords) for coords in coordsList)...)

		# Get whether we are interior (0), at boundary of vessel (1) or inside vessel (2)
		domain_info = get(domain_Dict,coord,0)

		if domain_info == 2 # Inside blood vessel we make no changes
			du.s[intChromArray...,coord...] .= 0.0
			du.E[coord...] = 0.0
			continue
		end

		# Gets indicies used to calculate spatial impact (assume periodic)
		cardinal_point_info = get_cardinal_gridpoints([elem for elem in coord],Np)

		# Gets the points to be used for spatial grid resolution
		E_cardinal = [E[info...] for info in cardinal_point_info]
		E_minus,E_plus = E_cardinal[end-1:-2:1],E_cardinal[end:-2:1]
		E_focal = E[coord...]

		if domain_info == 1	# boundary of vessel
			normal_vector = boundary_Dict[coord] 	# Grab normal vector
			for i = 1 : length(dp)
				if normal_vector[i] < 0 # The rightward point is in the vessel so we replace it with the B.C.
					E_plus[i] = E_minus[i] - 2.0*dp[i]*k*(E_focal - E_vessel)
				elseif normal_vector[i] > 0 # The leftward point is in the vessel
					E_minus[i] =  E_plus[i] - 2.0*dp[i]*k*(E_focal - E_vessel)
				end
			end
		end

		# Iterate over each compartment
		for focal in Iterators.product((1:length(chr) for chr in chromArray)...)

			# focal cells copy number state
			focalCN = collect(chromArray[k][focal[k]] for k in 1:nChrom)

			# All possible chromosome states
			parentCNList = calculateParents(focalCN,minChrom,maxChrom,stepsize)

			# Initialize inflow to 0
			inflow = 0.0

			# copy state to local variables for brevity
			s_focal = s[focal...,coord...]
			s_cardinal = [s[focal...,info...] for info in cardinal_point_info]
			s_minus,s_plus = s_cardinal[end-1:-2:1],s_cardinal[end:-2:1]

			# Checking the boundary
			if domain_info == 1	# boundary of vessel
				normal_vector = boundary_Dict[coord] 	# Grab normal vector
				for i = 1 : length(dp)
					if normal_vector[i] < 0 # The rightward point is in the vessel so we replace it with the B.C.
						s_plus[i] = s_minus[i] + χ/Γ*s_focal*(chemotaxis_form(E_plus[i],Ξ) -
						chemotaxis_form(E_minus[i],Ξ))
					elseif normal_vector[i] > 0 # The leftward point is in the vessel
						s_minus[i] = s_plus[i] - χ/Γ*s_focal*(chemotaxis_form(E_plus[i],Ξ) -
						chemotaxis_form(E_minus[i],Ξ))
					end
				end
			end

			energy_constraint = energy_constrained_birth(E_focal,ϕ)

			if max_birthRate*energy_constraint > log(2)/max_cell_cycle_duration
				# Inflow from the parents
				for parentCN in parentCNList

					s_parent = s[Int.(parentCN)...,coord...]

					# Get max birth rate
					birthRate = focal_birthRate[Int.(parentCN)...]

					if s_parent == 0.0 || birthRate == 0.0
						continue
					end

					# If birth rate is below a threshold determined by amount spent in G1 (default is 100 hours)
					# then we assume that birth does not occur...
					energy_modified_birth_rate = birthRate*energy_constraint
					if energy_modified_birth_rate > log(2)/max_cell_cycle_duration
						# Get flow rate from parentCN -> focalCN
						flowRate = q(parentCN,focalCN,misRate,nChrom)
						inflow += energy_modified_birth_rate*flowRate*s_parent
					else
						# println("At energy level $(E_focal), birth rate of $(parentCN) is too low...")
					end

				end

				if s_focal > 0
					# Grab the focal cells (max?) birth rate
					birthRate = focal_birthRate[Int.(focalCN)...]
					energy_modified_birth_rate = birthRate*energy_constraint

					if energy_modified_birth_rate > log(2)/max_cell_cycle_duration
						# Add it to the inflow to the focal compartment
						inflow += energy_modified_birth_rate*(1.0 - 2.0*misRate)*s_focal
					else
						# println("At energy level $(E_focal), birth rate of $(focalCN) is too low...")
					end
				end

			end

			# Death of the focal compartment
			outflow = deathRate*s_focal

			# u_xx -> (u[i+1] + u[i-1] - 2u[i])/dx^2
			# random migration
			diffusion = Γ*(sum((s_plus[i] + s_minus[i])/dp[i]^2 for i in 1:length(dp))
			- 2.0*sum_dp_invsq*s_focal)

			# chemotaxis
			chemotaxis = χ*(
				s_focal*(
				sum(
					(chemotaxis_form(E_plus[i],Ξ) + chemotaxis_form(E_minus[i],Ξ))/dp[i]^2
					for i in 1 : length(dp)
					)
				- 2.0*sum_dp_invsq*chemotaxis_form(E_focal,Ξ)) +
				sum(
					(s_plus[i] - s_minus[i])*(chemotaxis_form(E_plus[i],Ξ) - chemotaxis_form(E_minus[i],Ξ))/dp[i]^2
					for i in 1 : length(dp)))/4.0

			# Update the RHS dependent on where we are in the domain
			du.s[focal...,coord...] = inflow - outflow + diffusion - chemotaxis

			# dflow += diffusion
			# cflow += chemotaxis

		end

		# Diffusion of energy molecule through the tissue
		diffusion = Γₑ*(sum((E_plus[i] + E_minus[i])/dp[i]^2 for i in 1:length(dp)) - 2.0*sum_dp_invsq*E_focal)

		# Consumption by the cells at the grid point designated by coord.
		consumption = δ*E_focal*sum(s[intChromArray...,coord...])

		# update the energy compartment
		du.E[coord...] = diffusion - consumption

		# dflow += diffusion

	end

	# @show dflow,cflow

end

function get_cardinal_gridpoints(coordinate::Vector{Int},Np::Vector{Int})

	@assert length(coordinate) == length(Np) "Dimensions must match"
	dim = length(coordinate)

	[ (tmp = copy(coordinate); tmp[idx] = mod(v-1,Np[idx])+1; tmp)
		for idx in 1:dim for v in coordinate[idx]-1:coordinate[idx]+1
			if v != coordinate[idx]]
end

function calculateParents(offspring::Vector{T}, minChrom::Int,
	maxChrom::Int,stepChrom::Int) where T<:Real

	[ (x = copy(offspring); x[idx] = v; x)
		for idx in 1 : length(offspring)
			for v in max(offspring[idx]-1,minChrom):stepChrom:min(offspring[idx]+1,maxChrom)
				if v != offspring[idx] ]

end

function energy_constrained_birth(energy::Number,ϕ::Number)

	return energy/(ϕ + energy)

end

function chemotaxis_form(energy::Number,Ξ::Number)

	return log(Ξ + energy)

end

function energy_consumption(energy::Number,δ::Number,u::Vector{T}) where T <: Number

	return δ*energy*sum(u)

end

function q(parent::Vector{T},offspring::Vector{T},misRate::Float64,nChrom::Int) where T <: Real

	# parentCN = collect(linspaces[k][parent[k]] for k in 1:nChrom)
	# offspringCN = collect(linspaces[k][offspring[k]] for k in 1:nChrom)

	parentCN = parent
	offspringCN = offspring

	CNdist = norm(parentCN .- offspringCN,1)

	# No mutation occurred
	if (CNdist == 0.0)
		val = 1.0 - misRate
	else
		# FIXME: Should be a function dependent on parent ploidy
		val = misRate/nChrom
	end

	return val

end

function simulate_wellmixed(params,X::AbstractArray,Y::AbstractVector)

	# Grab options that determine output and other parts outside the model
	@unpack (
			debugging,
			compartmentMinimum,
			replating,
			maxPop,
			progress_check,
			saveat,
			interpolation_order,
			) = params.Options

	# Grab the model parameters from the struct
	@unpack (
			stepsize,
			minChrom,
			maxChrom,
			deathRate,
			misRate,
			finalDay,
			startPop,
			starting_copy_number,
			max_cell_cycle_duration) = params.UniversalParameters

	# Polyharmonic interpolator
	interp = PolyharmonicInterpolation.PolyharmonicInterpolator(X,Y,interpolation_order)

	# Number of chromosomes (inferred from interp)
	nChrom = interp.dim

	# The allowable states
	chromArray = [minChrom:stepsize:maxChrom for _ in 1:nChrom]

	# Total number of states
	nComp = prod(length(arr) for arr in chromArray)

	# Time interval to run simulation and initial condition
	tspan = (0.0,finalDay)
	u0 = zeros(Int(maxChrom)*ones(Int,nChrom)...)
	u0[starting_copy_number...] = startPop

	##compute birth rates associated with each copy number state
	focal_birthRate = zeros(Int(maxChrom)*ones(Int,nChrom)...);
	for focal in Iterators.product((1:length(j) for j in chromArray)...)
		focalCN = collect(chromArray[k][focal[k]] for k in 1:nChrom);
		focal_birthRate[focal...] = max(PolyharmonicInterpolation.polyharmonicSpline(interp,focalCN')[1],0.0);
	end

	# If we are replating, we do so when the population is million-fold in size
	cb1 = nothing
	if replating
		# Replate if above maxPop
		condition = (u,t,integrator) -> sum(u) - maxPop

		# This replates so that u sums to 1 (as in the initial condition)
		affect!(integrator) = begin
			if debugging > 1
				println("Max population reached at t = $(integrator.t)...replating.")
			end
			integrator.u /= maxPop
		end

		# Create callback
		cb1 = ContinuousCallback(condition,affect!,nothing)

	end

	##### This is printed to terminal to tell you what time point you're at #####
	cb2 = nothing
	if progress_check
		print_time(integrator) = @printf("t = %.2f\n",integrator.t)
		cb2 = PeriodicCallback(print_time, 0.05,save_positions=(false,false))
	end
	cbset = CallbackSet(cb1,cb2)

	if debugging > 0
		println("Beginning simulation...")
	end

	# run simulation
	odePars = (focal_birthRate=focal_birthRate,
				nChrom=nChrom,
				minChrom = chromArray[1].offset,
				maxChrom = chromArray[1].len,
				stepChrom=Int(chromArray[1].step),
				misRate=misRate,
				deathRate=deathRate,
				)
	prob = ODEProblem(ploidy_wellmixed_model,u0,tspan,odePars)
	sol = solve(prob,Tsit5(),maxiters=1e5,abstol=1e-8,reltol=1e-5,saveat=saveat,callback=cbset)

	if debugging > 0
		println("Simulation complete. Collecting results...")
	end

	# convert to array for output
	soln = reshape(Array(sol),nComp,length(sol.t))

	# Grab compartment -> CN array
	cnArray = Array{Float64}(undef,nComp,nChrom)
	for (i,focal) in enumerate(
		Iterators.product((1:length(j) for j in chromArray)...))
			cnArray[i,:] = collect(chromArray[k][focal[k]] for k in 1:nChrom)
	end

	return soln,cnArray,sol.t

end

function simulate_spatial(params,X::AbstractArray,Y::AbstractVector)

	# Grab options that determine output and other parts outside the model
	@unpack (
			debugging,
			compartmentMinimum,
			replating,
			maxPop,
			progress_check,
			saveat,
			interpolation_order,
			) = params.Options

	# Grab the model parameters from the struct
	@unpack (
			stepsize,
			minChrom,
			maxChrom,
			deathRate,
			misRate,
			finalDay,
			startPop,
			starting_copy_number,
			max_cell_cycle_duration) = params.UniversalParameters

	# Grab the spatial parameters from the struct
	@unpack (
			Γ,
			Γₑ,
			ϕ,
			Ξ,
			χ,
			δ,
			Lp,
			Np,
			k,
			E_vessel) = params.SpatialParameters

	# convert max cell cycle to days (the time scale used)
	max_cell_cycle_duration/= 24.0

	#= Build domain matrix:

	0 : In domain
	1 : Boundary of blood vessel
	2 : Interior of blood vessel

	=#
	df = CSV.read(params.blood_vessel_file,DataFrame)
	# maxX,maxY=maximum(df[!,"Centroid X µm"]),maximum(df[!,"Centroid Y µm"])
	# minX,minY=minimum(df[!,"Centroid X µm"]),minimum(df[!,"Centroid Y µm"])
	# df[!,"Centroid X µm"] .= (df[!,"Centroid X µm"] .- minX)./(maxX .- minX)
	# df[!,"Centroid Y µm"] .= (df[!,"Centroid Y µm"] .- minY)./(maxY .- minY)

	x = range(0,Lp[1],length=Np[1])
	y = range(0,Lp[2],length=Np[2])
	dp = Lp./(Np .- 1)

	# df[!,"Centroid X µm"] .= x[searchsortedfirst.(Ref(x),df[!,"Centroid X µm"])]
	# df[!,"Centroid Y µm"] .= y[searchsortedfirst.(Ref(y),df[!,"Centroid Y µm"])]
	# df = combine(groupby(df,["Centroid X µm","Centroid Y µm"]),last)

	domain_Dict = Dict{Tuple{Int,Int},Int}() # 1 = boundary blood vessel, 2 = inside blood vessel

	for i = 1 : Np[1]
		xdist = (x[i] .- df[!,"Centroid X µm"]).^2
		for j = 1 : Np[2]
			ydist = (y[j] .- df[!,"Centroid Y µm"]).^2
			if minimum(xdist .+ ydist) < sum(dp.^2)
				domain_Dict[(i,j)]=1
			end
		end
	end

	# Fill dict with info about the normal vector
	boundary_Dict = Dict{Tuple{Int,Int},Any}()

	# But first we find whether the point is interior
	for k in keys(domain_Dict)
		i,j = k
		i_plus = (i == Np[1]) ? 1 : i+1
		i_minus = (i == 1) ? Np[1] : i-1
		j_plus = (j == Np[2]) ? 1 : j+1
		j_minus = (j == 1) ? Np[2] : j-1

		# Vessel interior point
		if issubset([(i,j_plus),(i,j_minus),(i_plus,j),(i_minus,j)],keys(domain_Dict))
			domain_Dict[k] = 2
		end
	end

	# If no cardinal direction has an interior point we will delete it from the dict
	for (k,v) in domain_Dict

		# We are at an interior point
		if v == 1
			i,j = k[1],k[2]
			i_plus = (i == Np[1]) ? 1 : i+1
			i_minus = (i == 1) ? Np[1] : i-1
			j_plus = (j == Np[2]) ? 1 : j+1
			j_minus = (j == 1) ? Np[2] : j-1

			cardinal_directions = [(i,j_plus),(i,j_minus),(i_plus,j),(i_minus,j)]

			if any(z->z==2,get(domain_Dict,k,0) for k in cardinal_directions)
				# Vessel boundary point
				normal = [0;0]
				normal[1] += get(domain_Dict,(i,j_plus),0) == 2 ? -1 : 0
				normal[1] += get(domain_Dict,(i,j_minus),0) == 2 ? 1 : 0
				normal[2] += get(domain_Dict,(i_plus,j),0) == 2 ? -1 : 0
				normal[2] += get(domain_Dict,(i_minus,j),0) == 2 ? 1 : 0
				# normal/=norm(normal)
				boundary_Dict[k] = normal
			else
				delete!(domain_Dict,k)
			end
		end
	end

	# Initialize energy to be zero outside the interior of the blood vessels
	E0 = zeros(Np...)
	[E0[k...]=E_vessel for (k,v) in domain_Dict if v==2]

	# Polyharmonic interpolator
	interp = PolyharmonicInterpolation.PolyharmonicInterpolator(X,Y,interpolation_order)

	# Number of chromosomes (inferred from interp)
	nChrom = interp.dim

	# The allowable states
	chromArray = [minChrom:stepsize:maxChrom for i in 1:nChrom]

	# Total number of copy number states
	nComp = prod(length(arr) for arr in chromArray)

	# Time interval to run simulation and initial condition
	tspan = (0.0,finalDay)
	s0 = zeros(Int(maxChrom)*ones(Int,nChrom)...,Np...)
	s0[starting_copy_number...,:,:] .= rand(Np...)*startPop 	# CN is randomly placed on the grid
	[s0[starting_copy_number...,k...] = 0 for k in keys(domain_Dict)]

	u0 = ComponentArray(s=s0,E=E0)

	# If we are replating, we do so when the population is million-fold in size
	cb1 = nothing
	if replating
		# Replate if above maxPop
		condition = (u,t,integrator) -> sum(u) - maxPop

		# This replates so that u sums to 1 (as in the initial condition)
		affect!(integrator) = begin
			if debugging > 1
				println("Max population reached at t = $(integrator.t)...replating.")
			end
			integrator.u /= maxPop
		end

		# Create callback
		cb1 = ContinuousCallback(condition,affect!,nothing)

	end

	##### This is printed to terminal to tell you what time point you're at ######
	##### Add the smallest nonzero compartment size (and index) to the progress check
	cb2 = nothing
	if progress_check
		print_time(integrator) = @printf("t = %.2f\n",integrator.t)
		cb2 = PeriodicCallback(print_time, 0.1,save_positions=(false,false))
	end
	cbset = CallbackSet(cb1,cb2)

	if debugging > 0
		println("Beginning simulation...")
	end

	##compute birth rates associated with each copy number state
	focal_birthRate = zeros(Int(maxChrom)*ones(Int,nChrom)...);
	for focal in Iterators.product((1:length(j) for j in chromArray)...)
		focalCN = collect(chromArray[k][focal[k]] for k in 1:nChrom);
		focal_birthRate[focal...] = max(PolyharmonicInterpolation.polyharmonicSpline(interp,focalCN')[1],0.0);
		println("CN = $focalCN, birth rate = $(focal_birthRate[focal...])")
	end

	max_birthRate = maximum(focal_birthRate)

	# run simulation
	# isoutofdomain = (u,p,t)->any(x->x<0,u)
	odePars = (focal_birthRate=focal_birthRate,
				nChrom=nChrom,
				chromArray=chromArray,
				misRate=misRate,
				deathRate=deathRate,
				Γ=Γ,Γₑ=Γₑ,ϕ=ϕ,Ξ=Ξ,χ=χ,δ=δ,Np=Np,dp=dp,k=k,E_vessel=E_vessel,
				domain_Dict=domain_Dict,
				boundary_Dict=boundary_Dict,
				max_cell_cycle_duration=max_cell_cycle_duration,
				max_birthRate=max_birthRate,
				debugging=debugging)
	prob = ODEProblem(ploidy_spatial_model,u0,tspan,odePars)
	sol = solve(prob,VCABM(),abstol=1e-8,reltol=1e-5,saveat=saveat,callback=cbset)

	# Maybe lsoda() instead of VCABM() ????
	#=

	t = range(tspan..., length=1000) [u.s for u in sol(t)]
	[u.s for u in sol.u]
	=#

	if debugging > 0
		println("Simulation complete. Collecting results...")
	end

	# Grab copy number state array -> CN array
	cnArray = Array{Int}(undef,nComp,nChrom)
	for (i,focal) in enumerate(
		Iterators.product((1:length(j) for j in chromArray)...))
			cnArray[i,:] = collect(chromArray[k][focal[k]] for k in 1:nChrom)
	end

	return sol, cnArray

end

#############################################
#											#
#											#
#					Main		 			#
#											#
#											#
#############################################

function runPloidyMovement(params,X::AbstractArray,Y::AbstractVector,spatial::Bool)

	if spatial
		println("Running the spatial model...")
		return simulate_spatial(params,X,Y)
	else
		println("Running the well-mixed model...")
		return simulate_wellmixed(params,X,Y)
	end

end

max_stars_repo_name: MathOnco/ploidyEvolution

max_issues_repo_name: MathOnco/ploidyEvolution

max_forks_repo_name: MathOnco/ploidyEvolution


text:
function simu_dat_gen(
  ;
  n_features=20,
  n_samples=10,
  n_sig_features=5,
  sig_intensity=10,
  noise_intensity=1,
  features_as_cols=true,
)
  if n_sig_features > n_features
    error("n_sig_features ($n_sig_features) > n_features ($n_features)")
  end

  mat = randn(n_features, n_samples)

  if features_as_cols
    return mat
  else
    return mat'
  end
end

# test

print("Hello!")

max_stars_repo_name: granatumx/gbox-py

max_issues_repo_name: lanagarmire/granatumx

max_forks_repo_name: lanagarmire/granatumx


text:
#    myforeachfield((col, val) -> (@inbounds col[I] = val), s, vals)

using StructArrays
import StructArrays: staticschema, map_params, astuple, _map_params
using Base: tuple_type_cons, tuple_type_head, tuple_type_tail, tail

function sstuple(::Type{<:NTuple{N, Any}}) where {N}
    ntuple(j->Symbol(j), N)
end

function sstuple(::Type{NT}) where {NT<:NamedTuple}
    _map_params(x->sstuple(staticschema(x)), NT)
end

function f(s, es, exprs)
    if typeof(s) <: Symbol
        push!(exprs, es)
        return
    end
    for key in keys(s)
        f(getproperty(s,key), vcat(es, key), exprs)
    end
end

using StaticArrays
c = [(a=SHermitianCompact(@SVector(rand(3))), b=(@SVector(rand(2)))) for i=1:5]
d = d = StructArray(c, unwrap = t -> t <: Union{SHermitianCompact,SVector,Tuple})
a = []
b = []
S = StructArrays.staticschema(typeof(d))
s = sstuple(S)
f(s, a, b)
@show b


text:
using Test
using Distributions
using StatsBase
using Random

using EnsembleKalmanProcesses.ParameterDistributionStorage

@testset "ParameterDistributionStorage" begin
    @testset "ParameterDistributionType" begin
        # Tests for the ParameterDistributionType
        d = Parameterized(Gamma(2.0, 0.8))
        @test d.distribution == Gamma(2.0,0.8)

        d = Samples([1 2 3; 4 5 6])
        @test d.distribution_samples == [1 2 3; 4 5 6]
        d = Samples([1 4; 2 5; 3 6]; params_are_columns=false)
        @test d.distribution_samples == [1 2 3; 4 5 6]
        d = Samples([1, 2, 3, 4, 5, 6])
        @test d.distribution_samples == [1 2 3 4 5 6]
        d = Samples([1, 2, 3, 4, 5, 6]; params_are_columns=false)
        @test d.distribution_samples == reshape([1, 2, 3, 4, 5, 6],:,1)

    end
    @testset "ConstraintType" begin
        # Tests for the ConstraintType
        # The predefined transforms
        c1 = bounded_below(0.2)
        @test isapprox(c1.constrained_to_unconstrained(1.0) - (log(1.0 - 0.2)), 0)
        @test isapprox(c1.unconstrained_to_constrained(0.0) - (exp(0.0) + 0.2), 0)

        c2 = bounded_above(0.2)
        @test isapprox(c2.constrained_to_unconstrained(-1.0) - (log(0.2 - -1.0)), 0)
        @test isapprox(c2.unconstrained_to_constrained(10.0) - (0.2 - exp(10.0)), 0)

        c3 = bounded(-0.1,0.2)
        @test isapprox(c3.constrained_to_unconstrained(0.0) - (log( (0.0 - -0.1) / (0.2 - 0.0))) , 0 )
        @test isapprox(c3.unconstrained_to_constrained(1.0) - ((0.2 * exp(1.0) + -0.1) / (exp(1.0) + 1)) , 0)
        @test_throws DomainError bounded(0.2,-0.1)

        #an example with user defined invertible transforms
        c_to_u = (x -> 3 * x + 14)
        u_to_c = (x -> (x - 14) / 3)

        c4 = Constraint(c_to_u, u_to_c)
        @test isapprox( c4.constrained_to_unconstrained(5.0) - c_to_u(5.0) , 0)
        @test isapprox( c4.unconstrained_to_constrained(5.0) - u_to_c(5.0) , 0)

    end

    @testset "ParameterDistribution" begin

        d = Parameterized(Normal(0,1))
        c = no_constraint()
        name = "unconstrained_normal"
        u = ParameterDistribution(d,c,name)
        @test u.distributions == [d]
        @test u.constraints == [c]
        @test u.names == [name]

        # Tests for the ParameterDistribution
        d = Parameterized(MvNormal(4,0.1))
        c = [no_constraint(),
             bounded_below(-1.0),
             bounded_above(0.4),
             bounded(-0.1,0.2)]

        name = "constrained_mvnormal"
        u = ParameterDistribution(d,c,name)
        @test u.distributions == [d]
        @test u.constraints == c
        @test u.names == [name]
        @test_throws DimensionMismatch ParameterDistribution(d,c[1:3],name)
        @test_throws DimensionMismatch ParameterDistribution(d,c, [name,"extra_name"])

        # Tests for the ParameterDistribution
        d1 = Parameterized(MvNormal(4,0.1))
        c1 = [no_constraint(),
              bounded_below(-1.0),
              bounded_above(0.4),
              bounded(-0.1,0.2)]
        name1 = "constrained_mvnormal"
        u1 = ParameterDistribution(d1,c1,name1)

        d2 = Samples([1.0 3.0 5.0 7.0; 9.0 11.0 13.0 15.0])
        c2 = [bounded(10,15),
              no_constraint()]
        name2 = "constrained_sampled"
        u2 = ParameterDistribution(d2,c2,name2)

        u = ParameterDistribution([d1,d2], [c1,c2], [name1,name2])
        @test u.distributions == [d1,d2]
        @test u.constraints == cat([c1,c2]...,dims=1)
        @test u.names == [name1,name2]
    end

    @testset "getter functions" begin
        # setup for the tests:
        d1 = Parameterized(MvNormal(4,0.1))
        c1 = [no_constraint(),
              bounded_below(-1.0),
              bounded_above(0.4),
              bounded(-0.1,0.2)]
        name1 = "constrained_mvnormal"
        u1 = ParameterDistribution(d1,c1,name1)

        d2 = Samples([1 2 3 4])
        c2 = [bounded(10,15)]
        name2 = "constrained_sampled"
        u2 = ParameterDistribution(d2,c2,name2)

        u = ParameterDistribution([d1,d2], [c1,c2], [name1,name2])

        # Test for get_dimension(s)
        @test get_total_dimension(u1) == 4
        @test get_total_dimension(u2) == 1
        @test get_total_dimension(u) == 5
        @test get_dimensions(u1) == [4]
        @test get_dimensions(u2) == [1]
        @test get_dimensions(u) == [4,1]

        # Tests for get_name
        @test get_name(u1) == [name1]
        @test get_name(u) == [name1, name2]

        # Tests for get_n_samples
        @test typeof(get_n_samples(u)[name1]) <: String
        @test get_n_samples(u)[name2] == 4

        # Tests for get_distribution
        @test get_distribution(d1) == MvNormal(4,0.1)
        @test get_distribution(u1)[name1] == MvNormal(4,0.1)
        @test typeof(get_distribution(d2)) == Array{Int64, 2}
        @test typeof(get_distribution(u2)[name2]) == Array{Int64, 2}

        d = get_distribution(u)
        @test d[name1] == MvNormal(4,0.1)
        @test typeof(d[name2]) == Array{Int64, 2}

        # Test for get_all_constraints
        @test get_all_constraints(u) == cat([c1,c2]...,dims=1)
    end

    @testset "statistics functions" begin

        # setup for the tests:
        d1 = Parameterized(MvNormal(4,0.1))
        c1 = [no_constraint(),
              bounded_below(-1.0),
              bounded_above(0.4),
              bounded(-0.1,0.2)]
        name1 = "constrained_mvnormal"
        u1 = ParameterDistribution(d1,c1,name1)

        d2 = Samples([1 2 3 4])
        c2 = [bounded(10,15)]
        name2 = "constrained_sampled"
        u2 = ParameterDistribution(d2,c2,name2)

        d3 = Parameterized(Beta(2,2))
        c3 = [no_constraint()]
        name3 = "unconstrained_beta"
        u3 = ParameterDistribution(d3,c3,name3)

        u = ParameterDistribution([d1,d2],[c1,c2],[name1,name2])

        d4 = Samples([1 2 3 4 5 6 7 8; 8 7 6 5 4 3 2 1])
        c4 = [no_constraint(),
              no_constraint()]
        name4 = "constrained_MVsampled"
        v = ParameterDistribution([d1,d2,d3,d4],[c1,c2,c3,c4],[name1,name2,name3,name4])

        # Tests for sample distribution
        seed=2020
        Random.seed!(seed)
        s1 = rand(MvNormal(4,0.1),1)
        Random.seed!(seed)
        @test sample_distribution(u1) == s1

        Random.seed!(seed)
        s1 = rand(MvNormal(4,0.1),3)
        Random.seed!(seed)
        @test sample_distribution(u1,3) == s1

        Random.seed!(seed)
        idx = StatsBase.sample(collect(1:size(d2.distribution_samples)[2]),1)
        s2 = d2.distribution_samples[:, idx]
        Random.seed!(seed)
        @test sample_distribution(u2) == s2

        Random.seed!(seed)
        idx = StatsBase.sample(collect(1:size(d2.distribution_samples)[2]), 3)
        s2 = d2.distribution_samples[:, idx]
        Random.seed!(seed)
        @test sample_distribution(u2,3) == s2

        Random.seed!(seed)
        s1 = sample_distribution(u1,3)
        s2 = sample_distribution(u2,3)
        Random.seed!(seed)
        s = sample_distribution(u,3)
        @test s == cat([s1,s2]...,dims=1)

        #Test for get_logpdf
        @test_throws ErrorException get_logpdf(u,zeros(get_total_dimension(u)))
        x_in_bd = [0.5]
        Random.seed!(seed)
        lpdf3 = logpdf.(Beta(2,2),x_in_bd)[1] #throws deprecated warning without "."
        Random.seed!(seed)
        @test isapprox(get_logpdf(u3,x_in_bd) - lpdf3 , 0.0; atol=1e-6)
        @test_throws DimensionMismatch get_logpdf(u3, [0.5,0.5])

        #Test for get_cov, get_var
        block_cov = cat([get_cov(d1),get_var(d2),get_var(d3),get_cov(d4)]..., dims=(1,2))
        @test isapprox(get_cov(v) - block_cov, zeros(get_total_dimension(v),get_total_dimension(v)); atol=1e-6)
        #Test for get_mean
        means = cat([get_mean(d1), get_mean(d2), get_mean(d3), get_mean(d4)]...,dims=1)
        @test isapprox(get_mean(v) - means, zeros(get_total_dimension(v)); atol=1e-6)

    end

    @testset "transform functions" begin
        #setup for the tests
        d1 = Parameterized(MvNormal(4,0.1))
        c1 = [no_constraint(),
              bounded_below(-1.0),
              bounded_above(0.4),
              bounded(-0.1,0.2)]
        name1 = "constrained_mvnormal"
        u1 = ParameterDistribution(d1,c1,name1)

        d2 = Samples([1.0 3.0 5.0 7.0; 9.0 11.0 13.0 15.0])
        c2 = [bounded(10,15),
              no_constraint()]
        name2 = "constrained_sampled"
        u2 = ParameterDistribution(d2,c2,name2)

        u = ParameterDistribution([d1,d2], [c1,c2], [name1,name2])

        x_unbd = rand(MvNormal(6,3), 1000)  #6 x 1000
        # Tests for transforms
        x_real_constrained1 = mapslices(x -> transform_unconstrained_to_constrained(u1,x), x_unbd[1:4,:]; dims=1)
        @test isapprox(x_unbd[1:4,:] - mapslices(x -> transform_constrained_to_unconstrained(u1,x), x_real_constrained1; dims=1), zeros(size(x_unbd[1:4,:])); atol=1e-8)
        x_real_constrained2 = mapslices(x -> transform_unconstrained_to_constrained(u2,x), x_unbd[5:6,:]; dims=1)
        @test isapprox(x_unbd[5:6,:] -  mapslices(x -> transform_constrained_to_unconstrained(u2,x), x_real_constrained2; dims=1), zeros(size(x_unbd[5:6,:])); atol=1e-8)

        x_real = mapslices(x -> transform_unconstrained_to_constrained(u,x), x_unbd; dims=1)
        x_unbd_tmp = mapslices(x -> transform_constrained_to_unconstrained(u,x), x_real; dims=1)
        @test isapprox(x_unbd - x_unbd_tmp,zeros(size(x_unbd)); atol=1e-8)

    end

end


text:
# To use this version of 'Gravityfield':
#push!(LOAD_PATH, joinpath(@__DIR__, "../src")) (or similar)
# using Gravityfield_unitless
#
# Possible TODO: Use Traits to distinguish position vectors.
# Currently,
module Gravityfield_unitless
export γ,  m_earth ,  r_earth ,  r_sun ,  d_sun_earth ,  d_moon_earth, ρ_earth
export m_sun ,  m_moon ,  yr ,  ρ_sun
export V, ρ, Pos, Npos
export nest, slice_1, slice_2, slice_3, slice_4, slice_5, slice_6
export Body, Bodies, Sphere
export earthlikeSphere, earthlikeSphere_r_known, earthlikeSphere_m_known
export bodies_3d, bodies_2d, bodies_1d
export rvec, uvec, len, g, gvec, gfieldvals
export xygrids, sqgrids
import Base.show

"Newton's gravity constant"
const γ=6.67428E-11 #m^3*kg^-1*s^-2
const m_earth = 5.97E24 # kg
const r_earth = 6300000. # m
const r_sun = 695510000. #m
const d_sun_earth = 1.496e11  # m
const d_moon_earth = 384400000. # m
const m_sun = 1.989E30 # kg
const m_moon = 7.34767309e22
const yr = 31556926 # s
const ρ_sun = 1410. # kg/m^3

"Volume of a sphere"
V(r) = 4π / 3 * r^3
"Density of a sphere"
ρ(;m=1., r=1.) = m / V(r)
const ρ_earth = ρ(m = m_earth, r = r_earth)
"Position vector"
const Pos = Vector{T} where T<:Number
Pos(p) = Vector(p)
import Base.show
function show(io::IO, ::MIME"text/plain", p::Pos{T}) where T# long form
    typ = typeof(p)
    ioc = IOContext(io, :typeinfo => typ)
    print(ioc, "Pos{", T, "}([")
    # avoid recursive call
    typ = typeof(p)
    ioc = IOContext(io, :typeinfo => typ)
    for i = 1:(length(p)-1)
        print(ioc, p[i], ", ")
    end
    print(ioc, p[length(p)], "])")
end
"Nested positions. Each element is a position vector; the container can be a vector or matrix"
const Npos = Union{Array{S, 2}, Array{S,1}} where S<:Pos
Npos(p::Npos) = Array(p)
function show(io::IO, np::Npos) # short form
    np_unitless = map(p->ustrip(p), np)
    typ = typeof(np_unitless)
    ioc = IOContext(io, :typeinfo => typ)
    show(ioc, np_unitless)
    printstyled(ioc, unit(eltype(np[1])); color=:cyan)
    nothing
end
function show(io::IO, ::MIME"text/plain", np::Npos) # long form
    etyp = eltype(np[1])
    np_unitless = map(p->ustrip(p), np)
    typul = eltype(np_unitless[1])
    if ndims(np) == 1
        r= size(np)[1]
        println(io, r, "-element Npos{Pos{", typul, "}}:")
    else
        r, c = size(np)
        println(io, r, "x", c, " Npos{Pos{", typul, "}}:")
    end
    typ = typeof(np_unitless)
    ioc = IOContext(io, :typeinfo => typ)
    show(ioc, np_unitless)
    printstyled(io, unit(eltype(etyp)); color=:cyan)
end
const AA = AbstractArray
function nest(x::AA, y::AA, z::Number)
    Npos(map(x, y) do x, y
            [x, y, z]
        end)
end
function nest(x::AA, y::Number, z::Number)
    Npos(map(x) do x
            [x, y, z]
        end)
end
function nest(x::Number, y::AA, z::Number)
    Npos(map(y) do y
            [x, y, z]
        end)
end
function nest(x::Number, y::Number, z::AA)
    Npos(map(z) do z
            [x, y, z]
        end)
end

function nest(x::AA, y::Number)
    Npos(map(x) do x
            [x, y]
        end)
end
function nest(x::Number, y::AA)
    Npos(map(y) do y
            [x, y]
        end)
end
function nest(x::AA, y::AA, z::AA, α::AA, β::AA, γ::AA)
    Npos(map(x, y, z, α, β, γ) do x, y, z, α, β, γ
            [x, y, z, α, β, γ]
        end)
end
function nest(x::AA, y::AA, z::AA, α::AA, β::AA)
    Npos(map(x, y, z, α, β) do x, y, z, α, β
            [x, y, z, α, β]
        end)
end
function nest(x::AA, y::AA, z::AA, α::AA)
    Npos(map(x, y, z, α, β) do x, y, z, α
            [x, y, z, α]
        end)
end
function nest(x::AA, y::AA, z::AA)
    Npos(map(x, y, z) do x, y, z
            [x, y, z]
        end)
end
function nest(x::AA, y::AA)
    Npos(map(x, y) do x, y
            [x, y]
        end)
end
function nest(x::AA)
    Npos(map(x) do x
            [x]
        end)
end
"Slice out 1st component from a nested vector or matrix"
slice_1(np::Npos) = map(p-> p[1], np)
"Slice out 2nd component from a nested vector or matrix"
slice_2(np::Npos) = map(p-> p[2], np)
"Slice out 3rd component from a nested vector or matrix"
slice_3(np::Npos) = map(p-> p[3], np)
"Slice out 4th component from a nested vector or matrix"
slice_4(np::Npos) = map(p-> p[4], np)
"Slice out 5th component from a nested vector or matrix"
slice_5(np::Npos) = map(p-> p[5], np)
"Slice out sixth component from a nested vector or matrix"
slice_6(np::Npos) = map(p-> p[6], np)

abstract type Body end
"""
Constructor defaults to a spherical earth. Alternatively, use
```
    earthlikeSphere(p = my_p, m = my_p)
    earthlikeSphere(p = my_p, r = my_r)
```
"""
Base.@kwdef struct Sphere <: Body
    p::Pos = [0., 0., 0.]
    m::Number = m_earth
    r::Number = r_earth
    Sphere(p, m, r) = m > zero(m) && r > zero(r) ? new(p, m, r) : error("no negative m or r")
end

function show(io::IO, b::Sphere) # short form
    ioc = IOContext(io, :typeinfo => typeof(b.p))
    print(ioc, "Sphere(p = ", b.p, ", m = ", b.m, ", r = ", b.r, ")")
end

function show(io::IO, mime::MIME"text/plain", b::Sphere) # Long form
    print(io, "Sphere(p = ")
    # perhaps unnecessary splitting...
    show(io, mime, b.p)
    print(io, ", m = ")
    show(io, mime, b.m)
    print(io, ", r = ")
    show(io, mime,  b.r)
    print(io, ")")
end
"Alternative constructors of Sphere for when mass or radius is unknown"
earthlikeSphere_m_known(;p = [0, 0, 0], m = m_earth) = Sphere(p, m , (3m / (4π * ρ_earth))^(1/3))
earthlikeSphere_r_known(;p = [0, 0, 0], r = r_earth) = Sphere(p, ρ_earth * V(r), r)

const Bodies = Vector{T} where T <:Body
Bodies(bds::Bodies) = Vector(bds)
function bodies_3d()
    earth = Sphere()
    moon = earthlikeSphere_m_known(p = [d_moon_earth, 0m, 0m], m = m_moon)
    sun = Sphere([-d_sun_earth, 0, 0], m_sun, (3m_sun / (4π * ρ_sun))^(1/3))
    Bodies([sun, earth, moon])
end
function bodies_2d()
    earth = Sphere(p = [0,0]m)
    moon = earthlikeSphere_m_known(p = [d_moon_earth, 0m], m = m_moon)
    sun = Sphere([-d_sun_earth, 0], m_sun, (3m_sun / (4π * ρ_sun))^(1/3))
    Bodies([sun, earth, moon])
end
function bodies_1d()
    earth = Sphere(p = [0]m)
    moon = earthlikeSphere_m_known(p = [d_moon_earth], m = m_moon)
    sun = Sphere([-d_sun_earth], m_sun, (3m_sun / (4π * ρ_sun))^(1/3))
    Bodies([sun, earth, moon])
end

import Base.show_circular
import Base._truncate_at_width_or_chars
function show(io::IO, ::MIME"text/plain", bds::Bodies) # long form.
    # show more descriptively, with one line per body. Based on Dict.
    recur_io = IOContext(io, :SHOWN_SET => bds)
    limit::Bool = get(io, :limit, false)
    if !haskey(io, :compact)
        recur_io = IOContext(recur_io, :compact => true)
    end
    summary(io, bds)
    isempty(bds) && return
    print(io, ":")
    show_circular(io, bds) && return
    if limit
        sz = displaysize(io)
        rows, cols = sz[1] - 3, sz[2]
        rows < 2   && (print(io, " …"); return)
        cols < 12  && (cols = 12) # Minimum widths of 2 for key, 4 for value
        cols -= 6 # Subtract the widths of prefix "  " separator " => "
        rows -= 1 # Subtract the summary

        # determine max key width to align the output, caching the strings
        ps = Vector{AbstractString}(undef, min(rows, length(bds)))
        ms = Vector{AbstractString}(undef, min(rows, length(bds)))
        rs = Vector{AbstractString}(undef, min(rows, length(bds)))
        plen = 0
        mlen = 0
        rlen = 0
        for (i, b) in enumerate(bds)
            i > rows && break
            st = sprint(show, b.p, context=recur_io, sizehint=0)
            ps[i] = "p = " * st
            st = sprint(show,  b.m, context=recur_io, sizehint=0)
            ms[i] = " m = " * st
            st = sprint(show, b.r, context=recur_io, sizehint=0)
            rs[i] = " r = " * st
            plen = clamp(length(ps[i]), plen, cols)
            mlen = clamp(length(ms[i]), mlen, cols)
            rlen = clamp(length(rs[i]), rlen, cols)
        end
        if plen > max(div(cols, 2), cols - mlen - rlen)
            plen = max(cld(cols, 3), cols - mlen - rlen)
        end
    else
        rows = cols = typemax(Int)
    end

    for (i, bd) in enumerate(bds)
        print(io, "\n  ")
        i == rows < length(bds) && (print(io, rpad("⋮", plen), ", ⋮"); break)
        if limit
            pos = rpad(_truncate_at_width_or_chars(ps[i], plen, "\r\n"), plen)
        else
            st = sprint(show, bd.p, context=recur_io, sizehint=0)
            pos = "p = " * st
        end
        print(recur_io, pos)
        if limit
            val = rpad(_truncate_at_width_or_chars(ms[i], cols - plen, "\r\n"), mlen)
            print(io, val)
            val = _truncate_at_width_or_chars(rs[i], cols - plen - mlen, "\r\n")
            print(io, val)
        else
            st = " m = " * sprint(show, bd.m, context=recur_io, sizehint=0) *
                " r = " * sprint(show, bd.r, context=recur_io, sizehint=0)
            print(recur_io, st)
        end
    end
end

"Vector from position (of) a to position (of) b."
rvec(a::T, b::T) where T<: Body = b.p - a.p
rvec(a::Pos, b::Body) = b.p - a
rvec(a::Body, b::Pos) = b - a.p
rvec(a::T, b::T) where T<: Pos = b - a
"Unity vector of vector"
uvec(p::Pos) = len(p) > zero(eltype(p)) ? p / len(p) : zero(p)
"Unity vector from position (of) a to position (of) b."
uvec(a::T, b::T) where T<: Body = uvec(b.p - a.p)
uvec(a::Pos, b::Body) = uvec(b.p - a)
uvec(a::Body, b::Pos) = uvec(b - a.p)
uvec(a::T, b::T) where T<: Pos = uvec(b - a)

"Euclidean vector length"
len(p) =  sqrt(sum(p.*p))
"Euclidean length between positions (of) a and b"
len(a, b) = len(rvec(a, b))

"Gravitational acceleration magnitude"
g(m::Number, r::Number) = r > zero(r) ? γ * m / r^2 : zero(γ) / one(r)^2  # Test thoroughly for overflow

"Gravity vector on thisbody or position due to otherbody or bodies"
gvec(p::Pos, body::Body) = uvec(p, body) * g(body.m, len(p, body))
gvec(thisbody::Body, otherbody::Body) = gvec(thisbody.p, otherbody)
function gvec(p::Pos, sphere::Sphere)
    l = len(p, sphere)
    if l > sphere.r || l == 0.
        return uvec(p, sphere) * g(sphere.m, l)
    else
        # Inside sphere
        return uvec(p, sphere) * g(sphere.m, sphere.r) * l / sphere.r
    end
end
function gvec(p::Pos, bodies::Bodies)
    gvec_one = b -> gvec(p, b)
    sum(gvec_one.(bodies))
end
gvec(body::Body, bodies::Bodies) = gvec(body.p, bodies)
gvec(vp::Npos, bodies::Bodies) = map(p -> gvec(p, bodies), vp)

"Gravity field strength [length / time^2], nested position vectors"
gfieldvals(vp::Npos, bodies::Bodies) = map(p -> len(gvec(p, bodies)), vp)

"Return matrices mat_x, mat_y, size n * n"
function xygrids(;xmin = -3., xmax = 3., ymin, ymax, n = 100)
    mat_x = repeat(range(xmin, xmax, length = n)', n, 1)
    mat_y = repeat(range(ymin, ymax, length = n), 1, n)
    mat_x, mat_y
end
"100x100 grids around origo"
sqgrids(mi, ma) = xygrids( xmin = mi,xmax = ma, ymin = -(ma-mi)/2, ymax = (ma-mi)/2, n = 100)

nothing
end

max_stars_repo_name: hustf/Gravityfield

max_issues_repo_name: hustf/Gravityfield

max_forks_repo_name: hustf/Gravityfield


text:
type Object
	enabled::Bool
	obj::Any
end

type Group
	enabled::Bool
	objs::Dict{Symbol,Any}
end

# compiled after start, for faster loop?
#RenderProcess()

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
using Plots,Statistics, LaTeXStrings,JLD
dim=200
function neighbor_checking(network_,dim,i,j)
    neighbors_list=[]
    if i+1!=dim+1
        push!(neighbors_list,network_[i+1,j])
    end
    if i-1!=0
        push!(neighbors_list,network_[i-1,j])
    end
    if j+1!=dim+1
        push!(neighbors_list,network_[i,j+1])
    end
    if j-1!=0
        push!(neighbors_list, network_[i,j-1])
    end
    return neighbors_list
end
function on_or_blocked(network_,dim,i,j,p,S)
    if p>rand()
        network_[i,j]=1
        S+=1
    else
        network_[i,j]=-1
    end
    return network_,S
end

function operation(dim,p)
    network_=zeros(Int,dim,dim)
    network_[rand(1:dim),rand(1:dim)]=1 #the first entry
    S=1
    for step in 1:(dim^2)
        operation_list=[]
        for i in 1:dim
            for j in 1:dim
                if network_[i,j]==0 && (1 in neighbor_checking(network_,dim,i,j))
                    push!(operation_list,[i,j])
                end
            end
        end
        if length(operation_list)>0
            for i in operation_list
                network_,S=on_or_blocked(network_,dim,i[1],i[2],p,S)
            end
        elseif length(operation_list)==0
            break
        end
    end
    return network_,S
end
function RadiusOfGyration(network_,dim)
    i_list=[]
    j_list=[]
    for index in findall(x->x==1, network_)
        push!(i_list, index[1])
        push!(j_list, index[2])
    end
    i_com=mean(i_list)
    j_com=mean(j_list)
    fraction_list=[]
    for i in 1:dim
        for j in 1:dim
            if network_[i,j]==1
                push!(fraction_list, ((i-i_com)^2 + (j-j_com)^2))
            end
        end
    end
    RadiusOfGyration=sqrt(mean(fraction_list))
    return RadiusOfGyration
end
probability=[0.5, 0.55,0.59]
SavedData_xi=[]
SavedData_S=[]

for p in probability
    xi=[]
    S_list=[]
    for run_num in 1:100
        network_,S=operation(dim,p)
        push!(S_list,S)
        push!(xi,RadiusOfGyration(network_,dim))
    end
    push!(SavedData_xi, xi)
    push!(SavedData_S, S_list)
end
save("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter4\\4.7\\ClusterGrowth_S_xi.jld",
"data_xi",SavedData_xi,
"data_S", SavedData_S)
#plots:
load("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter4\\4.7\\ClusterGrowth_S_xi.jld")
xdata=[]
ydata=[]
for i in 1:3
    push!(ydata,log(mean(SavedData_S[i])))
    push!(xdata,log(mean(SavedData_xi[i])))
end
scatter(xdata, ydata,color=:orange, legend=false, dpi=400,
 xlabel=L"Log\ (\xi)", ylabel=L"log\ (S)",title=L"log\ (S) \_ Log\ (\xi)")
plot!(xdata, ydata,legend=false, dpi=400)
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter4\\Fig\\4.7_log(s)_log(xi).png")


text:
# This file contains functions for the Revenue Efficiency DEA model
"""
    RevenueDEAModel
An data structure representing a revenue DEA model.
"""
struct RevenueDEAModel <: AbstractRevenueDEAModel
    n::Int64
    m::Int64
    s::Int64
    rts::Symbol
    disposX::Symbol
    dmunames::Union{Vector{String},Nothing}
    eff::Vector
    lambda::SparseMatrixCSC{Float64, Int64}
    techeff::Vector
    alloceff::Vector
    Xtarget::Matrix
    Ytarget::Matrix
end

"""
    dearevenue(X, Y, P)
Compute revenue efficiency using data envelopment analysis for
inputs `X`, outputs `Y` and price of outputs `P`.

# Optional Arguments
- `rts=:VRS`: chooses variable returns to scale. For constant returns to scale choose `:CRS`.
- `dispos=:Strong`: chooses strong disposability of inputs. For weak disposability choose `:Weak`.
- `names`: a vector of strings with the names of the decision making units.

# Examples
```jldoctest
julia> X = [5 3; 2 4; 4 2; 4 8; 7 9.0];

julia> Y = [7 4; 10 8; 8 10; 5 4; 3 6.0];

julia> P = [3 2; 3 2; 3 2; 3 2; 3 2.0];

julia> dearevenue(X, Y, P)
Revenue DEA Model
DMUs = 5; Inputs = 2; Outputs = 2
Orientation = Output; Returns to Scale = VRS
──────────────────────────────────
    Revenue  Technical  Allocative
──────────────────────────────────
1  0.644444   0.777778    0.828571
2  1.0        1.0         1.0
3  1.0        1.0         1.0
4  0.5        0.5         1.0
5  0.456522   0.6         0.76087
──────────────────────────────────
```
"""
function dearevenue(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector},
    P::Union{Matrix,Vector}; rts::Symbol = :VRS, dispos::Symbol = :Strong,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)::RevenueDEAModel

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    np, sp = size(P, 1), size(P, 2)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if np != ny
        throw(DimensionMismatch("number of rows in P and Y ($np, $ny) are not equal"));
    end
    if sp != s
        throw(DimensionMismatch("number of columns in P and Y ($sp, $s) are not equal"));
    end

    if dispos != :Strong && dispos != :Weak
        throw(ArgumentError("`disposY` must be :Strong or :Weak"));
    end

    # Default optimizer
    if optimizer === nothing
        optimizer = DEAOptimizer(:LP)
    end

    # Get maximum revenue targets and lambdas
    n = nx

    Xtarget = X[:,:]
    Ytarget, rlambdaeff = deamaxrevenue(X, Y, P, rts = rts, dispos = dispos, optimizer = optimizer)

    # Revenue, technical and allocative efficiency
    refficiency  = vec( sum(P .* Y, dims = 2) ./ sum(P .* Ytarget, dims = 2) )
    techefficiency = 1 ./ efficiency(dea(X, Y, orient = :Output, rts = rts, slack = false, disposX = dispos, optimizer = optimizer))
    allocefficiency = refficiency ./ techefficiency
    return RevenueDEAModel(n, m, s, rts, dispos, names, refficiency, rlambdaeff, techefficiency, allocefficiency, Xtarget, Ytarget)

end

ismonetary(model::RevenueDEAModel)::Bool = false;

function Base.show(io::IO, x::RevenueDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    disposX = x.disposX
    dmunames = names(x)

    eff = efficiency(x)
    techeff = efficiency(x, :Technical)
    alloceff = efficiency(x, :Allocative)

    if !compact
        print(io, "Revenue DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Orientation = Output")
        print(io, "; Returns to Scale = ", string(x.rts))
        print(io, "\n")
        if disposX == :Weak print(io, "Weak disposability of inputs \n") end

        show(io, CoefTable(hcat(eff, techeff, alloceff), ["Revenue", "Technical", "Allocative"], dmunames))
    end

end


text:
#= Possible components for the gui
    - window displaying the simulation
    - ⏵/⏸ button for the simulation
    - slider for variables (temperature, adhesiveness, etc)
    - display for much of MH steps
    - count for number of cells
=#

# using CellularPotts
using GLMakie
using Colors

function CellGUI(CPM)

    #Create a blank scene with some padding
    outer_padding = 30
    scene, layout = layoutscene(outer_padding, resolution = (700, 700), backgroundcolor = RGB(0.98, 0.98, 0.98))

    #--------Simulation Screen--------
    axSim = layout[1, 1] = Axis(scene, title = "Simulation")

    #--------Buttons--------
    #Currently 2 buttons: a play/pause and a stop button
    #Place the buttons below the simulation and align to the left
    layout[2,1] = axButtons = GridLayout(tellwidth = false,halign = :left)
    buttonsLabels = ["▶","■","Divide!"]
    #Loop through and assign button labels, width, and color
    axButtons[1,1:length(buttonsLabels)] = [Button(
        scene,
        label = lab,
        width = 70,
        buttoncolor = :grey) for lab in buttonsLabels]

    #Set the buttons to individual variables
    playpause, stop, cellDivideButton = contents(axButtons)

    #--------Sliders--------
    #The first slider controls the inverse temperature (higher temperatures ⟹ accepting higher energy proposals)
    #The second slider changes volume constraint (lower values ⟹ cells more easily drift from desired volume)
    # slideLabels = ["Temperature (β):","Volume\nRestriction (λᵥ)","Memory\nRestriction (λₐ)","Memory\nLength (Ma)"]
    # slideRanges = [1:0.01:4, 1:10, 1.0:0.1:10.0, 10:40]
    slideLabels = ["Temperature (β):","Volume\nRestriction (λᵥ)"]
    slideRanges = [10:30, 30:70]

    #Put the sliders to the right of the simulation (tellheight=false means the slider and simulation heights can differ)
    Sublayout = GridLayout(valign = :top,tellheight=false)
    layout[1, 2] = Sublayout

    #Loop through sliders and assign labels and ranges
    sliders = [labelslider!(scene,
                            sl,
                            sr;
                            format = x -> "$(x)") for (sl, sr) in zip(slideLabels, slideRanges)]

    #Vertically stack the sliders
    sliderAxes = layout[1:length(slideLabels),2] = getfield.(sliders, :layout)
    Sublayout[:v] = sliderAxes

    #--------Update Simulation Screen with Model--------
    timestep = Node(1) #A Node is a variable being observed by the simulation
    frameskip = 10_000 #The MHStep will often fail b/c of bad random choice to change grid, so nothing will update

    #Everytime time is updated, run MHStep until the next frameskip
    heatmap_node = @lift begin
        currentTime = $timestep
        for t in 1:frameskip
            MHStep!(CPM)
        end
        CPM.cellMembership
    end

    #Create a heatmap on the simulation axis (layout 1,1)
    heatmap!(axSim,
            heatmap_node,
            show_axis = false,
            colormap = :Purples) #:Greys_3
    tightlimits!.(axSim)
    hidedecorations!.(axSim) #removes axis numbers

    #Generate all the adjacent squares in the grid
    n = CPM.graph.dimension[1]
    edgeConnectors = Edge2Grid(n)

    #Generate all of the edge Connections by putting a point on each cell corner
    horizontal = [Point2f0(x, y) => Point2f0(x+1, y) for x in 0:n-1, y in 0:n]
    vertical = [Point2f0(x, y) => Point2f0(x, y+1) for y in 0:n-1, x in 0:n]
    points = vcat(horizontal[:],vertical[:])

    #Determine the transparency of the linesegments
    gridflip = rotl90(CPM.cellMembership) #https://github.com/JuliaPlots/Makie.jl/issues/205

    #Cell borders are outlined in black
    black = RGBA{Float64}(0.0,0.0,0.0,1.0);
    clear = RGBA{Float64}(0.0,0.0,0.0,0.0);

    #Loop through all the grid connected and assign the correct color
    currentEdgeColors = [gridflip[edges[1]]==gridflip[edges[2]] ? clear : black for edges in edgeConnectors];

    #For each time update, recolor all of the edges
    lineColors_node = @lift begin
        currentTime = $timestep

        gridflip = rotl90(CPM.cellMembership)

        for (i,edges) in enumerate(edgeConnectors)
            currentEdgeColors[i] = gridflip[edges[1]]==gridflip[edges[2]] ? clear : black
        end

        currentEdgeColors
    end

    #Plot all of the line segments onto the simulation
    linesegments!(
        axSim,
        points,
        color = lineColors_node,
        linewidth = 2
    )

    #If the play/pause button is clicked, change the label
    lift(playpause.clicks) do clicks
        if isodd(clicks)
            playpause.label = "𝅛𝅛"
        else
            playpause.label = "▶"
        end
    end

    #If the play/pause button is clicked, change the label
    # lift(cellDivideButton.clicks) do clicks
    #     CellDivide!(CPM,rand(1:CPM.σ))
    # end

    #can this be a loop?
        #Temperature slider update
        lift(sliders[1][:slider].value) do val
            CPM.temperature = val
        end

        #Volume constraint slider update
        lift(sliders[2][:slider].value) do val
           VP.λᵥ = val #skip the medium (with index zero)
        end

        # #Memory Restriction
        # lift(sliders[3][:slider].value) do val
        #     CPM.pm.λact = val #skip the medium (with index zero)
        # end

        # lift(sliders[4][:slider].value) do val
        #     CPM.pm.Ma = val #skip the medium (with index zero)
        # end

    #Active cell movement (patrol)
    # if CPM.pm.on #if patrol movement is on
    #     #do I need this?
    #     heatmap_Gm = @lift begin
    #         currentTime = $timestep
    #         CPM.pm.Gm
    #     end

    #     #transparent colors
    #     colmap = RGBA.(colormap("Reds"),0.5)

    #     #Create a heatmap on the simulation axis (layout 1,1)
    #     heatmap!(axSim,
    #             heatmap_Gm,
    #             show_axis = false,
    #             colormap = colmap) #:Greys_3
    #     tightlimits!.(axSim)
    #     hidedecorations!.(axSim) #removes axis numbers
    # end

    #Create a window
    display(scene)

    #Keep running the simulation until stopped
    runsim = true
    stop.clicks[] = 0
    while runsim

        #Is the pause button pushed?
        if playpause.label[] == "▶"
            timestep[] += 1
        end

        #Has the stop button been pushed?

        if  stop.clicks[] == 1
            runsim = false
            GLMakie.destroy!(GLMakie.global_gl_screen()) #close the window
        end

        sleep(eps())
        #println(CPM.H)

    end
end

#Create a new simulation
# CPM = CellPotts(n=100,σ=1,Vd=[200],patrol=false)
# CPM = CellPotts(n=150,σ=500,Vd=rand(35:45,500),patrol=false)

# for i=1:10_000_000
#     MHStep!(CPM)
#     i % 10_000 == 0 ? println(i) : nothing
# end

CellGUI(CPM)

#heatmap(CPM.grid)
#CellDivide!(CPM,1)


text:
# used when thinking about the `alpha` exponent for xQ
using PyPlot
using LaTeXStrings

x = collect(linspace(0.,1.,1000))
y(r) = x.^r
rs = [1/4,1/3,1/2,2/3,3/4,1.]

fig = plt[:figure](figsize=(5.0,4.0))
ax = fig[:subplots](1,1)
fs = 10

for r in rs
    lw = 1
    if r == 1/2
        lw = 3
    end
    ax[:plot](x,y(r),label="$(@sprintf("r=%.2f",r))",lw = lw)
end

ax[:set_title](L"$y=x^r$ for different $r$",fontsize=fs)
ax[:set_xlabel](L"x",fontsize=fs)
ax[:set_ylabel](L"y=x^r",fontsize=fs)
ax[:legend]()

savefig("figs/power_comparison.pdf",dpi=300,transparent=true)

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
using TextStylometry
using Base.Test

@testset "I/O tests" begin include("corpus_tests.jl") end
@testset "measures tests" begin include("measures_tests.jl") end
@testset "bootstrap tests" begin include("bootstrap_tests.jl") end

using Lint
lintpkg("TextStylometry")


text:
function read_graphs(n, file)
    G1 = SimpleGraph(n)
    G2 = SimpleGraph(n)
    const_G2 = 0
    for line in eachline(file)
        #println(line)
        if line == "graph2"
            const_G2 += 1
        elseif line == "matrice log likelihood ratio"
            const_G2 += 1
        end
        if const_G2 == 0
            u = split(line, " : ")
            if length(u) >1
                add_edge!(G1, parse(Int, u[1]),parse(Int, u[2]))
                #println(parse(Int, u[1]), "|",  parse(Int, u[2]))
            end
        elseif const_G2 == 1
            u = split(line, " : ")
            if length(u) >1
                add_edge!(G2, parse(Int, u[1]),parse(Int, u[2]))
                #println(parse(Int, u[1]), "|",  parse(Int, u[2]))
            end
        end
    end
    return G1, G2
end

function Pair_IO(G1,G2,n,λ,s)
    N1 = Network(G1)
    N2 = Network(G2)
    #neig1 = create_neig_r(G1)
    #neig2 = create_neig_r(G2)
    deg1 = unique(sort(N1.deg))
    deg2 = unique(sort(N2.deg))
    min_deg = min(deg1[end],deg2[end])
    max_deg = max(deg1[end],deg2[end])

    dic_comb = Dict([((d,k),collect(combinations(range(1,d),k))) for d in 0:max_deg for k in 0:d])
    dic_perm = Dict([(k,collect(permutations(range(1,k)))) for k in 0:min_deg])
    return Pair_ER(N1,N2,SimpleGraph(n),n,λ,s,1:n,1:n,zeros(Int64,n,2),dic_comb,dic_perm)
end

function create_matrix_check(data)
    n = length(data)
    M_data = zeros(n,n)
    for (i,l) in enumerate(data)
        for (j,e) in enumerate(l)
            M_data[i,j] = e
        end
    end
    M_data
end

function check_data(PG,M_data,n_iter)
    n,_ = size(M_data)
    mp_fin = run_bp(PG,n_iter)
    M_loglr = create_matrix_lr(PG,mp_fin)
    check = true
    for i in 1:n
        for j in 1:n
            if (M_data[i,j] - M_loglr[i,j]) > 0.00001
                #println(i,"|",j,"|", M_data[i,j],"|",M_loglr[i,j])
                check = false
            end
        end
    end
    check
end


text:
struct LUFactorization{P} <: AbstractLinearAlgorithm
    pivot::P
end

function LUFactorization()
    pivot = @static if VERSION < v"1.7beta"
        Val(true)
    else
        RowMaximum()
    end
    LUFactorization(pivot)
end

function SciMLBase.solve(cache::LinearCache, alg::LUFactorization)
    cache.A isa Union{AbstractMatrix,AbstractDiffEqOperator} ||
        error("LU is not defined for $(typeof(prob.A))")
    cache = set_cacheval(cache, lu!(cache.A, alg.pivot))
    ldiv!(cache.cacheval, cache.b)
end

struct QRFactorization{P} <: AbstractLinearAlgorithm
    pivot::P
    blocksize::Int
end

function QRFactorization()
    pivot = @static if VERSION < v"1.7beta"
        Val(false)
    else
        NoPivot()
    end
    QRFactorization(pivot, 16)
end

function SciMLBase.solve(cache::LinearCache, alg::QRFactorization)
    cache.A isa Union{AbstractMatrix,AbstractDiffEqOperator} ||
        error("QR is not defined for $(typeof(prob.A))")
    cache = set_cacheval(
        cache,
        qr!(cache.A.A, alg.pivot; blocksize = alg.blocksize),
    )
    ldiv!(cache.cacheval, cache.b)
end

struct SVDFactorization{A} <: AbstractLinearAlgorithm
    full::Bool
    alg::A
end

SVDFactorization() = SVDFactorization(false, LinearAlgebra.DivideAndConquer())

function SciMLBase.solve(cache::LinearCache, alg::SVDFactorization)
    cache.A isa Union{AbstractMatrix,AbstractDiffEqOperator} ||
        error("SVD is not defined for $(typeof(prob.A))")
    cache = set_cacheval(cache, svd!(cache.A; full = alg.full, alg = alg.alg))
    ldiv!(cache.cacheval, cache.b)
end