text:
using RadiationDetectorSignals

function CalculateRiseTime(t, v::Array, low::Real, high::Real)
    tLow = 0
    tHigh = 0

    vMax = v[end]
    for i in 1:length(v)
        if v[i] >= low * vMax
            tLow = t[i]
            break
        end
    end
    for i in length(v):-1:1
        if v[i] <= high * vMax
            tHigh = t[i]
            break
        end
    end
    return tHigh-tLow
end

function CalculateRiseTime(w::RadiationDetectorSignals.RDWaveform, low::Real, high::Real)
    return CalculateRiseTime(w.time, abs.(w.value), low, high)
end

function CollectRiseTimes(events, low = 0.1, high = 0.9; contact = 1)
    riseTimes = Vector{Unitful.Time}()
    for i in 1:length(events)
        push!(riseTimes, CalculateRiseTime(events[i].waveforms[contact], low, high))
    end
    return riseTimes
end

function CollectRiseTimesFiltered(events, low = 0.1, high = 0.9; dt = 1e-9, contact = 1)
    riseTimes = Vector{Unitful.Time}()
    for i in 1:length(events)
        push!(riseTimes, CalculateRiseTime(events[i].waveforms[contact].time, CustomFilter(events[i].waveforms[contact].value, dt), low, high))
    end
    return riseTimes
end

max_stars_repo_name: leenderthayen/NDP

max_issues_repo_name: leenderthayen/NDP

max_forks_repo_name: leenderthayen/NDP


text:
function get_absolute_accuracy(obj::BracketedRealFieldUnivariateSolver)
    return jcall(obj, "getAbsoluteAccuracy", RealFieldElement, ())
end

function get_evaluations(obj::BracketedRealFieldUnivariateSolver)
    return jcall(obj, "getEvaluations", jint, ())
end

function get_function_value_accuracy(obj::BracketedRealFieldUnivariateSolver)
    return jcall(obj, "getFunctionValueAccuracy", RealFieldElement, ())
end

function get_max_evaluations(obj::BracketedRealFieldUnivariateSolver)
    return jcall(obj, "getMaxEvaluations", jint, ())
end

function get_relative_accuracy(obj::BracketedRealFieldUnivariateSolver)
    return jcall(obj, "getRelativeAccuracy", RealFieldElement, ())
end

function solve(obj::BracketedRealFieldUnivariateSolver, arg0::jint, arg1::RealFieldUnivariateFunction, arg2::RealFieldElement, arg3::RealFieldElement, arg4::AllowedSolution)
    return jcall(obj, "solve", RealFieldElement, (jint, RealFieldUnivariateFunction, RealFieldElement, RealFieldElement, AllowedSolution), arg0, arg1, arg2, arg3, arg4)
end

function solve(obj::BracketedRealFieldUnivariateSolver, arg0::jint, arg1::RealFieldUnivariateFunction, arg2::RealFieldElement, arg3::RealFieldElement, arg4::RealFieldElement, arg5::AllowedSolution)
    return jcall(obj, "solve", RealFieldElement, (jint, RealFieldUnivariateFunction, RealFieldElement, RealFieldElement, RealFieldElement, AllowedSolution), arg0, arg1, arg2, arg3, arg4, arg5)
end

function solve_interval(obj::BracketedRealFieldUnivariateSolver, arg0::jint, arg1::RealFieldUnivariateFunction, arg2::RealFieldElement, arg3::RealFieldElement)
    return jcall(obj, "solveInterval", BracketedRealFieldUnivariateSolver_Interval, (jint, RealFieldUnivariateFunction, RealFieldElement, RealFieldElement), arg0, arg1, arg2, arg3)
end

function solve_interval(obj::BracketedRealFieldUnivariateSolver, arg0::jint, arg1::RealFieldUnivariateFunction, arg2::RealFieldElement, arg3::RealFieldElement, arg4::RealFieldElement)
    return jcall(obj, "solveInterval", BracketedRealFieldUnivariateSolver_Interval, (jint, RealFieldUnivariateFunction, RealFieldElement, RealFieldElement, RealFieldElement), arg0, arg1, arg2, arg3, arg4)
end

ext: jl


text:
Random.seed!(192837465)
import Statistics.mean

function _getSignLinearValues(x0::Float64, n::Int64)
  μφ = 0.5-x0
  μg = x0^2-x0+0.5 # g = |φ|
  μgSq = x0^2-x0+1/3
  v = μφ/μg
  estimatorVariance = μg^2*v*(1-v)+v/n*(μgSq-μg^2)
  return μφ, μg, v, estimatorVariance
end

function _signLinearTest(x0::Float64, δ::Float64, c::Float64, n::Int64, m::Int64)
  μ() = rand()
  φ(x) = x-x0
  μφ, μg, v, estimatorVariance = _getSignLinearValues(x0, n)

  vs = Vector{Float64}(undef, m)
  flips = Vector{Int64}(undef, m)
  calls = Vector{Int64}(undef, m)
  for i in 1:m
    vs[i], flips[i], calls[i] = BernoulliFactory.signedEstimate(μ, φ, c, δ, n)
  end

  @test all(vs.>=0)
  @test all(vs.<=c)
  @test mean(vs.!= 0.0) ≈ v atol=0.01
  @test abs(mean(vs) - μφ)/sqrt(estimatorVariance/m) < 3

  for i in 1:m
    vs[i] > 0.0 && (calls[i] -= n)
  end
  @test mean(calls./flips) ≈ c/μg atol=0.1
end

function _signLinearTest(x0::Float64, a::Float64, b::Float64, δ::Float64,
  c::Float64, n::Int64, m::Int64)
  μ() = rand()
  φ(x) = x-x0
  _, μg, v, estimatorVariance = _getSignLinearValues(x0+b, n)
  μφ, _, _, _ = _getSignLinearValues(x0, n)

  vs = Vector{Float64}(undef, m)
  flips = Vector{Int64}(undef, m)
  calls = Vector{Int64}(undef, m)
  for i in 1:m
    vs[i], flips[i], calls[i] = BernoulliFactory.signedEstimate(μ, φ, a, b,
      δ, c, n)
  end

  @test all(vs.>=b)
  @test all(vs.<=max(2*b-a, c))
  @test mean(vs.!= b) ≈ v atol=0.01
  @test abs(mean(vs) - μφ)/sqrt(estimatorVariance/m) < 3

  for i in 1:m
    vs[i] > 0.0 && (calls[i] -= n)
  end
  @test mean(calls./flips) ≈ max(b-a,c-b)/μg atol=0.1
end

_signLinearTest(0.2, 0.1, 0.8, 3, ntrials)
_signLinearTest(0.4, 0.1, 1.0, 3, ntrials)
_signLinearTest(0.2, -0.2, 0.1, 0.2, 0.8, 5, ntrials)


text:
function babylonian_spiral(n)
    x = [0, 0]
    y = [0, 1]

    return x, y
end

max_stars_repo_name: basimr/lovelace

max_issues_repo_name: basimr/project-lovelace

max_forks_repo_name: basimr/project-lovelace


text:
## Exercise 4-5
## Make a copy of square and change the name to polygon. Add another parameter named n and modify the body so it draws an n-sided regular polygon.

## TIP
## The exterior angles of an n-sided regular polygon are 360/n degrees.
using ThinkJulia

println("Ans: ")

function polygon(turtle::Turtle, nsides::Int, len::Int = 100)
    angle = Int(360 / nsides)
    println("Angle: " * string(angle))

    for i in 1:nsides
        forward(turtle, len)
        turn(turtle, -angle)
    end
end

# (a) Try 1
# @svg begin
#     polygon(Turtle(), 5)
# end

# (b) Try 2
@svg begin
    polygon(Turtle(), 9)
end

println("End.")


text:
function init_constraint_struct(set::BoolSet{F1,F2}, internals) where {F1,F2}
    f = MOIU.eachscalar(internals.fct)

    lhs_fct = f[1:set.lhs_dimension]
    rhs_fct = f[end-set.rhs_dimension+1:end]

    if F1 <: MOI.ScalarAffineFunction
        lhs_fct = get_saf(lhs_fct)
    end
    if F2 <: MOI.ScalarAffineFunction
        rhs_fct = get_saf(rhs_fct)
    end

    if F1 <: MOI.VectorOfVariables
        lhs_fct = get_vov(lhs_fct)
    end
    if F2 <: MOI.VectorOfVariables
        rhs_fct = get_vov(rhs_fct)
    end

    lhs = get_constraint(lhs_fct, set.lhs_set)
    rhs = get_constraint(rhs_fct, set.rhs_set)

    return bool_constraint(set, internals, lhs, rhs)
end

function bool_constraint(::AndSet, internals, lhs, rhs)
    AndConstraint(
        internals,
        lhs,
        rhs
    )
end

function bool_constraint(::OrSet, internals, lhs, rhs)
    OrConstraint(
        internals,
        lhs,
        rhs
    )
end

"""
    activate_lhs!(com, constraint::BoolConstraint)

Activate the lhs constraint of `constraint` when not activated yet.
Saves at which stage it was activated.
"""
function activate_lhs!(com, constraint::BoolConstraint)
    lhs = constraint.lhs
    if !constraint.lhs_activated && lhs.impl.activate
        !activate_constraint!(com, lhs, lhs.fct, lhs.set) && return false
        constraint.lhs_activated = true
        constraint.lhs_activated_in_backtrack_idx = com.c_backtrack_idx
    end
    return true
end

"""
    activate_rhs!(com, constraint::BoolConstraint)

Activate the rhs constraint of `constraint` when not activated yet.
Saves at which stage it was activated.
"""
function activate_rhs!(com, constraint::BoolConstraint)
    rhs = constraint.rhs
    if !constraint.rhs_activated && rhs.impl.activate
        !activate_constraint!(com, rhs, rhs.fct, rhs.set) && return false
        constraint.rhs_activated = true
        constraint.rhs_activated_in_backtrack_idx = com.c_backtrack_idx
    end
    return true
end

function single_reverse_pruning_constraint!(
    com::CoM,
    constraint::BoolConstraint,
    fct::VAF{T},
    set::BoolSet,
    var::Variable,
    backtrack_idx::Int,
) where {
    T<:Real,
}
    for inner_constraint in (constraint.lhs, constraint.rhs)
        # the variable must be part of the inner constraint
        # Todo: Speed up the `in` here
        if inner_constraint.impl.single_reverse_pruning && var.idx in inner_constraint.indices
            single_reverse_pruning_constraint!(
                com,
                inner_constraint,
                inner_constraint.fct,
                inner_constraint.set,
                var,
                backtrack_idx,
            )
        end
    end
end

function reverse_pruning_constraint!(
    com::CoM,
    constraint::BoolConstraint,
    fct::VAF{T},
    set::BoolSet,
    backtrack_id::Int,
) where {
    T<:Real,
}
    # check if inner constraint should be deactived again
    if constraint.lhs_activated && backtrack_id == constraint.lhs_activated_in_backtrack_idx
        constraint.lhs_activated = false
        constraint.lhs_activated_in_backtrack_idx = 0
    end
    if constraint.rhs_activated && backtrack_id == constraint.rhs_activated_in_backtrack_idx
        constraint.rhs_activated = false
        constraint.rhs_activated_in_backtrack_idx = 0
    end

    for inner_constraint in (constraint.lhs, constraint.rhs)
        if inner_constraint.impl.reverse_pruning
            reverse_pruning_constraint!(
                com,
                inner_constraint,
                inner_constraint.fct,
                inner_constraint.set,
                backtrack_id,
            )
        end
    end
end

function restore_pruning_constraint!(
    com::CoM,
    constraint::BoolConstraint,
    fct::VAF{T},
    set::BoolSet,
    prune_steps::Union{Int,Vector{Int}},
) where {
    T<:Real,
}
    for inner_constraint in (constraint.lhs, constraint.rhs)
        if inner_constraint.impl.restore_pruning
            restore_pruning_constraint!(
                com,
                inner_constraint,
                inner_constraint.fct,
                inner_constraint.set,
                prune_steps,
            )
        end
    end
end

function finished_pruning_constraint!(
    com::CS.CoM,
    constraint::BoolConstraint,
    fct::VAF{T},
    set::BoolSet,
) where {
    T<:Real,
}
    for inner_constraint in (constraint.lhs, constraint.rhs)
        if inner_constraint.impl.finished_pruning
            finished_pruning_constraint!(
                com,
                inner_constraint,
                inner_constraint.fct,
                inner_constraint.set,
            )
        end
    end
end


text:
using Polynomials
using PolynomialEquations
using Test
using LinearAlgebra

@testset "Scaling of a polynomial" begin
    @test begin
        a = Polynomial(1:5,:s)
        ρ=10
        as = scale(a,ρ)
        r = Polynomial([1,20,300,4000,50000],:s)
        isequal(as,r)
    end
end

@testset "Shifting a Laurent polynomial" begin
    @test begin
        p = LaurentPolynomial([1,2,3,4],-2:1,:z)
        s = shift(p)
        r = LaurentPolynomial([1,2,3,4],-1:2,:z)
        isequal(s,r)
    end
    @test begin
        p = LaurentPolynomial([1,2,3,4],-2:1,:z)
        s = shift(p,-3)
        r = LaurentPolynomial([1,2,3,4],-5:-2,:z)
        isequal(s,r)
    end
end

@testset "Conjugate of a polynomial with respect to the imaginary axis" begin
    @test begin
        a = Polynomial(1:5,:s)
        p = cconj(a)
        r = Polynomial([1,-2,3,-4,5],:s)
        isequal(p,r)
    end
    @test begin
        a = Polynomial([1+1im, 2+2im, 3+3im, 4+4im, 5+5im],:s)
        p = cconj(a)
        r = Polynomial([1-1im, -(2-2im), 3-3im, -(4-4im), 5-5im],:s)
        isequal(p,r)
    end
end

@testset "Conjugate reciprocal polynomial" begin
    @test begin
        a = Polynomial(1:5,:s)
        p = conjreciprocal(a)
        r = Polynomial([5,4,3,2,1],:s)
        isequal(p,r)
    end
    @test begin
        c = Polynomial([1+1im, 2+2im, 3+3im, 4+4im, 5+5im],:s)
        p = conjreciprocal(c)
        r = Polynomial([5-5im, 4-4im, 3-3im, 2-2im, 1-1im],:s)
        isequal(p,r)
    end
end

@testset "Conjugate of a polynomial with respect to the unit circle" begin
    @test begin
        a = LaurentPolynomial(1:5,:z)
        d = dconj(a)
        r = LaurentPolynomial(5:-1:1,-4:0,:z)
        isequal(d,r)
    end
    @test begin
        a = LaurentPolynomial([1+1im, 2+2im, 3+3im, 4+4im, 5+5im],:z)
        d = dconj(a)
        r = LaurentPolynomial([5-5im, 4-4im, 3-3im, 2-2im, 1-1im],-4:0,:z)
        isequal(d,r)
    end
end

@testset "Building a lower triangular banded Toeplitz matrix for a polynomial" begin
    @test begin
        a = Polynomial([1.0,2.0,3.0,4.0],:s)
        A = ltbtmatrix(a,3)
        R = [1.0  0    0;
             2.0  1.0  0;
             3.0  2.0  1.0;
             4.0  3.0  2.0;
             0    4.0  3.0;
             0    0    4.0]
        isequal(A,R)
    end
    @test begin
        a = Polynomial([1,2,3,4],:s)
        A = ltbtmatrix(a,3)
        R = [1  0  0;
             2  1  0;
             3  2  1;
             4  3  2;
             0  4  3;
             0  0  4]
        isequal(A,R)
    end
end

@testset "Building a lower triangular banded Toeplitz matrix for a Laurent polynomial" begin
    @test begin
        a = LaurentPolynomial([1,2,3,4,5],-2:2,:z)
        A = ltbtmatrix(a,3)
        R = [1  0  0;
             2  1  0;
             3  2  1;
             4  3  2;
             5  4  3;
             0  5  4;
             0  0  5]
        isequal(A,R)
    end
end

@testset "Building a Sylvester resultant matrix" begin
    @test begin
        a = Polynomial([1,2,3],:s)
        b = Polynomial([4,5],:s)
        S = sylvestermatrix(a,b)
        R = [1  4  0;
             2  5  4;
             3  0  5]
        isequal(S,R)
    end
end

@testset "Building a Sylvester matrix" begin
    @test begin
        a = Polynomial([1.0,2.0,3.0],:s)
        b = Polynomial([4.0,5.0],:s)
        S = sylvestermatrix(a,b,(1,3))
        R = [1.0  4.0  0.0  0.0;
             2.0  5.0  4.0  0.0;
             3.0  0.0  5.0  4.0;
             0.0  0.0  0.0  5.0]
        isequal(S,R)
    end
    @test begin
        a = Polynomial([1.0,2.0,3.0],:s)
        b = Polynomial([4.0,5.0],:s)
        S = sylvestermatrix(a,b,(2,3))
        R = [1.0  0.0  4.0  0.0  0.0;
             2.0  1.0  5.0  4.0  0.0;
             3.0  2.0  0.0  5.0  4.0;
             0.0  3.0  0.0  0.0  5.0]
        isequal(S,R)
    end
end

@testset "Testing coprimeness of two polynomials" begin
    @test begin
        a = Polynomial([1,2,3],:s);
        b = Polynomial([1,2],:s);
        iscoprime(a,b)
    end
    @test begin
        a = Polynomial([1,2,3],:s);
        b = a*Polynomial([1,2],:s);
        !iscoprime(a,b)
    end
end

@testset "Testing Schur (discrete-time) stability of a polynomial" begin
    @test begin
        a = fromroots([-1/2,0, 2/3,-1/3im, 1/3im])
        isschurstable(a)
    end
    @test begin
        a = fromroots([2,0, 2/3,-3im, 3im])
        !isschurstable(a)
    end
end

@testset "Testing Hurwitz (continuous-time) stability of a polynomial" begin
    @test begin
        a = fromroots([-1, -2, -3+4im, -3-4im])
        ishurwitzstable(a)
    end
    @test begin
        a = fromroots([-1, -2, 3+4im, 3-4im])
        !ishurwitzstable(a)
    end
end

@testset "Solving the ax=b equations" begin
    @test begin
        a = Polynomial([1,2,3],:s);
        b = a*Polynomial([1,2],:s);
        x = axb(a,b)                # it should find the solution
        a*x≈b
    end
    @test begin
        a = Polynomial([1,2,3],:s);
        b = a*Polynomial([4,5,6],:s);
        x = axb(a,b)                # it should find the solution
        a*x≈b
    end
    @test begin
        a = Polynomial([1,2],:s)
        b = Polynomial([1,2,3],:s);
        x = axb(a,b)                # it should detect there is no solution
        isnothing(x)
    end
end

@testset "Solving the ax+by=0 equation" begin
    @test begin
        a = Polynomial([1,2,3],:s)
        b = Polynomial([4,5],:s)
        c = Polynomial([1,1],:s)
        a = a*c
        b = b*c
        x, y = axby0(a,b)
        norm(a*x+b*y)<1e-8
    end
end

@testset "Solving the ax+by=c equation" begin
    @test begin
        a = Polynomial([1,2,3],:s)
        b = Polynomial([4,5],:s)
        c = Polynomial([6],:s)
        x,y = axbyc(a,b,c)
        a*x+b*y≈c
    end
    @test begin
        a = Polynomial([1,2,3],:s)
        b = Polynomial([4,5],:s)
        c = Polynomial([6,7,8,9],:s)
        x,y = axbyc(a,b,c,deg=:miny)
        a*x+b*y≈c
    end
    @test begin
        a = Polynomial([1,2,3],:s)
        b = Polynomial([4,5],:s)
        c = Polynomial([6,7,8,9],:s)
        x,y = axbyc(a,b,c,deg=:minx)
        a*x+b*y≈c
    end
    @test begin
        a = Polynomial([1,2,3],:s)*Polynomial([1,1],:s)
        b = Polynomial([4,5],:s)*Polynomial([1,1],:s)
        c = Polynomial([6,7,8,9],:s)*Polynomial([1,1],:s)
        x,y = axbyc(a,b,c,deg=:miny)
        a*x+b*y≈c
    end
    @test begin
        a = Polynomial([1,2,3],:s)*Polynomial([1,1],:s)
        b = Polynomial([4,5],:s)*Polynomial([1,1],:s)
        c = Polynomial([6,7,8,9],:s)
        x,y = axbyc(a,b,c,deg=:miny)
        isnothing(x)
    end
end

@testset "Solving the ax̃+ãx=b+b̃ equation" begin
    @test begin
        a = Polynomial([-0.12, -0.29, 1],:s)
        b = Polynomial([1.86, -0.34, -1.14, -0.21, 1.19, -1.12],:s)
        x = axaxbb(a,b)
        cconj(a)*x+a*cconj(x)≈b+cconj(b)
    end
end

@testset "Solving the ax̃+b̃y=c+d̃ equation" begin
    @test begin
        a = Polynomial([-0.12, -0.29, 1],:s)
        b = Polynomial([1.86, -0.34, -1.14, -0.21, 1.19, -1.12],:s)
        c = Polynomial([1.2, 3.4, 5.6],:s)
        d = Polynomial([3.4, 5.6, 6.7, 7.8, 8.9, 9.1, 1.2],:s)
        x,y = axbycd(a,b,c,d)
        a*cconj(x)+cconj(b)*y≈c+cconj(d)
    end
end

@testset "Solving the equation ax+y=c while minimizing ||y||₁" begin
    @test begin
        a = Polynomial([-0.12, -0.29, 0.71, 2.79, 2.92, 1],:d)
        c = Polynomial([1.86, -0.34, -1.14, -0.21, 1.19, -1.12],:d)
        x, y, val = axycminl1(a,c)
        a*x+y≈c
    end
end

@testset "Computing a spectral factor of a polynomial" begin
    @test begin
        b = Polynomial([1.0,2.0,3.0],:s)*Polynomial([1,-1],:s)
        a = cconj(b)*b
        x = spectralfactor(a)
        s = Polynomial([1,3,5,3],:s)
        x ≈ s
    end
end

@testset "Computing a spectral factor of a Laurent polynomial" begin
    @test begin
        b = convert(LaurentPolynomial,fromroots([3,-1/2],var=:z))
        a = b*dconj(b)
        x = spectralfactor(a)
        s = LaurentPolynomial([-1/2,1/2,3],0:2,:z)
        x ≈ s
    end
end


text:
using MinimaxCalibratedEBayes
using RCall
using Plots
using StatsPlots
using LaTeXStrings
using ExponentialFamilies
using MosekTools
using Random

const MCEB = MinimaxCalibratedEBayes
pgfplotsx()

# # Prostate dáta analysis

# ## Preliminaries

# Let us make some analysis choices from the beginning

# ### Empirical Bayes Targets
# Our goal is to estimate the following two targets:

target_grid = -3:0.2:3
post_means = PosteriorMean.(StandardNormalSample.(target_grid))
lfsrs =  LFSR.(StandardNormalSample.(target_grid))

# ### Empirical Bayes Prior Class
# For the prior class we consider a mixture of Gaussians
# that have standard deviation 0.2 and with mixing measure
# supported on ... Note that we also specify the solver
# to be used for solving the convex programming problems;
# here we use Mosek.
gcal = GaussianMixturePriorClass(0.2, -3:0.02:3, Mosek.Optimizer)

# ### Discretization of marginal observations
marginal_grid = -6:0.02:6

# ## Load the data
prostz_file = MinimaxCalibratedEBayes.prostate_data_file
R"load($prostz_file)"
@rget prostz;
Zs = StandardNormalSample.(prostz)

# We will implement the preprocessing steps manually here.
# Afterwards, we will use wrapper that automate all of these steps.

# First we split our data into two folds
Random.seed!(200)
Zs_train, Zs_test, _, _ = MCEB._split_data(Zs)

# let us quickly check the range of values in the folds
extrema(response(Zs_train)), extrema(response(Zs_test))

# Now we use the first fold to estimate an $L_\infty$ neighborhood
# of the marginal density. We also fit a logspline model of the
# prior density to get our pilot estimators.

# ### Neighborhood construction:

infinity_band_options = KDEInfinityBandOptions(a_min=-6.0, a_max=6.0)
fkde_train = fit(infinity_band_options, response.(Zs_train))

# Let us visualize the result by showing a histogram of the datapoints
# the estimated density function and the $L_\infty$ bands.
prostate_marginal_plot = histogram(response.(Zs_train), bins=50, normalize=true,
          label="", alpha=0.6, fillcolor=:lightgray)
plot!(fkde_train, label="Kernel estimate \n  (and bands)",
      legend=:topright)

# ### G-model fitting

# We choose the following class
exp_family = ContinuousExponentialFamilyModel(Uniform(-3.6,3.6),
                         collect(-3.6:0.1:3.6); df=5, scale=true)
exp_family_solver = MCEB.ExponentialFamilyDeconvolutionMLE(cefm = exp_family, c0=0.001)

# Let us fit to the training data
exp_family_fit = fit(exp_family_solver, Zs_train)

# Let us look at the pilot estimates we got
lfsr_logspline_estimate = estimate.(lfsrs , exp_family_fit)
post_mean_logspline_estimate = estimate.(post_means , exp_family_fit)

expfamily_fits_plot = plot(plot(target_grid, lfsr_logspline_estimate),
                           plot(target_grid, post_mean_logspline_estimate))

# discretize our test samples to implement minimax calibration methodology
Zs_test_discr = DiscretizedStandardNormalSamples(Zs_test, marginal_grid)
Zs_test_discr = set_neighborhood(Zs_test_discr, fkde_train)

# add all of the above in one struct
mceb_setup = MinimaxCalibratorSetup(Zs_train = Zs_train,
                                    Zs_test = Zs_test,
						            prior_class = gcal,
						            fkde_train = fkde_train,
						            Zs_test_discr = Zs_test_discr,
						            pilot_method = exp_family_fit,
									delta_tuner =  MCEB.HalfCIWidth)#MCEB.RMSE)

post_means_fits = fit.(mceb_setup, post_means)
lfsrs_fits = fit.(mceb_setup, lfsrs)

gr()
prostate_postmeans_plot = plot(post_means_fits, ylims=(-2.5,2.5))
prostate_lfsr_plot = plot(lfsrs_fits, label=nothing)
plot(prostate_lfsr_plot, prostate_postmeans_plot, layout=(1,2))

# # The impact of Neighborhoods: Moving to opportunity

# We next seek to repeat the analysis in the second dataset. Here we will also
# see how most of the steps manually conducted for the prostate dataset
# can be automated.

# ## Loading the dataset
# Let us start by loading the dataset and filtering out missing observations.

using CSV
nbhood_file = MinimaxCalibratedEBayes.chetty_hendren_file
nbhoods_df = CSV.read(nbhood_file, types=Dict(3=>String));
missing_dx = [!ismissing(x) for x in nbhoods_df[!,:p25_coef]]
β_hat = [nbhoods_df[i,:p25_coef] for i in findall(missing_dx)]
se_hat = [nbhoods_df[i,:p25_se] for i in findall(missing_dx)]
nbhood_Zs = StandardNormalSample.(β_hat./se_hat);
length(nbhood_Zs)

# ## Setup the minimax calibrator
# The above steps repeat all the setup we conducted for the prostate data
# to compute the `mceb_setup` object.
Random.seed!(34)
mceb_options = MinimaxCalibratorOptions(prior_class = gcal,
                               marginal_grid = marginal_grid,
							   pilot_options = exp_family_solver)
nbhood_mceb_setup = fit(mceb_options, nbhood_Zs)

extrema(response.(nbhood_mceb_setup.Zs_train)),extrema(response.(nbhood_mceb_setup.Zs_test))

GR.settextfontprec(232, 3)
gr()
nbhood_post_means_fits = fit.(nbhood_mceb_setup, post_means)
nbhood_lfsrs_fits = fit.(nbhood_mceb_setup, lfsrs)

prostate_postmeans_plot = plot(post_means_fits, ylims=(-3.2,3.2),
                            fg_legend=:transparent, bg_legend=:transparent,
							legend = (0.35,0.84))
prostate_lfsrs_plot = plot(lfsrs_fits, label=nothing)

nbhood_postmeans_plot = plot(nbhood_post_means_fits,  ylims=(-3.2,3.2), label=nothing)
nbhood_lfsrs_plot = plot(nbhood_lfsrs_fits,  label=nothing)

prostate_marginal_plot = histogram(response.(Zs_train), bins=20, normalize=true,
          label="", alpha=0.4, fillcolor=:lightgray)
plot!(prostate_marginal_plot, fkde_train, fillalpha=0.45, #L"\bar{f}(x)\pm c_m",
      linewidth=1, label=nothing, #L"\bar{f}(x)\pm c_m",
	  legend=:topleft,  bg_legend=:transparent,
	  fg_legend=:transparent)

nbhood_marginal_plot = histogram(response.(nbhood_mceb_setup.Zs_train), bins=20, normalize=true,
	            label=nothing, alpha=0.4, fillcolor=:lightgray)
plot!(nbhood_marginal_plot, nbhood_mceb_setup.fkde_train, fillalpha=0.45,
                  label=nothing, linewidth=1)

plot()
using Plots.Measures
dataset_plot = plot(prostate_marginal_plot, prostate_postmeans_plot, prostate_lfsrs_plot,
	 nbhood_marginal_plot, nbhood_postmeans_plot, nbhood_lfsrs_plot,
	 title = ["a) Prostate" "b)" "c)" " d) Neighborhoods" "e)" "f)"],
	 layout= (2,3), size=(750,500))
dataset_plot

savefig(dataset_plot, "dataset_plots.pdf")
#savefig(dataset_plot, "mceb_datasets.pdf")


text:
function convert2d(loc)
    if loc == 1 || loc == 2 || loc == 3
        return (1,loc)
    elseif loc == 4 || loc == 5 || loc == 6
        return (2, loc-3)
    elseif loc == 7 || loc == 8 || loc == 9
        return (3, loc-6)
    end
end

function insertchar(char, r, c)
    board[r][c] = char
end

isspacefree(r, c) = board[r][c] == " "

function printboard(board)
    println("╔═══╦═══╦═══╗")
    println("║ "*board[1][1]*" ║ "*board[1][2]*" ║ "*board[1][3]*" ║")
    println("╠═══╬═══╬═══╣")
    println("║ "*board[2][1]*" ║ "*board[2][2]*" ║ "*board[2][3]*" ║")
    println("╠═══╬═══╬═══╣")
    println("║ "*board[3][1]*" ║ "*board[3][2]*" ║ "*board[3][3]*" ║")
    println("╚═══╩═══╩═══╝")
end

function countspaces(a)
    counter = 0
    for element in a
        if element == " "
            counter = counter + 1
        end
    end
    return counter
end

function isboardfull(board)
    if countspaces(board[1]) > 0 || countspaces(board[2]) > 0 || countspaces(board[3]) > 0
        return false
    else
        return true
    end
end

function player1smove()
    whethermove=true
    while whethermove == true
        println(" ")
        println("Player 1 Select position to place X (1-9)")
        movelocation = readline()
        try
            move = parse(Int64,movelocation)
            if move > 0 && move < 10
                if isspacefree(convert2d(move)[1], convert2d(move)[2]) == true
                    whethermove = false
                    return (convert2d(move)[1], convert2d(move)[2])

                else
                    println("This postion is already filled!")
                end
            else
                println("Please type a number between 1-9")
            end
        catch
            println("Please type a number!")
        end
    end
end

function player2smove()
    whethermove=true
    while whethermove == true
        println(" ")
        println("Player 2 Select position to place O (1-9)")
        movelocation = readline()
        try
            move = parse(Int64,movelocation)
            if move > 0 && move < 10
                if isspacefree(convert2d(move)[1], convert2d(move)[2]) == true
                    whethermove = false
                    return (convert2d(move)[1], convert2d(move)[2])
                else
                    println("This postion is already filled!")
                end
            else
                println("Please type a number between 1-9")
            end
        catch
            println("Please type a number!")
        end
    end
end

randommove(list) = rand(list, 1)

function pms(r)
    pm = []
    for i in [x for (x, char) in enumerate(board[r]) if char == " "]
        push!(pm, i+3*(r-1))
    end
    pm
end

function computersmove()
    possiblemoves = append!(append!(pms(1), pms(2)), pms(3)) #list of possible moves
    move = 0

    #see is there is a winning move first
    for element in ["O","X"]
        for i in possiblemoves
            #println(i)
            boardCopy = deepcopy(board)
            boardCopy[(convert2d(i)[1])][(convert2d(i)[2])] = element
            if winner(boardCopy, element) == true
                move = i
                return (convert2d(move[1])[1], convert2d(move[1])[2])
            end
        end
    end

    #see if there is a corner
    cornersavailable = []
    for i in possiblemoves
        if i in [1,3,7,9]
            push!(cornersavailable,i)
        end
    end

    if length(cornersavailable) > 0
        move = randommove(cornersavailable)
        return (convert2d(move[1])[1], convert2d(move[1])[2])
    end

    #see if there is centre
    if 5 in possiblemoves
        move = 5
        return (convert2d(move[1])[1], convert2d(move[1])[2])
    end

    #see if there is an edge
    edgesavailable = []
    for i in possiblemoves
        if i in [2,4,6,8]
            push!(edgesavailable,i)
        end
    end

    if length(edgesavailable) > 0
        move = randommove(edgesavailable)
        return (convert2d(move[1])[1], convert2d(move[1])[2])
    else
        return ()
    end
end

function winner(board,char)
    if board[1][1] == char && board[1][2] == char && board[1][3] == char
        return true
    elseif board[2][1] == char && board[2][2] == char && board[2][3] == char
        return true
    elseif board[3][1] == char && board[3][2] == char && board[3][3] == char
        return true
    elseif board[1][1] == char && board[2][1] == char && board[3][1] == char
        return true
    elseif board[1][2] == char && board[2][2] == char && board[3][2] == char
        return true
    elseif board[1][3] == char && board[2][3] == char && board[3][3] == char
        return true
    elseif board[1][1] == char && board[2][2] == char && board[3][3] == char
        return true
    elseif board[1][3] == char && board[2][2] == char && board[3][1] == char
        return true
    else
        return false
    end
end

function main_ai()
    while !(isboardfull(board))
        if !(winner(board, "O"))
            move = player1smove()
            insertchar("X", move[1], move[2])
            printboard(board)
        else
            println("O's win this time...")
            playagain()
        end
        if !(winner(board, "X"))
            move = computersmove()
            if length(move) == 0
                println(" ")
                println("Game is a Tie! No more spaces left to move.") #a tie must end on the ninth move, after the 4th computersmove (8th move in the game)
                playagain()
            else
                insertchar("O", move[1], move[2])
                println(" ")
                println("Computer placed an O in position ", string((move[2])+3*((move[1])-1)), ":")
                printboard(board)

            end
        else
            println("X's win, good job!")
            playagain()
        end
    end
end

function main_2p()
    while !(isboardfull(board))
        if !(winner(board, "O"))
            move = player1smove()
            insertchar("X", move[1], move[2])
            printboard(board)
        else
            println("Player 2 / O's win, good job!")
            playagain()
        end

        if isboardfull(board)
            println(" ")
            println("Game is a Tie! No more spaces left to move.")
            playagain()
        end

        if !(winner(board, "X"))
            move = player2smove()
            insertchar("O", move[1], move[2])
            printboard(board)
        else
            println("Player 1 / X's win, good job!")
            playagain()

        end
    end
end

function landing()
    global board = [[" " for i in 1:3] for k in 1:3]
    println("Welcome to Tic Tac Toe, to win complete a straight line of your letter (Diagonal, Horizontal, Vertical). The board has positions 1-9 starting at the top left.")
    printboard(board)
    version = lowercase(string(input("Please enter the number of players: 1 or 2\n")))
    while !(version == "1" || version == "2")
        version = lowercase(string(input("Please enter the number of players: 1 or 2\n")))
    end
    if version == "1"
        main_ai()
    elseif version == "2"
        main_2p()
    end
end

function playagain()
    playagain = lowercase(string(input("Play again? Y/N\n")))
    while !(playagain == "y" || playagain == "n")
        playagain = lowercase(string(input("Play again? Y/N\n")))
    end
    if playagain == "y"
        landing()
    elseif playagain == "n"
        exit()
    end
end

landing()

max_stars_repo_name: fungktt/JuliaGames

max_issues_repo_name: fungktt/JuliaGames

max_forks_repo_name: fungktt/JuliaGames


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -9.734460520574642 <= q <= 2.116911729869532)

                     add_NL_constraint(m, :(log(1 + exp(-0.8046949419530325 + -0.2278494590955038*log(1 + exp(0.1770285029406451 + -0.5067338360802731*log(1 + exp(-0.7998422425711689 + 0.5883237869246916*log(1 + exp(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5]))) + 0.47235910506966317*log(1 + exp(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))))) + -0.09889031570146667*log(1 + exp(-0.7622628423532047 + 0.6518912481623378*log(1 + exp(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5]))) + 0.4006903724184485*log(1 + exp(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))))))) + -0.18111143847071443*log(1 + exp(-0.12446319977672715 + 0.7018706953017659*log(1 + exp(-0.7998422425711689 + 0.5883237869246916*log(1 + exp(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5]))) + 0.47235910506966317*log(1 + exp(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))))) + 0.5275006439747312*log(1 + exp(-0.7622628423532047 + 0.6518912481623378*log(1 + exp(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5]))) + 0.4006903724184485*log(1 + exp(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))))))))) + log(1 + exp(-0.4025778812055263 + 0.08853333939955377*log(1 + exp(0.1770285029406451 + -0.5067338360802731*log(1 + exp(-0.7998422425711689 + 0.5883237869246916*log(1 + exp(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5]))) + 0.47235910506966317*log(1 + exp(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))))) + -0.09889031570146667*log(1 + exp(-0.7622628423532047 + 0.6518912481623378*log(1 + exp(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5]))) + 0.4006903724184485*log(1 + exp(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))))))) + -0.4269192713141794*log(1 + exp(-0.12446319977672715 + 0.7018706953017659*log(1 + exp(-0.7998422425711689 + 0.5883237869246916*log(1 + exp(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5]))) + 0.47235910506966317*log(1 + exp(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))))) + 0.5275006439747312*log(1 + exp(-0.7622628423532047 + 0.6518912481623378*log(1 + exp(-0.6653700708895021 + 0.6185402094345611*$(x[1]) + 0.05069108137655398*$(x[2]) + 0.9466113725109135*$(x[3]) + 0.06805104680691931*$(x[4]) + -0.30203159362192356*$(x[5]))) + 0.4006903724184485*log(1 + exp(-0.04665405044600712 + -0.8840756438370558*$(x[1]) + -0.20571795637475754*$(x[2]) + 0.8733859112268667*$(x[3]) + -0.5065045210605232*$(x[4]) + -0.7647800565294718*$(x[5]))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
export dmrg

function dmrg(H::MPO,
              psi0::MPS,
              sweeps::Sweeps;
              kwargs...)::Tuple{Float64,MPS}

  which_factorization::String = get(kwargs,:which_factorization,"automatic")

  psi = copy(psi0)
  N = length(psi)

  PH = ProjMPO(H)
  position!(PH,psi0,1)
  energy = 0.0

  quiet::Bool = get(kwargs,:quiet,false)
  for sw=1:nsweep(sweeps)
    sw_time = @elapsed begin

    @debug begin
      reset!(timer)
      println("\n\n~~~ This is sweep $sw ~~~\n\n")
    end

    for (b,ha) in sweepnext(N)
      position!(PH,psi,b)

      phi = psi[b]*psi[b+1]

      energy,phi = davidson(PH,phi;kwargs...)

      dir = ha==1 ? "fromleft" : "fromright"
      replaceBond!(psi,b,phi;
                   maxdim=maxdim(sweeps,sw),
                   mindim=mindim(sweeps,sw),
                   cutoff=cutoff(sweeps,sw),
                   dir=dir,
                   which_factorization=which_factorization)
    end
    end
    if !quiet
      @printf("After sweep %d energy=%.12f maxLinkDim=%d time=%.3f\n",sw,energy,maxLinkDim(psi),sw_time)
    end
  end
  return (energy,psi)
end

@doc """
dmrg(H::MPO,psi0::MPS,sweeps::Sweeps;kwargs...)::Tuple{Float64,MPS}

Optimize a matrix product state (MPS) to be the eigenvector
of the Hermitian matrix product operator (MPO) H with minimal
eigenvalue using the density matrix renormalization group
(DMRG) algorithm.

Inputs:
* `H::MPO` - a Hermitian MPO
* `psi0::MPS` - MPS used to initialize the optimization
* `sweeps::Sweeps` - `Sweeps` object used to control the algorithm

Returns:
* `energy::Float64` - eigenvalue of the optimized MPS
* `psi::MPS` - optimized MPS
""" dmrg


text:
"""
    SingularQuadratureRule{D,Q,S} <: AbstractQuadratureRule{D}

A quadrature rule over `D` intended to integrate functions which are singular at a point
`s ∈ D`.

A singular quadrature is rule is composed of a *regular* quadrature rule (e.g.
`GaussLegendre`) and a [`AbstractSingularityHandler`](@ref) to transform the regular
quadrature. The regular quadrature rule generates nodes and weights on the
`domain(sing_handler)`, and those are mapped into an appropriate quadrature over
`D = range(sing_handler)` using the singularity handler.

Besides the methods `(q::AbstractQuadratureRule)()` and
`(q::AbstractQuadratureRule)(el::AbstractElement)` described in
[`AbstractQuadratureRule`](@ref), singular quadrature objects also support
`(q::SingularQuadratureRule)(el::AbstractElement,s)`, which returns nodes `x`
and weights `w` for integrating a function over `el` with a possible
(integrable) singularity at parametric location `s ∈ D`. In practice this
typically means that the generated singular quadrature will accumulate points
near `s`, and the exact distribution of how the points are accumulated depends
on the `singularity_handler` employed.
"""
struct SingularQuadratureRule{D,Q,S} <: AbstractQuadratureRule{D}
    qrule::Q
    singularity_handler::S
    function SingularQuadratureRule(q::Q,s::S) where {Q,S}
        @assert domain(s) == domain(q) "domain of quadrature must coincide with domain of singularity handler"
        D = range(s) |> typeof
        new{D,Q,S}(q,s)
    end
end

# getters
qrule(q::SingularQuadratureRule) = q.qrule
singularity_handler(q::SingularQuadratureRule) = q.singularity_handler

@generated function (qs::SingularQuadratureRule{D,Q,S})() where {D,Q,S}
    qstd  = Q()
    shand = S()
    x̂,ŵ   = qstd() # reference quadrature
    # use the change of variables cov to modify reference quadrature
    x    = map(x->shand(x),x̂)
    w    = map(x̂,ŵ) do x̂,ŵ
        jac = jacobian(shand,x̂)
        μ   = abs(det(jac))
        μ*prod(ŵ)
    end
    return :($x,$w)
end

"""
    singular_quadrature(q::SingularQuadratureRule,s)

Use `q` to produce nodes `x` and weigths `w` for integrating a function over
`domain(q)` with a possible singularity at `s ∈ domain(q)`.

The actual implementation depends on the type of `q`.
"""
function singular_quadrature end

function singular_quadrature(q::SingularQuadratureRule{ReferenceLine},s)
    # remember that s is typically a static vector of length 1
    @assert 0 < s[1] < 1
    # qstd  = qrule(q)
    # shand = singularity_handler(q)
    # x̂,ŵ   = qstd()
    # x1    = map(x̂) do x
    #     x*(1-s) + s
    # end
    # x2    = map(x̂) do x
    #     x*(1-s) + s
    # end
    x̂,ŵ = q()
    # left domain
    x1  = @. s[1] * (1-x̂)
    w1  = @. ŵ*s[1]
    # right domain
    x2  = @. s[1]  + x̂*(1-s[1])
    w2  = @. ŵ*(1-s[1])
    # split the domain into two
    # l1    = line(s[1],0)
    # l2    = line(s[1],1)
    # # apply the quadrature to each segment
    # x1,w1 = q(l1)
    # x2,w2 = q(l2)
    # combine the nodes and weights
    return vcat(x1,x2), vcat(w1,w2)
end

"""
    singular_quadrature(k,q::SingularQuadratureRule,s)

Return nodes and weights to integrate a function over `domain(q)` with a
factored weight `k`.
"""
function singular_quadrature(k,q::SingularQuadratureRule,s)
    x,w = singular_quadrature(q,s)
    T = Base.promote_op(k,eltype(x))
    assert_concrete_type(T)
    w   = map(zip(x,w)) do (x,w)
        k(x)*w
    end
    return x,w
end

"""
    singular_weights(k,ui,q::SingularQuadratureRule,s)

Return weights to integrate a function over `domain(q)` with value known only at
nodes `ui` and a factored weight `k`.
"""
function singular_weights(k,ui,q::SingularQuadratureRule,s)
    x,w = singular_quadrature(k,q,s)
    x   = map(x->x[1],x)
    wlag = barycentric_lagrange_weights(ui)
    L    = barycentric_lagrange_matrix(ui,x,wlag)
    return transpose(L)*w
end

function singular_weights(k,qstd::AbstractQuadratureRule,q::SingularQuadratureRule,s)
    ui  = qnodes(qstd)
    ui  = map(x->x[1],ui)
    x,w = singular_quadrature(k,q,s)
    x   = map(x->x[1],x)
    wlag = barycentric_lagrange_weights(qstd)
    L    = barycentric_lagrange_matrix(ui,x,wlag)
    return transpose(L)*w
end

@generated function barycentric_lagrange_weights(Q::AbstractQuadratureRule)
    q  = Q()
    xs = qnodes(q) # static vector
    N  = length(xs)
    T  = eltype(xs) |> eltype
    x  = map(x->x[1],xs)
    # x  = Vector{T}(reinterpret(T,xs))
    w  =  barycentric_lagrange_weights(x)
    ws = SVector{N}(w)
    return :($ws)
end

"""
    singular_weights(q::SingularQuadratureRule,xi,k,s)

Return the weights to integrate `∫k(x)f(x)dx ≈ sum(f.(xi) .* w)` where a
singularity can be efficiently handled at location `s`.
"""
function singular_weights(qsin::SingularQuadratureRule,xi,k,s)
    x,w = qsin(k,s)
    wlag = barycentric_lagrange_weights(xi)
    L   = barycentric_lagrange_matrix(xi,x,wlag)
    return transpose(L)*w
end

# function singular_weights(qsin::SingularQuadratureRule,xi,K,s,el)
#     k   = (v) -> begin
#         jac = jacobian(el,v)
#         μ   = transpose(jac)*jac |> det |> sqrt
#         K(el(v))*μ
#     end
#     singular_weights(qsin,xi,k,s)
# end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
using Statistics

#=
Grading policy

Your overall score in this class will be a weighted average of your scores for the different components, with the following weights:

16% for the lecture exercises (divided equally among the 16 out of 19 lectures)
18% for the homeworks (divided equally among 6 (out of 7) homeworks)
2% for the Project 0
36% for the Projects (divided equally among 4 (out of 5)
12% for the Midterm exam (timed)
16% for the final exam (timed)

To earn a verified certificate for this course, you will need to obtain an overall score of 60% or more of the maximum possible overall score.

Lecture Exercises, Problem Sets, and Projects

    The lowest 3 scores among the 19 lectures will be dropped, so only 16 out of 19 lectures will count .

    The lowest 1 scores among the 7 homeworks will be dropped, so only 6 out of 7 homeworks will count .

    The lowest 1 score among the 5 projects (excluding project 0) will be dropped, so only 4 out of 5 projects will count .
=#

lectures_grades  = [0/11,0/18,0/12,0/11,0/9,0/13,0/9,0/14,0/5,0/18,0/18,0/11,0/24,0/10,0/24,0/20,0/20,0/9,0/8]
homeworks_grades = [0/61,0/28,0/12,0/25,0/26,0/32,0/36]
projects_grades  = [0/39,0/72,0/39,0/14,0/14] # exclude Project 0 here
project0_grade   = 0
midterm_grade    = 0/38
final_grade      = 0
lectures_grades_retained=sort(lectures_grades)[4:end]
homeworks_grades_retained=sort(homeworks_grades)[2:end]
projects_grades_retained=sort(projects_grades)[2:end]

grades = [mean(lectures_grades_retained),mean(homeworks_grades_retained),project0_grade,mean(projects_grades_retained),midterm_grade,final_grade]

w_lectures  = 0.16
w_homeworks = 0.18
w_project0  = 0.02
w_projects  = 0.36
w_midterm   = 0.12
w_final     = 0.16
weigths     = [w_lectures,w_homeworks,w_project0,w_projects,w_midterm,w_final]

overallGrade = grades' * weigths


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct FiveMassShaft <: Shaft
        H::Float64
        H_hp::Float64
        H_ip::Float64
        H_lp::Float64
        H_ex::Float64
        D::Float64
        D_hp::Float64
        D_ip::Float64
        D_lp::Float64
        D_ex::Float64
        D_12::Float64
        D_23::Float64
        D_34::Float64
        D_45::Float64
        K_hp::Float64
        K_ip::Float64
        K_lp::Float64
        K_ex::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int64
        internal::InfrastructureSystemsInternal
    end

Parameters of 5 mass-spring shaft model.
 It contains a High-Pressure (HP) steam turbine, Intermediate-Pressure (IP)
 steam turbine, Low-Pressure (LP) steam turbine, the Rotor and an Exciter (EX) mover.

# Arguments
- `H::Float64`: Rotor inertia constant in MWs/MVA, validation range: (0, nothing)
- `H_hp::Float64`: High pressure turbine inertia constant in MWs/MVA, validation range: (0, nothing)
- `H_ip::Float64`: Intermediate pressure turbine inertia constant in MWs/MVA, validation range: (0, nothing)
- `H_lp::Float64`: Low pressure turbine inertia constant in MWs/MVA, validation range: (0, nothing)
- `H_ex::Float64`:  Exciter inertia constant in MWs/MVA, validation range: (0, nothing)
- `D::Float64`: Rotor natural damping in pu, validation range: (0, nothing)
- `D_hp::Float64`: High pressure turbine natural damping in pu, validation range: (0, nothing)
- `D_ip::Float64`: Intermediate pressure turbine natural damping in pu, validation range: (0, nothing)
- `D_lp::Float64`: Low pressure turbine natural damping in pu, validation range: (0, nothing)
- `D_ex::Float64`: Exciter natural damping in pu, validation range: (0, nothing)
- `D_12::Float64`: High-Intermediate pressure turbine damping, validation range: (0, nothing)
- `D_23::Float64`: Intermediate-Low pressure turbine damping, validation range: (0, nothing)
- `D_34::Float64`: Low pressure turbine-Rotor damping, validation range: (0, nothing)
- `D_45::Float64`: Rotor-Exciter damping, validation range: (0, nothing)
- `K_hp::Float64`: High pressure turbine angle coefficient, validation range: (0, nothing)
- `K_ip::Float64`: Intermediate pressure turbine angle coefficient, validation range: (0, nothing)
- `K_lp::Float64`: Low pressure turbine angle coefficient, validation range: (0, nothing)
- `K_ex::Float64`: Exciter angle coefficient, validation range: (0, nothing)
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`
- `n_states::Int64`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct FiveMassShaft <: Shaft
    "Rotor inertia constant in MWs/MVA"
    H::Float64
    "High pressure turbine inertia constant in MWs/MVA"
    H_hp::Float64
    "Intermediate pressure turbine inertia constant in MWs/MVA"
    H_ip::Float64
    "Low pressure turbine inertia constant in MWs/MVA"
    H_lp::Float64
    " Exciter inertia constant in MWs/MVA"
    H_ex::Float64
    "Rotor natural damping in pu"
    D::Float64
    "High pressure turbine natural damping in pu"
    D_hp::Float64
    "Intermediate pressure turbine natural damping in pu"
    D_ip::Float64
    "Low pressure turbine natural damping in pu"
    D_lp::Float64
    "Exciter natural damping in pu"
    D_ex::Float64
    "High-Intermediate pressure turbine damping"
    D_12::Float64
    "Intermediate-Low pressure turbine damping"
    D_23::Float64
    "Low pressure turbine-Rotor damping"
    D_34::Float64
    "Rotor-Exciter damping"
    D_45::Float64
    "High pressure turbine angle coefficient"
    K_hp::Float64
    "Intermediate pressure turbine angle coefficient"
    K_ip::Float64
    "Low pressure turbine angle coefficient"
    K_lp::Float64
    "Exciter angle coefficient"
    K_ex::Float64
    ext::Dict{String, Any}
    states::Vector{Symbol}
    n_states::Int64
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function FiveMassShaft(H, H_hp, H_ip, H_lp, H_ex, D, D_hp, D_ip, D_lp, D_ex, D_12, D_23, D_34, D_45, K_hp, K_ip, K_lp, K_ex, ext=Dict{String, Any}(), )
    FiveMassShaft(H, H_hp, H_ip, H_lp, H_ex, D, D_hp, D_ip, D_lp, D_ex, D_12, D_23, D_34, D_45, K_hp, K_ip, K_lp, K_ex, ext, [:δ, :ω, :δ_hp, :ω_hp, :δ_ip, :ω_ip, :δ_lp, :ω_lp, :δ_ex, :ω_ex], 10, InfrastructureSystemsInternal(), )
end

function FiveMassShaft(; H, H_hp, H_ip, H_lp, H_ex, D, D_hp, D_ip, D_lp, D_ex, D_12, D_23, D_34, D_45, K_hp, K_ip, K_lp, K_ex, ext=Dict{String, Any}(), )
    FiveMassShaft(H, H_hp, H_ip, H_lp, H_ex, D, D_hp, D_ip, D_lp, D_ex, D_12, D_23, D_34, D_45, K_hp, K_ip, K_lp, K_ex, ext, )
end

# Constructor for demo purposes; non-functional.
function FiveMassShaft(::Nothing)
    FiveMassShaft(;
        H=0,
        H_hp=0,
        H_ip=0,
        H_lp=0,
        H_ex=0,
        D=0,
        D_hp=0,
        D_ip=0,
        D_lp=0,
        D_ex=0,
        D_12=0,
        D_23=0,
        D_34=0,
        D_45=0,
        K_hp=0,
        K_ip=0,
        K_lp=0,
        K_ex=0,
        ext=Dict{String, Any}(),
    )
end

"""Get FiveMassShaft H."""
get_H(value::FiveMassShaft) = value.H
"""Get FiveMassShaft H_hp."""
get_H_hp(value::FiveMassShaft) = value.H_hp
"""Get FiveMassShaft H_ip."""
get_H_ip(value::FiveMassShaft) = value.H_ip
"""Get FiveMassShaft H_lp."""
get_H_lp(value::FiveMassShaft) = value.H_lp
"""Get FiveMassShaft H_ex."""
get_H_ex(value::FiveMassShaft) = value.H_ex
"""Get FiveMassShaft D."""
get_D(value::FiveMassShaft) = value.D
"""Get FiveMassShaft D_hp."""
get_D_hp(value::FiveMassShaft) = value.D_hp
"""Get FiveMassShaft D_ip."""
get_D_ip(value::FiveMassShaft) = value.D_ip
"""Get FiveMassShaft D_lp."""
get_D_lp(value::FiveMassShaft) = value.D_lp
"""Get FiveMassShaft D_ex."""
get_D_ex(value::FiveMassShaft) = value.D_ex
"""Get FiveMassShaft D_12."""
get_D_12(value::FiveMassShaft) = value.D_12
"""Get FiveMassShaft D_23."""
get_D_23(value::FiveMassShaft) = value.D_23
"""Get FiveMassShaft D_34."""
get_D_34(value::FiveMassShaft) = value.D_34
"""Get FiveMassShaft D_45."""
get_D_45(value::FiveMassShaft) = value.D_45
"""Get FiveMassShaft K_hp."""
get_K_hp(value::FiveMassShaft) = value.K_hp
"""Get FiveMassShaft K_ip."""
get_K_ip(value::FiveMassShaft) = value.K_ip
"""Get FiveMassShaft K_lp."""
get_K_lp(value::FiveMassShaft) = value.K_lp
"""Get FiveMassShaft K_ex."""
get_K_ex(value::FiveMassShaft) = value.K_ex
"""Get FiveMassShaft ext."""
get_ext(value::FiveMassShaft) = value.ext
"""Get FiveMassShaft states."""
get_states(value::FiveMassShaft) = value.states
"""Get FiveMassShaft n_states."""
get_n_states(value::FiveMassShaft) = value.n_states
"""Get FiveMassShaft internal."""
get_internal(value::FiveMassShaft) = value.internal

"""Set FiveMassShaft H."""
set_H!(value::FiveMassShaft, val::Float64) = value.H = val
"""Set FiveMassShaft H_hp."""
set_H_hp!(value::FiveMassShaft, val::Float64) = value.H_hp = val
"""Set FiveMassShaft H_ip."""
set_H_ip!(value::FiveMassShaft, val::Float64) = value.H_ip = val
"""Set FiveMassShaft H_lp."""
set_H_lp!(value::FiveMassShaft, val::Float64) = value.H_lp = val
"""Set FiveMassShaft H_ex."""
set_H_ex!(value::FiveMassShaft, val::Float64) = value.H_ex = val
"""Set FiveMassShaft D."""
set_D!(value::FiveMassShaft, val::Float64) = value.D = val
"""Set FiveMassShaft D_hp."""
set_D_hp!(value::FiveMassShaft, val::Float64) = value.D_hp = val
"""Set FiveMassShaft D_ip."""
set_D_ip!(value::FiveMassShaft, val::Float64) = value.D_ip = val
"""Set FiveMassShaft D_lp."""
set_D_lp!(value::FiveMassShaft, val::Float64) = value.D_lp = val
"""Set FiveMassShaft D_ex."""
set_D_ex!(value::FiveMassShaft, val::Float64) = value.D_ex = val
"""Set FiveMassShaft D_12."""
set_D_12!(value::FiveMassShaft, val::Float64) = value.D_12 = val
"""Set FiveMassShaft D_23."""
set_D_23!(value::FiveMassShaft, val::Float64) = value.D_23 = val
"""Set FiveMassShaft D_34."""
set_D_34!(value::FiveMassShaft, val::Float64) = value.D_34 = val
"""Set FiveMassShaft D_45."""
set_D_45!(value::FiveMassShaft, val::Float64) = value.D_45 = val
"""Set FiveMassShaft K_hp."""
set_K_hp!(value::FiveMassShaft, val::Float64) = value.K_hp = val
"""Set FiveMassShaft K_ip."""
set_K_ip!(value::FiveMassShaft, val::Float64) = value.K_ip = val
"""Set FiveMassShaft K_lp."""
set_K_lp!(value::FiveMassShaft, val::Float64) = value.K_lp = val
"""Set FiveMassShaft K_ex."""
set_K_ex!(value::FiveMassShaft, val::Float64) = value.K_ex = val
"""Set FiveMassShaft ext."""
set_ext!(value::FiveMassShaft, val::Dict{String, Any}) = value.ext = val
"""Set FiveMassShaft states."""
set_states!(value::FiveMassShaft, val::Vector{Symbol}) = value.states = val
"""Set FiveMassShaft n_states."""
set_n_states!(value::FiveMassShaft, val::Int64) = value.n_states = val
"""Set FiveMassShaft internal."""
set_internal!(value::FiveMassShaft, val::InfrastructureSystemsInternal) = value.internal = val