text:
# DESCRIPTION
# This file contains the helicopter parameters.. Always re-run DI_design.jl
# after modifying the aircraft parameters.

# ------------------------- MISCELLANEOUS CONSTANTS ----------------------------

# conversion constants
D2R=pi/180
R2D=180/pi
g=32.17
FtLb_s2Hp=1/550
fps2kts=1/1.69
kts2fps=1.69
# air density [slug/ft³]
rhoSLSTD=0.002378
RHO=rhoSLSTD
# speed of sound [ft/s]
VSOUND=1135.0
# acceleration of gravity [ft/sec^2]
G=32.17

# -------------------------------- MAIN ROTOR ----------------------------------

# radius [ft]
R=26.83
# number of blade segments
NSEG=10
# number of blades
NB=4
# angular speed [rad/s]
OMEGA=27.0
# hinge offset [ft]
HOFFSET=1.25
# spar length (distance from hinge to start of rotor blade) [ft]
SPAR=2.25
# lift slope of airfoil [1/rad]
LIFTSLOPE=6.5
# zero lift angle of attack
ALPHA0=-2*pi/180.0
# profile drag coefficient
CD0=0.008
# hinge offset
e=HOFFSET/R
eprime=SPAR/R
# blade element locations
rseg=zeros(NSEG)
delseg=zeros(NSEG)
rseg[1]=sqrt((1-(e+eprime)^2)/(2*NSEG)+(e+eprime)^2)-e
delseg[1]=sqrt((rseg[1]+e)^2+(1-(e+eprime)^2)/(2*NSEG))-
          sqrt((rseg[1]+e)^2-(1-(e+eprime)^2)/(2*NSEG))
for iseg=2:NSEG
    rseg[iseg]=sqrt((1-(e+eprime)^2)/(NSEG)+(e+rseg[iseg-1])^2)-e
    delseg[iseg]=sqrt((rseg[iseg]+e)^2+(1-(e+eprime)^2)/(2*NSEG))-
                 sqrt((rseg[iseg]+e)^2-(1-(e+eprime)^2)/(2*NSEG))
end
RSEG = rseg
DELSEG = delseg
# dynamic twist mode shape
DynTwistMode = zeros(NSEG)
for iseg=1:NSEG
    DynTwistMode[iseg] = 0.28+(0.72*sin(0.5*pi*rseg[iseg]))
end
TauDT=0.05
# rigid blade prescribed twist
TWISTTABR = collect(0.0:0.05:1.0)
TWISTTABTHET= [0.0,0.0,0.0,0.0,-0.15,-0.95,-1.80,-2.75,-3.55,-4.40,-5.30,-6.15,
               -7.10,-7.90,-8.80,-9.65,-10.30,-10.75,-12.30,-13.10,-10.9]
# first mass moment of blade [sl-ft]
MBETA=86.7
# moment of intertia of blade about flap hinge [sl-ft²]
IBETA=1512.6
# blade weight [lbs]
WBLADE=256.9
# lag damper geometry [in]
ALD=0.227
BLD=3.242
CLD=12.04
DLD=10.0102
RLD=6.898
THETALDGEO=0.0
# shaft incidence
ITHSH=-3.0
IPHSH=0.0
Tshaft=[cos(ITHSH*D2R)                 0              -sin(ITHSH*D2R)
        sin(ITHSH*D2R)*sin(IPHSH*D2R)  cos(IPHSH*D2R)  cos(ITHSH*D2R)*sin(IPHSH*D2R)
        sin(ITHSH*D2R)*cos(IPHSH*D2R) -sin(IPHSH*D2R)  cos(ITHSH*D2R)*cos(IPHSH*D2R)]
# swashplate phasing angle [rad]
DELSP=-9.7*D2R
# δ₃ angle [rad]
DELTA3=0.0
# lag damping parameter [lbs/(ft/s)]
CLAG=1500.0
# lag damper stiffness [ft-lbs/rad]
KLAG=15060.0
# lag angle for 0 spring force [rad]
ZETA0=7*pi/180
# notional flap damper and spring
CFLAP=0.0
KFLAP=0.0
# blade tip chord [ft]
CHORDT=1.73
# blade root chord [ft]
CHORDR=1.73
# chord distribution
SEGOUT = zeros(NSEG)
SEGIN = zeros(NSEG)
CHORD = zeros(NSEG)
for iseg=1:NSEG
    SEGOUT[iseg]=sqrt(iseg*(1-(e + eprime)^2)/NSEG+(e + eprime)^2)
    SEGIN[iseg]=sqrt((iseg-1)*(1-(e + eprime)^2)/NSEG+(e + eprime)^2)
    CHORD[iseg] = (SEGOUT[iseg] + SEGIN[iseg] - 2*(e + eprime))*
                  (0.5/(1-e-eprime))*(CHORDT-CHORDR)+CHORDR;
end
# NACA 0012 airfloil data
fileIn = matopen("naca0012.mat")
ALTAB=read(fileIn, "ALTAB")
RETAB=read(fileIn, "RETAB")
CLTAB=read(fileIn, "CLTAB")
CDTAB=read(fileIn, "CDTAB")
close(fileIn)
MRAOATAB=ALTAB
MRCLTAB=CLTAB[:,5]
MRCDTAB=CDTAB[:,5]
# GenHel rotor airfoil data
fileIn = matopen("GHRotorData.mat")
CDR8UT=read(fileIn, "CDR8UT")
CDR8BT=read(fileIn, "CDR8BT")
CLR0BT=read(fileIn, "CLR0BT")
CLR8UT=read(fileIn, "CLR8UT")
CLR8BT=read(fileIn, "CLR8BT")
CDR0UT=read(fileIn, "CDR0UT")
CDR0BT=read(fileIn, "CDR0BT")
CLR0UT=read(fileIn, "CLR0UT")
close(fileIn)
CLR0UTAB=CLR0UT
AOAUTAB=collect(-180.0:2:180.0)
CLR0BTAB=reshape(CLR0BT,33,11)'
AOABTAB=collect(-32.0:2:32.0)
MACHTAB=collect(0.0:0.1:1.0)
CDR0UTAB=CDR0UT
CDR0BTAB=reshape(CDR0BT,33,11)'
ACL1 = 11.0
ACL2 = 172.0
ACL3 = -5.0
ACL4 = -172.0
DCDMR = 0.002
# wake curvature
KR = -1.0
# schedule from 1.5 in hover to 20, then 0 40 above
KRsched=[1.5, 1.5, 0.0, 0.0]
KRVsched=[0.0, 20.0, 40.0, 160.0]
# main rotor location (station, water line, butt line) [in]
FSMR=341.25
WLMR=315.0
BLMR=0.0

# ----------------------------- MASS PROPERTIES --------------------------------

# CG location (station, water line, butt line) [in]
FSCG=358.0
WLCG=248.2
BLCG=0.0
# aircraft gross weight (rotor blades not included) [lb]
WEIGHT=16000.0
WEIGHTNR = WEIGHT-NB*WBLADE;
MASS=WEIGHTNR/G;
# aircraft moment of inertia [lb-ft²]
IX=4659.0
IY=38512.0
IZ=36796.0
IXZ=1882.0

# --------------------------------- FUSELAGE -----------------------------------

# aerodynamic properties
fileIn = matopen("H60fuselage.mat")
# angle of attack
FUSEAOA=read(fileIn, "FUSEAOA")
# sideslip angle
FUSEBETA=read(fileIn, "FUSEBETA")
# sideslip angle absolute value
FUSEABETA=read(fileIn, "FUSEABETA")
# fuselage drag due to angle of attack [ft²]
FUSEDA=read(fileIn, "FUSEDA")
# fuselage drag due to sideslip [ft²]
FUSEDB=read(fileIn, "FUSEDB")
# fuselage siude force due to sideslip [ft²]
FUSEYB=read(fileIn, "FUSEYB")
# fuselage lift due to angle of attack [ft²]
FUSELA=read(fileIn, "FUSELA")
# fuselage lift due to sideslip [ft²]
FUSELB=read(fileIn, "FUSELB")
# fuselage roll moment due to sideslip [ft³]
FUSERB=read(fileIn, "FUSERB")
# fuselage pitch moment due to angle of attack [ft³]
FUSEMA=read(fileIn, "FUSEMA")
# fuselage pitch moment due to sideslip [ft³]
FUSEMB=read(fileIn, "FUSEMB")
# fuselage yaw moment due to sideslip [ft³]
FUSENB=read(fileIn, "FUSENB")
close(fileIn)
# location of fuselage aerodynamic center
FSfus=345.5
BLfus=0.0
WLfus=234.0
# fuselage CG location (without rotor) (station, water line, butt line) [in]
FSCGB=((WEIGHT*FSCG)-(NB*WBLADE*FSMR))/WEIGHTNR
WLCGB=((WEIGHT*WLCG)-(NB*WBLADE*WLMR))/WEIGHTNR
BLCGB=((WEIGHT*BLCG)-(NB*WBLADE*BLMR))/WEIGHTNR

# -------------------------------- TAIL ROTOR ----------------------------------

# location (station, water line, butt line) [in]
FSTR=732.0
BLTR=-14.0
WLTR=324.7
# cant angle
CantTR=20.0
# nominal RPS
OmegaTR=124.62
# chord
CHRDTR=0.81
# radius
RTR=5.5
# tip loss factor
BTLTR=0.92
# lift slope [1/rad]
a0TR=5.73
# moment of inertia about flapping hinge
IbTR=3.10
# twist
twistTR=-17.2*pi/180
# pitch bias
BIASTR = 5.0
# drag parameters
D0TR =  0.0087
D1TR = -0.0216
D2TR = 0.4
# TR airframe drag coefficient
CDTR=0.0
# tan(δ₃)
TD3TR=0.7002075
# change in coning with thrust
DELTTR = 0.001455
# blockage parameters with vertical tail
BVTTR  = 0.9
BVTTR1 = 1.0
VBVTTR = 30.0
# tail rotor gearing
GEARTR=OmegaTR/OMEGA

# ------------------- HORIZONTAL AND VERTICAL STABILIZERS ----------------------

# HORIZONTAL (idx=1) AND VERTICAL (idx=2) TAIL PROPERTIES
# tail area [ft²]
STAIL=zeros(2)
STAIL[1]=45.0
STAIL[2]=32.3
# cant of the tail (0 for horizontal, -90 for vertical)
PHITAIL=zeros(2)
PHITAIL[1]=0.0
PHITAIL[2]=-90.0
# pitch incidence of tail
ITAIL=zeros(2)
ITAIL[1]=0.0
ITAIL[2]=0.0
# location of tail aero center (station, water line, butt line) [in]
FSTAIL=zeros(2)
FSTAIL[1]=700.0
FSTAIL[2]=695.0
BLTAIL=zeros(2)
BLTAIL[1]=0.0
BLTAIL[2]=0.0
WLTAIL=zeros(2)
WLTAIL[1]=244.0
WLTAIL[2]=273.0
# aerodynamic data
ALTAIL=zeros(2,25)
ALTAIL[1,:]=[collect(-90.0:10:-40.0);collect(-30.0:5:30.0);collect(40.0:10:90.0)]
ALTAIL[2,:]=[collect(-90.0:10:-40.0);collect(-30.0:5:30.0);collect(40.0:10:90.0)]
CLTAIL=zeros(2,25)
CLTAIL[1,:]=[0, -0.294, -0.558, -0.745, -0.847, -0.847, -0.970, -1.030, -1.030,
             -0.930, -0.710, -0.356, 0, 0.356, 0.71, 0.93, 1.03, 1.03, 0.97,
              0.847, 0.847, 0.745, 0.558, 0.294, 0]
CLTAIL[2,:]=[0, -0.12, -0.28, -0.46, -0.66, -0.88, -1.0, -1.0, -0.93, -0.73,
             -0.5, -0.28, -0.06, 0.16, 0.38, 0.61, 0.82, 0.89, 0.89, 0.80, 0.63,
             0.48, 0.32, 0.17, 0.0]
CDTAIL=zeros(2,25)
CDTAIL[1,:]=[1.2000, 1.1610, 1.0500, 0.8880, 0.7020, 0.5310, 0.4300, 0.3700,
             0.3600, 0.1900, 0.0400, 0.0220, 0.01, 0.022, 0.04, 0.19, 0.36,
             0.37, 0.43, 0.531, 0.702, 0.888, 1.05, 1.161, 1.2]
CDTAIL[2,:]=[1.1, 1.025, 0.965, 0.875, 0.745, 0.575, 0.36, 0.265, 0.174, 0.118,
             0.066, 0.033, 0.018, 0.021, 0.044, 0.092, 0.162, 0.248, 0.355,
             0.58, 0.75, 0.875, 0.965, 1.02, 1.08]
# used for elevator control surface (not used here with all moving stab)
CLELEV=zeros(2)
CLELEV[1]=0.0
CLELEV[2]=0.0
# set at non-zero value to fix Horizontal stabilator at fixed position
# with value of 0 will use H60 schedule
STABSET = 0

# ------------------------------ CONTROL MIXING --------------------------------

# conversion from stick percentages to stick inches
perc2in = [0.1 0   0   0
           0   0.1 0   0
           0   0   0.1 0
           0   0   0   5.38/100]
# mixer gain matrix
LNKGAIN = [0.2062 0      0      0
           0      0.2172 0      0
           0      0      0.2025 0
           0      0      0      0.3780]
# convert stick inputs to servo commands
MIXGAIN = [0.6800  0      0.9930   0
           0       1.2020 0.9050  -0.4288
           0      -1.2020 1.2990   0.4288
           0       0      -0.8554  1.6043]
# gain for conversion servo to blade deflections
SWASHGAIN = [11.3413 -5.6706 -5.6706
              0      -5.6706  5.6706
              0       3.5034  3.5034]
TRGAIN = -9.142
# bias for conversion servo to blade deflections
SWASHBIAS = [-7.9738, 9.7280, 10.1111]
TRBIAS = 21.98

# ------------------------- AERODYNAMIC INTERACTION ----------------------------

# EK tables from GenHel
a1f_tab=[-6, 0, 6]
chi_tab=collect(0.0:10.0:100.0)
fileIn = matopen("AeroInteraction.mat")
FKXH1T=read(fileIn, "FKXH1T")
FKXWFT=read(fileIn, "FKXWFT")
FKZH1T=read(fileIn, "FKZH1T")
FKZWFT=read(fileIn, "FKZWFT")
close(fileIn)
ekxf_tab=reshape(FKXWFT,11,3)'
ekzf_tab=reshape(FKZWFT,11,3)'
ekxt_tab=reshape(FKXH1T,11,3)'
ekzt_tab=reshape(FKZH1T,11,3)'
alpha_tab=[-30, -25, -20, 10, 15, 20, 25, 30]*1.0
qlossht_tab=[1.0, 0.875, 0.76, 0.76, 0.82, 0.91, 1.0, 1.0]
psivt_tab=collect(-30.0:5:30.0)
qlossvt_tab=[1.0, 0.88, 0.79, 0.72, 0.66, 0.64, 0.62, 0.64, 0.66, 0.72, 0.79,
             0.88, 1.0]
alphaeps_tab=[collect(-90.0:10:-40.0);collect(-30.0:5:30.0);collect(40.0:10:90.0)]
eps_tab=[0, 0.25, 0.70, 1.2, 1.6, 1.9, 1.8, 1.4, 1.1, 0.8, 0.55, 0.5, 0.45, 0.4,
         0.38, 0.33, 0.19, -0.12, -0.4, -0.7, -0.75, -0.65, -0.45, -0.15, 0]
psisig_tab=[collect(-90.0:10:-40.0);collect(-30.0:5:30.0);collect(40.0:10:90.0)]
sig_tab=[0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 2.0, 2.0, 1.9, 1.8, 1.5, 0.88, 0.0,
         -0.88, -1.5, -1.8, -1.9, -2.0, -2.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0]

# --------------------- LINEARIZATION AND TRIM PARAMETERS ----------------------

# number of fuselage states
NFSTATES = 12
# number of rotor states
NRSTATES = 24
# number of tail rotor states
NTRSTATES = 1
# total number of states
NSTATES = NFSTATES+NRSTATES+NTRSTATES
# number of air frame states
NAFSTATES = NFSTATES+NRSTATES+NTRSTATES
# number of controls
NCTRLS = 4
# fuselage states index vector
IDXF=collect(1:NFSTATES)
# rotor states index vector
IDXR=collect(NFSTATES+1:NFSTATES+NRSTATES)
# tail rotor states index vector
IDXTR=collect(NFSTATES+NRSTATES+1:NFSTATES+NRSTATES+NTRSTATES)
# state scale factors
XSCALE=[1.0; 1.0; 1.0; (pi/180)*ones(6); 1.0; 1.0; 1.0; (pi/180)*ones(16);
        0.2/(OMEGA*R)*ones(3); pi/180; 1.0; 0.01; 0.01; 0.01;
        0.2/(OmegaTR*RTR)]
# perturbations for trim and linearization
DELCLIN=[1.0, 1.0, 1.0, 1.0]
DELXLIN=0.01*XSCALE
# define Trim Variables
TRIMVARS=[collect(1:8);IDXR[[collect(1:19);collect(21:24)]];IDXTR;
          collect(NSTATES+1:NSTATES+NCTRLS)]
TRIMTARG=[collect(1:12);IDXR[[collect(1:19);collect(21:24)]];IDXTR]
# number of azimuth locations to average trim over
NAZTRIM=10
# index of rotor azimuth (needed for azimuthal averaging in trimmer)
IDXAZ=32

max_stars_repo_name: umbertosaetti/J-GenHel

max_issues_repo_name: umbertosaetti/J-GenHel

max_forks_repo_name: umbertosaetti/J-GenHel


text:
function FieldQRDecomposition(arg0::FieldMatrix)
    return FieldQRDecomposition((FieldMatrix,), arg0)
end

function FieldQRDecomposition(arg0::FieldMatrix, arg1::RealFieldElement)
    return FieldQRDecomposition((FieldMatrix, RealFieldElement), arg0, arg1)
end

function get_h(obj::FieldQRDecomposition)
    return jcall(obj, "getH", FieldMatrix, ())
end

function get_q(obj::FieldQRDecomposition)
    return jcall(obj, "getQ", FieldMatrix, ())
end

function get_qt(obj::FieldQRDecomposition)
    return jcall(obj, "getQT", FieldMatrix, ())
end

function get_r(obj::FieldQRDecomposition)
    return jcall(obj, "getR", FieldMatrix, ())
end

function get_solver(obj::FieldQRDecomposition)
    return jcall(obj, "getSolver", FieldDecompositionSolver, ())
end

ext: jl


text:
using Raylib_jll

struct RayColor
    r::Cuchar
    g::Cuchar
    b::Cuchar
    a::Cuchar
end

const RAYWHITE = RayColor(245, 245, 245, 255)
const LIGHTGRAY = RayColor(200, 200, 200, 255)

function main()

    screenWidth = 800
    screenHeight = 450

    @ccall libraylib.InitWindow(screenWidth::Cint, screenHeight::Cint,
                             "raylib [core] example - basic window"::Cstring)::Cvoid

    @ccall libraylib.SetTargetFPS(60::Cint)::Cvoid

    while iszero(@ccall libraylib.WindowShouldClose()::Cint)
        @ccall libraylib.BeginDrawing()::Cvoid
        @ccall libraylib.ClearBackground(RAYWHITE::RayColor)::Cvoid
        @ccall libraylib.DrawText("Congrats! You created your first window!"::Cstring,
                               190::Cint, 200::Cint, 20::Cint, LIGHTGRAY::RayColor)::Cvoid
        @ccall libraylib.EndDrawing()::Cvoid
    end

    @ccall libraylib.CloseWindow()::Cvoid
end


text:
function solve_if_new!(Θ, last_Θ, U, X, V, conv_info, alg=Stocas.Poly())
    if Θ ≢ last_Θ
        copy!(last_Θ, Θ)
        copy!(U.Θ, Θ)
        update_utilities!(U)
        conv_info.iter += solve!(U, X, V, Stocas.Poly())
        true
    end
    false
end

function ll(Θ, last_Θ, U, X, D, E, conv_info)
    solve_if_new!(Θ, last_Θ, U, X, E.V, conv_info, Stocas.Poly())
    ll(U, D, E)
end
function ll(U, D, E)
    t_ll = 0.
    for ia = 1:length(E.P)
        map!(log, E.logP[ia], U.P[ia])
    end
    @inbounds for ia = 1:length(U.P), (x, a) in zip(D.x, D.a)
        t_ll -= E.logP[ia][x]*(a == ia)
    end
    t_ll/length(D.a)
end

function ∇ll!(grad, score, Θ, last_Θ, U, X, D, E, conv_info)
    solve_if_new!(Θ, last_Θ, U, X, E.V, conv_info, Stocas.Poly())
    ∇ll!(grad, score, U, X, D, E)
end
function ∇ll!(grad, score, U, X, D, E::NFXP)
    sc_logit!(score, U, X, D, E)
    grad[:] = mean(score,1)
end

function sc_logit!(score, U, X, D, E::NFXP)
    E.dVu .= (I-E.V.βFP)\sum(U.P[a].*U.Z[a] for a = 1:length(E.P))
    score_model = [U.Z[a] + U.β*(X.F[a]*E.dVu) - E.dVu for a = 1:length(E.P)]
    @inbounds for j = 1:size(score, 2)
        for i = 1:D.nobs
            score[i, j] = -score_model[D.a[i]][D.x[i],j]
        end
    end
end

function ∇²ll!(hessian, scores, Θ, last_Θ, U, X, D, E, conv_info)
    if solve_if_new!(Θ, last_Θ, U, X, E.V, conv_info, Stocas.Poly())
        sc_logit!(scores, U, X, D, E)
    end
    hessian .= scores'*scores/size(scores, 1)
end


text:
# File: fJBTest.jl
# Author:
# Date: 29-June-2020

# JARQUE-BERA TEST FOR NORMALITY
function fJBTest(x,sig,dof)
# -------------------------------------------------------------------
# x:         returns vector
# sig        significance level e.g. 95#  (0.95)
# dof        degrees of freedom e.g. 2
# jb:        Jarque-Bera value
# critval    inverse Chi-square at required significance level
# pval       Chi-square at Jarque-Bera value
# -------------------------------------------------------------------

    n = size(x,1); # # data points
    S = skewness(x,0); # sample skew
    K = kurtosis(x,0)-3; # sample excess kurtosis
    jb = (n/6)*(S^2+(K^2/4)); # Jarque-Bera value
    critval = chi2inv(sig,dof); # inverse Chi-square value
    pval = chi2pdf(jb,dof); # Chi-square at Jarque-Bera value

    return [jb, critval, pval]
end


text:
struct UniformScalingMap{T} <: LinearMap{T} # T will be determined from maps to which this is added
    λ::T
    M::Int
end
UniformScalingMap(λ::T, M::Int, N::Int) where {T} =
    (M == N ? UniformScalingMap(λ, M) : error("UniformScalingMap needs to be square"))
UniformScalingMap(λ::T, sz::Tuple{Int, Int}) where {T} =
    (sz[1] == sz[2] ? UniformScalingMap(λ, sz[1]) : error("UniformScalingMap needs to be square"))

# properties
Base.size(A::UniformScalingMap) = (A.M, A.M)
Base.isreal(A::UniformScalingMap) = isreal(A.λ)
LinearAlgebra.issymmetric(::UniformScalingMap) = true
LinearAlgebra.ishermitian(A::UniformScalingMap) = isreal(A)
LinearAlgebra.isposdef(A::UniformScalingMap) = isposdef(A.λ)

# special transposition behavior
LinearAlgebra.transpose(A::UniformScalingMap) = A
LinearAlgebra.adjoint(A::UniformScalingMap)   = UniformScalingMap(conj(A.λ), size(A))

# multiplication with vector
# call of LinearAlgebra.generic_mul! since order of arguments in mul! in stdlib/LinearAlgebra/src/generic.jl
# TODO: either leave it as is or use mul! (and lower bound on version) once fixed in LinearAlgebra
A_mul_B!(y::AbstractVector, A::UniformScalingMap, x::AbstractVector) =
    (length(x) == length(y) == A.M ? LinearAlgebra.generic_mul!(y, A.λ, x) : throw(DimensionMismatch("A_mul_B!")))
Base.:(*)(A::UniformScalingMap, x::AbstractVector) = A.λ * x

At_mul_B!(y::AbstractVector, A::UniformScalingMap, x::AbstractVector) =
    (length(x) == length(y) == A.M ? LinearAlgebra.generic_mul!(y, A.λ, x) : throw(DimensionMismatch("At_mul_B!")))

Ac_mul_B!(y::AbstractVector, A::UniformScalingMap, x::AbstractVector) =
    (length(x) == length(y) == A.M ? LinearAlgebra.generic_mul!(y, conj(A.λ), x) : throw(DimensionMismatch("Ac_mul_B!")))

# combine LinearMap and UniformScaling objects in linear combinations
Base.:(+)(A1::LinearMap, A2::UniformScaling{T}) where {T} = A1 + UniformScalingMap{T}(A2[1,1], size(A1, 1))
Base.:(+)(A1::UniformScaling{T}, A2::LinearMap) where {T} = UniformScalingMap{T}(A1[1,1], size(A2, 1)) + A2
Base.:(-)(A1::LinearMap, A2::UniformScaling{T}) where {T} = A1 - UniformScalingMap{T}(A2[1,1], size(A1, 1))
Base.:(-)(A1::UniformScaling{T}, A2::LinearMap) where {T} = UniformScalingMap{T}(A1[1,1], size(A2, 1)) - A2


text:
function DoubleDihedraFieldOfView(arg0::Vector3D, arg1::Vector3D, arg2::jdouble, arg3::Vector3D, arg4::jdouble, arg5::jdouble)
    return DoubleDihedraFieldOfView((Vector3D, Vector3D, jdouble, Vector3D, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5)
end

function get_footprint(obj::PolygonalFieldOfView, arg0::Transform, arg1::OneAxisEllipsoid, arg2::jdouble)
    return jcall(obj, "getFootprint", List, (Transform, OneAxisEllipsoid, jdouble), arg0, arg1, arg2)
end

function get_zone(obj::PolygonalFieldOfView)
    return jcall(obj, "getZone", SphericalPolygonsSet, ())
end

function offset_from_boundary(obj::PolygonalFieldOfView, arg0::Vector3D, arg1::jdouble, arg2::VisibilityTrigger)
    return jcall(obj, "offsetFromBoundary", jdouble, (Vector3D, jdouble, VisibilityTrigger), arg0, arg1, arg2)
end

function project_to_boundary(obj::PolygonalFieldOfView, arg0::Vector3D)
    return jcall(obj, "projectToBoundary", Vector3D, (Vector3D,), arg0)
end


text:
testOnAllContexts("GA") do ctx
    println("\tsolves a simple problem")
    testSimpleProblem(
        comparer ->
        begin
            ga = GA(ctx, 60, comparer, dna -> dna)
            @test ga.dnaSize == 60
            @test ga.populationSize == 100
            @test ga.mutationChance == 0.02f0
            @test ga.mutationStrength == 0.05f0
            @test ga.selectionStdDev == 0.25f0
            @test ga.keepElitesRate == 0.05f0
            ga
        end)

    println("\tsolves a problem with synthesizer")
    testSimpleProblem(
        ctx,
        (comparer, synth) ->
        begin
            ga = GA(
                ctx,
                dnaSize!(synth),
                comparer,
                dna -> decode(synth, dna),
                mutationChance = 0.05f0,
                selectionStdDev = 0.3f0)
            @test ga.dnaSize == dnaSize!(synth)
            @test ga.populationSize == 100
            @test ga.mutationChance == 0.05f0
            @test ga.mutationStrength == 0.05f0
            @test ga.selectionStdDev == 0.3f0
            @test ga.keepElitesRate == 0.05f0
            ga
        end)
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
using Plots, LinearAlgebra, Random, Statistics
pyplot()

# compute gradient
function df(X::Matrix{Float64}, y::Vector{Float64}, theta::Vector{Float64})
    return X' * (X * theta - y)
end

# gradient descent
function lr_gd(X::Matrix{Float64}, y::Vector{Float64}, alpha::Float64, num_iterations::Int64)
    (num_samples, num_features) = size(X)
    theta = zeros(num_features)
    for i = 1:num_iterations
        theta += df(X, y, theta)
    end
    return theta
end

# stochastic gradient descent, batch size of 1
function lr_sgd(X::Matrix{Float64}, y::Vector{Float64}, alpha::Float64, batch_size::Int64, num_iterations::Int64)
    (num_samples, num_features) = size(X)
    theta = zeros(num_features)
    for i = 1:num_iterations
        # use num_samples minibatches of size 1 to compute gradient to update theta
        # using k:k in order to keep type as a matrix
        theta += mean(df(X[k:k,:], y[k:k], theta) for k in (rand(1:num_samples) for b = 1:batch_size));
    end
    return theta
end
function main()
    num_features = 1024;
    dataset_sizes = 2 .^ [7,8,9,10,11,12,13,14,15,16];
    num_iterations = 1024;
    alpha = 0.001;
    # we will do this for the vanilla sgd case, without varying minibatch size
    batch_size = 1;
    gd_times = Float64[];
    sgd_times = Float64[];
    for num_samples in dataset_sizes
        X = randn(num_samples, num_features);
        y = randn(num_samples);
        gd_time = (@elapsed lr_gd(X, y, alpha, num_iterations)) / num_iterations;
        sgd_time = (@elapsed lr_sgd(X, y, alpha, batch_size, num_iterations)) / num_iterations;
        push!(gd_times, gd_time);
        push!(sgd_times, sgd_time);
        println("dataset size $num_samples --- GD time = $gd_time, SGD time = $sgd_time");
    end
    plot(dataset_sizes, gd_times, label="Gradient Descent (GD)", color=:blue,
        markershape=:circle, xaxis=:log, yaxis=:log);
    plot!(dataset_sizes, sgd_times, label="Stochastic Gradient Descent (SGD, with minibatch size=$batch_size)",
        markershape=:rect, color=:orange, marker=:, xaxis=:log, yaxis=:log,
        xticks = 10 .^ [2, 3, 4, 5, 6],
        xlabel="Number of data points", ylabel="Time per iteration (seconds)",
        title="Comparing Runtime per Iteration for GD and SGD");
    savefig("timesgd.pdf")
end
main()


text:
#####
# printing
#####

function env_traits()
    [eval(x) for x in RLBase.ENV_API if endswith(String(x), "Style")]
end

Base.show(io::IO, t::MIME"text/plain", env::AbstractEnv) =
    show(io, MIME"text/markdown"(), env)

function Base.show(io::IO, t::MIME"text/markdown", env::AbstractEnv)
    s = """
    # $(nameof(env))

    ## Traits
    | Trait Type | Value |
    |:---------- | ----- |
    $(join(["|$(string(f))|$(f(env))|" for f in env_traits()], "\n"))

    ## Is Environment Terminated?
    $(is_terminated(env) ? "Yes" : "No")

    """

    if get(io, :is_show_state_space, true)
        s *= """
        ## State Space
        `$(state_space(env))`

        """
    end

    if get(io, :is_show_action_space, true)
        s *= """
        ## Action Space
        `$(action_space(env))`

        """
    end

    if NumAgentStyle(env) !== SINGLE_AGENT
        s *= """
        ## Players
        $(join(["- `$p`" for p in players(env)], "\n"))

        ## Current Player
        `$(current_player(env))`
        """
    end

    if get(io, :is_show_state, true)
        s *= """
        ## Current State

        ```
        $(state(env))
        ```
        """
    end

    show(io, t, Markdown.parse(s))
end

#####
# testing
#####

using Test

"""
Call this function after writing your customized environment to make sure that
all the necessary interfaces are implemented correctly and consistently.
"""
function test_interfaces!(env)
    rng = Random.MersenneTwister(666)

    @info "testing $(nameof(env)), you need to manually check these traits to make sure they are implemented correctly!" NumAgentStyle(
        env,
    ) DynamicStyle(env) ActionStyle(env) InformationStyle(env) StateStyle(env) RewardStyle(
        env,
    ) UtilityStyle(env) ChanceStyle(env)

    @testset "copy" begin
        X = copy(env)
        Y = copy(env)
        reset!(X)
        reset!(Y)

        if ChanceStyle(Y) ∉ (DETERMINISTIC, EXPLICIT_STOCHASTIC)
            s = 888
            Random.seed!(Y, s)
            Random.seed!(X, s)
        end

        @test Y !== X

        @test state(Y) == state(X)
        @test action_space(Y) == action_space(X)
        @test reward(Y) == reward(X)
        @test is_terminated(Y) == is_terminated(X)

        while !is_terminated(Y)
            A, A′ = legal_action_space(X), legal_action_space(Y)
            @test A == A′
            a = rand(rng, A)
            Y(a)
            X(a)
            @test state(Y) == state(X)
            @test reward(Y) == reward(X)
            @test is_terminated(Y) == is_terminated(X)
        end
    end

    @testset "SingleAgent" begin
        if NumAgentStyle(env) === SINGLE_AGENT
            reset!(env)
            total_reward = 0.0
            while !is_terminated(env)
                if StateStyle(env) isa Tuple
                    for ss in StateStyle(env)
                        @test state(env, ss) ∈ state_space(env, ss)
                    end
                end

                A = legal_action_space(env)
                if ActionStyle(env) === MINIMAL_ACTION_SET
                    action_space(env) == legal_action_space
                elseif ActionStyle(env) === FULL_ACTION_SET
                    # @test legal_action_space(env) == action_space(env)[legal_action_space_mask(env)]
                else
                    @error "TODO:"
                end

                a = rand(rng, A)
                @test a ∈ action_space(env)
                env(a)
                @test state(env) ∈ state_space(env)

                total_reward += reward(env)
            end
            r = reward(env)  # make sure we can still get the reward no matter what the RewardStyle of the env is.
            if RewardStyle(env) === TERMINAL_REWARD
                @test total_reward == r
            end
        end
    end

    @testset "MultiAgent" begin
        if NumAgentStyle(env) isa MultiAgent
            reset!(env)
            rewards = [0.0 for p in players(env)]
            while !is_terminated(env)
                if InformationStyle(env) === PERFECT_INFORMATION
                    for p in players(env)
                        @test state(env) == state(env, p)
                    end
                end
                a = rand(rng, legal_action_space(env))
                env(a)
                for (i, p) in enumerate(players(env))
                    @test state(env, p) ∈ state_space(env, p)
                    rewards[i] += reward(env, p)
                end
            end
            # even the game is already terminated
            # make sure each player can still get some necessary info
            for p in players(env)
                state(env, p)
                reward(env, p)
                @test is_terminated(env)
                # we don't need to check legal_action_space when the game is
                # already over
                # @test isempty(legal_action_space(env, p))
            end
            if RewardStyle(env) === TERMINAL_REWARD
                for (p, r) in zip(players(env), rewards)
                    @test r == reward(env, p)
                end
            end

            if UtilityStyle(env) === ZERO_SUM
                @test sum(rewards) == 0
            elseif UtilityStyle(env) == IDENTICAL_UTILITY
                @test all(rewards[1] .== rewards)
            end
        end
    end

    reset!(env)
end

function test_runnable!(env, n = 1000; rng = Random.GLOBAL_RNG)
    @testset "random policy with $(nameof(env))" begin
        reset!(env)
        for _ in 1:n
            A = legal_action_space(env)
            a = rand(rng, A)
            @test a in A

            if ActionStyle(env) === EXPLICIT_STOCHASTIC &&
               current_player(env) == chance_player(env)
                @test isapprox(sum(prob(env)), 1)
            end

            S = state_space(env)
            s = state(env)
            @test s in S
            env(a)
            if is_terminated(env)
                reset!(env)
            end
        end
        reset!(env)
    end
end


text:
export preprocessLDLt, preprocessLDLt2
function preprocessLDLt(H ,g, params :: Tparams, n1, n2)
    # printstyled("dans preprocessLDLt ⤈ \n", color = :red)
    T = eltype(H)
    # @show T
    # @show H[1, 1]
    # @show H[1, 2]
    # @show H[2, 1]
    # @show H[2, 2]
    # @show g[1]
    # @show g[2]
    n, = size(g)
    # @show n
    global L = Matrix{T}(undef, n, n)
    global D = Matrix{T}(undef, n, n)
    global pp = Vector{Int}(undef, n)
    global ρ = nothing
    global ncomp = nothing

    # @show H

    try
        (L, D, pp, ρ, ncomp) = ldlt_symm(H, 'r')
    catch
 	println("*******   Problem in LDLt")
        res = PDataLDLt()
        res.OK = false
        return res
    end

    # printstyled("dans preprocessLDLt apres ldlt_symm eltype(L) = $(eltype(L)) \n", color = :cyan)
    # @show L
    # @show D
    # @show L[1, 1]
    # @show L[2, 1]
    # @show L[1, 2]
    # @show L[2, 2]
    # @show D[1, 1]
    # @show D[2, 1]
    # @show D[1, 2]
    # @show D[2, 2]
    # @show pp
    # @show ρ
    # @show ncomp

    if true in isnan.(D)
 	println("*******   Problem in D from LDLt: NaN")
        println(" cond (H) = $(cond(H))")
        res = PDataLDLt()
        res.OK = false
        return res
    end

    # @show L
    # @show D
    # @show pp

    # D = Hermitian(D)
    # @show D
    X = eigen(SymTridiagonal(D))
    Δ = X.values
    # @show Δ[1]
    # @show Δ[2]
    Q =  X.vectors
    # @show Q[1, 1]
    # @show Q[2, 1]
    # @show Q[1, 2]
    # @show Q[2, 2]
    l_m, = findmin(Δ)
    # @show l_m
    ĝ = L \ (g[pp])
    # @show ĝ[1]
    # @show ĝ[2]
    # @show eltype(ĝ)
    # @show ĝ
    g̃ = Q' * ĝ
    # @show g̃[1]
    # @show g̃[2]
    # @show eltype(g̃)
    # @show g̃

    n_g = norm(g)
    # @show n_g
    λ =  max(-l_m, 0.0)
    # @show eltype(λ)
    # @show eltype(λ)
    # @show λ
    # printstyled("avant de sortir de PreProcessLDLt ↑ \n", color = :red)
    return  PDataLDLt(L, D, pp, Δ, Q, g̃, λ, true, true)
end

function preprocessLDLt2(H ,g, params :: Tparams, n1, n2)
    # printstyled("dans preprocessLDLt \n", color = :red)
    T = eltype(H)
    # @show H
    # @show g
    n, = size(g)
    # @show n
    global L = Matrix{T}(undef, n, n)
    global D = Matrix{T}(undef, n, n)
    global pp = Vector{Int}(undef, n)
    global ρ = nothing
    global ncomp = nothing
    global S = nothing

    # @show H
    # printstyled("ici ↩ \n", color = :red)

    try
        # (L, D, pp, ρ, ncomp) = ldlt_symm(H, 'r')
        # @show H
        S = bunchkaufman(Symmetric(H, :L), true, check = false)
    catch
 	println("*******   Problem in LDLt")
        res = PDataLDLt()
        res.OK = false
        return res
    end
    L  = S.L
    D  = S.D
    pp = S.p

    # printstyled("dans preprocessLDLt apres ldlt_symm eltype(L) = $(eltype(L)) \n", color = :cyan)
    # @show L
    # @show D
    # @show pp
    # @show ρ
    # @show ncomp

    if true in isnan.(D)
 	println("*******   Problem in D from LDLt: NaN")
        println(" cond (H) = $(cond(H))")
        res = PDataLDLt()
        res.OK = false
        return res
    end

    # @show L
    # @show D
    # @show pp

    D = Hermitian(D)
    X = eigen(Matrix(D))
    Δ = X.values
    # @show Δ
    Q =  X.vectors
    l_m, = findmin(Δ)
    # @show l_m
    ĝ = L \ (g[pp])
    # @show ĝ[1]
    # @show ĝ[2]
    # @show eltype(ĝ)
    g̃ = Q' * ĝ
    # @show g̃[1]
    # @show g̃[2]
    # @show eltype(g̃)

    n_g = norm(g)
    # @show n_g
    λ =  max(-l_m, 0.0)
    # @show λ
    # printstyled("avant de sortir de PreProcessLDLt ↑ \n", color = :red)
    return  PDataLDLt(L, D, pp, Δ, Q, g̃, λ, true, true)
end

function AInv(X :: PDataLDLt, d̃ ::  Array{T,1}) where T
    d̂ = X.Q * d̃
    u = X.L' \ d̂
    # @show u[invperm(X.pp)][1]
    # @show u[invperm(X.pp)][2]
    # printstyled("avant la sortie de AInv ⟰ \n", color = :red)
    return u[invperm(X.pp)]
end

function reconstructH(X :: PDataLDLt)
    A = X.L * X.D * X.L'
    return A[X.pp.X.pp]
end


text:
using Julog
using Test

include("parse.jl")
include("unify.jl")
include("builtins.jl")
include("ext_builtins.jl")
include("math_ops.jl")
include("lists.jl")
include("custom_funcs.jl")
include("prolog.jl")
include("conversions.jl")
include("lineage.jl")
include("naturals.jl")
include("permutations.jl")


text:
using Test, Sockets, LanguageServer, CSTParser, SymbolServer, SymbolServer.Pkg, StaticLint, LanguageServer.JSON, LanguageServer.JSONRPC
using LanguageServer:
    Document,
    TextDocument,
    apply_text_edits,
    applytextdocumentchanges,
    get_text_document,
    set_text_document!,
    get_text,
    get_offset,
    get_line_offsets,
    get_position_from_offset,
    get_open_in_editor,
    set_open_in_editor,
    is_workspace_file,
    get_uri
using LanguageServer.URIs2
const LS = LanguageServer
const Range = LanguageServer.Range

# TODO Replace this with a proper mock endpoint
JSONRPC.send(::Nothing, ::Any, ::Any) = nothing
function settestdoc(text)
    empty!(server._documents)
    LanguageServer.textDocument_didOpen_notification(LanguageServer.DidOpenTextDocumentParams(LanguageServer.TextDocumentItem(uri"untitled:testdoc", "julia", 0, text)), server, nothing)

    doc = LanguageServer.getdocument(server, uri"untitled:testdoc")
    LanguageServer.parse_all(doc, server)
    doc
end

function on_all_docs(server, f)
    for doc in values(server._documents)
        f(doc)
    end
end

function on_all_offsets(doc, f)
    offset = 1
    while offset <= lastindex(get_text(doc))
        f(doc, offset)
        offset = nextind(get_text(doc), offset)
    end
end

@testset "LanguageServer" begin
    @testset "document" begin
        include("test_document.jl")
    end
    @testset "communication" begin
        include("test_communication.jl")
    end
    @testset "intellisense" begin
        include("test_intellisense.jl")

        server = LanguageServerInstance(IOBuffer(), IOBuffer(), dirname(Pkg.Types.Context().env.project_file), first(Base.DEPOT_PATH))
        server.runlinter = true
        server.jr_endpoint = nothing
        LanguageServer.initialize_request(init_request, server, nothing)

        @testset "completions" begin
            include("requests/completions.jl")
        end
        @testset "actions" begin
            include("requests/actions.jl")
        end
        @testset "features" begin
            include("requests/features.jl")
        end
        @testset "hover" begin
            include("requests/hover.jl")
        end
        @testset "textdocument" begin
            include("requests/textdocument.jl")
        end
        @testset "misc" begin
            include("requests/misc.jl")
        end
        @testset "brute force tests" begin
            @info "Self-parse test"
            if get(ENV, "CI", false) != false
                @info "skipping brute-force tests on CI"
            else
                # run tests against each position in each document
                empty!(server._documents)
                LanguageServer.load_folder(dirname(String(first(methods(LanguageServer.eval)).file)), server)
                on_all_docs(server, doc -> begin
                    @info "Testing LS functionality at all offsets" file=get_uri(doc)
                    on_all_offsets(doc, function (doc, offset)
                        tdi = LanguageServer.TextDocumentIdentifier(get_uri(doc))
                        pos = LanguageServer.Position(LanguageServer.get_position_from_offset(doc, offset)...)
                        @test LanguageServer.get_offset(doc, LanguageServer.get_position_from_offset(doc, offset)...) == offset
                        LanguageServer.textDocument_completion_request(LanguageServer.CompletionParams(tdi, pos, missing), server, server.jr_endpoint)
                        LanguageServer.textDocument_hover_request(LanguageServer.TextDocumentPositionParams(tdi, pos), server, server.jr_endpoint)
                        LanguageServer.textDocument_signatureHelp_request(LanguageServer.TextDocumentPositionParams(tdi, pos), server, server.jr_endpoint)
                        LanguageServer.textDocument_definition_request(LanguageServer.TextDocumentPositionParams(tdi, pos), server, server.jr_endpoint)
                        LanguageServer.textDocument_references_request(LanguageServer.ReferenceParams(tdi, pos, missing, missing, LanguageServer.ReferenceContext(true)), server, server.jr_endpoint)
                        LanguageServer.textDocument_rename_request(LanguageServer.RenameParams(tdi, pos, missing, "newname"), server, server.jr_endpoint)
                    end)
                end)

                on_all_docs(server, doc -> begin
                    symbols=length(LanguageServer.textDocument_documentSymbol_request(LanguageServer.DocumentSymbolParams(LanguageServer.TextDocumentIdentifier(get_uri(doc)),missing, missing), server, server.jr_endpoint))
                    @info "Found $symbols symbols" file=get_uri(doc)
                end)

                LanguageServer.workspace_symbol_request(LanguageServer.WorkspaceSymbolParams("", missing, missing), server, server.jr_endpoint)
            end
        end
    end
    @testset "edit" begin
        include("test_edit.jl")
    end
    @testset "paths" begin
        include("test_paths.jl")
    end

    include("test_uris2.jl")

    @testset "misc" begin
        include("test_misc.jl")
    end
end


text:
# Empirical estimation of CDF and PDF

## Empirical CDF

struct ECDF{T <: AbstractVector{<:Real}, W <: AbstractWeights{<:Real}}
    sorted_values::T
    weights::W
end

function (ecdf::ECDF)(x::Real)
    n = searchsortedlast(ecdf.sorted_values, x)
    evenweights = isempty(ecdf.weights)
    weightsum = evenweights ? length(ecdf.sorted_values) : sum(ecdf.weights)
    partialsum = evenweights ? n : sum(view(ecdf.weights, 1:n))
    partialsum / weightsum
end

function (ecdf::ECDF)(v::RealVector)
    evenweights = isempty(ecdf.weights)
    weightsum = evenweights ? length(ecdf.sorted_values) : sum(ecdf.weights)
    ord = sortperm(v)
    m = length(v)
    r = similar(ecdf.sorted_values, m)
    r0 = zero(weightsum)
    i = 1
    for (j, x) in enumerate(ecdf.sorted_values)
        while i <= m && x > v[ord[i]]
            r[ord[i]] = r0
            i += 1
        end
        r0 += evenweights ? 1 : ecdf.weights[j]
        if i > m
            break
        end
    end
    while i <= m
        r[ord[i]] = weightsum
        i += 1
    end
    return r / weightsum
end

"""
    ecdf(X; weights::AbstractWeights)

Return an empirical cumulative distribution function (ECDF) based on a vector of samples
given in `X`. Optionally providing `weights` returns a weighted ECDF.

Note: this function that returns a callable composite type, which can then be applied to
evaluate CDF values on other samples.

`extrema`, `minimum`, and `maximum` are supported to for obtaining the range over which
function is inside the interval ``(0,1)``; the function is defined for the whole real line.
"""
function ecdf(X::RealVector; weights::AbstractVector{<:Real}=Weights(Float64[]))
    isempty(weights) || length(X) == length(weights) || throw(ArgumentError("data and weight vectors must be the same size," *
        "got $(length(X)) and $(length(weights))"))
    ord = sortperm(X)
    ECDF(X[ord], isempty(weights) ? weights : Weights(weights[ord]))
end

minimum(ecdf::ECDF) = first(ecdf.sorted_values)

maximum(ecdf::ECDF) = last(ecdf.sorted_values)

extrema(ecdf::ECDF) = (minimum(ecdf), maximum(ecdf))


text:
# coding:utf-8

doc =
"""
usage:
  cabocha_parser("クロールで泳いでいる少女を見た")
"""

""" cabocha_parser """
function cabocha_parser(sentence, option)
  return readall(pipeline(`echo $sentence`, `cabocha -$option`))
end

function cabocha_parser(sentence)
  sentence = replace(sentence, "\n", "")
  result = readall(pipeline(`echo $sentence`, `cabocha -f1`))
  d = Dict(:chunk_id => Array{Int, 1}(), :chunk => Array{UTF8String,1}(), :tok_surface => Array{Array{UTF8String,1},1}(),
           :tok_feature => Array{Array{UTF8String,1},1}(), :link => Array{Int, 1}(), :head => Array{Int, 1}(), :func => Array{Int, 1}())
  for line in split(result, "\n")
    line == "EOS" && break
    if ismatch(r"^\* \d+", line)
     tmp = split(line)
     head, func = map(parse, split(tmp[4], "/"))
     push!(d[:chunk_id], parse(tmp[2])+1)
     push!(d[:chunk], "")
     push!(d[:link], parse(tmp[3][1:end-1])+1)
     push!(d[:head], head+1)
     push!(d[:func], func+1)
     push!(d[:tok_surface], Array{UTF8String,1}())
     push!(d[:tok_feature], Array{UTF8String,1}())
    else
     surface, feature = split(line, "\t")
     d[:chunk][end] = d[:chunk][end]*surface
      push!(d[:tok_surface][end], surface)
      push!(d[:tok_feature][end], feature)
    end
  end
  return d
end


text:
"""
    MaternKernel(; ν::Real=1.5)

Matérn kernel of order `ν`.

# Definition

For inputs ``x, x' \\in \\mathbb{R}^d``, the Matérn kernel of order ``\\nu > 0`` is
defined as
```math
k(x,x';\\nu) = \\frac{2^{1-\\nu}}{\\Gamma(\\nu)}\\big(\\sqrt{2\\nu}\\|x-x'\\|_2\\big) K_\\nu\\big(\\sqrt{2\\nu}\\|x-x'\\|_2\\big),
```
where ``\\Gamma`` is the Gamma function and ``K_{\\nu}`` is the modified Bessel function of
the second kind of order ``\\nu``.

A Gaussian process with a Matérn kernel is ``\\lceil \\nu \\rceil - 1``-times
differentiable in the mean-square sense.

See also: [`Matern12Kernel`](@ref), [`Matern32Kernel`](@ref), [`Matern52Kernel`](@ref)
"""
struct MaternKernel{Tν<:Real} <: SimpleKernel
    ν::Vector{Tν}
    function MaternKernel(; nu::Real=1.5, ν::Real=nu)
        @check_args(MaternKernel, ν, ν > zero(ν), "ν > 0")
        return new{typeof(ν)}([ν])
    end
end

@functor MaternKernel

@inline function kappa(κ::MaternKernel, d::Real)
    result = _matern(first(κ.ν), d)
    return ifelse(iszero(d), one(result), result)
end

function _matern(ν::Real, d::Real)
    y = sqrt(2ν) * d
    return exp((one(d) - ν) * logtwo - loggamma(ν) + ν * log(y) + log(besselk(ν, y)))
end

metric(::MaternKernel) = Euclidean()

Base.show(io::IO, κ::MaternKernel) = print(io, "Matern Kernel (ν = ", first(κ.ν), ")")

## Matern12Kernel = ExponentialKernel aliased in exponential.jl

"""
    Matern32Kernel()

Matérn kernel of order ``3/2``.

# Definition

For inputs ``x, x' \\in \\mathbb{R}^d``, the Matérn kernel of order ``3/2`` is given by
```math
k(x, x') = \\big(1 + \\sqrt{3} \\|x - x'\\|_2 \\big) \\exp\\big(- \\sqrt{3}\\|x - x'\\|_2\\big).
```

See also: [`MaternKernel`](@ref)
"""
struct Matern32Kernel <: SimpleKernel end

kappa(::Matern32Kernel, d::Real) = (1 + sqrt(3) * d) * exp(-sqrt(3) * d)

metric(::Matern32Kernel) = Euclidean()

Base.show(io::IO, ::Matern32Kernel) = print(io, "Matern 3/2 Kernel")

"""
    Matern52Kernel()

Matérn kernel of order ``5/2``.

# Definition

For inputs ``x, x' \\in \\mathbb{R}^d``, the Matérn kernel of order ``5/2`` is given by
```math
k(x, x') = \\bigg(1 + \\sqrt{5} \\|x - x'\\|_2 + \\frac{5}{3}\\|x - x'\\|_2^2\\bigg)
           \\exp\\big(- \\sqrt{5}\\|x - x'\\|_2\\big).
```

See also: [`MaternKernel`](@ref)
"""
struct Matern52Kernel <: SimpleKernel end

kappa(::Matern52Kernel, d::Real) = (1 + sqrt(5) * d + 5 * d^2 / 3) * exp(-sqrt(5) * d)

metric(::Matern52Kernel) = Euclidean()

Base.show(io::IO, ::Matern52Kernel) = print(io, "Matern 5/2 Kernel")


text:
## dcf.jl  Tests related to decision cost functions

using Test
import CSV
import GZip

x = GZip.open("ru.2009.table.gz") do fd
    CSV.read(fd, delim='\t', truestrings=["TRUE"], falsestrings=["FALSE"])
end
tnt = TNT(x)
tar, non = tnt.tar, tnt.non
ruc = roc(tnt, collapse=false)
r = roc(tnt)

## scalar
d = DCF(0.01, 1, 10)            # nist 2008
@test oeff(d) ≈ 1/9.9
@test peff(d) ≈ 1/10.9
@test plo(d) ≈ -log(9.9)
## print
show(stdout, MIME"text/plain"(), d)
println([d,d])

lo = plo(d)
b = 0.03291547867387523
@test ber(x, lo) ≈ b
@test ber(tnt, lo) ≈ b
@test ber(tar, non, lo) ≈ b
@test ber(ruc, lo) ≈ b

m1 = 0.03173675448189074
@test minber(x, lo) ≈ m1
@test minber(tnt, lo) ≈ m1
@test minber(tar, non, lo) ≈ m1
@test minber(r, lo) ≈ m1

## set a global DCF as well
setdcf(cmiss=10)                # evalita 2009

for norm in (false, true)
    c = 0.035877871754524004 * 10^norm
    @test dcf(x, d=d, norm=norm) ≈ c
    @test dcf(tnt, d=d, norm=norm) ≈ c
    @test dcf(tar, non, d=d, norm=norm) ≈ c
    @test dcf(ruc, d=d, norm=norm) ≈ c

    m2 = 0.034593062385260914 * 10^norm
    @test mindcf(x, d=d, norm=norm) ≈ m2
    @test mindcf(tnt, d=d, norm=norm) ≈ m2
    @test mindcf(tar, non, d=d, norm=norm) ≈ m2
    @test mindcf(r, d=d, norm=norm) ≈ m2

    ## global DCF
    c = 0.23094942466561763084 * 2^norm
    @test dcf(x, norm=norm) ≈ c
    @test dcf(tnt, norm=norm) ≈ c
    @test dcf(tar, non, norm=norm) ≈ c
    @test dcf(ruc, norm=norm) ≈ c

    m3 = 0.21843024685287174003 * 2^norm
    @test mindcf(x, norm=norm) ≈ m3
    @test mindcf(tnt, norm=norm) ≈ m3
    @test mindcf(tar, non,norm=norm) ≈ m3
    @test mindcf(ruc, norm=norm) ≈ m3
end

## array
lo = collect(-7:0.01:7)
dc = DCF(sigmoid.(lo), 1, 1)

@test ber(x, lo) ≈ ber(ruc, lo) atol=1e-5
@test ber(tnt, lo) == ber(tar, non, lo)

@test minber(x, lo) == minber(tnt, lo) == minber(tar, non, lo) == minber(r, lo)
@test all(minber(r, lo) .≤ ber(ruc, lo))

for norm in (false, true)
    @test dcf(x, d=dc, norm=norm) ≈ dcf(ruc, d=dc, norm=norm)
    @test dcf(x, d=dc, norm=norm) == dcf(tnt, d=dc, norm=norm) == dcf(tar, non, d=dc, norm=norm)
    @test mindcf(x, d=dc, norm=norm) == mindcf(tnt, d=dc, norm=norm) == mindcf(tar, non, d=dc, norm=norm) == mindcf(r, d=dc, norm=norm)
    @test all(mindcf(r, d=dc, norm=norm) .≤ dcf(ruc, d=dc, norm=norm))
end


text:
function multiplyexistingvar(v::Vector{PolyVar{C}}, x::PolyVar{C}, i::Int) where {C}
    updatez = z -> begin
        newz = copy(z)
        newz[i] += 1
        newz
    end
    # /!\ v not copied for efficiency, do not mess up with vars
    v, updatez
end
function insertvar(v::Vector{PolyVar{C}}, x::PolyVar{C}, i::Int) where {C}
    n = length(v)
    I = 1:i-1
    J = i:n
    K = J.+1
    w = Vector{PolyVar{C}}(undef, n+1)
    w[I] = v[I]
    w[i] = x
    w[K] = v[J]
    updatez = z -> begin
        newz = Vector{Int}(undef, n+1)
        newz[I] = z[I]
        newz[i] = 1
        newz[K] = z[J]
        newz
    end
    w, updatez
end

include("cmult.jl")
include("ncmult.jl")

MP.multconstant(α, x::Monomial)   = Term(α, x)
MP.mapcoefficientsnz(f::Function, p::Polynomial) = Polynomial(f.(p.a), p.x)

# I do not want to cast x to TermContainer because that would force the promotion of eltype(q) with Int
function Base.:(*)(x::DMonomialLike, p::Polynomial)
    # /!\ No copy of a is done
    Polynomial(p.a, x*p.x)
end
function Base.:(*)(x::DMonomialLike{false}, p::Polynomial)
    # /!\ No copy of a is done
    # Order may change, e.g. y * (x + y) = y^2 + yx
    Polynomial(monovec(p.a, [x*m for m in p.x])...)
end
function Base.:(*)(p::Polynomial, x::DMonomialLike)
    # /!\ No copy of a is done
    Polynomial(p.a, p.x*x)
end

function _term_poly_mult(t::Term{C, S}, p::Polynomial{C, T}, op::Function) where {C, S, T}
    U = Base.promote_op(op, S, T)
    if iszero(t)
        zero(Polynomial{C, U})
    else
        n = nterms(p)
        allvars, maps = mergevars([t.x.vars, p.x.vars])
        nv = length(allvars)
        # Necessary to annotate the type in case it is empty
        Z = Vector{Int}[zeros(Int, nv) for i in 1:n]
        for i in 1:n
            Z[i][maps[1]] = t.x.z
            Z[i][maps[2]] += p.x.Z[i]
        end
        Polynomial(op.(t.α, p.a), MonomialVector(allvars, Z))
    end
end
Base.:(*)(p::Polynomial, t::Term) = _term_poly_mult(t, p, (α, β) -> β * α)
Base.:(*)(t::Term, p::Polynomial) = _term_poly_mult(t, p, *)
_sumprod(a, b) = a * b + a * b
function Base.:(*)(p::Polynomial{C, S}, q::Polynomial{C, T}) where {C, S, T}
    U = Base.promote_op(_sumprod, S, T)
    if iszero(p) || iszero(q)
        zero(Polynomial{C, U})
    else
        samevars = _vars(p) == _vars(q)
        if samevars
            allvars = _vars(p)
        else
            allvars, maps = mergevars([_vars(p), _vars(q)])
        end
        N = length(p)*length(q)
        Z = Vector{Vector{Int}}(undef, N)
        a = Vector{U}(undef, N)
        i = 0
        for u in p
            for v in q
                if samevars
                    z = u.x.z + v.x.z
                else
                    z = zeros(Int, length(allvars))
                    z[maps[1]] += u.x.z
                    z[maps[2]] += v.x.z
                end
                i += 1
                Z[i] = z
                a[i] = u.α * v.α
            end
        end
        polynomialclean(allvars, a, Z)
    end
end


text:
using FearlessPoker
using FactCheck
using StatsBase

facts("Hand Parsing") do
  @fact Hand(13757003685786313) --> Hand("8d,Tc,Td,Th,Ts,Qh,Qs")
  @fact Hand(405331374) --> Hand("2c,2d,2s,3h,4s,6h,Td")
  @fact Hand(8039256302) --> Hand("2c,2s,3c,5c,6s,7d,Tc")
  @fact Hand(690690) --> Hand("2c,2d,2h,2s,3c,3d,4c")
end

# write your own tests here
facts("Hand evaluation") do
  @fact evaluate(Hand("Ah,Ad,Ac,As,Ts,8d,2c")) --> 814.1
  @fact evaluate(Hand("Ah,Ad,Ac,As,7s,8d,2c")) --> 814.08
  @fact evaluate(Hand("Ah,Kh,Th,Jh,Qh,9d,Td")) --> 914
  @fact evaluate(Hand("Ah,Ad,Ac,Td,Tc,2c,8s")) --> 714.10
  @fact evaluate(Hand("Ah,Kh,8h,Jh,Qh,9d,Td")) --> 614
  @fact evaluate(Hand("Ah,Jd,Qc,Ks,9d,Tc,2h")) --> 514
  @fact evaluate(Hand("Ah,Ad,2s,2c,9d,Td,Tc")) --> 314.1002
  @fact evaluate(Hand("Kc,Kd,9h,7c,5c,Tc,3h")) --> 213.100907
  @fact evaluate(Hand("Ah,Jc,9s,7s,7d,6h,4d")) --> 207.141109
  @fact evaluate(Hand("7s,Kc,4d,5c,8s,Jd,9s")) --> 113.11090807
  @fact evaluate(Hand("Qh,9h,4d,9c,Td,Jc,3d")) --> 209.121110
  @fact evaluate(Hand("2c,2d,5s,6h,Qc,Qs,Ad")) --> 312.0214
  @fact evaluate(Hand("2c,2d,2s,3h,4s,6h,Td")) --> 402.1006
  @fact evaluate(Hand("2c,2s,3c,5c,6s,7d,Tc")) --> 202.100706
end

facts("Evaluate Straights") do
  @fact FearlessPoker.findstraight(Hand("Ah,2d,3c,4d,5s,9d,4h")) --> Hand("Ah,2d,3c,4d,4h,5s")
  @fact FearlessPoker.findstraight(Hand("Th,Jh,Qc,Ks,As,2d,3d")) --> Hand("Th,Jh,Qc,Ks,As")
  @fact FearlessPoker.findstraight(Hand("Ah,2s,3s,6d,7d,9s")) --> nothing
  @fact FearlessPoker.findstraight(Hand("2s,3s,4d,5d,6c,7c,8h")) --> Hand("4d,5d,6c,7c,8h")
  @fact FearlessPoker.findstraight(Hand("Ts,Jc,Qd,Ks,Ah")) --> Hand("Ts,Jc,Qd,Ks,Ah")
  @fact_throws AssertionError FearlessPoker.findstraight(Hand("2s,3s,4d"))
end

facts("Evaluate Flush") do
  @fact FearlessPoker.findflush(Hand("2d,3d,4d,7d,9d,Ah,Ks")) --> Hand("2d,3d,4d,7d,9d")
  @fact FearlessPoker.findflush(Hand("2d,3s,4d,5h,As,Kc")) --> nothing
end

facts("Mapping Hands to an Integer Id") do
  @fact convert(Integer, Hand("2d,3d,4d,7d,9d,Ah,Ks")) --> convert(Integer, Hand("2d,3d,9d,Ah,Ks,4d,7d"))

  @fact length(unique(map(h -> convert(Integer, FearlessPoker.Hand(h)), combinations(DECK, 4)))) --> length(combinations(1:52, 4))
end


text:
module MPIExecutor

using Serialization
using MPI

import Base: size, get!

include("RemoteFunction.jl")
include("worker.jl")

macro prof(name::Symbol, ex)
	if true
			quote
					local elapsedtime = time_ns()
					local val = $(esc(ex))
					elapsedtime = time_ns() - elapsedtime
					println($(String(name))," ", elapsedtime/1e9)
					val
			end
	else
			quote
				$(esc(ex))
			end
	end
end

export MPIPoolExecutor, shutdown!, @remote,
    submit!, run!, run_until!, then!, finally!,
    fulfill!, whenall!, get!,
    run_broadcast!, main_worker, @everywhere

abstract type BaseExecutor end

include("Future.jl")

struct WorkUnit
    f::Function
    args::Tuple
    fut::Future

    function WorkUnit(f, args, fut)
        new(f, args, fut)
    end
end

mutable struct MPIPoolExecutor <: BaseExecutor
    slaves::Array{Int64,1}
    idle::Array{Int64,1}
    comm::MPI.Comm
    identifier::Int64
    tracker::Int64
    runnable::Array{WorkUnit,1}
    running::Dict{Int64, WorkUnit}
    io::IOBuffer

    function MPIPoolExecutor(comm::MPI.Comm=MPI.COMM_WORLD)
      worker_count = MPI.Comm_size(comm) - 1
      slaves = Int64[1:worker_count;]
      new(slaves, copy(slaves), comm, 0, 0, WorkUnit[], Dict{Int64, WorkUnit}(), IOBuffer())
    end
end

Base.size(pool::MPIPoolExecutor) = isempty(pool.slaves) ? 1 : length(pool.slaves)

function shutdown!(pool::MPIPoolExecutor)
    while !all_idle(pool)
        wait_any!(pool)
    end

    io = pool.io
    seek(io, 0)

    buf = MPI.Buffer(io.data, io.size, MPI.Datatype(UInt8))
    for worker in pool.slaves
        MPI.Send(buf, worker, 2, pool.comm)
    end
end

function MPIPoolExecutor(f::Function, ::Nothing=nothing)
    if !MPI.Initialized()
        MPI.Init()
    end

    @assert !MPI.Finalized()
    MPIPoolExecutor(f, MPI.COMM_WORLD)
end

function MPIPoolExecutor(f::Function, worker_count::Int64, comm=MPI.COMM_WORLD)
    if !MPI.Initialized()
        MPI.Init()
    end

    @assert !MPI.Finalized()

    if MPI.Comm_size(comm) <= worker_count
        additional_workers = worker_count - MPI.Comm_size(comm) + 1
        intercomm = MPI.Comm_spawn("julia", ["-e", "import MPIExecutor; MPIExecutor.main_worker()"], additional_workers, comm)
        comm = MPI.Intercomm_merge(intercomm, false)
        try
            @assert MPI.Comm_size(comm) == worker_count + 1
            MPIPoolExecutor(f, comm)
        finally
            # free the merged intercomm https://github.com/open-mpi/ompi/issues/8426
            MPI.free(comm)
        end
    elseif MPI.Comm_size(comm) == worker_count + 1
        MPIPoolExecutor(f, comm)
    else
        rank = MPI.Comm_rank(comm)
        color = (rank <= worker_count) ? 0 : 1
        subset_comm = MPI.Comm_split(comm, color, rank)

        ret = try
            if color == 0
                MPIPoolExecutor(f, subset_comm)
            else
                nothing
            end
        finally
            MPI.free(subset_comm)
        end

        MPI.Barrier(comm)
        ret
    end
end

function MPIPoolExecutor(f::Function, comm::MPI.Comm)
    if MPI.Comm_rank(comm) == 0
        pool = MPIPoolExecutor(comm)

        try
            f(pool)
        finally
            shutdown!(pool)
        end
    else
        main_worker(comm)
    end
end

function is_anon_function(f::Function)
    t = typeof(f)
    tn = t.name
    if isdefined(tn, :mt)
        name = tn.mt.name
        mod = tn.module
        return mod === Main && # only Main
            t.super === Function && # only Functions
            unsafe_load(Base.unsafe_convert(Ptr{UInt8}, tn.name)) == UInt8('#') && # hidden type
            (!isdefined(mod, name) || t != typeof(getfield(mod, name))) # XXX: 95% accurate test for this being an inner function
            # TODO: more accurate test? (tn.name !== "#" name)
    end
    return false
end

register!(pool::MPIPoolExecutor, x::Function...) = map(f -> register!(pool, f), x)

function send_to_workers(pool::MPIPoolExecutor, tag, args...)
    io = pool.io
    seek(io, 0)

    for x in args
      serialize(io, x)
    end

    buf = MPI.Buffer(io.data, io.size, MPI.Datatype(UInt8))
    for worker in pool.slaves
        MPI.Send(buf, worker, tag, pool.comm)
    end
end

function register!(pool::MPIPoolExecutor, f::Function)
    @assert is_anon_function(f)
    rid = (pool.identifier += 1)

    send_to_workers(pool, 1, rid, f)
    RemoteFunction(rid, f)
end

macro everywhere(pool, ex::Expr)
  quote
    eval!($(esc(pool)), $(QuoteNode(ex)))
  end
end

function eval!(pool::MPIPoolExecutor, ex::Expr)
    if isempty(pool.slaves)
        Main.eval(ex)
    else
        send_to_workers(pool, 0, ex)
    end
end

function submit!(pool::MPIPoolExecutor, f::Function, args...)
    t = WorkUnit(f, args, Future(pool))
    push!(pool.runnable, t)
    t.fut
end

function all_idle(pool::MPIPoolExecutor)
    length(pool.idle) == length(pool.slaves)
end

function run!(pool::MPIPoolExecutor)
    if isempty(pool.slaves) && !isempty(pool.runnable)
        todo = pop!(pool.runnable)
        fulfill!(todo.fut, todo.f(todo.args...))
    else
        while !isempty(pool.runnable) && !isempty(pool.idle)
            todo = pop!(pool.runnable)
            dispatch!(pool, todo, pop!(pool.idle))
        end

        receive_any!(pool)
    end

    isempty(pool.runnable) && all_idle(pool)
end

function run_until!(pool::MPIPoolExecutor)
    run_until!(pool, () -> false)
end

function run_until!(pool::MPIPoolExecutor, pull::Function)
    if isempty(pool.slaves)
        # master-only mode
        if isempty(pool.runnable) && ! pull()
            return nothing
        end

        @assert !isempty(pool.runnable)
        todo = pop!(pool.runnable)
        fulfill!(todo.fut, todo.f(todo.args...))
    else
        while !isempty(pool.runnable) && !isempty(pool.idle)
            todo = pop!(pool.runnable)
            dispatch!(pool, todo, pop!(pool.idle))
        end

        while !isempty(pool.idle) && pull()
            @assert !isempty(pool.runnable)
            todo = pop!(pool.runnable)
            dispatch!(pool, todo, pop!(pool.idle))
        end

        if !all_idle(pool)
            wait_any!(pool)
        end
    end
end

function handle_recv!(pool::MPIPoolExecutor, s::MPI.Status)
    received_from = MPI.Get_source(s)
    count = MPI.Get_count(s, UInt8)

    io = pool.io
    Base.ensureroom(io, count)

    buf = MPI.Buffer(io.data, count, MPI.Datatype(UInt8))
    MPI.Recv!(buf, received_from, 0, pool.comm)
    io.size = count
    seek(io, 0)

    tracker_id = deserialize(io)::Int64
    push!(pool.idle, received_from)
    fulfill!(pool.running[tracker_id].fut, deserialize(io)::Any)
end

function receive_any!(pool::MPIPoolExecutor)
    result, s = MPI.Iprobe(MPI.MPI_ANY_SOURCE, 0, pool.comm)
    if result
      handle_recv!(pool, s)
    else
      nothing
    end
end

function wait_any!(pool::MPIPoolExecutor)
    s = MPI.Probe(MPI.MPI_ANY_SOURCE, 0, pool.comm)
    handle_recv!(pool, s)
end

function dispatch!(pool::MPIPoolExecutor, work::WorkUnit, worker)
    io = pool.io
    seek(io, 0)

    tracker_id = (pool.tracker += 1)
    pool.running[tracker_id] = work
    serialize(io, work.f)
    serialize(io, tracker_id)
    serialize(io, work.args)

    buf = MPI.Buffer(io.data, io.size, MPI.Datatype(UInt8))
    MPI.Send(buf, worker, 3, pool.comm)
end

function run_broadcast!(pool::MPIPoolExecutor, f::Function, args...)
    @assert all_idle(pool)

    if !isempty(pool.slaves)
      all_futs = Future[]

      for worker in pool.idle
          t = WorkUnit(f, args, Future(pool))
          dispatch!(pool, t, worker)
          push!(all_futs, t.fut)
      end

      pool.idle = []

      whenall!(all_futs)
    else
      Future(pool, Some(f(args...)))
    end
end

end # module


text:
# Inefficient: creates copies, no in place
function mergesorted(sortedarray1::Vector{T}, sortedarray2::Vector{T})::Vector{T} where T <: Real
	n = length(sortedarray1)
	m = length(sortedarray2)
	result = Vector{Int}(undef, n + m)
	i = 1
	j = 1

	while i ≤ n && j ≤ m
		if sortedarray1[i] ≤ sortedarray2[j]
			result[i + j - 1] = sortedarray1[i]
			i += 1
		else
			result[i + j - 1] = sortedarray2[j]
			j += 1
		end
	end

	if i > n
		result[i + j - 1 : end] = sortedarray2[j:end]
	else
		result[i + j - 1 : end] = sortedarray1[i:end]
	end

	return result
end

function mergesort(array::Vector{T})::Vector{T} where T <: Real
	# Split
	len = length(array)
	len == 1 && (return array)

	# Sort
	sortedleft  = mergesort(array[1:(len ÷ 2)])
	sortedright = mergesort(array[(len ÷ 2 + 1):end])

	# Combine
	return mergesorted(sortedleft, sortedright)
end

println(mergesort([6, 4, 1, 9, 8, 2, 6, 4]))

max_stars_repo_name: SudipSinha/edu

max_issues_repo_name: SudipSinha/edu

max_forks_repo_name: SudipSinha/edu


text:
module PCA

    include("LinearRegression.jl")
    using LinearAlgebra

    function helper_findPCAs(
        X    :: AbstractArray{T, 2};
        num  :: Integer = 1,  # number of eignevectors wish to return
    ) where T <: AbstractFloat

        # dim means the spatial dimension of X
        # N means the number of sampling (time, for example)
        dim, N = size(X)

        Σ = Symmetric(X * X')

        println("Dimension of Σ : ", size(Σ))
        println("Any missing data? ", any(isnan.(Σ)))

        if any(isnan.(Σ))
            throw(ErrorException("Data contains NaN!"))
        end

        @time F = eigen(Σ)

        order = collect(1:dim)
        sort!(order; by = i -> F.values[i], rev=true)

        evs = zeros(T, dim, num)
        for i = 1:num
            evs[:, i] = normalize(F.vectors[:, order[i]])
        end

        return evs

    end

    # This function removes seasonality and detrend the input data
    function findPCAs(
        raw_data,
        EOF_idx;
        modes::Integer
    )

        data = copy(raw_data)

        data_size = size(data)
        field_size = data_size[1:end-1]
        N  = reduce(*, field_size)
        Nt = data_size[end]
        flat_data = reshape(data, N, Nt)
        flat_EOF_idx = reshape(EOF_idx, :)

        # construct mapping of PCA
        grid_numbering = reshape( collect(1:N), field_size... )
        EOF_mapping = grid_numbering[EOF_idx]  # The i-th reduced grid idx onto the original ?-th grid

        EOF_pts = sum(EOF_idx)
        EOF_input = zeros(Float64, EOF_pts, Nt)
        t_vec = collect(eltype(EOF_input), 1:Nt)

        for t = 1:Nt
            _data = view(flat_data, :, t)
            EOF_input[:, t] = _data[flat_EOF_idx]
        end

        for i=1:EOF_pts
            EOF_input[i, :] = detrend(t_vec, view(EOF_input, i, :), order=2)
        end

        # remove seasonality
        for i=1:EOF_pts
            vw = view(EOF_input, i, :)
            rmSeasonality!(vw; period=12)
        end

        # Doing PCAs
        println("Solving for PCA...")
        eigen_vectors = helper_findPCAs(EOF_input, num=modes)

        PCAs = zeros(Float64, N, modes)
        PCAs .= NaN

        for i = 1:EOF_pts
            PCAs[EOF_mapping[i], :] = eigen_vectors[i, :]
        end

        # Project anomalies onto PCAs
        PCAs_ts = (transpose(eigen_vectors) * EOF_input) # dimension = (modes, Nt)
#        for m = 1:modes
#            PCAs_ts[m, :] ./= std(view(PCAs_ts, m, :))
#        end

        return reshape(PCAs, field_size..., modes), PCAs_ts
    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
struct DisplaySimulator
    display::Union{AbstractDisplay, Nothing}
    render_kwargs
    max_fps::Float64
    predisplay::Function
    extra_initial::Bool
    extra_final::Bool
    stepsim::StepSimulator
end

"""
    DisplaySimulator(;kwargs...)

Create a simulator that displays each step of a simulation.

Given a POMDP or MDP model `m`, this simulator roughly works like

    for step in stepthrough(m, ...)
        display(render(m, step))
    end

# Keyword Arguments
- `display::AbstractDisplay`: the display to use for the first argument to the `display` function. If this is `nothing`, `display(...)` will be called without an `AbstractDisplay` argument.
- `render_kwargs::NamedTuple`: keyword arguments for `POMDPModelTools.render(...)`
- `max_fps::Number=10`: maximum number of frames to be displayed per second - `sleep` will be used to skip extra time, so this is not designed for high precision
- `predisplay::Function`: function to call before every call to `display(...)`. The only argument to this function will be the display (if it is specified) or `nothing`
- `extra_initial::Bool=false`: if `true`, display an extra step at the beginning with only elements `t`, `sp`, and `bp` for POMDPs (this can be useful to see the initial state if `render` displays only `sp` and not `s`).
- `extra_final`::Bool=true`: if `true`, display an extra step at the end with only elements `t`, `done`, `s`, and `b` for POMDPs (this can be useful to see the final state if `render` displays only `s` and not `sp`).
- `max_steps::Integer`: maximum number of steps to run for
- `spec::NTuple{Symbol}`: specification of what step elements to display (see `eachstep`)
- `rng::AbstractRNG`: random number generator

See the POMDPSimulators documentation for more tips about using specific displays.
"""
function DisplaySimulator(;display=nothing,
                           render_kwargs=NamedTuple(),
                           max_fps=10,
                           predisplay=(d)->nothing,
                           extra_initial=false,
                           extra_final=true,
                           max_steps=nothing,
                           spec=CompleteSpec(),
                           rng=Random.GLOBAL_RNG
                         )
    stepsim = StepSimulator(rng, max_steps, spec)
    return DisplaySimulator(display,
                            render_kwargs,
                            max_fps,
                            predisplay,
                            extra_initial,
                            extra_final,
                            stepsim)
end

function simulate(sim::DisplaySimulator, m, args...)
    rsum = 0.0
    disc = 1.0
    dt = 1/sim.max_fps
    tm = time()
    isinitial = true
    last = NamedTuple() # for extra_final

    for step in simulate(sim.stepsim, m, args...)
        if isinitial && sim.extra_initial
            isinitial = false
            istep = initialstep(m, step)
            vis = render(m, istep; sim.render_kwargs...)
            perform_display(sim, vis)
            sleep_until(tm += dt)
        end

        vis = render(m, step; sim.render_kwargs...)
        perform_display(sim, vis)
        rsum += disc*get(step, :r, missing)
        disc *= discount(m)
        sleep_until(tm += dt)

        last = step # save for extra final
    end

    if sim.extra_final
        fstep = finalstep(m, last)
        vis = render(m, fstep; sim.render_kwargs...)
        perform_display(sim, vis)
    end

    if ismissing(rsum)
        return nothing
    else
        return rsum
    end
end

sleep_until(t) = sleep(max(t-time(), 0.0))

initialstep(m::MDP, step) = (t=0, sp=get(step, :s, missing))
initialstep(m::POMDP, step) = (t=0,
                               sp=get(step, :s, missing),
                               bp=get(step, :b, missing))
finalstep(m::MDP, last) = (done=true,
                           t=get(last, :t, missing) + 1,
                           s=get(last, :sp, missing))
finalstep(m::POMDP, last) = (done=true,
                             t=get(last, :t, missing) + 1,
                             s=get(last, :sp, missing),
                             b=get(last, :bp, missing))

function perform_display(sim::DisplaySimulator, vis)
    sim.predisplay(sim.display)
    if sim.display===nothing
        display(vis)
    else
        display(sim.display, vis)
    end
end


text:
function pulled_back_integrand(op::HH3DSingleLayerFDBIO,
    test_local_space::LagrangeRefSpace{<:Any,0},
    trial_local_space::LagrangeRefSpace{<:Any,0},
    test_chart, trial_chart)

    (u,v) -> begin

        x = neighborhood(test_chart,u)
        y = neighborhood(trial_chart,v)

        f = test_local_space(x)
        g = trial_local_space(y)

        j = jacobian(x) * jacobian(y)

        α = op.alpha
        γ = op.gamma
        R = norm(cartesian(x)-cartesian(y))
        G = exp(-γ*R)/(4*π*R)

        αjG = α*G*j

        SA[f[1].value * αjG * g[1].value]
    end
end

function pulled_back_integrand(op::HH3DHyperSingularFDBIO,
    test_local_space::LagrangeRefSpace{<:Any,1},
    trial_local_space::LagrangeRefSpace{<:Any,1},
    test_chart, trial_chart)

    (u,v) -> begin

        x = neighborhood(test_chart,u)
        y = neighborhood(trial_chart,v)

        nx = normal(x)
        ny = normal(y)

        f = test_local_space(x)
        g = trial_local_space(y)

        j = jacobian(x) * jacobian(y)

        α = op.alpha
        β = op.beta
        γ = op.gamma
        R = norm(cartesian(x)-cartesian(y))
        G = exp(-γ*R)/(4*π*R)

        αjG = ny*α*G*j
        βjG = β*G*j

        A = SA[ αjG*g[1].value, αjG*g[2].value, αjG*g[3].value]
        B = SA[ βjG*g[1].curl, βjG*g[2].curl, βjG*g[3].curl ]

        SMatrix{3,3}((
            dot(nx*f[1].value, A[1]) + dot(f[1].curl, B[1]),
            dot(nx*f[2].value, A[1]) + dot(f[2].curl, B[1]),
            dot(nx*f[3].value, A[1]) + dot(f[3].curl, B[1]),
            dot(nx*f[1].value, A[2]) + dot(f[1].curl, B[2]),
            dot(nx*f[2].value, A[2]) + dot(f[2].curl, B[2]),
            dot(nx*f[3].value, A[2]) + dot(f[3].curl, B[2]),
            dot(nx*f[1].value, A[3]) + dot(f[1].curl, B[3]),
            dot(nx*f[2].value, A[3]) + dot(f[2].curl, B[3]),
            dot(nx*f[3].value, A[3]) + dot(f[3].curl, B[3])))
    end
end

function momintegrals!(op::HH3DSingleLayerFDBIO,
    test_local_space::LagrangeRefSpace{<:Any,0}, trial_local_space::LagrangeRefSpace{<:Any,0},
    test_triangular_element, trial_triangular_element, out, strat::SauterSchwabStrategy)

    I, J, K, L = SauterSchwabQuadrature.reorder(
        test_triangular_element.vertices,
        trial_triangular_element.vertices, strat)

    test_triangular_element  = simplex(
        test_triangular_element.vertices[I[1]],
        test_triangular_element.vertices[I[2]],
        test_triangular_element.vertices[I[3]])

    trial_triangular_element = simplex(
        trial_triangular_element.vertices[J[1]],
        trial_triangular_element.vertices[J[2]],
        trial_triangular_element.vertices[J[3]])

    # igd = MWSL3DIntegrand(test_triangular_element, trial_triangular_element,
    #     op, test_local_space, trial_local_space)
    igd = pulled_back_integrand(op, test_local_space, trial_local_space,
        test_triangular_element, trial_triangular_element)
    G = SauterSchwabQuadrature.sauterschwab_parameterized(igd, strat)
    out[1,1] += G[1,1]

    nothing
end

function momintegrals!(op::HH3DHyperSingularFDBIO,
    test_local_space::LagrangeRefSpace{<:Any,1}, trial_local_space::LagrangeRefSpace{<:Any,1},
    test_triangular_element, trial_triangular_element, out, strat::SauterSchwabStrategy)

    I, J, K, L = SauterSchwabQuadrature.reorder(
        test_triangular_element.vertices,
        trial_triangular_element.vertices, strat)

    test_triangular_element  = simplex(
        test_triangular_element.vertices[I[1]],
        test_triangular_element.vertices[I[2]],
        test_triangular_element.vertices[I[3]])

    trial_triangular_element = simplex(
        trial_triangular_element.vertices[J[1]],
        trial_triangular_element.vertices[J[2]],
        trial_triangular_element.vertices[J[3]])

    test_sign = Combinatorics.levicivita(I)
    trial_sign = Combinatorics.levicivita(J)
    σ = test_sign * trial_sign

    igd = pulled_back_integrand(op, test_local_space, trial_local_space,
        test_triangular_element, trial_triangular_element)
    G = SauterSchwabQuadrature.sauterschwab_parameterized(igd, strat)
    for j ∈ 1:3, i ∈ 1:3
        out[i,j] += G[K[i],L[j]] * σ
    end

    nothing
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule CImGui_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("CImGui")
JLLWrappers.@generate_main_file("CImGui", UUID("7dd61d3b-0da5-5c94-bbf9-a0296c6e3925"))
end  # module CImGui_jll


text:
using Test
import MathOptInterface
const MOI = MathOptInterface

mul2(x) = x * 2
div2(x) = div(x, 2)
d = MOI.Utilities.DenseDict{Int, Float64}(div2, mul2, 3)

d[4] = 0.25
@test !haskey(d, 2)
@test haskey(d, 4)
@test !haskey(d, 6)
@test length(d) == 1
@test collect(d) == [4 => 0.25]
@test d[4] == 0.25

d[2] = 1.5
@test haskey(d, 2)
@test haskey(d, 4)
@test !haskey(d, 6)
@test length(d) == 2
@test collect(d) == [2 => 1.5, 4 => 0.25]
@test d[2] == 1.5
@test d[4] == 0.25

d[6] = 0.75
@test haskey(d, 2)
@test haskey(d, 4)
@test haskey(d, 6)
@test length(d) == 3
@test collect(d) == [2 => 1.5, 4 => 0.25, 6 => 0.75]
@test d[2] == 1.5
@test d[4] == 0.25
@test d[6] == 0.75


text:
# Triangle, square, pentagonal, hexagonal, heptagonal, and octagonal numbers are
# all figurate (polygonal) numbers and are generated by the following formulae:
#
# Triangle    P3,n=n(n+1)/2    1, 3, 6, 10, 15, ...
# Square      P4,n=n2          1, 4, 9, 16, 25, ...
# Pentagonal  P5,n=n(3n−1)/2   1, 5, 12, 22, 35, ...
# Hexagonal   P6,n=n(2n−1)     1, 6, 15, 28, 45, ...
# Heptagonal  P7,n=n(5n−3)/2   1, 7, 18, 34, 55, ...
# Octagonal   P8,n=n(3n−2)     1, 8, 21, 40, 65, ...
#
# The ordered set of three 4-digit numbers: 8128, 2882, 8281, has three
# interesting properties.
# 1. The set is cyclic, in that the last two digits of each number is the first
# two digits of the next number (including the last number with the first).
# 2. Each polygonal type: triangle (P3,127=8128), square (P4,91=8281), and
# pentagonal (P5,44=2882), is represented by a different number in the set.
# 4. This is the only set of 4-digit numbers with this property.
#
# Find the sum of the only ordered set of six cyclic 4-digit numbers for which
# each polygonal type: triangle, square, pentagonal, hexagonal, heptagonal, and
# octagonal, is represented by a different number in the set.

using ProjectEulerSolutions

# Solve for n given x and s: P(s,n) = x.
function polygonal_inv(s::Integer, x::Integer)::Float32
    return (sqrt(8*(s-2)*x + (s-4)^2) + (s-4))/(2*(s-2))
end

# Implements polygonal number P(s,n) = x, where s is the number of sides.
function polygonal_number(s::Integer, n::Integer)::Integer
   return div((s - 2) * n * n - n * (s - 4), 2)
end

# Helper recursive function to find the numbers that match the beginning and
# ending pairs.
function find_loop!(numbers::Array{Integer,1}, used_polys::Array{Integer,1}, end_pair::Integer, polys::Array{Array{Integer,1},1})::Bool
    d = 10^(ndigits(numbers[1])-2)
    if length(numbers) == 6 && end_pair == div(numbers[1], d)
        return true
    end
    ks = filter(x -> !(x in used_polys), 2:length(polys))
    for k in ks
        push!(used_polys, k)
        for n in polys[k]
            if div(n, d) == end_pair
                push!(numbers, n)
                found = find_loop!(numbers, used_polys, n % 100, polys)
                if found
                    return true
                end
                pop!(numbers)
            end
        end
        pop!(used_polys)
    end
    return false
end

# This is a verbose solution primarily due to manipulating the data into the
# correct data structures.  Only cycle through 4 digit polygonal numbers, then
# recursively loop through collections, remembering which polygonal numbers
# have already been used.  There is probably a graph-based solution which find
# paths through a network (6-partite graph?).
function p061solution(ndig_start::Integer=4, ndig_end::Integer=4)::Integer

    min_n = 10^(ndig_start-1)
    max_n = 10^ndig_end-1

    polys = Array{Array{Integer,1},1}()
    # Initialize data structure, making use of vectorized instructions
    for s = 3:8
        min_num = Integer(ceil(polygonal_inv(s, min_n)))
        max_num = Integer(floor(polygonal_inv(s, max_n)))

        # Calculate all 4 digit polygonal numbers for s
        poly_s = polygonal_number.(s, min_num:max_num)

        # Only keep numbers whose lower halves have two digits
        lower_s = poly_s .% 100
        inds = lower_s .> 9
        push!(polys, poly_s[inds])
    end

    # Loop through potential starting points, starting with s = 3.  The order
    # doesn't matter since all polygonal values are eventually used.
    used_polys = Array{Integer,1}([1])
    numbers = Array{Integer,1}()
    for n in polys[1]
        push!(numbers, n)
        if find_loop!(numbers, used_polys, n % 100, polys)
             return sum(numbers)
        end
        pop!(numbers)
    end

    return -1
end

p061 = Problems.Problem(p061solution)

Problems.benchmark(p061, 4, 4)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project