text:
using ColorLab,LinearAlgebra,YAML,Plots,NeuroAnalysis

"Get digital RGB color spectral measured from a specific display"
function RGBSpectral(measurement)
    C = map(i->parse.(Float64,split(i)),measurement["Color"])
    λ = measurement["WL"]
    I = measurement["Spectral"]
    return C,λ,I
end
colorstring(c::Vector) = join(string.(c)," ")

## Generate all color related data for a display
cd(@__DIR__)
displayname = "ROGPG279Q"
resultdir = "./$displayname"
isdir(resultdir) || mkpath(resultdir)
colordatapath = joinpath(resultdir,"colordata.yaml")
figfmt = [".png",".svg"]

config = YAML.load_file("../Configuration/CommandConfig_SNL-C.yaml")
C,λ,I = RGBSpectral(config["Display"][displayname]["SpectralMeasurement"])
RGBToLMS,LMSToRGB = RGBLMSMatrix(C,λ,I,observer=2)
RGBToXYZ,XYZToRGB = RGBXYZMatrix(C,λ,I,observer=2)

## Maximum Cone Isolating RGBs through a RGB color
th = [0.5,0.5,0.5]
# each column of LMSToRGB is the cone isolating RGB direction
ConeIsoRGBVec = trimatrix(LMSToRGB)
# Scale RGB direction into Unit Cube
ConeIsoRGBVec./=maximum(abs.(ConeIsoRGBVec),dims=1)
# since through color is the center of RGB cube, the line intersects at two symmatric points on the faces of unit cube
minc_lms = clamp.(quavectors(th .- 0.5ConeIsoRGBVec),0,1);maxc_lms = clamp.(quavectors(th .+ 0.5ConeIsoRGBVec),0,1)
mcc = contrast_michelson.(RGBToLMS*maxc_lms,RGBToLMS*minc_lms)

## Cone Isolating RGBs with same pooled michelson cone contrast
# since all the maximum Cone Isolating color pairs are symmatric around `th` color, i.e. `th = (maxc+minc)/2`,
# then the weber contrast using `maxc and th` is equivalent to michelson contrast using `maxc and minc`.
# Here we use Weber Cone Contrast to scale to the minimum michelson contrast.
bg = RGBToLMS*[th;1]
LMSToContrast,ContrastToLMS = LMSContrastMatrix(bg)
maxcc = LMSToContrast*RGBToLMS*maxc_lms
mincc = LMSToContrast*RGBToLMS*minc_lms
pcc = [norm(maxcc[1:3,i]) for i in 1:3]
ccf = minimum(pcc)./pcc
# since S cone contrast is much larger than L and M, scale S to L and M will significently decrease effectiveness
# of the S Cone Isolating color, so here we don't scale S Cone
ccf[3]=1
mmaxc_lms = clamp.(LMSToRGB*ContrastToLMS*quavectors(ccf'.*trivectors(maxcc)),0,1)
mminc_lms = clamp.(LMSToRGB*ContrastToLMS*quavectors(ccf'.*trivectors(mincc)),0,1)
mmcc = contrast_michelson.(RGBToLMS*mmaxc_lms,RGBToLMS*mminc_lms)

## DKL Isolating RGBs through a background RGB color
LMSToDKL,DKLToLMS = LMSDKLMatrix(bg,isnorm=true)
# each column of DKLToLMS is the DKL Isolating LMS direction, then it's converted to RGB direction
DKLIsoRGBVec = trimatrix(LMSToRGB*DKLToLMS)
# Scale RGB direction into Unit Cube
DKLIsoRGBVec./=maximum(abs.(DKLIsoRGBVec),dims=1)
# since through color is the center of RGB cube, the line intersects at two symmatric points on the faces of unit cube
minc_dkl = clamp.(quavectors(th .- 0.5DKLIsoRGBVec),0,1);maxc_dkl = clamp.(quavectors(th .+ 0.5DKLIsoRGBVec),0,1)

## DKL Isoluminance plane
DKLToRGB = LMSToRGB*DKLToLMS
RGBToDKL = LMSToDKL*RGBToLMS
lum = 0
hueangle_dkl_ilp = [0,30,60,75,80,83,90,95,101,105,110,120,150]
append!(hueangle_dkl_ilp,hueangle_dkl_ilp .+ 180)
# rotate l+ direction around l+m axis within the `lum` Isoluminance plane
DKLIsoLumRGBVec = trivectors(hcat(map(i->DKLToRGB*RotateXYZMatrix(deg2rad(i),dims=1)*[0,1,0,0], hueangle_dkl_ilp)...))
# find Intersections of Isoluminance directions with faces of unit RGB cube
anglec = hcat(map(i->intersectlineunitorigincube((DKLToRGB*[lum,0,0,1])[1:3],DKLIsoLumRGBVec[:,i]),1:size(DKLIsoLumRGBVec,2))...)
hue_dkl_ilp = clamp.(quavectors(anglec),0,1)
wp_dkl_ilp = repeat([th;1],inner=(1,size(hue_dkl_ilp,2)))

# Plot DKL Hues
IsoLumc = [RGB(hue_dkl_ilp[1:3,i]...) for i in 1:size(hue_dkl_ilp,2)]
IsoLumdkl = RGBToDKL*hue_dkl_ilp

title="DKL_$(length(hueangle_dkl_ilp))Hue_L$lum"
p=plot()
foreach(i->plot!(p,[0,IsoLumdkl[2,i]],[0,IsoLumdkl[3,i]],color=RGBA(0.5,0.5,0.5,0.5)),1:size(hue_dkl_ilp,2))
plot!(p,IsoLumdkl[2,:],IsoLumdkl[3,:],aspectratio=:equal,color=IsoLumc,lw=1.5,markersize=6,marker=:circle,
markerstrokewidth=0,legend=false,xlabel="L-M",ylabel="S-(L+M)",title=title)
p
foreach(ext->savefig(joinpath(resultdir,"$title$ext")),figfmt)

## HSL equal angular distance hue[0:30:330] and equal energy white with matched luminance in CIE [x,y,Y] coordinates
# Name:     R:1             Y:3             G:5                             B:9                             WP
hslhues =  [0.63    0.54    0.42    0.34    0.3     0.27    0.22    0.17    0.15    0.2     0.32    0.5     0.33;
            0.34    0.41    0.5     0.57    0.6     0.5     0.33    0.15    0.07    0.1     0.16    0.27    0.33;
            1.0     1.0     1.0     1.0     1.0     1.0     1.0     1.0     1.0     1.0     1.0     1.0     1.0]

# Computed RGBAs
hslhues[3,:] .= 14.5
hues_hsl = XYZToRGB*quavectors(xyY2XYZ(hslhues))
hues_hsl = clamp.(hues_hsl,0,1)

# Mannual adjusted and tested RGBAs
hues_xyY = [0.6345  0.5461  0.4214  0.3396  0.3074  0.2728  0.2251  0.1742  0.1524  0.2052  0.3225  0.5086  0.3294;
            0.3458  0.4147  0.5152  0.5719  0.6035  0.5029  0.3332  0.1502  0.0691  0.1103  0.1652  0.2724  0.3292;
            14.410  14.540  14.520  14.340  14.430  14.540  14.480  14.510  14.420  14.400  14.300  14.480  14.45]

hues_hsl = [0.370   0.225   0.087   0.028   0.006   0.000   0.000   0.000   0.000   0.120   0.280   0.350   0.098;
            0.010   0.051   0.089   0.108   0.115   0.113   0.103   0.078   0.025   0.030   0.012   0.010   0.078;
            0.000   0.000   0.001   0.003   0.001   0.030   0.110   0.410   1.000   0.550   0.280   0.067   0.083;
            1.0     1.0     1.0     1.0     1.0     1.0     1.0     1.0     1.0     1.0     1.0     1.0     1.0]

hueangle_hsl = 0:30:330
hue_hsl = hues_hsl[:,1:end-1]
wp_hsl = repeat(hues_hsl[:,end],inner=(1,size(hue_hsl,2)))

# Plot HSL Hues
IsoLumc = [RGB(hue_hsl[1:3,i]...) for i in 1:size(hue_hsl,2)]
IsoLumhsl = [cos.(deg2rad.(hueangle_hsl)) sin.(deg2rad.(hueangle_hsl))]'

title="HSL_$(length(hueangle_hsl))Hue_Ym"
p=plot()
foreach(i->plot!(p,[0,IsoLumhsl[1,i]],[0,IsoLumhsl[2,i]],color=RGBA(0.5,0.5,0.5,0.5)),1:size(hue_hsl,2))
plot!(p,IsoLumhsl[1,:],IsoLumhsl[2,:],aspectratio=:equal,color=IsoLumc,lw=1.5,markersize=9,marker=:circle,
markerstrokewidth=0,legend=false,xlabel="S",ylabel="S",title=title)
p
foreach(ext->savefig(joinpath(resultdir,"$title$ext")),figfmt)

## Save color data
colordata = Dict{String,Any}("LMS_X"=>colorstring.([minc_lms[:,1],maxc_lms[:,1]]),
                            "LMS_Y"=>colorstring.([minc_lms[:,2],maxc_lms[:,2]]),
                            "LMS_Z"=>colorstring.([minc_lms[:,3],maxc_lms[:,3]]),
                            "LMS_X_WP"=>colorstring.([[th;1],[th;1]]),
                            "LMS_Y_WP"=>colorstring.([[th;1],[th;1]]),
                            "LMS_Z_WP"=>colorstring.([[th;1],[th;1]]),
                            "LMS_XYZ_MichelsonContrast" => diag(mcc),
                            "LMS_Xmcc"=>colorstring.([mminc_lms[:,1],mmaxc_lms[:,1]]),
                            "LMS_Ymcc"=>colorstring.([mminc_lms[:,2],mmaxc_lms[:,2]]),
                            "LMS_Zmcc"=>colorstring.([mminc_lms[:,3],mmaxc_lms[:,3]]),
                            "LMS_Xmcc_WP"=>colorstring.([[th;1],[th;1]]),
                            "LMS_Ymcc_WP"=>colorstring.([[th;1],[th;1]]),
                            "LMS_Zmcc_WP"=>colorstring.([[th;1],[th;1]]),
                            "LMS_XYZmcc_MichelsonContrast" => diag(mmcc),
                            "DKL_X"=>colorstring.([minc_dkl[:,1],maxc_dkl[:,1]]),
                            "DKL_Y"=>colorstring.([minc_dkl[:,2],maxc_dkl[:,2]]),
                            "DKL_Z"=>colorstring.([minc_dkl[:,3],maxc_dkl[:,3]]),
                            "DKL_X_WP"=>colorstring.([[th;1],[th;1]]),
                            "DKL_Y_WP"=>colorstring.([[th;1],[th;1]]),
                            "DKL_Z_WP"=>colorstring.([[th;1],[th;1]]),
                            "HSL_HueYm_Angle" => hueangle_hsl,
                            "HSL_HueYm" => colorstring.([hue_hsl[:,i] for i in 1:size(hue_hsl,2)]),
                            "HSL_HueYm_WP" => colorstring.([wp_hsl[:,i] for i in 1:size(wp_hsl,2)]),
                            "HSL_RGYm" => colorstring.([hue_hsl[:,i] for i in (1,5)]),
                            "HSL_RGYm_WP" => colorstring.([wp_hsl[:,i] for i in (1,5)]),
                            "HSL_YBYm" => colorstring.([hue_hsl[:,i] for i in (3,9)]),
                            "HSL_YBYm_WP" => colorstring.([wp_hsl[:,i] for i in (3,9)]),
                            "HSL_RBYm" => colorstring.([hue_hsl[:,i] for i in (1,9)]),
                            "HSL_RBYm_WP" => colorstring.([wp_hsl[:,i] for i in (1,9)]),
                            "HSL_YGYm" => colorstring.([hue_hsl[:,i] for i in (3,5)]),
                            "HSL_YGYm_WP" => colorstring.([wp_hsl[:,i] for i in (3,5)]),
                            "DKL_HueL0_Angle"=>hueangle_dkl_ilp,
                            "DKL_HueL0"=>colorstring.([hue_dkl_ilp[:,i] for i in 1:size(hue_dkl_ilp,2)]),
                            "DKL_HueL0_WP"=>colorstring.([wp_dkl_ilp[:,i] for i in 1:size(wp_dkl_ilp,2)]),
                            "RGBToLMS" => vec(RGBToLMS),
                            "LMSToRGB" => vec(LMSToRGB),
                            "RGBToXYZ" => vec(RGBToXYZ),
                            "XYZToRGB" => vec(XYZToRGB),
                            "LMSToContrast" => vec(LMSToContrast),
                            "ContrastToLMS" => vec(ContrastToLMS),
                            "LMSToDKL" => vec(LMSToDKL),
                            "DKLToLMS" => vec(DKLToLMS),
                            "DKLToRGB" => vec(DKLToRGB),
                            "RGBToDKL" => vec(RGBToDKL))
YAML.write_file(colordatapath,colordata)

a=rand(3,10)
a_xyz = trivectors(RGBToXYZ*quavectors(a))

a_xyz = [19.01	57.06	3.53	19.01;
        20.00	43.06	6.56	20.00;
        21.78   31.96   2.14    21.78]

camview = cam16view(W=[100,100.00,100],La=200.31,Yb=20)

cam = XYZ2CAM16(a_xyz;camview...)
ucs = CAM16UCS(cam.J,cam.M,cam.h)

a=[1,2,3]

string(a)

rcam = CAM16UCSinv(ucs)

rxyz = CAM162XYZ(cam.J,cam.M,cam.h;camview...)

rxyz = CAM162XYZ(rcam...;camview...)

cam.J

rad2deg(cam.h[2])

RGBToXYZ*[0.5,0.5,0.5,1]

cam16uhs = hcat([[j,m,h] for h in deg2rad.(cam16uniquehue[1,1:end-1]), m in 25:25, j in 15:15]...)
#rcam = CAM16UCSinv(cam16uhs,form=:polar)
#uhs_rgb = XYZToRGB*quavectors(CAM162XYZ(rcam...;camview...))
uhs_rgb = XYZToRGB*quavectors(CAM162XYZ(cam16uhs[1,:],cam16uhs[2,:],cam16uhs[3,:];camview...))

uhs_rgb = quavectors(desaturate2gamut!(trivectors(uhs_rgb)))
IsoLumc = [RGB(uhs_rgb[1:3,i]...) for i in 1:size(uhs_rgb,2)]
IsoLumdkl = RGBToDKL*uhs_rgb
p=plot(IsoLumdkl[2,:],IsoLumdkl[3,:],aspectratio=:equal,color=IsoLumc,lw=1.5,markersize=4.5,marker=:circle,markerstrokewidth=0,legend=false,xlabel="L-M",ylabel="S-(L+M)")
foreach(i->plot!(p,[0,IsoLumdkl[2,i]],[0,IsoLumdkl[3,i]],color=:gray,linestyle=:dot),1:size(hue_dkl_ilp,2))
p

## HSL hues in DKL
hues_hsl=[0.3137 0 0 1;
           0.1922 0.03137 0 1;
           0.07059 0.07059 0 1;
           0.01569 0.08627 0 1;
           0 0.0902 0 1;
           0 0.08627 0.01961 1;
           0 0.07843 0.09412 1;
           0 0.0549 0.3608 1;
           0 0 0.9843 1;
           0.1255 0 0.5686 1;
           0.2314 0 0.2353 1;
           0.2902 0 0.05882 1;
           0.07843 0.05882 0.06275 1]'
huewp_hsl_dkl = trivectors(RGBToDKL*hues_hsl)

plot(huewp_hsl_dkl[2,1:end-1],huewp_hsl_dkl[3,1:end-1],markersize=9,marker=:circle,color=cm_hsl.colors[1:30:331],xlabel="L-M",ylabel="S-(L+M)",frame=:origin,
annotations=[(huewp_hsl_dkl[2,i]+0.05,huewp_hsl_dkl[3,i],text(hueangle_hsl[i],7,:gray5,:bottom,:left)) for i in eachindex(hueangle_hsl)])
scatter!(huewp_hsl_dkl[2,end:end],huewp_hsl_dkl[3,end:end],markersize=9,marker=:circle,color=:gray,leg=false)
foreach(ext->savefig(joinpath(resultdir,"HSL_$(length(hueangle_hsl))Hue_Ym_KL$ext")),figfmt)

plot(huewp_hsl_dkl[2,1:end-1],huewp_hsl_dkl[1,1:end-1],markersize=9,marker=:circle,color=cm_hsl.colors[1:30:331],xlabel="L-M",ylabel="Lum",frame=:origin,
annotations=[(huewp_hsl_dkl[2,i]+0.05,huewp_hsl_dkl[1,i],text(hueangle_hsl[i],7,:gray5,:bottom,:left)) for i in eachindex(hueangle_hsl)])
scatter!(huewp_hsl_dkl[2,end:end],huewp_hsl_dkl[1,end:end],markersize=9,marker=:circle,color=:gray,leg=false)
foreach(ext->savefig(joinpath(resultdir,"HSL_$(length(hueangle_hsl))Hue_Ym_KD$ext")),figfmt)

## Cone Isolating colormap
L➕ = RGBA(maxc_lms[:,1]...)
L➖ = RGBA(minc_lms[:,1]...)
M➕ = RGBA(maxc_lms[:,2]...)
M➖ = RGBA(minc_lms[:,2]...)
S➕ = RGBA(maxc_lms[:,3]...)
S➖ = RGBA(minc_lms[:,3]...)

cs = range(L➖,L➕,length=360)
cm_liso = (colors=cs,notes="The L Cone Isolating colors linearly generated between[min,max], constrained by a `$displayname` LCD Display.")

cs = range(M➖,M➕,length=360)
cm_miso = (colors=cs,notes="The M Cone Isolating colors linearly generated between[min,max], constrained by a `$displayname` LCD Display.")

cs = range(S➖,S➕,length=360)
cm_siso = (colors=cs,notes="The S Cone Isolating colors linearly generated between[min,max], constrained by a `$displayname` LCD Display.")

cs = range(L➖,RGBA(1,1,1,1.0),L➕,length=360)
cm_lisow = (colors=cs,notes="The L Cone Isolating colors linearly generated between[min,white,max], constrained by a `$displayname` LCD Display.")

cs = range(M➖,RGBA(1,1,1,1.0),M➕,length=360)
cm_misow = (colors=cs,notes="The M Cone Isolating colors linearly generated between[min,white,max], constrained by a `$displayname` LCD Display.")

cs = range(S➖,RGBA(1,1,1,1.0),S➕,length=360)
cm_sisow = (colors=cs,notes="The S Cone Isolating colors linearly generated between[min,white,max], constrained by a `$displayname` LCD Display.")

cs = range(L➖,RGBA(0,0,0,1.0),L➕,length=360)
cm_lisob = (colors=cs,notes="The L Cone Isolating colors linearly generated between[min,black,max], constrained by a `$displayname` LCD Display.")

cs = range(M➖,RGBA(0,0,0,1.0),M➕,length=360)
cm_misob = (colors=cs,notes="The M Cone Isolating colors linearly generated between[min,black,max], constrained by a `$displayname` LCD Display.")

cs = range(S➖,RGBA(0,0,0,1.0),S➕,length=360)
cm_sisob = (colors=cs,notes="The S Cone Isolating colors linearly generated between[min,black,max], constrained by a `$displayname` LCD Display.")

## DKL Axis Isolating colormap
Lum➕ = RGBA(maxc_dkl[:,1]...)
Lum➖ = RGBA(minc_dkl[:,1]...)
L_M➕ = RGBA(maxc_dkl[:,2]...)
L_M➖ = RGBA(minc_dkl[:,2]...)
S_LM➕ = RGBA(maxc_dkl[:,3]...)
S_LM➖ = RGBA(minc_dkl[:,3]...)

cs = range(Lum➖,Lum➕,length=360)
cm_lumiso = (colors=cs,notes="The DKL Lum Axis colors linearly generated between[min,max], constrained by a `$displayname` LCD Display.")

cs = range(L_M➖,L_M➕,length=360)
cm_lmiso = (colors=cs,notes="The DKL L-M Axis colors linearly generated between[min,max], constrained by a `$displayname` LCD Display.")

cs = range(S_LM➖,S_LM➕,length=360)
cm_slmiso = (colors=cs,notes="The DKL S-(L+M) Axis colors linearly generated between[min,max], constrained by a `$displayname` LCD Display.")

cs = range(L_M➖,RGBA(1,1,1,1.0),L_M➕,length=360)
cm_lmisow = (colors=cs,notes="The DKL L-M Axis colors linearly generated between[min,white,max], constrained by a `$displayname` LCD Display.")

cs = range(S_LM➖,RGBA(1,1,1,1.0),S_LM➕,length=360)
cm_slmisow = (colors=cs,notes="The DKL S-(L+M) Axis colors linearly generated between[min,white,max], constrained by a `$displayname` LCD Display.")

cs = range(L_M➖,RGBA(0,0,0,1.0),L_M➕,length=360)
cm_lmisob = (colors=cs,notes="The DKL L-M Axis colors linearly generated between[min,black,max], constrained by a `$displayname` LCD Display.")

cs = range(S_LM➖,RGBA(0,0,0,1.0),S_LM➕,length=360)
cm_slmisob = (colors=cs,notes="The DKL S-(L+M) Axis colors linearly generated between[min,black,max], constrained by a `$displayname` LCD Display.")

## DKL IsoLum Plane colormap
angles = 0:360;lum=0
DKLIsoLumRGBVec = trivectors(hcat(map(i->DKLToRGB*RotateXYZMatrix(deg2rad(i),dims=1)*[0,1,0,0], angles)...))
anglec = hcat(map(i->intersectlineunitorigincube((DKLToRGB*[lum,0,0,1])[1:3],DKLIsoLumRGBVec[:,i]),1:size(DKLIsoLumRGBVec,2))...)
dkl_ilp = clamp.(quavectors(anglec),0,1)

cs = [RGBA(dkl_ilp[:,i]...) for i in 1:size(dkl_ilp,2)]
cm_lumiso_plane = (colors=cs,notes="The exact DKL max_cone_contrast colors angled[0,360] in lum=$lum plane, constrained by a `$displayname` LCD Display.")
plotcolormap(cm_lumiso_plane.colors,xlabel="L-M",ylabel="S-(L+M)")

# Linear colormap generated from DKL L-M and S-(L+M) axis colors
cs = range(L_M➕,S_LM➕,L_M➖,S_LM➖,L_M➕,length=365)
cm_lumiso_planel = (colors=cs,notes="The DKL colors linearly generated from max_cone_contrast colors of L-M and S-(L+M) axis at lum=$lum plane, constrained by a `$displayname` LCD Display.")
plotcolormap(cm_lumiso_planel.colors,xlabel="L-M",ylabel="S-(L+M)")

## DKL IsoLM Plane colormap
angles = 0:360;lm=0
DKLIsoLMRGBVec = trivectors(hcat(map(i->DKLToRGB*RotateXYZMatrix(deg2rad(i),dims=2)*[0,0,1,0], angles)...))
anglec = hcat(map(i->intersectlineunitorigincube((DKLToRGB*[0,lm,0,1])[1:3],DKLIsoLMRGBVec[:,i]),1:size(DKLIsoLMRGBVec,2))...)
dkl_ilmp = clamp.(quavectors(anglec),0,1)

cs = [RGBA(dkl_ilmp[:,i]...) for i in 1:size(dkl_ilmp,2)]
cm_lmiso_plane = (colors=cs,notes="The exact DKL max_cone_contrast colors angled[0,360] in lm=$lm plane, constrained by a `$displayname` LCD Display.")
plotcolormap(cm_lmiso_plane.colors,xlabel="S-(L+M)",ylabel="Lum")

# Linear colormap generated from DKL S-(L+M) and Lum axis colors
cs = range(S_LM➕,Lum➕,S_LM➖,Lum➖,S_LM➕,length=365)
cm_lmiso_planel = (colors=cs,notes="The DKL colors linearly generated from max_cone_contrast colors of S-(L+M) and Lum axis at lm=$lm plane, constrained by a `$displayname` LCD Display.")
plotcolormap(cm_lmiso_planel.colors,xlabel="S-(L+M)",ylabel="Lum")

## DKL IsoSLM Plane colormap
angles = 0:360;slm=0
DKLIsoSLMRGBVec = trivectors(hcat(map(i->DKLToRGB*RotateXYZMatrix(deg2rad(i),dims=3)*[0,-1,0,0], angles)...))
anglec = hcat(map(i->intersectlineunitorigincube((DKLToRGB*[0,0,slm,1])[1:3],DKLIsoSLMRGBVec[:,i]),1:size(DKLIsoSLMRGBVec,2))...)
dkl_islmp = clamp.(quavectors(anglec),0,1)

cs = [RGBA(dkl_islmp[:,i]...) for i in 1:size(dkl_islmp,2)]
cm_slmiso_plane = (colors=cs,notes="The exact DKL max_cone_contrast colors angled[0,360] in slm=$slm plane, constrained by a `$displayname` LCD Display.")
plotcolormap(cm_slmiso_plane.colors,xlabel="M-L",ylabel="Lum")

# Linear colormap generated from DKL M-L and Lum axis colors
cs = range(L_M➖,Lum➕,L_M➕,Lum➖,L_M➖,length=365)
cm_slmiso_planel = (colors=cs,notes="The DKL colors linearly generated from max_cone_contrast colors of M-L and Lum axis at slm=$slm plane, constrained by a `$displayname` LCD Display.")
plotcolormap(cm_slmiso_planel.colors,xlabel="M-L",ylabel="Lum")

## HSL Hues colormap
l=0.4
cs = map(i->RGBA(HSLA(i,1,l,1)),0:360)
cm_hsl_hue = (colors=cs, notes="The HSL max_saturated hues angled[0,360] at L=$l.")
plotcolormap(cm_hsl_hue.colors,xlabel="S",ylabel="S")

## Save Color Maps
save(joinpath(resultdir,"colormaps.jld2"),"lms_mccliso",cm_liso,"lms_mccmiso",cm_miso,"lms_mccsiso",cm_siso,
    "lms_mcclisow",cm_lisow,"lms_mccmisow",cm_misow,"lms_mccsisow",cm_sisow,
    "lms_mcclisob",cm_lisob,"lms_mccmisob",cm_misob,"lms_mccsisob",cm_sisob,
    "dkl_mcclumiso",cm_lumiso,"dkl_mcclmiso",cm_lmiso,"dkl_mccslmiso",cm_slmiso,
    "dkl_mcclmisow",cm_lmisow,"dkl_mccslmisow",cm_slmisow,
    "dkl_mcclmisob",cm_lmisob,"dkl_mccslmisob",cm_slmisob,
    "dkl_mcchue_l$lum",cm_lumiso_plane,"lidkl_mcchue_l$lum",cm_lumiso_planel,
    "dkl_mcchue_lm$lm",cm_lmiso_plane,"lidkl_mcchue_lm$lm",cm_lmiso_planel,
    "dkl_mcchue_slm$slm",cm_slmiso_plane,"lidkl_mcchue_slm$slm",cm_slmiso_planel,"hsl_mshue_l$l",cm_hsl_hue)

max_stars_repo_name: VLAB-ExperimentSystem/VLab

max_issues_repo_name: VLAB-ExperimentSystem/VLab

max_forks_repo_name: VLAB-ExperimentSystem/VLab


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

JLLWrappers.@generate_main_file_header("Gloo")
JLLWrappers.@generate_main_file("Gloo", UUID("5909ec93-c24e-5240-9f48-bc73cd7db194"))
end  # module Gloo_jll

ext: jl


text:
abstract type Tri end
abstract type Tet end

"""
$(TYPEDEF)

Data structure for holding geometry data.
"""
struct eMesh{T1<:Union{Nothing,Tri},T2<:Union{Nothing,Tet}}
    point::Vector{SVector{3,Float64}}
    tri::Union{Nothing,Vector{SVector{3,Int64}}}
    tet::Union{Nothing,Vector{SVector{4,Int64}}}
    ϵ::Union{Nothing,Vector{Float64}}
    function eMesh( point::Vector{SVector{3,Float64}},
                    tri::Union{Nothing,Vector{SVector{3,Int64}}},
                    tet::Union{Nothing,Vector{SVector{4,Int64}}}=nothing,
                    ϵ::Union{Nothing,Vector{Float64}}=nothing)

        T1_ = ifelse(tri == nothing, Nothing, Tri)
        T2_ = ifelse(tet == nothing, Nothing, Tet)
        if T2_ == Tet
            @assert(isa(ϵ, Vector{Float64}))
            @assert(length(ϵ) == length(point), "length(ϵ) = $(length(ϵ)) but length(point) = $(length(point))")
            (length(ϵ) != 0) && @assert(0.0 < maximum(ϵ), "normalized penetration extent must be non-negative")
            (length(ϵ) != 0) && @assert(minimum(ϵ) == 0.0, "normalized penetration extent must be zero on the surface of the volume mesh")
            for k = 1:length(tet)
                (0.0 < volume(point[tet[k]])) || error("inverted tetrahedron")
            end
        else
            @assert(ϵ == nothing)
        end
        (T1_ == T2_ == Nothing) && error("a whole lot of nothing")
        return new{T1_,T2_}(deepcopy(point), deepcopy(tri), deepcopy(tet), deepcopy(ϵ))
    end
    function eMesh(hm::HomogenousMesh, tet::Union{Nothing,Vector{SVector{4,Int64}}}=nothing,
            ϵ::Union{Nothing,Vector{Float64}}=nothing)

        point = [SVector{3,Float64}(k) for k = hm.vertices]
        tri = [SVector{3,Int64}(k) for k = hm.faces]
        return eMesh(point, tri, tet, ϵ)
    end
    eMesh{Tri,Nothing}() = eMesh(Vector{SVector{3,Float64}}(), Vector{SVector{3,Int64}}(), nothing, nothing)
    eMesh{Nothing,Tet}() = eMesh(Vector{SVector{3,Float64}}(), nothing, Vector{SVector{4,Int64}}(), Vector{Float64}())
    eMesh{Tri,Tet}()     = eMesh(Vector{SVector{3,Float64}}(), Vector{SVector{3,Int64}}(), Vector{SVector{4,Int64}}(), Vector{Float64}())
end

"""
$(SIGNATURES)

Converts an `eMesh` to a mesh that contains only tetrahedrons.
You need to do this before adding an `eMesh` as contact geometry.
"""
as_tet_eMesh(e_mesh::eMesh{Tri,Tet}) = eMesh(e_mesh.point, nothing, e_mesh.tet, e_mesh.ϵ)
as_tet_eMesh(e_mesh::eMesh{Nothing,Tet}) = deepcopy(e_mesh)

"""
$(SIGNATURES)

Converts an `eMesh` to a mesh that contains only triangles.
You need to do this before adding an `eMesh` as contact geometry.
"""
as_tri_eMesh(e_mesh::eMesh{Tri,Tet}, is_repair::Bool=true) = eMesh(e_mesh.point, e_mesh.tri, nothing, nothing)
as_tri_eMesh(e_mesh::eMesh{Tri,Nothing}, is_repair::Bool=true) = deepcopy(e_mesh)
function as_tri_eMesh(eM::eMesh{Nothing,Tet}, is_repair::Bool=true)
	i3 = Vector{SVector{3,Int64}}()
	for k = 1:n_tet(eM)
		i_tet = eM.tet[k]
		ϵ_tet = eM.ϵ[i_tet]
		i_tet_new = sort_so_big_ϵ_last(ϵ_tet, i_tet)
		push!(i3, i_tet_new[1:3])
	end
	eM_tri = eMesh(eM.point, i3)
	if is_repair
		mesh_repair!(eM_tri)
	end
	return eM_tri
end

vertex_pos_for_tri_ind(eM::eMesh{Tri,T2}, k::Int64) where {T2} = eM.point[eM.tri[k]]
vertex_pos_for_tet_ind(eM::eMesh{T1,Tet}, k::Int64) where {T1} = eM.point[eM.tet[k]]

n_point(eM::eMesh) = length(eM.point)
n_tri(eM::eMesh) = length(eM.tri)
n_tet(eM::eMesh) = length(eM.tet)
get_tri(eM::eMesh) = eM.tri
get_tet(eM::eMesh) = eM.tet
get_point(eM::eMesh) = eM.point

function Base.empty!(e_mesh::eMesh{T1,T2}) where {T1,T2}
    empty!(e_mesh.point)
    (T1 == Nothing) || empty!(e_mesh.tri)
    (T2 == Nothing) || empty!(e_mesh.tet)
    (T2 == Nothing) || empty!(e_mesh.ϵ)
    return nothing
end

function Base.isempty(e_mesh::eMesh{T1,T2}) where {T1,T2}
    is_emp = isempty(e_mesh.point)
    if T1 == Tri
        is_emp = is_emp || isempty(e_mesh.tri)
    end
    if T2 == Tet
        is_emp = is_emp || isempty(e_mesh.tet)
        is_emp = is_emp || isempty(e_mesh.ϵ)
    end
    return is_emp
end

function Base.append!(eM_1::eMesh{T1,T2}, eM_2::eMesh{T1,T2}) where {T1,T2}
    n_1 = n_point(eM_1)
    if T1 != Nothing
        for k = eM_2.tri
            push!(eM_1.tri, k .+ n_1)
        end
    end
    if T2 != Nothing
        for k = eM_2.tet
            push!(eM_1.tet, k .+ n_1)
        end
        append!(eM_1.ϵ, eM_2.ϵ)
    end
    append!(eM_1.point, eM_2.point)
    return nothing
end

### VERIFICATION
verify_mesh(eM::eMesh{Tri,Nothing}) = verify_mesh_tri(eM)
verify_mesh(eM::eMesh{Nothing,Tet}) = verify_mesh_tet(eM)
function verify_mesh(eM::eMesh{Tri,Tet})
	eM_tri = as_tri_eMesh(eM)
    eM_tet = as_tet_eMesh(eM)
	verify_mesh_tri(eM_tri)
	verify_mesh_tet(eM_tet)
    eM_tri_2 = as_tri_eMesh(eM_tet)
    (area(eM_tri) ≈ area(eM_tri_2)) || error("surfaces where penetration extent is 0.0 is not the same as the tri surface")
	return nothing
end

function verify_mesh_tet(eM::eMesh{Nothing,Tet})
    length_tet = n_tet(eM)
    length_point = n_point(eM)
    length_ϵ = length(eM.ϵ)
    for k = 1:length_tet
        iΔ = eM.tet[k]
        point_Δ = eM.point[iΔ]
        (0.0 < volume(point_Δ)) || error("tet with negative volume")
    end
    (length_ϵ == length_point) || error("number of points not equal to number or ϵ")
    if !isempty(eM.ϵ)
        (minimum(eM.ϵ) < 0.0) && error("normalized penetration extent must be non-negative")
        (minimum(eM.ϵ) == 0.0) || error("normalized penetration extent must be 0.0 somewhere")
    end
    return nothing
end

function verify_mesh_tri(eM::eMesh{Tri,Nothing})
    for k = 1:n_tri(eM)
        iΔ = eM.tri[k]
        point_Δ = eM.point[iΔ]
    end
    return nothing
end

### MESH MANIPULATION
"""
$(SIGNATURES)

Creates a transformation matrix by creating a `basic_dh` from `extra_arguments` and then transforms the `eMesh`.
See [`basic_dh`](@ref) for how to create common transformation types.
Reflections are not allowed.
"""
function transform!(e_mesh::eMesh{T1,T2}, extra_args...) where {T1,T2}
    dh = basic_dh(extra_args...)
    point = e_mesh.point
    for k = 1:n_point(e_mesh)
        point[k] = dh_vector_mul(dh, point[k])
    end
end

struct PointObj
    p::SVector{3,Float64}
    w::Float64
    b::Bool
    k::Int64
end

function crop_mesh(eM::eMesh{Tri,Nothing}, plane::SMatrix{1,4,Float64,4})
    function add_cropped_triangle!(eM_new::eMesh{Tri,Nothing}, iΔ::SVector{3,Int64})
        function plane_dot_2(eM_new::eMesh{Tri,Nothing}, plane::SMatrix{1,4,Float64,4}, k::Int64)
            p = get_point(eM_new)[k]
            w = dot(plane[:], onePad(p))
            b = -1.0e-12 < w
            return PointObj(p, w, b, k)
        end

        o1 = plane_dot_2(eM_new, plane, iΔ[1])
        o2 = plane_dot_2(eM_new, plane, iΔ[2])
        o3 = plane_dot_2(eM_new, plane, iΔ[3])
        s_bool = o1.b + o2.b + o3.b
        (s_bool == 3) && (push!(eM_new.tri, iΔ))  # all entire triangle
        ((s_bool == 0) || (s_bool == 3)) && (return nothing)  # exit function
        (o1.b == o2.b) && ((o1, o2, o3) = (o2, o3, o1))  # put 1 and 2 on opposite sides
        (o2.b == o3.b) && ((o1, o2, o3) = (o3, o1, o2))  # put 2 and 3 on opposite sides
        (o1.b == o2.b) || (o2.b == o3.b) && error("something is wrong")  # sanity check
        push!(eM_new.point, weightPoly(o1.p, o2.p, o1.w, o2.w))
        i12 = length(get_point(eM_new))
        push!(eM_new.point, weightPoly(o2.p, o3.p, o2.w, o3.w))
        i23 = length(get_point(eM_new))
        if o2.b  # add top of triangle
            push!(eM_new.tri, SVector{3,Int64}(i12, o2.k, i23))
        else  # add bottom quadralateral
            push!(eM_new.tri, SVector{3,Int64}(o1.k, i12, i23))
            push!(eM_new.tri, SVector{3,Int64}(o1.k, i23, o3.k))
        end
    end

    eM_new = eMesh{Tri,Nothing}()
    append!(eM_new.point, eM.point)
    for iΔ = eM.tri
        add_cropped_triangle!(eM_new, iΔ)
    end
    mesh_repair!(eM_new)
    return eM_new
end

function invert!(eM::eMesh{Tri,Nothing})
    for k = 1:n_tri(eM)
        tri_k = eM.tri[k]
        eM.tri[k] = SVector{3,Float64}(tri_k[3], tri_k[2], tri_k[1])
    end
end

### MESH REPAIR
function mesh_repair!(e_mesh::eMesh{T1,T2}) where {T1,T2}
    mesh_remove_unused_points!(e_mesh)
    mesh_inplace_rekey!(e_mesh)
    mesh_remove_unused_points!(e_mesh)
    return delete_triangles!(e_mesh)
end

function remove_degenerate!(eM::eMesh)
	area_or_vol(v::SVector{3,SVector{3,Float64}}) = area(v)
	area_or_vol(v::SVector{4,SVector{3,Float64}}) = volume(v)

	tol = 1.0e-6
	for vec_ind = (eM.tet, eM.tri)
		if vec_ind != nothing
			v = [area_or_vol(eM.point[vec_ind[k]]) for k = 1:length(vec_ind)]
			max_v = maximum(v)
			i_bad = findall(v .< max_v * tol)
			deleteat!(vec_ind, i_bad)
		end
	end
end

rekey!(v::Nothing, i::Vector{Int64}) = nothing
rekey!(v::Vector{SVector{N,Int64}}, i::Vector{Int64}) where {N} = replace!(x -> i[x], v)

function mesh_inplace_rekey(e_mesh::eMesh{T1,T2}) where {T1,T2}
    e_mesh_copy = deepcopy(e_mesh)
    mesh_inplace_rekey!(e_mesh_copy)
    return e_mesh_copy
end

function mesh_inplace_rekey!(e_mesh::eMesh{T1,T2}) where {T1,T2}
    # Expresses the mesh using the minumum number of points

    function shortest_side(e_mesh::eMesh)
        shortest_side(v::Nothing) = Inf
        shortest_side(ind::Vector{SVector{N,Int64}}) where {N} =  minimum([shortest_side(e_mesh.point[k]) for k = ind])
        function shortest_side(v::SVector{N,SVector{3,Float64}}) where {N}
            d_min = Inf
            for k = 1:N
                for kk = 1:(k-1)
                    d_min = min(d_min, norm(v[k] - v[kk]))
                end
            end
            return d_min
        end

        return min(shortest_side(e_mesh.tri), shortest_side(e_mesh.tet))
    end

    function inplace_rekey(point::Vector{SVector{3,Float64}}, min_side_length::Float64)
        balltree = BallTree(point; reorder = false)
        idxs = inrange(balltree, point, min_side_length * 0.499)  # Assume that all points closer to each other than hald a side length are duplicates
        return first.(sort!.(idxs))
    end

    if !isempty(e_mesh)  # reducing over an empty collection is not allowed
        min_side_length = shortest_side(e_mesh)
        new_key = inplace_rekey(e_mesh.point, min_side_length)
        rekey!(e_mesh.tri, new_key)
        rekey!(e_mesh.tet, new_key)
        mesh_remove_unused_points!(e_mesh)
        return new_key
    end
    return nothing
end

function mesh_remove_unused_points!(e_mesh::eMesh{T1,T2}) where {T1,T2}
    truth_vector = falses(n_point(e_mesh))
    (T1 != Nothing) && (for k = e_mesh.tri; truth_vector[k] = true; end)
    (T2 != Nothing) && (for k = e_mesh.tet; truth_vector[k] = true; end)
    new_key = cumsum(truth_vector)
    rekey!(e_mesh.tri, new_key)
    rekey!(e_mesh.tet, new_key)
    ind_truth = findall(truth_vector)
    point_new = e_mesh.point[ind_truth]
    empty!(e_mesh.point)
    append!(e_mesh.point, point_new)
    if T2 != Nothing
        ϵ_new = e_mesh.ϵ[ind_truth]
        empty!(e_mesh.ϵ)
        append!(e_mesh.ϵ, ϵ_new)
    end
    return nothing
end

delete_triangles!(e_mesh::eMesh{Nothing,T2}) where {T2} = -9999
function delete_triangles!(e_mesh::eMesh{Tri,T2}) where {T2}
	# TODO: fix this to only delete triangles for connected meshes

    ### make the first index lowest
    for k = 1:n_tri(e_mesh)
        iΔ = e_mesh.tri[k]
        min_index = minimum(iΔ)
        (min_index == iΔ[1]) || (iΔ = SVector(iΔ[3], iΔ[1], iΔ[2]))
        (min_index == iΔ[1]) || (iΔ = SVector(iΔ[3], iΔ[1], iΔ[2]))
    end

    ### create a dictionary of key repition counts
    key_type = Vector{Tuple{Int64,SVector{3,Int64}}}
    dict_ind = Dict{SVector{3,Int64},key_type}()
    for k = 1:n_tri(e_mesh)
        iΔ = e_mesh.tri[k]
        sort_iΔ = sort(iΔ)
        !haskey(dict_ind, sort_iΔ) && (dict_ind[sort_iΔ] = key_type() )
        push!(dict_ind[sort_iΔ], (k, iΔ))
    end

    ### check for issues is triangle pairs
    i_delete = Vector{Int64}()
    for (key_k, val_k) = dict_ind
        length_val_k = length(val_k)
        if length_val_k == 2
            val_k1 = val_k[1]
            val_k2 = val_k[2]
            (val_k1[2] == val_k2[2]) && error("non-opposing triangles")
            push!(i_delete, val_k1[1], val_k2[1])
        elseif 3 <= length_val_k
            error("something is wrong")
        end
    end

    ### Actually delete the opposing triangles
    i_delete = sort(i_delete)
    deleteat!(e_mesh.tri, i_delete)
	return length(i_delete)
end

function sub_div_mesh(eM_ico::eMesh{Tri,T2}, n_div::Int64) where {T2}
    n_end(n::Int64) = div((n + 1) * n, 2)
    n_start(n::Int64) = 1 + n_end(n - 1)

    function sub_div_triangle(p::SVector{3,SVector{3,Float64}}, n_div::Int64)
        function sub_div_triangle_vert_index(n_div::Int64)
            i_tri = Vector{SVector{3,Int64}}()
            for k = 1:n_div
                for kk = 0:(k - 1)
                    i1 = n_start(k) + kk
                    i2 = i1 + k
                    i3 = i2 + 1
                    push!(i_tri, SVector{3,Int64}(i1, i2, i3))
                end
                for kk = 0:(k - 2)
                    i1 = n_start(k) + kk
                    i2 = i1 + k + 1
                    i3 = i2 - k
                    push!(i_tri, SVector{3,Int64}(i1, i2, i3))
                end
            end
            return i_tri
        end

        function get_new_point(n_vert::Int64, n_div::Int64, p::SVector{3,SVector{3,Float64}})
            function find_layer_1(i_point::Int64)
                i_end_layer = 1
                i_layer = 1
                while i_end_layer < i_point
                    i_layer += 1
                    i_end_layer += i_layer
                end
                norm_extent_1 = ifelse(i_point == 1, 0.0, (i_end_layer - i_point) / (i_layer - 1) )
                return i_layer, norm_extent_1
            end

            i_layer, norm_extent_1 = find_layer_1(n_vert)
            ϕ_1 = (n_div - i_layer + 1) / n_div
            ϕ_2 = (1 - ϕ_1) * norm_extent_1
            ϕ_3 = 1 - ϕ_1 - ϕ_2
            ϕ = SVector{3,Float64}(ϕ_1, ϕ_2, ϕ_3)
            return sum(p .* ϕ)
        end

        point = Vector{SVector{3,Float64}}()
        i_tri_div = sub_div_triangle_vert_index(n_div)
        for k = 1:n_end(n_div + 1)
            push!(point, get_new_point(k, n_div, p))
        end
        return eMesh(point, i_tri_div, nothing, nothing)
    end

    eM_ico_div = eMesh{Tri,Nothing}()
    for k = 1:n_tri(eM_ico)
        p = vertex_pos_for_tri_ind(eM_ico, k)
        append!(eM_ico_div, sub_div_triangle(p, n_div))
    end
    mesh_repair!(eM_ico_div)
    return eM_ico_div
end

### BASIC SHAPES
"""
$(SIGNATURES)

Outputs an `eMesh` for a half-plane.
"""
function eMesh_half_plane(plane_w::Float64=1.0, is_include_vis_sides::Bool=false)
    v1, v2, v3 = [SVector{3,Float64}(cos(theta), sin(theta), 0.0) for theta = (0.0, 2*pi/3, 4*pi/3)]
    v4 = SVector{3,Float64}(0.0, 0.0, -1.0) * plane_w
    point = [v1, v2, v3, v4]
    if is_include_vis_sides
        tri = [SVector{3,Int64}(1,2,3), SVector{3,Int64}(1,3,4),  SVector{3,Int64}(1,4,2), SVector{3,Int64}(2,4,3)]
    else
        tri = [SVector{3,Int64}(1,2,3)]
    end
    tet = [SVector{4,Int64}(4,1,2,3)]
    ϵ = [0.0, 0.0, 0.0, plane_w]
    return eMesh(point, tri, tet, ϵ)
end

"""
$(SIGNATURES)

Outputs an `eMesh` for a sphere. Larger values of `n_div` create finer discretizations.
"""
function eMesh_sphere(rad::Union{Float64,SVector{3,Float64}}=1.0, n_div::Int64=4)
    function make_icosahedron()
        φ = Base.MathConstants.golden

        v = Vector{SVector{3,Float64}}()
        for s_1 = (-1.0, +1.0)
            for s_2 = (-1.0, +1.0)
                push!(v, SVector{3,Float64}(    0.0,     s_1, φ * s_2) )
                push!(v, SVector{3,Float64}(    s_1, φ * s_2,     0.0) )
                push!(v, SVector{3,Float64}(φ * s_2,     0.0,     s_1) )
            end
        end

        n_ico_vert = 12
        d = zeros(n_ico_vert, n_ico_vert)
        for k = 1:n_ico_vert
            for kk = 1:n_ico_vert
                d[k, kk] = norm(v[k] - v[kk])
            end
        end

        v_face_vert = Vector{SVector{3,Int64}}()
        b = d .== 2.0
        for i1 = 1:n_ico_vert
            for i2 = 1:n_ico_vert
                for i3 = 1:n_ico_vert
                    if (i1 < i2 < i3) && b[i1, i2] && b[i2, i3] && b[i1, i3]
                        p1 = v[i1]
                        p2 = v[i2]
                        p3 = v[i3]
                        n = triangleNormal(p1, p2, p3)
                        c = normalize(centroid(p1, p2, p3))
                        i_face = ifelse(n ≈ c, SVector{3,Int64}(i1, i2, i3), SVector{3,Int64}(i1, i3, i2))
                        push!(v_face_vert, i_face)
                    end
                end
            end
        end

        push!(v, SVector{3,Float64}(0,0,0))
        v_tet = Vector{SVector{4,Int64}}()
        for k = 1:length(v_face_vert)
            push!(v_tet, SVector{4,Int64}(13, v_face_vert[k]...))
        end
        ϵ = vcat(zeros(n_ico_vert), 1.0)
        return eMesh(v, v_face_vert, v_tet, ϵ)
    end

    function volumize_about(eM::eMesh{Tri,Nothing})
        i_tet = Vector{SVector{4,Int64}}()
        n_vert = n_point(eM)
        n_center = n_vert + 1
        for k = 1:n_tri(eM)
            push!(i_tet, SVector{4,Int64}(n_center, eM.tri[k]...))
        end
        ϵ = vcat(zeros(n_vert), 1.0)
        point = deepcopy(eM.point)
        push!(point, zeros(SVector{3,Float64}))
        return eMesh(point, eM.tri, i_tet, ϵ)
    end

    function project_to_sphere!(eM::eMesh)
        for k = 1:n_point(eM)
            eM.point[k] = normalize(eM.point[k])
        end
        return nothing
    end

    eM_ico = make_icosahedron()
    eM_ico_div = sub_div_mesh(eM_ico, n_div)
    project_to_sphere!(eM_ico_div)

    rad = ones(SVector{3,Float64}) .* rad
    transform!(eM_ico_div, Diagonal(rad))

    return volumize_about(eM_ico_div)
end

function output_box_ind()
    oriented_box_faces = (
        SVector{4, Int64}(1,3,5,7),  # -x
        SVector{4, Int64}(2,6,4,8),  # +x
        SVector{4, Int64}(1,5,2,6),  # -y
        SVector{4, Int64}(3,4,7,8),  # +y
        SVector{4, Int64}(1,2,3,4),  # -z
        SVector{4, Int64}(5,7,6,8),  # +z
    )
    ϵ = zeros(Float64,8)
    push!(ϵ, 1.0)
    tri = Vector{SVector{3,Int64}}()
    tet = Vector{SVector{4,Int64}}()
    for k = 1:6
        bf_k = oriented_box_faces[k]
        push!(tri, bf_k[SVector{3,Int64}(1,3,4)])
        push!(tri, bf_k[SVector{3,Int64}(1,4,2)])
    end
    for k = 1:length(tri)
        tri_k = tri[k]
        push!(tet, SVector{4,Int64}(9, tri_k[1], tri_k[2], tri_k[3]))
    end
    return tri, tet, ϵ
end

"""
$(SIGNATURES)

Outputs an `eMesh` for a box.
"""
function eMesh_box(r::Union{Float64,SVector{3,Float64}}=1.0, c::SVector{3,Float64}=zeros(SVector{3,Float64}))
    point = [
        SVector{3,Float64}(-1,-1,-1),
        SVector{3,Float64}(+1,-1,-1),
        SVector{3,Float64}(-1,+1,-1),
        SVector{3,Float64}(+1,+1,-1),
        SVector{3,Float64}(-1,-1,+1),
        SVector{3,Float64}(+1,-1,+1),
        SVector{3,Float64}(-1,+1,+1),
        SVector{3,Float64}(+1,+1,+1),
        SVector{3,Float64}( 0, 0, 0),
    ]
    tri, tet, ϵ = output_box_ind()
    e_mesh = eMesh(point, tri, tet, ϵ)
    r = r .* ones(SVector{3,Float64})
    transform!(e_mesh, SMatrix{3,3,Float64,9}(r[1], 0, 0, 0, r[2], 0, 0, 0, r[3]))
    transform!(e_mesh, c)
    return e_mesh
end

const SF3 = SVector{3,Float64}
const SI3 = SVector{3,Int64}
const SI4 = SVector{4,Int64}

function create_2D_circle(rad::Float64=1.0; n::Int64=12)
    points_2D = Vector{SF3}()
    tri_2D = Vector{SI3}()
    i_center = n + 1
    θ_range = LinRange(0.0, 2 * π, n + 1)[2:end]
    for (k, θ) = enumerate(θ_range)
        s, c = sincos(θ)
        push!(points_2D, rad * SF3(c, s, 0.0))
        push!(tri_2D, SI3(k, mod1(k + 1, n), i_center))
    end
    push!(points_2D, SF3(0, 0, 0))
    return eMesh(points_2D, tri_2D)
end

function smallest_first(v::SI4)
	_, i = findmin(v)
	return SI4(v[i], v[mod1(i + 1, 4)], v[mod1(i + 2, 4)], v[mod1(i + 3, 4)])
end

function add_stuff_for_extrude!(tri::Vector{SI3}, tet::Vector{SI4}, k::Int64, i_tri_k::SI3, n_✴_2D::Int64)
	b1, b2, b3 = i_tri_k
	t4, t5, t6 = i_tri_k .+ n_✴_2D
	i_center = k + n_✴_2D * 2
	oriented_slice_faces = (
        SI4(b1,b2,t5,t4),
        SI4(b2,b3,t6,t5),
        SI4(b1,t4,t6,b3)
    )
    tri_add = Vector{SI3}()
    push!(tri_add, SI3(b1,b3,b2))
    push!(tri_add, SI3(t4,t5,t6))
    for bf_k = oriented_slice_faces
		bf_k = smallest_first(bf_k)
		push!(tri_add, bf_k[SI3(1,2,3)])
        push!(tri_add, bf_k[SI3(1,3,4)])
    end
    for tri_k = tri_add
        push!(tri, tri_k)
        push!(tet, SI4(i_center, tri_k[1], tri_k[2], tri_k[3]))
    end
end

function verify_trianges_have_same_normal(surf::eMesh{Tri,Nothing})
	n̂ = [triangleNormal(vertex_pos_for_tri_ind(surf, k)) for k = 1:n_tri(surf)]
	n̂_mean = n̂[1]
	for n̂_k = n̂
		(n̂_k ≈ n̂_mean) || error("All triangles must have the same normal.")
	end
	return n̂_mean
end

"""
$(SIGNATURES)

Extrudes a planar `eMesh` of type `eMesh{Tri,Nothing}` into a type `eMesh{Tri,Tet}` in the direction of the plane
normal. This direction is detected automatically.
The function errors if the points do not lie on a plane.
"""
function extrude_mesh(surf::eMesh{Tri,Nothing}, thick::Float64)
	n̂ = verify_trianges_have_same_normal(surf)

    tri_2D    = surf.tri
    points_2D = surf.point
	n_▲_2D = length(tri_2D )
	n_✴_2D = length(points_2D)
    point_lo = points_2D .- [n̂ * thick / 2]
    point_hi = points_2D .+ [n̂ * thick / 2]
    point_centroid = [centroid(points_2D[k]) for k = tri_2D]
    point = vcat(point_lo, point_hi, point_centroid)
    ϵ = vcat(zeros(2 * n_✴_2D), ones(n_▲_2D))
    tri = Vector{SI3}()
    tet = Vector{SI4}()
    for (k, i_tri_k) = enumerate(tri_2D)
		add_stuff_for_extrude!(tri, tet, k, i_tri_k, n_✴_2D)
    end
    eM_thick = eMesh(point, tri, tet, ϵ)
	mesh_repair!(eM_thick)
    return eM_thick
end

"""
$(SIGNATURES)

Outputs an `eMesh` for a cylinder.
"""
function eMesh_cylinder(rad::Float64=1.0, height::Float64=1.0; n::Int64=6)
	eM_circle = create_2D_circle(rad, n=n)
	return extrude_mesh(eM_circle, height)
end


text:
using SinkholeDynamics
using Documenter

makedocs(;
    modules=[SinkholeDynamics],
    authors="Patrick Eastham <peastham@math.fsu.edu> and contributors",
    repo="https://github.com/pseastham/SinkholeDynamics.jl/blob/{commit}{path}#L{line}",
    sitename="SinkholeDynamics.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://pseastham.github.io/SinkholeDynamics.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/pseastham/SinkholeDynamics.jl",
)


text:
"""
    make_cache(source, cacheidx)

Return a in-memory copy of `source` at observation indices `cacheidx`.
Defaults to `getobs(source, cacheidx)`.
"""
make_cache(source, cacheidx) = getobs(source, cacheidx)

"""
    CachedDataset(source, cachesize = numbobs(source))
    CachedDataset(source, cacheidx = 1:numbobs(source))
    CachedDataset(source, cacheidx, cache)

Wrap a `source` data container and cache `cachesize` samples in memory.
This can be useful for improving read speeds when `source` is a lazy data container,
but your system memory is large enough to store a sizeable chunk of it.

By default the observation indices `1:cachesize` are cached.
You can manually pass in a set of `cacheidx` as well.

See also [`make_cache`](@ref) for customizing the default cache creation for `source`.
"""
struct CachedDataset{T, S} <: AbstractDataContainer
    source::T
    cacheidx::Vector{Int}
    cache::S
end

CachedDataset(source, cachesize::Int) = CachedDataset(source, 1:cachesize)

CachedDataset(source, cacheidx::AbstractVector{<:Integer} = 1:numobs(source)) =
    CachedDataset(source, collect(cacheidx), make_cache(source, cacheidx))

function Base.getindex(dataset::CachedDataset, i::Integer)
    _i = findfirst(==(i), dataset.cacheidx)

    return isnothing(_i) ? getobs(dataset.source, i) : getobs(dataset.cache, _i)
end
Base.length(dataset::CachedDataset) = numobs(dataset.source)


text:
"""
    RProp_params
"""
struct RProp_params{T<:Real}
    eta_plus::T
    eta_minus::T
    Delta_min::T
    Delta_max::T
    Delta_0::T
end

"""
    train_with_RProp(X::Array{Float64,2}, y::Array{Float64,1}, cov_func_parameters::gaussian_kernel_hyperparameters, MaxIter::Int, noise::Float64, params::RProp_params)
Trains kriging hyperparameters with RProp.
"""
function calibrate_by_ML_with_Rprop(X::AbstractArray{T,2}, y::AbstractArray{R,1}, cov_func_parameters::gaussian_kernel_hyperparameters, MaxIter::Integer, noise::Real, params::RProp_params) where T<:Real where R<:Real
    # An implementation of this paper https://pdfs.semanticscholar.org/aa65/042ae494455a14811927eb0574871d276454.pdf
    iter = 0
    w_0 = cov_func_parameters.w_0
    w_i = cov_func_parameters.w_i
    threshold = params.Delta_min * params.eta_plus + 10*eps()
    old_0_sign = 1
    old_i_signs = ones(length(w_i))
    while iter < MaxIter
        marginal_likelihood, K, invK = marginal_likelihood_gaussian_derivatives(X, y, gaussian_kernel_hyperparameters(w_0, w_i), noise)
        Delta0 = params.Delta_0
        Deltas = Array{typeof(w_0),1}(undef, length(w_i)) .= params.Delta_0
        new_0_sign = sign(marginal_likelihood[1])
        new_i_signs = sign.(marginal_likelihood[2:length(marginal_likelihood)])
        w_0 = w_0 + new_0_sign*Delta0
        w_i = w_i .+ new_i_signs .* Deltas
        if iter > 0
            Delta0 = new_0_sign * old_0_sign > 0 ? Delta0*params.eta_plus : Delta0*params.eta_minus
            Delta0 = min(max(Delta0, params.Delta_min ), params.Delta_max )
            for i in 1:length(w_i)
                Deltas[i] = new_i_signs[i]*old_i_signs[i] > 0 ? Deltas[i]*params.eta_plus : Deltas[i] *params.eta_minus
                Deltas[i] = min(max(Deltas[i], params.Delta_min ), params.Delta_max )
            end
        end
        if all(Deltas .< threshold)
            return gaussian_kernel_hyperparameters(w_0, w_i)
        end
        old_0_sign  = new_0_sign
        old_i_signs = new_i_signs
        iter = iter + 1
    end
    return gaussian_kernel_hyperparameters(w_0, w_i)
end

"""
    sample(dim::Int, batch_size::Int, replace::Bool)
This does sampling with or without replacement.
"""
function sample(dim::Integer, batch_size::Integer)
    return randperm(dim)[1:batch_size]
end

"""
    calibrate_by_ML_with_SGD(X::Array{Float64,2}, y::Array{Float64,1}, cov_func_parameters::gaussian_kernel_hyperparameters, steps::Int, batch_size::Int, step_multiple::Float64, noise::Float64, seed::Int = 1988)
This trains a kriging model by using maximum likelihood with stochastic gradient descent.
"""
function calibrate_by_ML_with_SGD(X::AbstractArray{T,2}, y::AbstractArray{R,1}, cov_func_parameters::gaussian_kernel_hyperparameters, steps::Integer, batch_size::Integer, step_multiple::Real, noise::Real, seed::Integer = 1988) where T<:Real where R<:Real
    Random.seed!(seed)
    ow_0 = cov_func_parameters.w_0
    ow_i = cov_func_parameters.w_i
    ndims = length(ow_i)
    nobs = length(y)
    for s in 1:steps
        samples = sample(nobs,batch_size)
        XSample = X[samples, :]
        ySample = y[samples]
        marginal_likelihood, K, invK = marginal_likelihood_gaussian_derivatives(XSample, ySample, gaussian_kernel_hyperparameters(ow_0, ow_i), noise)
        normalised_grad =  sign.(marginal_likelihood) .* ( abs.(marginal_likelihood) ./ (abs.(marginal_likelihood) .+ maximum(abs.(marginal_likelihood))) )
        ow_0 = ow_0 .* (1 .+ normalised_grad[1] .* step_multiple)
        ow_i = ow_i .* (1 .+ normalised_grad[2:(ndims+1)] .* step_multiple)
    end
    return gaussian_kernel_hyperparameters(ow_0, ow_i)
end


text:
# main tests

function safe_mapslices(op, A, region)
    newregion = intersect(region, 1:ndims(A))
    return isempty(newregion) ? A : mapslices(op, A, newregion)
end
safe_sum{T}(A::Array{T}, region) = safe_mapslices(sum, A, region)
safe_prod{T}(A::Array{T}, region) = safe_mapslices(prod, A, region)
safe_maximum{T}(A::Array{T}, region) = safe_mapslices(maximum, A, region)
safe_minimum{T}(A::Array{T}, region) = safe_mapslices(minimum, A, region)

Areduc = rand(3, 4, 5, 6)
for region in {
        1, 2, 3, 4, 5, (1, 2), (1, 3), (1, 4), (2, 3), (2, 4), (3, 4),
    (1, 2, 3), (1, 3, 4), (2, 3, 4), (1, 2, 3, 4)}

    r = fill(NaN, Base.reduced_dims(size(Areduc), region))
    @test_approx_eq sum!(r, Areduc) safe_sum(Areduc, region)
    @test_approx_eq prod!(r, Areduc) safe_prod(Areduc, region)
    @test_approx_eq maximum!(r, Areduc) safe_maximum(Areduc, region)
    @test_approx_eq minimum!(r, Areduc) safe_minimum(Areduc, region)

    @test_approx_eq sum(Areduc, region) safe_sum(Areduc, region)
    @test_approx_eq prod(Areduc, region) safe_prod(Areduc, region)
    @test_approx_eq maximum(Areduc, region) safe_maximum(Areduc, region)
    @test_approx_eq minimum(Areduc, region) safe_minimum(Areduc, region)
end

@test reducedim((a,b) -> a|b, [true false; false false], 1, false) == [true false]
R = reducedim((a,b) -> a+b, [1 2; 3 4], 2, 0.0)
@test eltype(R) == Float64
@test_approx_eq R [3,7]
@test reducedim((a,b) -> a+b, [1 2; 3 4], 1, 0) == [4 6]

# inferred return types
rt = Base.return_types(reducedim, (Function, Array{Float64, 3}, Int, Float64))
@test length(rt) == 1 && rt[1] == Array{Float64, 3}

max_stars_repo_name: tobig/julia

max_issues_repo_name: tobig/julia

max_forks_repo_name: tobig/julia


text:
using MPIMeasurements
using HDF5

filename = "systemMatrix.h5"

# define Grid
shp = [9,1,1]
fov = [90.0,1.0,1.0]Unitful.mm
ctr = [156.0,-11.2,71.0]Unitful.mm
positions = RegularGridPositions(shp,fov,ctr)

scanner = MPIScanner("HeadScanner.toml")
robot = getRobot(scanner)
daq = getDAQ(scanner)
safety = getSafety(scanner)
su = getSurveillanceUnit(scanner)

params = toDict(daq.params)

params["studyName"]="TestTobi"
params["studyDescription"]="A very cool measurement"
params["scannerOperator"]="Tobi"
params["dfStrength"]=[5e-3]
params["acqNumAverages"]=10000

x = linspace(0,1,3)
params["acqFFValues"] = Float64[]
#params["acqFFValues"] = repeat( cat(1,x,reverse(x[2:end-1])),inner=5)
params["acqNumPeriodsPerFrame"]=1 #length(params["acqFFValues"])

currents = [10.0, 10.0]

data = measurementSystemMatrix(su, daq, robot, safety, positions,
                    currents, params,
                    controlPhase=true, waitTime = 0.5, voltToCurrent = 0.08547008547008547)

h5write(filename, "/matrix", data.signals)


text:
using Transpiler, GeometryTypes, StaticArrays

function vert_image(vertex, args)
    p, uv = vertex
    proj = args[1]
    position = proj * Vec{4, eltype(p)}(p[1], p[2], 0f0, 1f0)
    position, (uv,)
end
function frag_image(fragment_in, args)
    uv, = fragment_in
    image = args[2]
    color = image[uv]
    (color,)
end

vertex = Vec2f0(0), Vec2f0(0)
args = (eye(Mat4f0),)
vert_image(vertex, args)
src, typ = Transpiler.emit_fragment_shader(vert_image, (typeof(vertex), typeof(args)))
write(STDOUT, src)

using Transpiler: GLMethod
using Sugar
x = GLMethod((vert_image, Tuple{typeof(vertex), typeof(args)}))
rasterizer(
    window,
    vbo, args,
    vert_linesegments, frag_linesegments;
    geometryshader = geom_linesegments,
    primitive_in = :lines
)


text:
module ProgressMeter

using Printf: @printf, @sprintf

"""
Holds the five characters that will be used to generate the progress bar.
"""
mutable struct BarGlyphs
    leftend::Char
    fill::Char
    front::Union{Vector{Char}, Char}
    empty::Char
    rightend::Char
end

BarGlyphs() = BarGlyphs(
    '|',
    '█',
    Sys.iswindows() ? '█' : ['▏', '▎', '▍', '▌', '▋', '▊', '▉'],
    ' ',
    '|',
)

"""
    printprogress(io::IO, barglyphs::BarGlyphs, tfirst::Float64, desc, progress, eta_seconds::Real)

Print progress bar to `io`.

# Arguments
- `io::IO`
- `barglyphs::BarGlyphs`
- `tfirst::Float64`
- `desc`: description to be printed at left side of progress bar.
- `progress`: a number between 0 and 1 or `nothing`.
- `eta_seconds::Real`: ETA in seconds
"""
function printprogress(
    io::IO,
    barglyphs::BarGlyphs,
    tfirst::Float64,
    desc,
    progress,
    eta_seconds::Real,
)
    t = time()
    percentage_complete = 100.0 * (isnothing(progress) || isnan(progress) ? 0.0 : progress)

    #...length of percentage and ETA string with days is 29 characters
    barlen = max(0, displaysize(io)[2] - (length(desc) + 29))

    if !isnothing(progress) && progress >= 1
        bar = barstring(barlen, percentage_complete, barglyphs=barglyphs)
        dur = durationstring(t - tfirst)
        @printf io "%s%3u%%%s Time: %s" desc round(Int, percentage_complete) bar dur
        return
    end

    bar = barstring(barlen, percentage_complete, barglyphs=barglyphs)
    if 0 <= eta_seconds <= typemax(Int)
        eta_sec = round(Int, eta_seconds)
        eta = durationstring(eta_sec)
    else
        eta = "N/A"
    end
    @printf io "%s%3u%%%s  ETA: %s" desc round(Int, percentage_complete) bar eta
    return
end

function compute_front(barglyphs::BarGlyphs, frac_solid::AbstractFloat)
    barglyphs.front isa Char && return barglyphs.front
    idx = round(Int, frac_solid * (length(barglyphs.front) + 1))
    return idx > length(barglyphs.front) ? barglyphs.fill :
           idx == 0 ? barglyphs.empty :
           barglyphs.front[idx]
end

function barstring(barlen, percentage_complete; barglyphs)
    bar = ""
    if barlen>0
        if percentage_complete == 100 # if we're done, don't use the "front" character
            bar = string(barglyphs.leftend, repeat(string(barglyphs.fill), barlen), barglyphs.rightend)
        else
            n_bars = barlen * percentage_complete / 100
            nsolid = trunc(Int, n_bars)
            frac_solid = n_bars - nsolid
            nempty = barlen - nsolid - 1
            bar = string(barglyphs.leftend,
                         repeat(string(barglyphs.fill), max(0,nsolid)),
                         compute_front(barglyphs, frac_solid),
                         repeat(string(barglyphs.empty), max(0, nempty)),
                         barglyphs.rightend)
        end
    end
    bar
end

function durationstring(nsec)
    days = div(nsec, 60*60*24)
    r = nsec - 60*60*24*days
    hours = div(r,60*60)
    r = r - 60*60*hours
    minutes = div(r, 60)
    seconds = floor(r - 60*minutes)

    hhmmss = @sprintf "%u:%02u:%02u" hours minutes seconds
    if days>9
        return @sprintf "%.2f days" nsec/(60*60*24)
    elseif days>0
        return @sprintf "%u days, %s" days hhmmss
    end
    hhmmss
end

# issue #31: isnothing require Julia 1.1
# copy-over from
# https://github.com/JuliaLang/julia/blob/0413ef0e4de83b41b637ba02cc63314da45fe56b/base/some.jl
if !isdefined(Base, :isnothing)
    isnothing(::Any) = false
    isnothing(::Nothing) = true
end

end


text:
import LinearAlgebra: BlasInt
const LA = LinearAlgebra

@enum Factorization begin
    QR_FACT
    LU_FACT
end
struct MatrixWorkspace{T} <: AbstractMatrix{Complex{T}}
    A::Matrix{Complex{T}} # Matrix
    d::Vector{T} # Inverse of scaling factors
    # factorizations
    fact::Base.RefValue{Factorization}
    factorized::Base.RefValue{Bool}
    lu::Union{Nothing,LA.LU{Complex{T},Matrix{Complex{T}}}} # LU Factorization of D * J
    qr::LA.QRPivoted{Complex{T},Matrix{Complex{T}}} # Pivoted QR Factorization of D * J
    # qr worksapce
    qr_perm::Vector{Int} # work copy of the permutation in the qr
    qr_work::Vector{Complex{T}} # preallocated memory for pivoted qr
    qr_rwork::Vector{T} # preallocated memory for pivoted qr
    qr_ormqr_work::Vector{Complex{T}}
    # overdetermined
    qr_rhs::Vector{Complex{T}}
    # iterative refinement
    ir_r::Vector{Complex{T}}
    ir_δx::Vector{Complex{T}}
    residual_abs::Vector{T}
    # cond estimators
    lu_cond_work::Vector{Complex{T}}
    lu_cond_rwork::Vector{T}
    qr_cond_work::Vector{ComplexF64}
    qr_cond_rwork::Vector{T}
    # inf_norm_est
    inf_norm_est_work::Vector{Complex{T}}
    inf_norm_est_rwork::Vector{T}
    row_scaling::Vector{T}
end

MatrixWorkspace(Â::AbstractMatrix) =
    MatrixWorkspace(Â, promote_type(real(eltype(Â)), Float64))

function MatrixWorkspace(Â::AbstractMatrix, ::Type{T}) where {T}
    A = Matrix{Complex{T}}(Â)
    m, n = size(A)
    d = ones(m)
    fact = m == n ? LU_FACT : QR_FACT
    # lu
    lu = m == n ? LA.lu(A; check = false) : nothing
    # qr
    qr = LA.qr(A, Val(true))
    qr_perm = zeros(Int, length(qr.p))
    qr_work = Vector{Complex{T}}(undef, 1)
    qr_ormqr_work = Vector{Complex{T}}(undef, 1)
    qr_rwork = Vector{T}(undef, 2n)

    ir_r = Vector{Complex{T}}(undef, m)
    ir_δx = Vector{Complex{T}}(undef, n)
    residual_abs = Vector{T}(undef, m)
    qr_rhs = copy(ir_r)
    lu_cond_work = Vector{Complex{T}}(undef, 2n)
    lu_cond_rwork = Vector{T}(undef, 2n)
    qr_cond_work = Vector{Complex{T}}(undef, 2n)
    qr_cond_rwork = Vector{T}(undef, n)
    inf_norm_est_work = Vector{Complex{T}}(undef, n)
    inf_norm_est_rwork = Vector{T}(undef, n)
    row_scaling = ones(m)
    MatrixWorkspace(
        A,
        d,
        Ref(fact),
        Ref(true),
        lu,
        qr,
        qr_perm,
        qr_work,
        qr_rwork,
        qr_ormqr_work,
        qr_rhs,
        ir_r,
        ir_δx,
        residual_abs,
        lu_cond_work,
        lu_cond_rwork,
        qr_cond_work,
        qr_cond_rwork,
        inf_norm_est_work,
        inf_norm_est_rwork,
        row_scaling,
    )
end

Base.size(MW::MatrixWorkspace) = size(MW.A)

import Base: @propagate_inbounds
@propagate_inbounds Base.getindex(MW::MatrixWorkspace, i::Integer) = getindex(MW.A, i)
@propagate_inbounds Base.getindex(MW::MatrixWorkspace, i::Integer, j::Integer) =
    getindex(MW.A, i, j)
@propagate_inbounds Base.setindex!(MW::MatrixWorkspace, x, i::Integer) =
    setindex!(MW.A, x, i)
@propagate_inbounds Base.setindex!(MW::MatrixWorkspace, x, i::Integer, j::Integer) =
    setindex!(MW.A, x, i, j)
@propagate_inbounds Base.copyto!(MW::MatrixWorkspace, A::AbstractArray) = copyto!(MW.A, A)

"""
    updated!(MW::MatrixWorkspace)

Indicate that the matrix `MW` got updated.
"""
updated!(MW::MatrixWorkspace) = (MW.factorized[] = false; MW)

"""
    update!(MW::MatrixWorkspace, A::Matrix)

Update the matrix in `MW` with `A`.
"""
function update!(MW::MatrixWorkspace, A::Matrix)
    @boundscheck (size(MW.A) == size(A) || throw(ArgumentError("Matrix of invalid size.")))
    @inbounds copyto!(MW.A, A)
    updated!(MW)
end

"""
    factorization!(MW::MatrixWorkspace, fact::Factorization)

Update the used factorization in `MW` to fact.
"""
function factorization!(MW::MatrixWorkspace, fact::Factorization)
    if fact == LU_FACT && MW.lu === nothing
        throw(ArgumentError("Cannot set `LU_FACT` for non-square matrix."))
    end
    MW.factorized[] = MW.fact[] == fact
    MW.fact[] = fact
    MW
end

##################
# Factorizations #
##################

###########################
# CUSTOM LU Factorization #
###########################
# This is an adoption of LA.generic_lufact!
# with 3 changes:
# 1) For choosing the pivot we use abs2 instead of abs
# 2) Instead of using the robust complex division we
#    just use the naive division. This shouldn't
#    lead to problems since numerical diffuculties only
#    arise for very large or small exponents
# 3) We fold lu! and ldiv! into one routine
#    this has the effect that we do not need to allocate
#    the pivot vector anymore and also avoid the allocations
#    coming from the LU wrapper
function lu!(
    A::AbstractMatrix{T},
    b::Union{AbstractVector,Nothing} = nothing,
    ipiv::Union{Vector{I},Nothing} = nothing,
    ::Val{Pivot} = Val(true),
) where {T,I<:Integer,Pivot}
    m, n = size(A)
    minmn = min(m, n)
    # LU Factorization
    @inbounds begin
        for k = 1:minmn
            # find index max
            kp = k
            if Pivot
                amax = zero(real(T))
                for i = k:m
                    if T <: Complex
                        absi = abs2(A[i, k])
                    else
                        absi = abs(A[i, k])
                    end
                    if absi > amax
                        kp = i
                        amax = absi
                    end
                end
            end
            if !(ipiv isa Nothing)
                ipiv[k] = kp
            end
            if !iszero(A[kp, k])
                if k != kp
                    # Interchange
                    for i = 1:n
                        tmp = A[k, i]
                        A[k, i] = A[kp, i]
                        A[kp, i] = tmp
                    end
                    if !(b isa Nothing)
                        b[k], b[kp] = b[kp], b[k]
                    end
                end
                # Scale first column
                # Akkinv = @fastmath inv(A[k,k])
                Akkinv = inv(A[k, k])
                for i = k+1:m
                    A[i, k] *= Akkinv
                end
            end
            # Update the rest
            for j = k+1:n
                A_kj = A[k, j]
                for i = k+1:m
                    A[i, j] -= A[i, k] * A_kj
                end
            end
        end
    end
    A
end

#############
# Custom QR #
#############
"""
    geqp3!(A::AbstractMatrix{ComplexF64},
        jpvt::AbstractVector{BlasInt},
        tau::AbstractVector{ComplexF64},
        work::Vector{ComplexF64},
        rwork=Vector{Float64}(undef, 2*size(A,2)))

Version of `LAPACK.geqp3!` which also accepts the work buffers.
"""
function geqp3!(
    A::AbstractMatrix{ComplexF64},
    jpvt::AbstractVector{BlasInt},
    tau::AbstractVector{ComplexF64},
    work::Vector{ComplexF64},
    rwork = Vector{Float64}(undef, 2 * size(A, 2)),
)
    m, n = size(A)
    lda = stride(A, 2)
    lda == 0 && return nothing
    jpvt .= BlasInt(0)
    lwork = BlasInt(-1)
    info = Ref{BlasInt}()
    for i = 1:2
        ccall(
            (LA.LAPACK.@blasfunc(zgeqp3_), LA.LAPACK.liblapack),
            Cvoid,
            (
                Ref{BlasInt},
                Ref{BlasInt},
                Ptr{ComplexF64},
                Ref{BlasInt},
                Ptr{BlasInt},
                Ptr{ComplexF64},
                Ptr{ComplexF64},
                Ref{BlasInt},
                Ptr{Float64},
                Ptr{BlasInt},
            ),
            m,
            n,
            A,
            lda,
            jpvt,
            tau,
            work,
            lwork,
            rwork,
            info,
        )
        LA.LAPACK.chklapackerror(info[])
        if i == 1
            lwork = BlasInt(real(work[1]))
            resize!(work, lwork)
        end
    end
    nothing
end

function factorize!(WS::MatrixWorkspace)
    if WS.fact[] == LU_FACT
        @inbounds copyto!(WS.lu.factors, WS.A)
        lu!(WS.lu.factors, nothing, WS.lu.ipiv)
    elseif WS.fact[] == QR_FACT
        @inbounds copyto!(WS.qr.factors, WS.A)
        geqp3!(WS.qr.factors, WS.qr.jpvt, WS.qr.τ, WS.qr_work, WS.qr_rwork)
    end
    WS.factorized[] = true
    WS
end

##########
## ldiv ##
##########
@inline _ipiv!(A::LA.LU, b::AbstractVector) = apply_ipiv!(b, 1:length(A.ipiv), A.ipiv)
@inline _inverse_ipiv!(A::LA.LU, b::StridedVecOrMat) =
    apply_ipiv!(b, length(A.ipiv):-1:1, A.ipiv)
@inline function apply_ipiv!(b::AbstractVector, range::OrdinalRange, ipiv)
    @inbounds for i in range
        if i != ipiv[i]
            _swap_rows!(b, i, ipiv[i])
        end
    end
    b
end
@inline function _swap_rows!(B::StridedVector, i::Integer, j::Integer)
    B[i], B[j] = B[j], B[i]
    B
end

@inline function ldiv_upper!(
    A::AbstractMatrix,
    b::AbstractVector,
    x::AbstractVector = b;
    singular_exception::Bool = false,
)
    n = size(A, 2)
    @inbounds for j = n:-1:1
        singular_exception && iszero(A[j, j]) && throw(LA.SingularException(j))
        # xj = x[j] = (@fastmath A[j,j] \ b[j])
        xj = x[j] = A[j, j] \ b[j]
        for i = 1:(j-1)
            b[i] -= A[i, j] * xj
        end
    end
    b
end
@inline function ldiv_unit_lower!(
    A::AbstractMatrix,
    b::AbstractVector,
    x::AbstractVector = b,
)
    n = size(A, 2)
    @inbounds for j = 1:n
        xj = x[j] = b[j]
        for i = j+1:n
            b[i] -= A[i, j] * xj
        end
    end
    x
end

function lu_ldiv!(x, LU::LA.LU, b::AbstractVector; check::Bool = false)
    x === b || copyto!(x, b)
    _ipiv!(LU, x)
    ldiv_unit_lower!(LU.factors, x)
    ldiv_upper!(LU.factors, x; singular_exception = check)
    x
end

function ldiv_adj_unit_lower!(A::AbstractMatrix, b::AbstractVector, x::AbstractVector = b)
    n = size(A, 1)
    @inbounds for j = n:-1:1
        z = b[j]
        for i = n:-1:j+1
            z -= conj(A[i, j]) * x[i]
        end
        x[j] = z
    end
    x
end

function ldiv_adj_upper!(
    A::AbstractMatrix,
    b::AbstractVector,
    x::AbstractVector = b;
    singular_exception::Bool = false,
)
    n = size(A, 1)
    @inbounds for j = 1:n
        z = b[j]
        for i = 1:j-1
            z -= conj(A[i, j]) * x[i]
        end
        iszero(A[j, j]) && singular_exception && throw(SingularException(j))
        x[j] = @fastmath conj(A[j, j]) \ z
    end
    x
end

function lu_ldiv_adj!(x, LU::LA.LU, b::AbstractVector; check::Bool = false)
    x === b || copyto!(x, b)
    ldiv_adj_upper!(LU.factors, x; singular_exception = check)
    ldiv_adj_unit_lower!(LU.factors, x)
    _inverse_ipiv!(LU, x)
    x
end

"""
    ormqr!(side, trans, A, tau, C, work)

Version of `LAPACK.ormqr!` which also accepts the work buffer.
Computes `Q * C` (`trans = N`), `transpose(Q) * C` (`trans = T`), `adjoint(Q) * C`
(`trans = C`) for `side = L` or the equivalent right-sided multiplication
for `side = R` using `Q` from a `QR` factorization of `A` computed using
`geqrf!`. `C` is overwritten.
"""
function ormqr!(
    side::AbstractChar,
    trans::AbstractChar,
    A::Matrix{ComplexF64},
    tau::Vector{ComplexF64},
    C::Vector{ComplexF64},
    work::Vector{ComplexF64},
)
    m, n = (size(C, 1), 1)
    mA = size(A, 1)
    k = length(tau)
    lwork = BlasInt(-1)
    info = Ref{BlasInt}()
    for i = 1:2  # first call returns lwork as work[1]
        ccall(
            (LA.LAPACK.@blasfunc(zunmqr_), LA.LAPACK.liblapack),
            Cvoid,
            (
                Ref{UInt8},
                Ref{UInt8},
                Ref{BlasInt},
                Ref{BlasInt},
                Ref{BlasInt},
                Ptr{ComplexF64},
                Ref{BlasInt},
                Ptr{ComplexF64},
                Ptr{ComplexF64},
                Ref{BlasInt},
                Ptr{ComplexF64},
                Ref{BlasInt},
                Ptr{BlasInt},
            ),
            side,
            trans,
            m,
            n,
            k,
            A,
            max(1, stride(A, 2)),
            tau,
            C,
            max(1, stride(C, 2)),
            work,
            lwork,
            info,
        )
        LA.LAPACK.chklapackerror(info[])
        if i == 1
            lwork = BlasInt(real(work[1]))
            resize!(work, lwork)
        end
    end
    C
end

# This is an adoption of the Julia implementation of ldiv! for QRPivoted
function qr_ldiv!(
    x,
    A::LA.QRPivoted{ComplexF64},
    b::Vector{ComplexF64},
    perm::Vector{Int},
    ormqr_work::Vector{ComplexF64},
      #corank::Int,
)
    mA, nA = size(A.factors)
    nr = min(mA, nA)
    nrhs = length(b)
    nr == 0 && return x

    # The following is equivalent to LA.lmul!(LA.adjoint(A.Q), b) but
    # we can also pass the preallocated work vector
    ormqr!('L', 'C', A.factors, A.τ, b, ormqr_work)
    ldiv_upper!(A.factors, b; singular_exception = false)
    @inbounds for i = 1:nr
        x[i] = b[i]
        perm[i] = A.p[i]
    end
    Base.invpermute!!(x, perm)
    # The following would cover the case of rank deficient qr.
    # This is currently not used but we keep it here since it will be used later.
    # if corank == nr
    #     fill!(x, zero(eltype(x)))
    #     return x
    # end
    # rank = nr - corank
    # if rank ≠ nA || rank ≠ nr
    #     C, τ = LA.LAPACK.tzrzf!(A.factors[1:rank,:])
    #     LA.ldiv!(LA.UpperTriangular(C[1:rank,1:rank]),view(LA.lmul!(LA.adjoint(A.Q), b), 1:rank))
    #     for i in 1:rank
    #         x[i] = b[i]
    #     end
    #     for i in rank+1:nA
    #         x[i] = zero(eltype(x))
    #     end
    #     LA.LAPACK.ormrz!('L', eltype(x)<:Complex ? 'C' : 'T', C, τ, reshape(x, length(x), 1))
    #     invpermute!(x, A.p)
    # end
    return x
end

function LA.ldiv!(x::AbstractVector, WS::MatrixWorkspace, b::AbstractVector)
    WS.factorized[] || factorize!(WS)
    if WS.fact[] == LU_FACT
        lu = WS.lu
        if lu !== nothing
            lu_ldiv!(x, lu, b)
        end
    elseif WS.fact[] == QR_FACT
        copyto!(WS.qr_rhs, b)
    # if issquare(WS.A)
            # # we applied row scaling for square matrices
            # @inbounds for i in eachindex(Jac.b)
            #     Jac.b[i] = b[i] / Jac.D[i]
            # end
    #
        qr_ldiv!(x, WS.qr, WS.qr_rhs, WS.qr_perm, WS.qr_ormqr_work)
    end
    x
end

##########################
## Iterative refinement ##
##########################
"""
    residual!(u, A, x, b, ::Type{T}=eltype(u))

Compute the residual `Ax-b` in precision `T` and store in `u`.
"""
function residual!(
    u::AbstractVector{S},
    A::AbstractMatrix{S},
    x::AbstractVector{S},
    b::AbstractVector{S},
    ::Type{T} = eltype(u),
) where {S,T}
    @boundscheck size(A, 1) == length(b) && size(A, 2) == length(x)
    m, n = size(A)
    @inbounds for i = 1:m
        dot = A[i, 1] * convert(T, x[1])
        for j = 2:n
            dot = muladd(A[i, j], convert(T, x[j]), dot)
        end
        u[i] = dot - b[i]
    end
    u
end

"""
    iterative_refinement_step!([x,] MW::MatrixWorkspace, x̂, b, norm::AbstractNorm=InfNorm(), T=eltype(x̂))

Apply one step of iterative refinement where the residual is computed with precision `T`.
Stores the result in `x` and returns the norm of the update `δx`. If `x` is not provided
`x̂` is updated inplace.
"""
function iterative_refinement_step!(
    WS::MatrixWorkspace,
    x̂,
    b,
    norm::AbstractNorm = InfNorm(),
    ::Type{T} = eltype(x̂),
) where {T}
    iterative_refinement_step!(x̂, WS, x̂, b, norm, T)
end
function iterative_refinement_step!(
    x,
    WS::MatrixWorkspace,
    x̂,
    b,
    norm::AbstractNorm = InfNorm(),
    ::Type{T} = eltype(x̂),
) where {T}
    residual!(WS.ir_r, WS.A, x̂, b, T)
    δx = LA.ldiv!(WS.ir_δx, WS, WS.ir_r)
    for i in eachindex(x)
        x[i] = x̂[i] - δx[i]
    end
    return norm(δx)
end

#########################
## Condition estimates ##
#########################
"""
    rcond(WS::MatrixWorkspace)

Computes the reciprocal condition number of the matrix `A` in `WS` wrt the ∞-norm.
If `A` is a non-square matrix, the reciprocal condition number of `R` in the
QR decomposition is computed.
"""
function rcond(WS::MatrixWorkspace)
    WS.factorized[] || factorize!(WS)
    if WS.fact[] == LU_FACT
        lu = WS.lu
        if lu === nothing
            NaN
        else
            rcond(lu, LA.opnorm(WS.A, Inf), WS.lu_cond_work, WS.lu_cond_rwork)
        end
    else # WS.fact[] == QR_FACT
        rcond(WS.qr, WS.qr_cond_work, WS.qr_cond_rwork)
    end
end
function rcond(LU::LA.LU, anorm, work, rwork)
    A = LU.factors
    n = size(A, 1)
    rcond = Ref{Float64}()
    info = Ref{Int64}()
    lda = max(1, stride(A, 2))
    ccall(
        (LA.LAPACK.@blasfunc(zgecon_), LinearAlgebra.LAPACK.liblapack),
        Cvoid,
        (
            Ref{UInt8},
            Ref{Int64},
            Ptr{ComplexF64},
            Ref{Int64},
            Ref{Float64},
            Ref{Float64},
            Ptr{ComplexF64},
            Ptr{Float64},
            Ptr{Int64},
        ),
        'I',
        n,
        A,
        lda,
        anorm,
        rcond,
        work,
        rwork,
        info,
    )
    rcond[]
end

function rcond(qr::LA.QRPivoted, work, rwork)
    # use trcon to estimate cond of R
    m, n = size(qr.factors)
    rcond = Ref{Float64}(1)
    info = Ref{BlasInt}()
    ccall(
        (LA.LAPACK.@blasfunc(ztrcon_), LA.LAPACK.liblapack),
        Cvoid,
        (
            Ref{UInt8},
            Ref{UInt8},
            Ref{UInt8},
            Ref{BlasInt},
            Ptr{ComplexF64},
            Ref{BlasInt},
            Ref{Float64},
            Ptr{ComplexF64},
            Ptr{Float64},
            Ptr{BlasInt},
        ),
        'I',
        'U',
        'N',
        n,
        qr.factors,
        m,
        rcond,
        work,
        rwork,
        info,
    )
    if m == n
        # Q has condition number 1 in the 2 norm
        # Some experiments indicate that in the Inf norm this is roughly
        # 0.8√n
        rcond[] / (0.8 * √(n))
    else
        rcond[]
    end
end

"""
    inf_norm_est(WS::MatrixWorkspace, g::Union{Nothing,Vector{<:Real}}=nothing, d::Union{Nothing,Vector{<:Real}}=nothing)
    inf_norm_est(lu::LA.LU, g::Union{Nothing,Vector{<:Real}}=nothing, d::Union{Nothing,Vector{<:Real}}=nothing)
    inf_norm_est(lu::LA.LU, g::Union{Nothing,Vector{<:Real}}, d::Union{Nothing,Vector{<:Real}}, work{<:Complex}, rwork::Vector{<:Real})

Estimation of the infinity norm of `diag(d)⁻¹A⁻¹g` where `g` and `d` are optional positive vectors
and `A=LU`. If `g` is `nothing` the infinity norm of `A⁻¹` is estimated.
This uses the 1-norm lapack condition estimator described by Highahm in [^H88].

[^H88]: Higham, Nicholas J. "FORTRAN codes for estimating the one-norm of a real or complex matrix, with applications to condition estimation." ACM Transactions on Mathematical Software (TOMS) 14.4 (1988): 381-396.
"""
function inf_norm_est(
    lu::LA.LU,
    g::Union{Nothing,Vector{<:Real}} = nothing,
    d::Union{Nothing,Vector{<:Real}} = nothing,
)
    n = size(lu.factors, 2)
    work = Vector{eltype(lu.factors)}(undef, n)
    rwork = Vector{real(eltype(lu.factors))}(undef, n)
    inf_norm_est(lu, g, d, work, rwork)
end
function inf_norm_est(
    WS::MatrixWorkspace,
    g::Union{Nothing,Vector{<:Real}} = nothing,
    d::Union{Nothing,Vector{<:Real}} = nothing,
)
    WS.fact[] == LU_FACT || factorization!(WS, LU_FACT)
    WS.factorized[] || factorize!(WS)
    inf_norm_est(WS.lu, g, d, WS.inf_norm_est_work, WS.inf_norm_est_rwork)
end
function inf_norm_est(
    lu::LA.LU,
    g::Union{Nothing,Vector{<:Real}},
    d::Union{Nothing,Vector{<:Real}},
    work::Vector{<:Complex},
    rwork::Vector{<:Real},
)
    z = ξ = y = work
    x = rwork

    n = size(lu.factors, 1)
    x .= inv(n)
    if d !== nothing
        x ./= d
    end
    lu_ldiv_adj!(y, lu, x)
    if g !== nothing
        y .*= g
    end
    γ = sum(abs, y)
    y ./= abs.(y)
    if g !== nothing
        y .*= g
    end
    lu_ldiv!(z, lu, ξ)
    if d !== nothing
        x .= real.(z) ./ d
    else
        x .= real.(z)
    end
    k = 2
    while true
        j, max_xᵢ = 1, abs(x[1])
        for i = 2:n
            abs_xᵢ = abs(x[i])
            if abs_xᵢ > max_xᵢ
                j, max_xᵢ = i, abs_xᵢ
            end
        end
        x .= zero(eltype(x))
        x[j] = one(eltype(x))
        if d !== nothing
            x ./= d
        end
        lu_ldiv_adj!(y, lu, x)
        if g !== nothing
            y .*= g
        end
        γ̄ = γ
        γ = sum(abs, y)
        γ ≤ γ̄ && break
        ξ .= y ./ abs.(y)
        if g !== nothing
            ξ .*= g
        end
        lu_ldiv!(z, lu, ξ)
        if d !== nothing
            x .= real.(z) ./ d
        else
            x .= real.(z)
        end
        k += 1
        if x[j] == LA.norm(x, Inf) || k > 5
            break
        end
    end
    γ
end

"""
    row_scaling!(WS::MatrixWorkspace,
        norm::Union{InfNorm, WeightedNorm{<:InfNorm}},
        residual::Union{Nothing,AbstractVector{<:Real}} = nothing)

Compute a suitable scaling `D` of the rows of `WS.A` and store it in `WS.row_scaling`.
The returned scaling has the property that the from `norm` induced operator for `D⁻¹ * WS.A`
is 1.
The vector `residual` is a an optional vector containing the absolute values of the limiting
residual of Newton's method. The scaling is performaned in such a way that `norm(residual)`
is *not* increase. Unless the from `norm` induced operator norm is less than 1.
Then the residual is increased by the minimal amount necessary.
"""
function row_scaling!(
    WS::MatrixWorkspace{T},
    norm::Union{InfNorm,WeightedNorm{<:InfNorm}},
    r::Union{Nothing,AbstractVector{<:Real}} = nothing,
) where {T}
    row_scaling!(WS.row_scaling, WS, norm, r)
end
function row_scaling!(
    d::AbstractVector,
    WS::MatrixWorkspace{T},
    norm::Union{InfNorm,WeightedNorm{<:InfNorm}},
    r::Union{Nothing,AbstractVector{<:Real}} = nothing,
) where {T}
    d .= 1e-3
    m = length(d)
    if isa(norm, WeightedNorm)
        for j = 1:size(WS.A, 2), i = 1:m
            d[i] += abs(WS.A[i, j]) * norm[j]
        end
    else
        for j = 1:size(WS.A, 2), i = 1:m
            d[i] += abs(WS.A[i, j])
        end
    end

    if size(WS, 1) == 1
        d[1] = max(1, d[1])
        return d
    end

    if r !== nothing
        d̂ = min(maximum(d), 1.0)
        # compute the norm of r under the scaling of all rows whose norm is larger
        # than d̂, i.e. all rows that either decrease the value of |rᵢ| or
        # or increase it by the minimal abount to still have opnorm 1.
        norm_scaled_r = r[1] / max(d[1], d̂)
        for i = 2:m
            norm_scaled_r = max(norm_scaled_r, r[i] / max(d[i], d̂))
        end
        for i = 1:length(d)
            r̂ᵢ = r[i] / norm_scaled_r
            if isfinite(r̂ᵢ)
                d[i] = max(d[i], r̂ᵢ)
            end
        end
    else
        # Don't scale zero rows -> set dᵢ = 1 if too small
        d_max_eps = sqrt(eps())
        for i = 1:m
            if d[i] < d_max_eps
                d[i] = one(T)
            end
        end
    end
    d
end

"""
    scaled_cond!(WS::MatrixWorkspace, norm, residual)

Perform first a row scaling as described by [`row_scaling!`](@ref) and then
compute the condition number ``κ`` of the row-equilibrated matrix `WS.A`.
This is the best possible condition number under row scaling (e.g. [Higham02, Thm 7.5]).

[Higham02]: Higham, Nicholas J. Accuracy and stability of numerical algorithms. Vol. 80. Siam, 200
"""
function scaled_cond!(
    WS::MatrixWorkspace,
    norm::Union{<:WeightedNorm{<:InfNorm},InfNorm},
    r::Union{Nothing,AbstractVector{<:Real}} = nothing,
)
    d = row_scaling!(WS, norm, r)
    # row scaling makes ||A||_∞ = 1 so we don't this here
    if (norm isa WeightedNorm)
        cond_est = inf_norm_est(WS, d, weights(norm))
    else
        cond_est = inf_norm_est(WS, d)
    end
    cond_est
end

#####################
## JacobianMonitor ##
#####################
struct JacobianMonitor{T}
    J::MatrixWorkspace{T}
    cond::Base.RefValue{Float64}
    forward_err::Base.RefValue{Float64}
    # stats
    factorizations::Base.RefValue{Int}
    ldivs::Base.RefValue{Int}
end
function JacobianMonitor(A::AbstractMatrix)
    J = MatrixWorkspace(A)
    cond = Ref(1.0)
    forward_err = Ref(0.0)
    JacobianMonitor(J, cond, forward_err, Ref(0), Ref(0))
end

updated!(JM::JacobianMonitor) = updated!(JM.J)
jacobian(JM::JacobianMonitor) = JM.J

function Base.show(io::IO, JM::JacobianMonitor{T}) where {T}
    println(io, "JacobianMonitor{$T}:")
    println(io, " • cond → ", round(JM.cond[], sigdigits = 5))
    println(io, " • forward_err → ", round(JM.forward_err[], sigdigits = 5))
    println(io, " • # factorizations → ", JM.factorizations[])
    println(io, " • # ldivs → ", JM.ldivs[])
end

"""
    init!(JM::JacobianMonitor)

(Re-)initialize the `JacobianMonitor`.
"""
function init!(JM::JacobianMonitor; keep_stats::Bool = false)
    JM.cond[] = 1.0
    JM.forward_err[] = 0.0
    JM.factorizations[] = 0.0
    JM.ldivs[] = 0.0
    JM
end

"""
    forward_err(JM::JacobianMonitor)

Return the last estimate of the forward error as computed by [`forward_err`](@ref) or [`strong_forward_err`](@ref).
"""
forward_err(JM::JacobianMonitor) = JM.forward_err[]

"""
    forward_err!(x̂::AbstractVector, JM::JacobianMonitor, b::AbstractVector, norm::AbstractNorm)

Compute an estimate of the forward error ||x-x̂||/||x|| and update the value in `JM`.
"""
function forward_err!(
    JM::JacobianMonitor,
    x̂::AbstractVector,
    b::AbstractVector,
    norm::AbstractNorm,
    T = eltype(x̂),
)
    forward_err!(x̂, JM, x̂, b, norm, T)
end
function forward_err!(
    x::AbstractVector,
    JM::JacobianMonitor,
    x̂::AbstractVector,
    b::AbstractVector,
    norm::AbstractNorm,
    T = eltype(x̂),
)
    norm_x̂ = norm(x̂)
    JM.forward_err[] = iterative_refinement_step!(x, JM.J, x̂, b, norm, T) / norm_x̂
end

"""
    strong_forward_err!(x̂::AbstractVector, JM::JacobianMonitor, b::AbstractVector, norm::AbstractNorm)

Compute a more robust estimate of the forward error ||x-x̂||/||x|| using eq. (2.14) in [Demmel, Section 2.4.4]
and update the value in `JM`.
If `jacobian(JM).fact[] == QR_FACT` then this falls back to [`forward_err!`](@ref).

[Demmel]: Demmel, James W. Applied numerical linear algebra. Vol. 56. Siam, 1997.
"""
function strong_forward_err!(
    JM::JacobianMonitor,
    x̂::AbstractVector,
    b::AbstractVector,
    norm::InfNorm,
    T = eltype(x̂),
)
    strong_forward_err!(x̂, JM, x̂, b, norm, T)
end
function strong_forward_err!(
    x::AbstractVector,
    JM::JacobianMonitor,
    x̂::AbstractVector,
    b::AbstractVector,
    norm::Union{InfNorm,WeightedNorm{InfNorm}},
    T = eltype(x̂),
)
    if jacobian(JM).fact[] == QR_FACT
        return JM.forward_err[] = forward_err!(x, JM, x̂, b, norm)
    end
    norm_x̂ = norm(x̂)
    WS = jacobian(JM)
    residual!(WS.ir_r, WS.A, x̂, b, T)
    WS.residual_abs .= abs.(b)
    for j = 1:size(WS.A, 2)
        ax̂ⱼ = abs(x̂[j])
        for i = 1:size(WS.A, 1)
            WS.residual_abs[i] += abs(WS.A[i, j]) * ax̂ⱼ
        end
    end
    WS.residual_abs .*= eps()
    WS.residual_abs .+= abs.(WS.ir_r)

    if norm isa WeightedNorm
        JM.forward_err[] = inf_norm_est(WS, WS.residual_abs, weights(norm)) / norm_x̂
    else
        JM.forward_err[] = inf_norm_est(WS, WS.residual_abs) / norm_x̂
    end
end

"""
    cond!(JM::JacobianMonitor, norm, residual)

Compute the condition number ``κ`` of the Jacobian ``J`` with respect to the
infinity norm and store it. If ``J`` is a square matrix the condition number of the
row-equilibrated Jacobian is computed. See [`scaled_cond!`](@ref) for details.
"""
function cond!(
    JM::JacobianMonitor,
    norm::AbstractNorm = InfNorm(),
    r::Union{Nothing,AbstractVector{<:Real}} = nothing,
)
    if jacobian(JM).fact[] == LU_FACT
        JM.cond[] = scaled_cond!(jacobian(JM), norm, r)
    else
        # TODO: The weighted norm is NOT considered currently.
        JM.cond[] = inv(rcond(jacobian(JM)))
    end
end

"""
    apply_row_scaling!(r, JM::JacobianMonitor)

Apply the computed row scaling used in [`cond!`](@ref) resp. [`scaled_cond!`](@ref) to `r`.
"""
apply_row_scaling!(r::AbstractVector, JM::JacobianMonitor) = (r ./= JM.J.row_scaling; r)

"""
    cond(JM::JacobianMonitor)

Returns the with [`cond!`](@ref) computed condition number.
"""
LA.cond(JM::JacobianMonitor) = JM.cond[]

"""
    enum JacobianMonitorUpdates

## Cases

* `JAC_MONITOR_UPDATE_NOTHING`: Do nothing
* `JAC_MONITOR_UPDATE_FERR`: Update the forward error estimate
* `JAC_MONITOR_UPDATE_COND`: Update the condition number estimate
* `JAC_MONITOR_UPDATE_ALL`: Update the forward error and condition number estimate
"""
@enum JacobianMonitorUpdates begin
    JAC_MONITOR_UPDATE_NOTHING
    JAC_MONITOR_UPDATE_FERR
    JAC_MONITOR_UPDATE_COND
    JAC_MONITOR_UPDATE_ALL
end

"""
    ldiv!(x̂::AbstractVector,
                  JM::JacobianMonitor,
                  b::AbstractVector,
                  norm::AbstractNorm,
                  update::JacobianMonitorUpdates=JAC_MONITOR_UPDATE_NOTHING)

solve the linear system `jacobian(JM)x̂=b`. `update` controls the computation of additional
informations, see [`JacobianMonitorUpdates`](@ref).
"""
function LA.ldiv!(
    x̂::AbstractVector,
    JM::JacobianMonitor,
    b::AbstractVector,
    norm::AbstractNorm = InfNorm(),
    update::JacobianMonitorUpdates = JAC_MONITOR_UPDATE_NOTHING,
)
    # stats update
    JM.factorizations[] += !jacobian(JM).factorized[]
    JM.ldivs[] += 1

    LA.ldiv!(x̂, jacobian(JM), b)
    if update == JAC_MONITOR_UPDATE_FERR || update == JAC_MONITOR_UPDATE_ALL
        forward_err!(JM, x̂, b, norm)
    end
    if update == JAC_MONITOR_UPDATE_COND || update == JAC_MONITOR_UPDATE_ALL
        cond!(JM)
    end
    x̂
end

#########################
## Hermite Normal Form ##
#########################

"""
    hnf(A, T=elytpe(A))

Compute the hermite normal form `H` of `A` by overwriting `A` and the corresponding transformation
matrix `U` using precision T. This is `A*U == H` and `H` is a lower triangular matrix.

The implementation follows the algorithm described in [1].

[1] Kannan, Ravindran, and Achim Bachem. "Polynomial algorithms for computing the Smith and Hermite normal forms of an integer matrix." SIAM Journal on Computing 8.4 (1979): 499-507.
"""
function hnf(A, T = eltype(A))
    H = similar(A, T)
    H .= A
    U = similar(A, T)
    hnf!(H, U)
    H, U
end

# use checked arithmethic
⊡(x, y) = Base.checked_mul(x, y)
⊞(x, y) = Base.checked_add(x, y)

"""
    hnf!(H, U, A)

Inplace version of [hnf](@ref).

    hnf!(A, U)

Inplace version of [hnf](@ref) overwriting `A` with `H`.
"""
hnf!(H, U, A) = hnf!(copyto!(H, A), U)
function hnf!(A, U)
    n = size(A, 1)
    U .= 0
    @inbounds for i = 1:n
        U[i, i] = one(eltype(U))
    end
    @inbounds for i = 1:(n-1)
        ii = i ⊞ 1
        for j = 1:i
            if !iszero(A[j, j]) || !iszero(A[j, ii])
                # 4.1
                r, p, q = gcdx(A[j, j], A[j, ii])
                # 4.2
                d_j = -A[j, ii] ÷ r
                d_ii = A[j, j] ÷ r
                for k = 1:n
                    a_kj, a_kii = A[k, j], A[k, ii]
                    A[k, j] = a_kj ⊡ p ⊞ a_kii ⊡ q
                    A[k, ii] = a_kj ⊡ d_j ⊞ a_kii ⊡ d_ii

                    u_kj, u_kii = U[k, j], U[k, ii]
                    U[k, j] = u_kj ⊡ p ⊞ u_kii ⊡ q
                    U[k, ii] = u_kj ⊡ d_j ⊞ u_kii ⊡ d_ii
                end
            end
            # 4.3
            if j > 1
                reduce_off_diagonal!(A, U, j)
            end
        end
        #5
        reduce_off_diagonal!(A, U, ii)
    end
    # Special case for 1 × 1 matrices to guarantee that the diagonal is positive
    # This comes up in polyhedral homotopy for univariate polynomials
    if n == 1
        U[1, 1] = flipsign(U[1, 1], A[1, 1])
        A[1, 1] = flipsign(A[1, 1], A[1, 1])
    end

    nothing
end

@inline function reduce_off_diagonal!(A, U, k)
    n = size(A, 1)
    @inbounds if A[k, k] < 0
        for i = 1:n
            A[i, k] = -A[i, k]
            U[i, k] = -U[i, k]
        end
    end
    @inbounds for z = 1:(k-1)
        if !iszero(A[z, z])
            r = -ceil(eltype(A), A[k, z] / A[z, z])
            for i = 1:n
                U[i, z] = U[i, z] ⊞ r ⊡ U[i, k]
                A[i, z] = A[i, z] ⊞ r ⊡ A[i, k]
            end
        end
    end
    nothing
end


text:
using Distributions
using Plots
using ProgressMeter
using DifferentialEquations
using Statistics
using DelimitedFiles

function ODEmodel(du,u,params,t)

    (xC,xD) = u
    (alphaC,alphaD,sigma,beta,kappa,n,dC,dD,eta,Nt) = params

    # Compute expected public good benefit
    y = xC/(xC+xD)

    pC = beta/n*(1.0 + (n-1)*y)
    pD = beta/n*(n-1)*y

    (rC,rD) = nonlinearPublicGoodForm(pC,pD,params)

    du[1] = rC*xC*(1.0 - (xC+xD)/eta) - dC*xC
    du[2] = rD*xD*(1.0 - (xC+xD)/eta) - dD*xD

end

function cumsum(x)
    y = Array{Float64}(undef,length(x))
    for n = 1 : length(x)
        y[n] = sum(x[1:n])
    end

    y = y./y[end]

    return y
end

function computePublicGoodBenefit(c,d,n,beta)

    if c-1 >= 0 && d >= 0 && n > 1
        ############ A cooperator sees n-1 other cells: ############
        dist = Hypergeometric(c-1,d,n-1)
        cdfdist = cdf.(Ref(dist),support(dist))

        # Find the number of cooperators in the group
        roll = rand()
        cInGroup = findall(x->(x>roll),cdfdist)[1]-1

        if d < n-1         # Less in population than group size
            cInGroup += n-d-1
        end

        # Public good benefit
        pC = beta/n*(1.0 + cInGroup)

        # println(cInGroup)

    else
        pC = beta/n
    end
    if c >=0 && d-1 >= 0 && n > 1
        ############ A defector sees n-1 other cells: ############
        # Probability of picking individuals without replacement
        dist = Hypergeometric(c,d-1,n-1)
        cdfdist = cdf.(Ref(dist),support(dist))

        # Find the number of cooperators in the group
        roll = rand()
        cInGroup = findall(x->(x>roll),cdfdist)[1]-1

        if d < n-1         # Less in population than group size
            cInGroup += n-d-1
        end

        # Public good benefit
        pD = beta/n*cInGroup
    else
        pD = 0.0
    end

    # println(cInGroup)

    return (pC,pD)

end

function nonlinearPublicGoodForm(pC,pD,params)

    (alphaC,alphaD,sigma,beta,kappa,n,dC,dD,eta) = params

    # Compute rC and rD
    rC = alphaC*(1+exp(sigma))/(1+exp(sigma-pC)) - kappa
    rD = alphaD*(1+exp(sigma))/(1+exp(sigma-pD))

    return (rC,rD)

end

function transitionRates(x,params)

    # Grab populations and parameters
    (c, d) = x
    (alphaC,alphaD,sigma,beta,kappa,n,dC,dD,eta,Nt) = params

    count = 0
    # pCvec = Vector{Float64}
    # pDvec = Vector{Float64}
    # t_10  = Vector{Float64}

    # Initialize the rates array
    rates = Array{Float64}(undef,4)

    cVec = Array{Float64}(undef,Nt); cVec[1] = c
    dVec = Array{Float64}(undef,Nt); dVec[1] = d
    pCvec = Array{Float64}(undef,Nt-1);
    pDvec = Array{Float64}(undef,Nt-1);
    t = Array{Float64}(undef,Nt); t[1] = 0.0

    for i = 2 : Nt

        # Compute public good benefit
        (pC,pD) = computePublicGoodBenefit(c,d,n,beta)
        # Compute rC and rD
        (rC,rD) = nonlinearPublicGoodForm(pC,pD,params)

        # print(rC," ",rD,"\n")

        # Rate function (includes all the possible transitions)
        rates[1] = max(0,rC*c*(1.0 - (c + d)/eta))      # Birth of a cooperator
        rates[2] = max(0,rD*d*(1.0 - (c + d)/eta))      # Birth of a defector
        rates[3] = dC*c                                 # Death of a cooperator
        rates[4] = dD*d                                 # Death of a defector

        # Calculate the ratetotal this will be used to compute time till next event
        # and which event occurs
        ratetotal = sum(rates)

        # println(rates)

        # Calculate time till next event
        roll = rand()
        t[i] = t[i-1] - log(roll)/ratetotal

        # Cumulative sum of the rates
        cumrates = cumsum(rates)

        # Now we roll to see which event will occur
        roll = rand()
        # print(roll,"\n")
        index = findall(x->(x>roll),cumrates)[1]

        decreaseC = [3;5;6]     # Events which lead to decrease in cooperators
        decreaseD = [4;7;8]     # Events which lead to decrease in defectors

        # Change the population size
        if any(decreaseC .== index)        # C -> C - 1
            c -= 1
            # println("c dies\n")
        elseif any(decreaseD .== index)    # D -> D - 1
            d -= 1
            # println("d dies\n")
        elseif index == 1                           # C -> C + 1
            c += 1
            # println("c born\n")
        else
            d += 1
            # println("d born\n")
        end

        cVec[i] = c
        dVec[i] = d

        pCvec[i-1] = pC
        pDvec[i-1] = pD

        if c*d == 0
            t    = t[1:i]
            cVec = cVec[1:i]
            dVec = dVec[1:i]
            pCvec = pCvec[1:i]
            pDvec = pDvec[1:i]
            break
        end

        if i % 10000 == 0
            # println(pC," ", pD,"\n")
        end

    end

return (t,cVec,dVec,pCvec,pDvec)

end

# (C0,D0)
C0 = 50; D0 = 50;
x = (C0,D0)

betaval = 5.0; n = 15; sigma = 2.0;
Nt = 10^5
plotsteps = Int(max(Nt/1000,100))

# # (alphaC,alphaD,sigma,beta,kappa,n,dC,dD,eta,Nt)
# params = (1.0,1.0,2.0,betaval,0.5,n,0.1,0.1,1000.0,Nt)

Nruns = 200

savematrix = Array{Any}(undef,Nruns)

params = (1.0,1.0,sigma,betaval,0.5,n,0.1,0.1,1000.0,Nt)

@showprogress 1 "Computing..." for i = 1 : Nruns

    (t,cooperators,defectors,pCvec,pDvec) = transitionRates(x,params)

    t = t[1:plotsteps:end];
    cooperators = cooperators[1:plotsteps:end]
    defectors = defectors[1:plotsteps:end]
    savematrix[i] = [t cooperators defectors]

end

# writedlm(outfile,savematrix)

u0 = [float(C0),float(D0)]
tspan = (0.0,1000.0)

prob = ODEProblem(ODEmodel,u0,tspan,params)
sol = solve(prob,Tsit5())
plotting = false

Nt_ODE = length(sol.t)
xC      = Array{Float64}(undef,Nt_ODE)
xD      = Array{Float64}(undef,Nt_ODE)
y       = Array{Float64}(undef,Nt_ODE)
pC_ODE  = Array{Float64}(undef,Nt_ODE)
pD_ODE  = Array{Float64}(undef,Nt_ODE)
for i = 1 : Nt_ODE
    xC[i] = sol.u[i][1]
    xD[i] = sol.u[i][2]
    y[i]  = xC[i]/(xC[i]+xD[i])
    pC_ODE[i] = betaval/n*(1.0+y[i]*(n-1))
    pD_ODE[i] = betaval/n*y[i]*(n-1)
end

# writedlm(outfile,[sol.t xC])

max_stars_repo_name: MathOnco/nonlinearPGG

max_issues_repo_name: MathOnco/nonlinearPGG

max_forks_repo_name: MathOnco/nonlinearPGG


text:
export get_info, print_info, get_info_deep

# Distance evaluation & summarisation
# -----------------------------------

Base.show(io::IO, d::CompleteMatchingDistance) = print(io, typeof(d))

function (d::T where {T<:CompleteMatchingDistance})(
    X::Vector{S}, Y::Vector{S}
    ) where {S}
    C = get_cost_matrix_dynamic(d, X, Y)
    if size(C) == (0,0)
        println("Weird hungarian output")
        @show X, Y
    else
        out = hungarian(C)
        return out[2]
    end
end

function Base.show(io::IO, d_gen::General{T}) where {T<:CompleteMatchingDistance}
    b = IOBuffer()
    show(b, d_gen.d)
    d_str = String(take!(b))
    print(io, "General{$(d_str)}")
end

function (d::General{T} where {T<:CompleteMatchingDistance})(
    X::Vector{S}, Y::Vector{S}
    ) where {S}
    C = get_cost_matrix_fixed(d, X, Y)
    return hungarian(C)[2]
end

const MatchingBasedDistance = Union{CompleteMatchingDistance, General{T}} where {T<:CompleteMatchingDistance}

function get_info(
    d::T,
    X::Vector{S}, Y::Vector{S}
    ) where {T<:MatchingBasedDistance, S}

    C = get_cost_matrix_fixed(d, X, Y)
    assignment, cost = hungarian(C)
    indx, indy = (Int[], Int[])
    for i in eachindex(X)
        j = assignment[i]
        if j ≤ length(Y)
            push!(indx, i)
            push!(indy, j)
        end
    end
    return indx, indy
end

function print_info(
    d::T,
    X::Vector{S}, Y::Vector{S}
    ) where {T<:MatchingBasedDistance, S}

    title = "Matching Summary"
    println(title)
    println("-"^length(title) * "\n")
    indx, indy = get_info(d, X, Y)
    max_len = maximum(map(z->length(@sprintf("%s",z)), X[indx]))
    for (i,j) in zip(indx,indy)
        xtmp, ytmp = (X[i],Y[j])
        pad = max_len - length(@sprintf("%s", xtmp))
        println(" "^pad * "$xtmp", " → $ytmp")
    end
    println("\nUnmatched entries (first):")
    for i in eachindex(X)
        if i ∉ indx
            print("$(X[i]), ")
        end
    end
    println("\nUnmatched entries (second):")
    for i in eachindex(Y)
        if i ∉ indy
            print("$(Y[i]), ")
        end
    end
end

function get_info_deep(
    d::T,
    X::Vector{S}, Y::Vector{S}
    ) where {T<:MatchingBasedDistance, S}

    d_g = d.ground_dist

    indx, indy = get_info(d, X, Y)

    outx, outy = (
        [zeros(Bool,i) for i in length.(X)],
        [zeros(Bool,i) for i in length.(Y)]
    )

    for (i,j) in zip(indx, indy)
        outx[i], outy[j] = get_info(d_g, X[i], Y[j])
    end

    return outx, outy

end


text:
"""
    `generatePoint(dims::Int, range::AbstractFloat, offset::AbstractFloat = 0)`

Generates a random `dims`-dimensional point represented as an array.

d ∈ (offset - range / 2, offset + range / 2) for each dimension of the point.
"""
function generatePoint(dims::Int, range::AbstractFloat, offset::AbstractFloat = 0.0)
  return range * rand(dims) - range / 2 + offset
end

"""
Generates `n` random points without outputs. The rest of parameters are same as parameters for `generatePoint`.
"""
function generatePoints(n::Int, dims::Int, range::AbstractFloat, offset::AbstractFloat = 0.0)
  return [generatePoint(dims, range, offset) for i in 1:n]
end

"""
# Arguments
- `amount`: the amount of points to be generated,
- `polyDegree`: the degree of polynomial used for data generation,
- `polyDim`: the spatial dimension of each term in the polynomial used for data generation,
- `ptRange` and `ptOffset`: the same as `generatePoint`,
- `polyRange` and `polyOffset`: the same as `generatePoint`, but for each term in the polynomial used for data generation,
- `outputFile`: name of the output file
"""
function generateRandomData(amount::Int, polyDegree::Int, polyDim::Int, ptRange::Float64, ptOffset::Float64, polyRange::Float64, polyOffset::Float64, outputFile::String)
  vars, mons = generateMonomials(polyDegree, polyDim)
  coeffs = polyRange * rand(length(mons)) - (polyRange / 2 + polyOffset)
  poly = polynomial(coeffs, mons)
  res = generatePoints(amount, polyDegree, ptRange, ptOffset)
  for p in res
    push!(p, poly(vars => p))
  end
  writecsv(outputFile, res)
  println("Polynomial '$poly' was used for data generation")
  println("Done")
end