text:
function moments(data,rx,ry;x0=-2.2,x1=2.2,y0=-4.5,y1=4.5,v0=-400.0,v1=400.0,dx=1.0,dv=20.0,vmax=80.0,denoise=4.0,show_beam=true)
	figmom=Figure(resolution = (rx,ry))
    fv=findall(x -> (x .> v0) .& (x .< v1), data["V"])
    fx=findall(x -> (x .> x0) .& (x .< x1), data["X"])
    fy=findall(x -> (x .> y0) .& (x .< y1), data["Y"])

    cube=data["data"][fx,fy,fv]
	#cube=filter(!isnan, cube)
    cube[cube .< (denoise * data["rms"])] .= 0.0
    X=data["X"][fx]
    Y=data["Y"][fy]
    V=data["V"][fv]

    axmom0= Axis(figmom[2,1])
    axmom0.xticks = collect(-10:dx:10)
    axmom0.yticks = collect(-10:dx:10)
    axmom0.aspect=DataAspect()
    mom0=sum(cube,dims=3)[:,:,1]
    mom0[mom0 .== 0] .= NaN
    hm=heatmap!(axmom0,X,Y,mom0)
    Colorbar(figmom[1, 1], hm, vertical = false,label=L"Flux $[\mathrm{Jy\,km\,s^{1}\,beam^{-1}}]$")
    axmom0.xlabel="X [kpc]"
    axmom0.ylabel="Y [kpc]"

    axmom1= Axis(figmom[2,2])
    axmom1.xticks = collect(-10:dx:10)
    axmom1.yticks = collect(-10:dx:10)
    axmom1.aspect=DataAspect()
    mom0=sum(cube,dims=3)[:,:,1]
    mom1=[sum(V .* cube[ix,iy,:]) for ix in 1:size(cube,1), iy in 1:size(cube,2)] ./mom0
    hm=heatmap!(axmom1,X,Y,mom1,colorrange=[v0,v1],colormap=:seismic)
    cb=Colorbar(figmom[1, 2], hm, vertical = false,label=L"Mean Velocity $[\mathrm{km\,s^{1}}]$")
    cb.ticks = (v0:100:v1)
    axmom1.xlabel="X [kpc]"
    axmom1.ylabel="Y [kpc]"

    axmom2= Axis(figmom[2,3])
    axmom2.xticks = collect(-10:dx:10)
    axmom2.yticks = collect(-10:dx:10)
    axmom2.aspect=DataAspect()
    mom0=sum(cube,dims=3)[:,:,1]
    mom1=[sum(V .* cube[ix,iy,:]) for ix in 1:size(cube,1), iy in 1:size(cube,2)] ./mom0
    mom2=sqrt.([sum((V .- mom1[ix,iy]) .^2.0 .* cube[ix,iy,:]) for ix in 1:size(cube,1), iy in 1:size(cube,2)] ./mom0)
    hm=heatmap!(axmom2,X,Y,mom2,colorrange=[0.0,vmax])#,colormap=:seismic)
    cb=Colorbar(figmom[1, 3], hm, vertical = false,label=L"Velocity Dispersion $[\mathrm{km\,s^{1}}]$")
    cb.ticks = (0.0:20.0:vmax)
    axmom2.xlabel="X [kpc]"
    axmom2.ylabel="Y [kpc]"
    #rowsize!(figmom.layout, 1, 0.95*axmom2.scene.px_area[].widths[2])
	#colsize!(figmom.layout, 1, 0.95*axmom2.scene.px_area[].widths[1])
	if show_beam
		phis_beam=0.0:0.01:2*pi
		major_beam=maximum(data["beam"][1:2])/2.0
		minor_beam=minimum(data["beam"][1:2])/2.0
		X_beam=0.6
		Y_beam=0.25
		i_beam=acos(minor_beam/major_beam)
		pa_beam=deg2rad(data["beam"][3])
		R_beam=major_beam ./(1.0 .+tan.(i_beam)^2.0 .* cos.(phis_beam .-pa_beam).^2.0).^0.5
		x_beam=R_beam .*cos.(phis_beam) .+0.8 *x1
		y_beam=R_beam .*sin.(phis_beam) .+ (y0+abs(0.1*y0))
		lines!(axmom0,x_beam,y_beam, label=false, transparency = true,linewith=4.0, color=:blue)
	end

    return figmom
end

function plot_pvd(data,ang;cdata=nothing,clouds=nothing,disks=nothing,merged=false,slit=0.1,zoom_pvd=[-4.1,4.1,-370.0,370.0,1.0,50.0,150.0])
    θ=deg2rad(ang)
    ϕ=θ+pi/2
    figpvf=Figure(resolution=(1.9*600,650))
    axpvf=Axis(figpvf[1,1], xlabel=L"Projected radius [$\mathrm{kpc}$]",ylabel=L"Projected Velocity [$\mathrm{km\,s^{-1}}$]", labelsize=20)
    contourf!(axpvf, data["Y"], data["V"], pvd(ang,data;slit=slit), levels=collect(2:2:20), colormap=:roma,extendhigh=:red)
    if !isnothing(cdata)
        contour!(axpvf, cdata["Y"], cdata["V"], pvd(ang,cdata;slit=slit), levels=collect(2:2:20))
    end
    #heatmap!(axpvf,data["Y"],data["V"],pvd(ang,data;slit=slit))
    if !isnothing(clouds)
        X=merged ? clouds.Xc : clouds.Xp
        Y=merged ? clouds.Yc : clouds.Yp
        V=merged ? clouds.Vc : clouds.Vp
        S=merged ? clouds.Sc : clouds.Sp

        Φ=atan.(Y,X)
        Rp=sqrt.(Y .^2.0 .+X .^2.0 )
        R_project=F(θ,X,Y,Rp;slit=slit).* (-sign.(cos.(F(θ,X,Y,Φ;slit=slit) .- (θ-pi/2.0))))
        Vslit=F(θ,X,Y,V;slit=slit)
        dVslit=F(θ,X,Y,S;slit=slit)
        chm=GLMakie.errorbars!(R_project, Vslit,dVslit, markersize=10)
    end

    if !isnothing(disks)
        for (n,di) in enumerate(keys(disks))
			r_model_m=Vector{Float64}()
			v_model_m=Vector{Float64}()
			r_model=Vector{Float64}()
			v_model=Vector{Float64}()

			for (rd,i_d,phi_d,v_d) in zip(disks[di].R, disks[di].I, disks[di].PA,disks[di].V)
				Rd_sky=rd/sqrt(1.0+tan(i_d)^2.0*cos(ϕ-phi_d)^2.0)
				Vd_sky=v_d*sin(i_d)*sin(ϕ-phi_d)/sqrt(1.0 +tan(i_d)^2.0*cos(ϕ-phi_d)^2.0)
				push!(r_model_m,-Rd_sky)
				push!(v_model_m,-Vd_sky)
				push!(r_model,Rd_sky)
				push!(v_model,Vd_sky)
			end
			lines!(axpvf,r_model_m,v_model_m,color=Cycled(n),linewidth=4)
			lines!(axpvf,r_model,v_model,color=Cycled(n),linewidth=4,label=di)
		end
		axislegend(axpvf)
    end
    axpvf.xticks = -37:zoom_pvd[5]:37
    axpvf.yticks =	-500:zoom_pvd[6]:500
    xlims!(axpvf,zoom_pvd[1],zoom_pvd[2])
    ylims!(axpvf,zoom_pvd[3],zoom_pvd[4])
    return figpvf
end

function plot_disks(disks,P;sample=nothing)
	fig=Figure()
	axpa=Axis(fig[1, 1])
	axi=Axis(fig[2, 1])

	for di in keys(disks)
		lines!(axpa,disks[di].R,rad2deg.(disks[di].PA),label=di)
		lines!(axi,disks[di].R,rad2deg.(disks[di].I),label=di)
	end
	axislegend(axpa)
	axislegend(axi)
	if !isnothing(sample)
        P_tmp=copy(P)
		alpha_pvd=7.0/size(sample)[1]
		for ppi in 1:size(sample)[1]
			pp_i=sample[ppi,:]
			update_parameters!(P_tmp,pp_i)
			disks_tmp = make_disks(P.Rd,P_tmp)
			lines!(axpa,disks_tmp.R,rad2deg.(disks_tmp.PA), transparency=true, color=RGBA(0.0,0.0,0.0,alpha_pvd))
			lines!(axi,disks_tmp.R,rad2deg.(disks_tmp.I), transparency=true, color=RGBA(0.0,0.0,0.0,alpha_pvd))
		end
	end

	axpa.ylabel="Position Angle [deg]"
	axi.ylabel="Inclination [deg]"
	axi.xlabel="Radius [kpc]"
	axpa.xticks = 0.0:0.5:16.0
	axi.xticks = 0.0:0.5:16.0
	fig
end

function plot_disksv(disks;sample=nothing)
	fig=Figure()
	axpa=Axis(fig[1, 1])
	axi=Axis(fig[2, 1])

	for di in keys(disks)
		lines!(axpa,disks[di].R,disks[di].V,label=di)
	end
	axislegend(axpa)
	# if !isnothing(sample)
	# 	alpha_pvd=7.0/size(sample)[1]
	# 	for ppi in 1:size(sample)[1]
	# 		pp_i=sample[ppi,:]
	# 		update_parameters!(P,pp_i)
	# 		disks_tmp = make_disks(P.Rd,P)
	# 		lines!(axpa,disks_tmp.R,rad2deg.(disks_tmp.PA), transparency=true, color=RGBA(0.0,0.0,0.0,alpha_pvd))
	# 		lines!(axi,disks_tmp.R,rad2deg.(disks_tmp.I), transparency=true, color=RGBA(0.0,0.0,0.0,alpha_pvd))
	# 	end
	# end
	axpa.ylabel="Rotational Velocity [km/s]"
	axpa.xlabel="Radius [kpc]"
	axpa.xticks = 0.0:0.5:16.0
	fig
end

function plot_sky(clouds::Clouds,disks::Disks;msize=10,zoom_image=[-2.0,2.0,-5.0,5.0,1.0,1.0])
    imratio=(zoom_image[2]-zoom_image[1])/(zoom_image[4]-zoom_image[3])
    figsky=Figure(resolution=(imratio*800,820))
	#ai = show_data ? 2 : 1
	axsky = Axis(figsky[2, 1],xgridvisible=true,ygridvisible=true)
	axsky.aspect = DataAspect()
    cmap=((:blue, 0.95), (:red, 0.95))
    sc=scatter!(axsky,clouds.Xp,clouds.Yp,color=clouds.dV,colorrange=(-80,80),markersize=msize,marker=:diamond,transparency=true,colormap=:bluesreds,label="CO(2-1) residuals")
    Colorbar(figsky[1,1], sc,vertical=false,label=L"dV [$\mathrm{km\,s^{-1}}$]")

    phis=0.0:0.01:2*pi
	for (rdi,i_di,pa_di) in zip(disks.R,disks.I,disks.PA)
        Rd_sky=rdi ./(1.0 .+tan.(i_di)^2.0 .* cos.(phis .-pa_di).^2.0).^0.5
        xd=Rd_sky .*cos.(phis)
        yd=Rd_sky .*sin.(phis)
        lines!(axsky,xd,yd,transparency = true, color=RGBA(0.0,0.0,0.0,0.1))
	end
    axsky.xticks = -4:zoom_image[5]:4
	axsky.yticks =	-5:zoom_image[6]:5
	xlims!(axsky,zoom_image[1],zoom_image[2])
	ylims!(axsky,zoom_image[3],zoom_image[4])
	axsky.xlabel="X [kpc]"
	axsky.ylabel="Y [kpc]"
    return figsky
end

function plot_cube_model(data,model;disks=nothing,zoom_image=[-2.0,2.0,-5.0,5.0,1.0,1.0])
    imratio=(zoom_image[2]-zoom_image[1])/(zoom_image[4]-zoom_image[3])

    figd=Figure(resolution=(imratio*800,820))
    axd = Axis(figd[2, 1],xgridvisible=true,ygridvisible=true,ylabel="Y [kpc]",xlabel="X [kpc]")
    axd.aspect = DataAspect()

    momm=sum(model["data"],dims=3)[:,:,1]
    data_cube=copy(data["data"])
    momd=sum(data_cube,dims=3)[:,:,1]

    maxmask=findall(x->x < 4*data["rms"],maximum(data["data"],dims=3)[:,:,1])
    momd[maxmask] .= NaN

    rres=(momm .- momd) ./(data["rms"]*size(data["data"])[3])
    sc=GLMakie.contourf!(axd,model["X"],model["Y"], rres,levels=[-2,-1.5,-1.25,-1,-0.75,-0.5,-0.25,0.25,0.5,0.75,1,1.25,1.5,2],colormap=:bluesreds)
    GLMakie.Colorbar(figd[1,1],sc,label="(Model-Data)/(rms * channels)",vertical=false,fontsize=25)

    if !isnothing(disks)
        phis_=collect(0:0.01:2*pi)
        for (rdi,i_di,pa_di) in zip(disks.R,disks.I,disks.PA)
            Rd_sky=rdi ./(1.0 .+tan.(i_di)^2.0 .* cos.(phis_ .-pa_di).^2.0).^0.5
            xd=Rd_sky .*cos.(phis_)
            yd=Rd_sky .*sin.(phis_)
            lines!(axd,xd,yd, label=false, transparency = true, color=RGBA(0.0,0.0,0.0,0.06))
        end
    end

    axd.xticks = -4:zoom_image[5]:4
	axd.yticks =	-5:zoom_image[6]:5
	xlims!(axd,zoom_image[1],zoom_image[2])
	ylims!(axd,zoom_image[3],zoom_image[4])
	axd.xlabel="X [kpc]"
	axd.ylabel="Y [kpc]"
    figd
end

function cloud_fitting_diagnostics(data,clouds::Clouds,savename;sigma=4.0,pvds=[0.0,90.0],slit=0.1)
    mom0_data = sum(data["data"] .> sigma*data["rms"], dims = 3)[:, :, 1];

    Xp=clouds.Xp
    Yp=clouds.Yp
    Vp=clouds.Vp
    Ip=clouds.Ip
    Sp=clouds.Sp
    Rp=sqrt.(Xp .^2 .+ Yp .^2)
    Φp=atan.(Yp,Xp)

    figd=Figure(resolution = (1500,700))
    axRc= Axis(figd[1,1],Aspect=:equal)
    axxy= Axis(figd[1,2],Aspect=:equal)
    axS= Axis(figd[1,3],Aspect=:equal)

    heatmap!(axxy,data["X"],data["Y"],mom0_data)
    lines!(axRc,Rp[sortperm(Rp)],log10.(cumsum(Ip[sortperm(Rp)])))
    scatter!(axxy,Xp,Yp,color=Vp,colorrange=[-300.0,300.0],colormap=:seismic,markersize=3,trasnparency=true)

    scatter!(axS,Rp,Sp,markersize=2,trasnparency=true)
    if length(clouds.Xc)<1
       #save("$(savename)_fit.png",figd)
    end
    # for pvd in pvds
    #     save("$(savename)_pvd_$(pvd).png",plot_pvd(data,pvd;clouds=clouds,merged=false))
    # end

    if length(clouds.Xc)>1
        Xc=clouds.Xc
        Yc=clouds.Yc
        Vc=clouds.Vc
        Ic=clouds.Ic
        Sc=clouds.Sc
        Rc=sqrt.(Xc .^2 .+ Yc .^2)
        Φc=atan.(Yc,Xc)

        lines!(axRc,Rc[sortperm(Rc)],log10.(cumsum(Ic[sortperm(Rc)])))
        #scatter!(axxy,Xc,Yc,color=Vc,colorrange=[-300.0,300.0],colormap=:seismic,markersize=10,trasnparency=true)

        scatter!(axS,Rc,Sc,markersize=10,trasnparency=true)
        #save("$(savename)_fit.png",figd)

        figd=Figure(resolution = (1500,700))
        axxy1= Axis(figd[1,1],Aspect=:equal)
        axxy= Axis(figd[1,2],Aspect=:equal)
        scatter!(axxy1,Xc,Yc,color=Vc,colorrange=[-300.0,300.0],colormap=:seismic,markersize=10,trasnparency=true)
        for (n,i) in enumerate(unique(clouds.I))
            color=Cycled(n)
            xx=Xp[clouds.belongs .== i]
            yy=Yp[clouds.belongs .== i]
            scatter!(axxy,xx,yy,color=color,markersize=8)
        end
        scatter!(axxy,Xc,Yc,color=:black,marker='X',markersize=12)
        save("$(savename)_merging.png",figd)

        for pvdi in pvds
            θ=deg2rad(pvdi)
            figpvf=Figure(resolution=(1.9*1000,1050))
            axpvf=Axis(figpvf[1,1], xlabel=L"Projected radius [$\mathrm{kpc}$]",ylabel=L"Projected Velocity [$\mathrm{km\,s^{-1}}$]", labelsize=20)
            try
                contourf!(axpvf, data["Y"], data["V"], pvd(pvdi,data;slit=slit), levels=collect(2:2:20), colormap=:roma,extendhigh=:red)
            catch er
                @warn "Cannot make contour plots, trying heatmap"
                heatmap!(axpvf, data["Y"], data["V"], pvd(pvdi,data;slit=slit),colorrange=[2,20],colormap=:roma,extendhigh=:red)
            end
            for (n,i) in enumerate(unique(clouds.I))
                color=Cycled(n)
                xx=Xp[clouds.belongs .== i]
                yy=Yp[clouds.belongs .== i]
                vv=Vp[clouds.belongs .== i]

                rc=sqrt.(xx .^2.0 .+ yy .^2.0)
                phic=atan.(yy,xx)
                rr_p=F(θ,xx,yy,rc;slit=slit).* (-sign.(cos.(F(θ,xx,yy,phic;slit=slit) .- (θ-pi/2.0))))

                vv_p=F(θ,xx,yy,vv;slit=slit)

                scatter!(axpvf,rr_p,vv_p,color=color,markersize=3,transparency=true)

                xx=[Xp[i]]
                yy=[Yp[i]]
                vv=[Vp[i]]
                rc=sqrt.(xx .^2.0 .+ yy .^2.0)
                phic=atan.(yy,xx)
                rr_p=F(θ,xx,yy,rc;slit=slit).* (-sign.(cos.(F(θ,xx,yy,phic;slit=slit) .- (θ-pi/2.0))))

                vv_p=F(θ,xx,yy,vv;slit=slit)

                scatter!(axpvf,rr_p,vv_p,color=:black,marker=:rect,markersize=13,transparency=true)
            end
            xx=Xp
            yy=Yp
            vv=Vp
            rc=sqrt.(xx .^2.0 .+ yy .^2.0)
            phic=atan.(yy,xx)
            rr_p=F(θ,xx,yy,rc;slit=slit).* (-sign.(cos.(F(θ,xx,yy,phic;slit=slit) .- (θ-pi/2.0))))

            vv_p=F(θ,xx,yy,vv;slit=slit)

            scatter!(axpvf,rr_p,vv_p,color=:black,marker=:rect,markersize=2,transparency=true)

            #scatter!(axpvf,F(θ,Xc,Yc,Rc;slit=slit).* (-sign.(cos.(F(θ,Xc,Yc,Φc;slit=slit) .- (θ-pi/2.0)))),F(θ,Xc,Yc,Vc;slit=slit),color=:black,marker='X',markersize=12)

            #save("$(savename)_pvd_$(pvdi)_merged.png",figpvf)
        end

        # figb=Figure(resolution = (2000,1000))
        # axxyb= Axis(figb[1,1],Aspect=:equal)
        # axyvb = Axis(figb[1,2],Aspect=:equal)
        # scatter!(axxyb,Xp,Yp,color=clouds.belongs,makersize=7)
        # scatter!(axyvb,Yp,Vp,color=clouds.belongs,makersize=7)
        # for (n,i) in enumerate(clouds.I)
        #     scatter!(axxyb,[Xp[i]],[Yp[i]],color=:red,marker="$(n%9)",markersize=15)
        #     scatter!(axyvb,[Yp[i]],[Vp[i]],color=:red,marker="$(n%9)",markersize=15)
        # end
        # save("$(savename)_fit_merging.png",figb)
    end
    return figd,figpvf
end


text:
# Contact between an ellipsoid and a plane
# ========================================
# References
# ----------
# [1] P. P. Klein, On the Ellipsoid and Plane Intersection
#     Equation, Applied Mathematics, 2012, 3, 1634-1640
# -----------------------------------------------------------------------
# ellipsoid is defined by centre x, radii a, b, c, and orientation R;
# plane is defined by point and normal;
# -------------------------------------------------------------
# returned: (conpnt, depth, A, B), where:
# conpnt - contact point (or [0.0, 0.0, 0.0] if no contact)
# depth - positive penetration depth (or -1.0 if no contact)
# A, B - contact ellipse radii (or 0.0, 0.0 if no contact)
# _____________________________________________________________
function ellipsoid_plane_contact(x, a, b, c, R, point, normal)

  # these are 3-vectors
  x = vec(x);
  point = vec(point);
  normal = vec(normal);

  # we move the plane into the natural coordinate system of our input ellipsoid,
  # which is centered at x and oriented according to the R rotation

  point = R'*(point - x); # so that plane's point, relative to the origin of
                          # the ellipsoid is suitably shifted and rotated
                          # and this needs to be followed by

  normal = R' * normal; # to redefine the normal in the natural coordinate
                        # system of the ellipsoid

  C1 = [1.0/a 0.0 0.0; 0.0 1.0/b 0.0; 0.0 0.0 1.0/c]; # this is matrix D1 in [1]

  v = vec(rand(3,1));   # let's generate a random vector first
  r = cross(normal, v); # now r'*normal = 0 as required
  r = r / norm(r);      # lets normalize it (make it's length = 1)
  s = cross(normal, r); # now s'*normal = 0 and s'*r = 0 as required
  q = point;

  # now lets check for condition (7) to be fullfilled

  # by doing this transformation we are already satisfying condition (4-6)
  # as rr and ss are orthogonal to each other and to the normal vector

  # now lets choose omega and fulfill condition (7)

  # here I will define all the dot products needed
  D1 = dot((C1*r),(C1*s));
  D2 = dot((C1*r),(C1*r));
  D3 = dot((C1*s),(C1*s));

  if abs(D2-D3) < 1E-10
    omega = 0.25*pi;
  else
    omega = 0.5 * atan((2*D1)/(D2-D3));
  end

  rr = cos(omega)* r + sin(omega) * s;  # here we are transforming r and
  ss = -sin(omega) * r + cos(omega) * s; # s to fullfill condition (7) in [1]

  r = rr;
  s = ss;

  D2 = dot((C1*r),(C1*r));
  D3 = dot((C1*s),(C1*s));
  D4 = dot((C1*q),(C1*r));
  D5 = dot((C1*q),(C1*s));
  D6 = dot((C1*q),(C1*q));

  # now lets find the ellipse of the intersection
  d = D6 -((D4)^2/D2)-((D5)^2/D3); # this is just a value

  # semi-axes squared radii of the contact ellipse
  Aa = (1-d)/D2;
  Bb = (1-d)/D3;

  if Aa > 0. && Bb > 0. # there is contact
    # radii of the contact ellipse
    A  = sqrt(Aa);
    B  = sqrt(Bb);

    # now the centre of the contact ellipse
    t0 = -D4/D2;
    u0 = -D5/D3;
    x0 = q + t0*r + u0*s;

    # now lets find the contact point (in ellipsoid coordinates)
    alpha = sqrt(1.0/((x0[1]^2/a^2) + (x0[2]^2/b^2) + (x0[3]^2/c^2)));
    conpnt = alpha * x0;

    # and the penetration depth
    depth = dot(normal,  x0 - conpnt);

    # now contact point in input coordinates
    conpnt = R*conpnt + x;

  else # there is no contact
    conpnt = [0.0; 0.0; 0.0];
    depth = -1.0;
    A = 0.0;
    B = 0.0;
  end

  return conpnt, depth, A, B
end

ext: jl

max_stars_repo_name: tkoziara/edyna

max_issues_repo_name: tkoziara/edyna

max_forks_repo_name: tkoziara/edyna


text:
# ------------------------------------------------------------------------------
#                                   TESTS
# ------------------------------------------------------------------------------

using Base.Test
importall MKL

# ==============================================================================
# BLAS LEVEL 1
# ==============================================================================

## ?asum
fname = :mkl_asum
@eval begin
    res = $fname(4,Float32[1,2,3,4],1)
    @test  res == 10 && typeof(res) == Float32
    res = $fname(4,Float64[1,2,3,4],1)
    @test  res == 10 && typeof(res) == Float64
end
# ------------------------------------------------------------------------------

## ?axpy!
fname = :mkl_axpy!
@eval begin
    res = $fname(4,convert(Float32,2.0),Float32[1,2,3,4],1,Float32[-1,-1,-1,-1],1)
    @test  res == Float32[1,3,5,7] && typeof(res) == Vector{Float32}
    res = $fname(4,convert(Float64,2.0),Float64[1,2,3,4],1,Float64[-1,-1,-1,-1],1)
    @test  res == Float64[1,3,5,7] && typeof(res) == Vector{Float64}
end
# ------------------------------------------------------------------------------

## ?copy!
fname = :mkl_copy!
@eval begin
    res = $fname(4,Float32[1,2,3,4],1,Float32[-1,-1,-1,-1],1)
    @test  res == Float32[1,2,3,4] && typeof(res) == Vector{Float32}
    res = $fname(4,Float64[1,2,3,4],1,Float64[-1,-1,-1,-1],1)
    @test  res == Float64[1,2,3,4] && typeof(res) == Vector{Float64}
end
# ------------------------------------------------------------------------------

## ?dot
fname = :mkl_dot
@eval begin
    res = $fname(4,Float32[1,2,3,4],1,Float32[1,2,3,4],1)
    @test  res == 30 && typeof(res) == Float32
    res = $fname(4,Float64[1,2,3,4],1,Float64[1,2,3,4],1)
    @test  res == 30 && typeof(res) == Float64
end
# ------------------------------------------------------------------------------

## sdot
fname = :mkl_sdot
@eval begin
    res = $fname(4,Float32[1,2,3,4],1,Float32[1,2,3,4],1)
    @test  res == 30 && typeof(res) == Float64
end
# ------------------------------------------------------------------------------

## nrm2
fname = :mkl_nrm2
@eval begin
    res = $fname(4,Float32[1,2,3,4],1)
    @test  res == convert(Float32,sqrt(1+2^2+3^2+4^2)) && typeof(res) == Float32
    res = $fname(4,Float64[1,2,3,4],1)
    @test  res == convert(Float64,sqrt(1+2^2+3^2+4^2)) && typeof(res) == Float64
end
# ------------------------------------------------------------------------------

## scal!
fname = :mkl_scal!
@eval begin
    res = $fname(4,convert(Float32,2.0),Float32[1,2,3,4],1)
    @test  res == Float32[2,4,6,8] && typeof(res) == Vector{Float32}
    res = $fname(4,convert(Float64,2.0),Float64[1,2,3,4],1)
    @test  res == Float64[2,4,6,8] && typeof(res) == Vector{Float64}
end
# ------------------------------------------------------------------------------

# ==============================================================================
# VML
# ==============================================================================

## Add!
fname = :mkl_add!
@eval begin
    res = $fname(4,Float32[1,2,3,4],Float32[1,2,3,4], zeros(Float32,4))
    @test  res == Float32[2,4,6,8] && typeof(res) == Vector{Float32}
    res = $fname(4,Float64[1,2,3,4],Float64[1,2,3,4], zeros(Float64,4))
    @test  res == Float64[2,4,6,8] && typeof(res) == Vector{Float64}
end
# ------------------------------------------------------------------------------

# Sub!
fname = :mkl_sub!
@eval begin
    res = $fname(4,Float32[1,2,3,4],Float32[1,2,3,4], zeros(Float32,4))
    @test  res == Float32[0,0,0,0] && typeof(res) == Vector{Float32}
    res = $fname(4,Float64[1,2,3,4],Float64[1,2,3,4], zeros(Float64,4))
    @test  res == Float64[0,0,0,0] && typeof(res) == Vector{Float64}
end
# ------------------------------------------------------------------------------

## Sqr!
fname = :mkl_sqrt!
@eval begin
    res = $fname(4,Float32[1,4,9,16],Float32[1,1,1,1])
    @test  res == sqrt(Float32[1,4,9,16]) && typeof(res) == Vector{Float32}
    res = $fname(4,Float64[1,4,9,16],Float64[1,1,1,1])
    @test  res == sqrt(Float64[1,4,9,16]) && typeof(res) == Vector{Float64}
end

# ------------------------------------------------------------------------------

# Mul!
fname = :mkl_mul!
@eval begin
    res = $fname(4,Float32[1,2,3,4],2*Float32[1,2,3,4], zeros(Float32,4))
    @test  res == Float32[2,8,18,32] && typeof(res) == Vector{Float32}
    res = $fname(4,Float64[1,2,3,4],2*Float64[1,2,3,4], zeros(Float64,4))
    @test  res == Float64[2,8,18,32] && typeof(res) == Vector{Float64}
end
# ------------------------------------------------------------------------------

## Abs!
fname = :mkl_abs!
@eval begin
    res = $fname(4,Float32[-1,2,-3,4],zeros(Float32,4))
    @test  res == Float32[1,2,3,4] && typeof(res) == Vector{Float32}
    res = $fname(4,Float64[-1,2,-3,4],zeros(Float64,4))
    @test  res == Float64[1,2,3,4] && typeof(res) == Vector{Float64}
end


text:
@testset "consistent.jl" begin
    @testset "probabilities" begin
        probs = rand(100)
        c = @inferred(Consistent(probs))
        @test c isa Consistent{Float64,Bool}
        @test c.predictions == probs
        @test eltype(c) === Tuple{Float64,Bool}

        @test_throws ErrorException Consistent(Float64[])
        @test_throws ErrorException Consistent([0.5, prevfloat(0.0)])
        @test_throws ErrorException Consistent([0.1, nextfloat(1.0)])
    end

    @testset "categorical" begin
        predictions = map(1:100) do _
            prediction = rand(5)
            lmul!(inv(sum(prediction)), prediction)
            return prediction
        end
        c = @inferred(Consistent(predictions))
        @test c isa Consistent{Vector{Float64},Int}
        @test c.predictions == predictions
        @test eltype(c) === Tuple{Vector{Float64},Int}

        @test_throws ErrorException Consistent(Vector{Float64}[])
        @test_throws ErrorException Consistent([[0.0, 1.0], [0.1, 0.2]])
    end

    @testset "Distributions" begin
        predictions = map(Normal, randn(100), rand(100))
        c = @inferred(Consistent(predictions))
        @test c isa Consistent{Normal{Float64},Float64}
        @test c.predictions == predictions
        @test eltype(c) === Tuple{Normal{Float64},Float64}

        @test_throws ErrorException Consistent(Normal{Float64}[])
    end
end


text:
function numquad(integrand, test_local_space, trial_local_space,
    test_chart, trial_chart, out)

    outer_qps = quadpoints(test_chart, 10)
    inner_qps = quadpoints(trial_chart, 11)

    M = numfunctions(test_local_space)
    N = numfunctions(trial_local_space)

    for (x,wx) in outer_qps
        f = test_local_space(x)
        for (y,wy) in inner_qps
            g = trial_local_space(y)
            G = integrand(x,y)
            ds = wx * wy
            out += SMatrix{M,N}([dot(f[i][1], G*g[j][1]) for i ∈ 1:M, j ∈ 1:N])*ds
        end
    end

    out
end

function numquad_cf(integrand, test_local_space, trial_local_space,
    test_chart, trial_chart, out)

    outer_qps = quadpoints(test_chart, 10)

    inner_ksi = quadpoints(simplex(point(0),point(1)), 13)
    inner_eta = quadpoints(simplex(point(0),point(1)), 13)

    @assert sum(w for (x,w) in inner_ksi) ≈ 1
    @assert sum(w for (x,w) in inner_eta) ≈ 1

    @assert sum(w for (p,w) in outer_qps) ≈ volume(test_chart)
    @assert sum(w1*w2*(1-cartesian(p2)[1]) for (p1,w1) in inner_ksi for (p2,w2) in inner_eta) ≈ 0.5

    p1, p2, p3 = trial_chart.vertices

    M = numfunctions(test_local_space)
    N = numfunctions(trial_local_space)

    @assert volume(trial_chart) ≈ volume(test_chart)

    @assert M == N == 3

    check = 0.0

    for (x,wx) in outer_qps
        s1 = simplex(p3,cartesian(x),p2)
        s2 = simplex(p1,cartesian(x),p3)
        s3 = simplex(p2,cartesian(x),p1)

        @assert volume(s1) + volume(s2) + volume(s3) ≈ volume(trial_chart)

        qq = carttobary(trial_chart, cartesian(x))
        @assert 0 <= qq[1] <= 1
        @assert 0 <= qq[2] <= 1
        @assert 0 <= 1-qq[1]-qq[2] <= 1

        f = test_local_space(x)
        @assert length(f) == 3
        @assert f[1][1] isa SVector
        inner_check = 0.0
        for (ksi,wksi) in inner_ksi
            for (eta,weta) in inner_eta

                p = cartesian(ksi)[1]
                q = cartesian(eta)[1]

                @assert p isa Float64
                @assert q isa Float64

                u = p*(1-q)
                v = q
                wuv = wksi*weta*(1-cartesian(eta)[1])

                @assert 0 <= u <= 1
                @assert 0 <= v <= 1
                @assert 0 <= 1-u-v <= 1

                y1 = neighborhood(s1,(u,v))
                y2 = neighborhood(s2,(u,v))
                y3 = neighborhood(s3,(u,v))

                u1 = carttobary(s1, cartesian(y1))
                u2 = carttobary(s2, cartesian(y2))
                u3 = carttobary(s3, cartesian(y3))
                @assert 0 <= u1[1] <= 1
                @assert 0 <= u1[2] <= 1
                @assert 0 <= u2[1] <= 1
                @assert 0 <= u2[2] <= 1
                @assert 0 <= u3[1] <= 1
                @assert 0 <= u3[2] <= 1

                g1 = trial_local_space(neighborhood(trial_chart,carttobary(trial_chart,cartesian(y1))))
                g2 = trial_local_space(neighborhood(trial_chart,carttobary(trial_chart,cartesian(y2))))
                g3 = trial_local_space(neighborhood(trial_chart,carttobary(trial_chart,cartesian(y3))))

                @assert length(g1) == 3
                @assert length(g2) == 3
                @assert length(g3) == 3

                @assert g1[1][1] isa SVector
                @assert g2[1][1] isa SVector
                @assert g3[1][1] isa SVector

                @assert jacobian(y1) >= 0
                @assert jacobian(y2) >= 0
                @assert jacobian(y3) >= 0
                @assert jacobian(y1)+jacobian(y2)+jacobian(y3) ≈ 2*volume(test_chart)

                G1 = integrand(x,y1)*jacobian(y1)
                G2 = integrand(x,y2)*jacobian(y2)
                G3 = integrand(x,y3)*jacobian(y3)

                @assert wksi > 0
                @assert weta > 0
                @assert wx*wuv > 0

                check += wx*wuv
                inner_check += wuv

                out += wx*wuv*(SMatrix{M,N}([dot(f[i][1], (G1*g1[j][1]+G2*g2[j][1]+G3*g3[j][1])) for i=1:M, j=1:N]))
            end
        end
        @assert inner_check ≈ 0.5
    end

    @assert check ≈ volume(test_chart)*0.5
    out
end


text:
using DataStructures
import IterTools

# From Julia's v0.6 reslease notes:
# The Collections module has been removed, and all functions defined
# therein have been moved to the DataStructures package (#19800).
"""
`erdos_renyi_undirected`
========================

Generate an undirected Erdős-Rényi graph.
An undirected
Erdős-Rényi graph is generated by letting each
undirected edge between `n` nodes
be present with probability `p`.

There is another form of the call where the
input is the average degree of the generated
graph.

**The current implementation uses sprand, this may change in the future.**
**Do not depend on this routine for reliable output between versions.**

Input
-----
- `n`: the number of nodes
- `p`: the probability of an edge, or the average degree.
   if \$`p` >= 1\$, then \$`p`\$ is interpreted as an average degree
   instead of a probability. (There is no point in generating
   an Erdős-Rényi graph with probability \$`p`=1\$)
- `d`: The desired average degree, converted into a probabiltity
   via d/n

Functions
---------
- `erdos_renyi_undirected(n::Int,p::Float64)` specify the
  probability or average degree
- `erdos_renyi_undirected(n::Int,d::Int)` specify the average degree directly

Output
------
- A matrix network type for the Erdős-Rényi graph.

Example
-------
~~~~
# show the connected phase transition
n = 100
avgdegs = linspace(1.,2*log(n),100)
compsizes = map( (dbar) ->
        maximum(scomponents(erdos_renyi_undirected(n,dbar)).sizes),
    avgdegs )
using Plots
unicodeplots()
plot(avgdegs,compsizes,xaxis=("average degree"),yaxis=("largest component size"))
~~~~
"""
function erdos_renyi_undirected(n::Int, p::Float64)
    if p < 0 || p > n throw(DomainError(p)) end

    if p >= 1. # interpret as average degree
        p = p/n # convert to probability
    end
    if isnan(p)
        return _matrix_network_direct(spzeros(n,n),1)
    end
    A = sprand(n,n,p)
    if n > 0
        Aup = triu(A,1)
    else
        Aup = A
    end
    Asym = max.(Aup,Aup')
    return _matrix_network_direct(Asym,1)
end
erdos_renyi_undirected(n::Int, d::Int) = erdos_renyi_undirected(n, d/n)

erdős_rényi_undirected = erdos_renyi_undirected

"""
`erdos_renyi_directed`
========================

Generate an directed Erdős-Rényi graph. A directed
Erdős-Rényi graph is generated by letting each
directed edge between `n` nodes
be present with probability `p`.

There is another form of the call where the
input is the average degree of the generated
graph.

**The current implementation uses sprand, this may change in the future.**
**Do not depend on this routine for reliable output between versions.**

Input
-----
- `n`: the number of nodes
- `p`: the probability of an edge.
- `d`: the average degree

Functions
---------
- `erdos_renyi_directed(n::Int,p::Float64)` specify the probability
  (p < 1.) or the average degree (p >= 1.)
- `erdos_renyi_directed(n::Int,d::Int)` specify the average degree
  directly

Output
------
- A matrix network type for the Erdős-Rényi graph.
"""
function erdos_renyi_directed(n::Int, p::Float64)
    if p < 0 || p > n throw(DomainError(p)) end
    if p >= 1. # interpret as average degree
        p = p/n # convert to probability
    end
    if isnan(p)
        return _matrix_network_direct(spzeros(n,n))
    end

    A = sprand(n,n,p)

    # return _matrix_network_direct(A-spdiagm(diag(A),0)) # directions don't matter
    return _matrix_network_direct(A-spdiagm(0 => diag(A))) # directions don't matter
end

erdos_renyi_directed(n::Int, d::Int) = erdos_renyi_directed(n,d/n)

erdős_rényi_directed = erdos_renyi_directed

"""
`chung_lu_undirected`
========================

Generate an approximate undirected Chung-Lu graph. The approximation is because we draw
exactly |E| edges where each edge is sampled from the Chung-Lu model. But then
we discard duplicate edges and self-loops. So the new graph will always have fewer
edges than the input degree sequence.

**This will likely change in future versions and provide an exact Chung-Lu model.**

If the graph is

Usage
-----
- `chung_lu_undirected(d)`
- `chung_lu_undirected(d,nedges)`

Input
-----
- `d`: the degree sequence vector. This vector is non-negative and has the expected
degree for each vertex.

Output
------
- A MatrixNetwork for the undirected graph that results from the Chung-Lu sample.

Example
-------
~~~~
A = load_matrix_network("tapir")
d = vec(sum(A,1))
B = sparse(chung_lu_undirected(d))
nnz(A)
nnz(B)
~~~~

"""
function chung_lu_undirected end

function chung_lu_undirected(d::Vector{Int})
    chung_lu_undirected(d, floor(Int,sum(d)/2))
end

function chung_lu_undirected(d::Vector{Int}, nedges::Int)
    n = length(d)

    # TODO find some standard function for this
    for v in d
        if v < 0
            throw(DomainError(v))
        end
    end

    if nedges < 0 || nedges > div(n*(n-1),2)
        throw(ArgumentError("nedges $nedges is too large for $n node undirected graph"))
    end

    nodevec = zeros(Int,sum(d))
    curedge = 1
    for i=1:n
        for j=1:d[i]
            nodevec[curedge] = i
            curedge += 1
        end
    end

    ei, ej = unique_edge_sample_undirected(nodevec, nedges)
    A = sparse(ei,ej,1.,n,n)

    return _matrix_network_direct(A) # avoid the transpose
end

function unique_edge_sample_undirected(nodevec, nedges::Int)
    edges = Set{Tuple{Int,Int}}()
    sizehint!(edges, nedges)

    curedges = 0
    while curedges < nedges
        src = rand(nodevec)
        dst = rand(nodevec)
        if src == dst
            continue
        else
            # fix the order
            if src < dst
                dst, src = src, dst
            end
            if !((src,dst) in edges)
                push!(edges, (src,dst))
                curedges += 1
            end
        end
    end

    ei = zeros(Int,nedges*2)
    ej = zeros(Int,nedges*2)

    k = 1
    for edge in edges
        ei[k] = edge[1]
        ej[k] = edge[2]
        k+=1
        ei[k] = edge[2]
        ej[k] = edge[1]
        k+=1
    end

    return ei, ej
end

"""
Not public right now
"""
function _chung_lu_dense_undirected(d::Vector{Int})
    M = zeros(Int,n,n)
    idenom = 1/sum(d)
    for i=1:n
        for j=i+1:n
            if rand(Float64) < d[i]*d[j]*idenom
                M[i,j] = 1
            end
        end
    end
    M = M + M'
    return _matrix_network_direct(sparse(M))
end

"""
The internal Havel-Hakimi function has an optional store
behavior that saves the edges as they come out of the algorithm.
This enables us to generate a Havel Hakimi graph, which can
be useful.
"""
function _havel_hakimi(degs::Vector{Int}, store::Bool, ei::Vector{Int}, ej::Vector{Int})
    q = PriorityQueue{Int,Int}(Base.Order.Reverse)
    n = length(degs)
    effective_n = n
    degsum = 0
    dmax = 0
    for (i,d) in enumerate(degs)
        q[i] = d
        degsum += d
        dmax = max(d,dmax)
        if d < 0
            throw(ArgumentError("the degree sequence must be non-negative"))
        end
    end

    if mod(degsum,2) != 0; return false; end
    if n > 0 && dmax >= n; return false; end # n > 0 checks for the empty graph

    if store
        resize!(ei,degsum)
        resize!(ej,degsum)
    end

    dlist = Vector{Pair{Int,Int}}(undef,dmax)
    enum = 1

    while !isempty(q)
        vi,d = peek(q) # vi is the cur vertex, d is the cur deg
        dequeue!(q)    # remove it
        for n=1:d                  # make a list of each neighbor
            if isempty(q); return false; end
            dlist[n] = peek(q)
            dequeue!(q)
        end
        # now "add" an edge from vi->neighbor, and thus, decrease it's
        # degree when we re-add it
        for n=1:d
            if store
                ei[enum] = vi
                ej[enum] = dlist[n][1]
                enum += 1
                ei[enum] = dlist[n][1]
                ej[enum] = vi
                enum += 1
            end
            if dlist[n][2] < 1
                return false
            elseif dlist[n][2] == 1
                # don't both re-adding the vertex
            else
                q[dlist[n][1]] = dlist[n][2] - 1
            end
        end
    end
    return true
end

"""
`is_graphical_sequence`
=======================

Check whether or not a degree sequence is graphical,
which means that it is a valid degree sequence for
an undirected graph.

Note that this does not mean it is a valid degree
sequence for a connected undirected graph. So,
for instance,
`[1,1,1,1]`
is a valid degree sequence for two
disconnected edges

Usage
-----
`is_graphical_sequence(d)` returns true or false

Input
-----
- `d::Vector{Int}`:  a vector of integer valued degrees

Output
------
- a boolean that is true if the sequence is graphical
"""
function is_graphical_sequence(d::Vector{Int})
    return _havel_hakimi(d, false, Int[],Int[])
end

"""
`havel_hakimi_graph`
====================

Create a graph with a given degree sequence

Usage
-----
`A = havel_hakimi_graph(d)` returns an instance of the
a graph with degree sequence d or throws ArgumentError
if the degree sequence is not graphical.

Input
-----
- `d::Vector{Int}`:  a vector of integer valued degrees

Output
------
-`A`: a matrix network for the undirected graph that
results from the Havel-Hakimi procedure.
"""
function havel_hakimi_graph(d::Vector{Int})
    ei = Int[]
    ej = Int[]
    if _havel_hakimi(d, true, ei, ej) == false
        throw(ArgumentError("the degree sequence is not graphical"))
    end
    return MatrixNetwork(ei,ej,length(d))
end

"""
`preferential_attachment_graph`
===============================

Generate an instance of a preferential attachment
graph. This is an undirected graph that is generated
as follows:

* Start with a k0-node clique.
* Add n-k0 vertices where
     each vertex links to k nodes chosen
     based their degree (and repeats
     are allowed).

Functions
---------
The following functions are synonyms

- `preferential_attachment_graph`
- `pa_graph`

and

- `preferential_attachment_edges!`
- `pa_edges!`

The computational functions are

- `pa_graph(n,k,k0)` Generate a PA graph with a k0 node clique
  and n total nodes and k edges added per node. This returns
  a MatrixNetwork type

The edge functions are

- `pa_edges!(nnew,k,edges)` Add new edges to an
  existing set by adding `nnew` nodes to the set of edges
  where each node picks k edges based on the degrees. The
  new node ids are based on the largest entry in the edges
  array.
- `pa_edges!(nnew,k,edges,n0)` Generate a set of edges total`
  nodes to the set of edges where n0+1 is the starting index
  for the new set of nodes

Input
-----
- `n`: the number of nodes in the final graph
- `k`: the number of links picked by each node when it is added.
  The actual degree can be larger or smaller than this number because
  of links from other nodes or duplicates selected.
- `k0`: The number of nodes in the starting clique.
- `edges`: A list of edges to be manipulated in the process of
 generating new edges.

Output
------
- A matrix network type for the preferential attachment graph.
- `edges` An updated list of edges.

Example
-------
~~~~
pa_graph(100,5,2)
~~~~
"""
:preferential_attachment_graph, :pa_graph, :pa_edges!, :preferential_attachment_edges!

function preferential_attachment_graph(n::Int,k::Int,k0::Int)
    #n >= 0 || throw(ArgumentError(@sprintf("n=%i must be non-negative",n)))
    #k >= 1 || throw(ArgumentError(@sprintf("k=%i must be strictly positive",k0)))
    k0 >= 0 || throw(ArgumentError(@sprintf("k0=%i must be non-negative",n)))
    n >= k0 || throw(ArgumentError(@sprintf("n=%i must be >= k0=%i",n, k0)))
    #k >= 0 || throw(ArgumentError(@sprintf("k=%i must be non-negative",k)))
    edges = Vector{Tuple{Int,Int}}()
    # add the clique
    for i=1:k0
        for j=1:i-1
            push!(edges, (i,j))
            push!(edges, (j,i))
        end
    end
    return MatrixNetwork(preferential_attachment_edges!(n-k0,k,edges,k0),n)
end

function preferential_attachment_edges!(nnew::Int,k::Int,edges::Vector{Tuple{Int,Int}})
    if length(edges) == 0
        throw(ArgumentError("the list of initial edges must be non-empty"))
    end
    n0 = max(edges[1]...)
    for j=2:length(edges)
        n0 = max(n0,edges[j]...)
    end
    return preferential_attachment_edges!(nnew,k,edges,n0)
end

function preferential_attachment_edges!(
            nnew::Int,k::Int,edges::Vector{Tuple{Int,Int}},n0::Int)
    for iter=1:nnew
        i = n0+iter
        newedges = unique([rand(edges)[1] for j=1:k])
        for v in newedges
            push!(edges, (i, v))
            push!(edges, (v, i))
        end
    end
    return edges
end

pa_graph = preferential_attachment_graph
pa_edges! = preferential_attachment_edges!

"""
'generalized_preferential_attachment_graph'
===========================================

Generate an instance of a generalized preferential attachment graph which
follows the Avin,Lotker,Nahum,Peleg description. This is an undirected graph
that is generated as follows:

- Start with a k0-node clique
- Add n - k0 vertices where at each time step one of three events occurs: A new
node is added with probability p, a new edge between two existing nodes is added
with probability r, two new nodes with an edge between them is added with
probability 1 - p - r

Functions
---------
The following functions are synonyms

- 'generalized_preferential_attachment_graph'
- 'gpa_graph'

and

- 'generalized_preferential_attachment_edges!'
- 'gpa_edges!'

The computational functions are

- 'gpa_graph(n,p,r,k0)' Generate a GPA graph with a k0 clique and n total nodes.
    This returns a MatrixNetwork type
- 'gpa_graph(n,p,r,k0,Val{true})' Generate a GPA graph with a k0 clique and
    n total nodes, allowing self-loops. This returns a MatrixNetwork type

The edge functions are

-   'gpa_edges!(n,p,r,edges,n0)' Add new edges to an existing set, by taking
    n0 time steps. Edges are added in one of three ways: From a new node to
    an existing node with probability p, between two existing nodes with
    probability r, between two new nodes with probability 1-p-r
-   'gpa_edges!(n,p,r,edges,n0,Val{true})' Add new edges to an existing set, by
    taking n0 time steps. Edges are added in one of three ways: From a new node
    to an existing node with probability p, between two existing nodes with
    probability r (allowing self-loops), between two new nodes with probability
    1-p-r

Input
-----
- 'n': the number of nodes in the final graph.
- 'p': The probability of a node event, p must be a constant.
- 'r': The probability of an edge event, r must be a constant. p+r <=1
- 'k0': the number of nodes in the starting clique.
- 'Val{true}': Include this parameter if self-loops are allowed. Default is false
- 'edges': A list of edges to be manipulated in the process of generating
          new edges.

Output
------
- A matrix network type for the generalized preferential attachment graph.
- 'edges': An updated list of edges.

Example:
generalized_preferential_attachment_graph(100,1/3,1/2,2)

"""
:generalized_preferential_attachment_graph, :gpa_graph,
:generalized_preferential_attachment_edges!, :gpa_edges!

generalized_preferential_attachment_graph(n::Int,p::Float64,r::Float64,k0::Int) =
    generalized_preferential_attachment_graph(n,p,r,k0,Val{false})
generalized_preferential_attachment_edges!(n::Int,p::Float64,r::Float64,edges::Vector{Tuple{Int,Int}},n0::Int) =
    generalized_preferential_attachment_edges!(n,p,r,edges,n0,Val{false})

function generalized_preferential_attachment_graph(
    n::Int,p::Float64,r::Float64,k0::Int,::Type{Val{true}})
    k0 >= 0 || throw(ArgumentError(@sprintf("k0=%i must be non-negative",k0)))
    n >= k0 || throw(ArgumentError(@sprintf("n=%i must be >= k0=%i",n,k0)))
    0<=p<=1 || throw(ArgumentError(@sprintf("p=%0.3f must be between 0 and 1",p)))
    0<=r<=1 || throw(ArgumentError(@sprintf("r=%0.3f must be between 0 and 1",r)))
    p+r <= 1 || throw(ArgumentError(@sprintf("(p=%0.3f)+(r=%0.3f) must be <=1",p,r)))
    edges = Vector{Tuple{Int,Int}}()
    #add the clique
    for i = 1:k0
        for j = 1:i-1
            push!(edges,(i,j))
            push!(edges, (j,i))
        end
    end
    return MatrixNetwork(generalized_preferential_attachment_edges!(n,p,r,edges,k0,Val{true}),n)
end

function generalized_preferential_attachment_graph(
    n::Int,p::Float64,r::Float64,k0::Int,::Type{Val{false}})
    k0 >= 0 || throw(ArgumentError(@sprintf("k0=%i must be non-negative",k0)))
    n >= k0 || throw(ArgumentError(@sprintf("n=%i must be >= k0=%i",n,k0)))
    0<=p<=1 || throw(ArgumentError(@sprintf("p=%0.3f must be between 0 and 1",p)))
    0<=r<=1 || throw(ArgumentError(@sprintf("r=%0.3f must be between 0 and 1",r)))
    p+r <= 1 || throw(ArgumentError(@sprintf("(p=%0.3f)+(r=%0.3f) must be <=1",p,r)))
    edges = Vector{Tuple{Int,Int}}()
    #add the clique
    for i = 1:k0
        for j = 1:i-1
            push!(edges,(i,j))
            push!(edges, (j,i))
        end
    end
    return MatrixNetwork(generalized_preferential_attachment_edges!(n,p,r,edges,k0,Val{false}),n)
end

function generalized_preferential_attachment_edges!(
    n::Int,p::Float64,r::Float64,edges::Vector{Tuple{Int,Int}},n0::Int,::Type{Val{true}})
    i = n0
    while i < n
        #generate a random value between 0 and 1
        x = rand()
        if x < p #node event
            v = rand(edges)[1]
            push!(edges, (i+1,v[1]))
            push!(edges, (v[1], i+1))
            i = i+1;
        elseif x < p+r #edge event, self loops permitted
            v1 = rand(edges)[1]
            v2 = rand(edges)[1]
            push!(edges, (v1, v2))
            push!(edges, (v2, v1))
        else #component event
            if i+2 <= n #only allow this step if there is room for two more nodes
                push!(edges, (i+1, i+2))
                push!(edges, (i+2, i+1))
                i = i+2;
            end
        end
    end
    return edges
end

function _check_for_two_distinct_nodes(edges::Vector{Tuple{Int,Int}})
    length(edges) > 0 || throw(ArgumentError("requires at least one edge"))
    firstnode = edges[1][1]
    return any(IterTools.imap(x -> firstnode != x[1] || firstnode != x[2], edges))
end

function generalized_preferential_attachment_edges!(
    n::Int,p::Float64,r::Float64,edges::Vector{Tuple{Int,Int}},n0::Int,::Type{Val{false}})
    i = n0
    if i >= n
        return edges
    end

    if !_check_for_two_distinct_nodes(edges::Vector{Tuple{Int,Int}})
        throw(ArgumentError("The starting graph must have at least two distinct nodes"))
    end

    while i < n
        #generate a random value between 0 and 1
        x = rand()
        if x < p #node event
            v = rand(edges)[1]
            push!(edges, (i+1,v[1]))
            push!(edges, (v[1], i+1))
            i = i+1;
        elseif x < p+r #edge event, no self-loops permitted
            v1 = rand(edges)[1]
            v2 = rand(edges)[1]
            while (v1 == v2 && i != 1) #i != 1 because we want more than 1 node for this to work
                v1 = rand(edges)[1]
                v2 = rand(edges)[1]
            end
            push!(edges, (v1, v2))
            push!(edges, (v2, v1))
        else #component event
            if i+2 <= n #only allow this step if there is room for two more nodes
                push!(edges, (i+1, i+2))
                push!(edges, (i+2, i+1))
                i = i+2;
            end
        end
    end
    return edges
end

gpa_graph = generalized_preferential_attachment_graph
gpa_edges! = generalized_preferential_attachment_edges!

"""
`roach_graph`
=============

Generate a roach graph on 4n vertices which follows the Guattery-Miller
description. The roach graph has a body that consists of 2n vertices which
ard two n-vertex line-graphs that have been connected together. The
body has two antennae that result from adding an n-vertex line graph
to one vertex on each side.

    # the graph looks like
    #
    # (           top body               )   (       top antennae       )
    #
    # o - o - o - ... n vertices total - o - o - ... n vertices total - 0
    # |   |   |   ...    |   |   |   |   |
    # o - o - o - ... n vertices total - o - o - ... n vertices total - 0
    #
    # (         bottom body              )   (     bottom antennae      )
    #
    # there are 4n vertices and 2(2n-1) + n edges

Functions
---------
  * `roach_graph(n) -> A::MatrixNetwork`
  * `roach_graph(n, Val{true}) -> (A::MatrixNetwork, Matrix{Float64})` this also
    returns coordinates for the graph.

Example
-------
~~~~
A = sparse_transpose(roach_graph(3, Val{true})) # get back the matrix
L = spdiagm(vec(sum(A,2))) - A

#lams,vecs = eig(full(L))

~~~~

"""
function roach_graph end

roach_graph(n::Integer) = roach_graph(n::Integer, Val{false})

function roach_graph(n::Integer, ::Type{Val{false}})
    n >= 0 || throw(ArgumentError("n=$(n) must be larger than 0"))
    line1 = 1:(2n-1)
    line2 = 2:2n
    ei = [line1; line2; line1.+2n; line2.+2n; 1:n; 2n+1:2n+n]
    ej = [line2; line1; line2.+2n; line1.+2n; 2n+1:2n+n; 1:n]
    return _matrix_network_direct(sparse(ei,ej,1,4n,4n))
end

# get coordinates too
function roach_graph(n::Integer, ::Type{Val{true}})
    A = roach_graph(n, Val{false})
    xy = [-n:-1 -ones(n); 1:n -ones(n); -n:-1 ones(n); 1:n ones(n)]
    return A,xy
end

"""
`lollipop_graph`
================

Generate a lollipop graph, which consists of a clique with a line tail, so
it looks like a lollipop.

Functions
---------
* `lollipop_graph(n)` generate the graph with an n-node tail and n-node clique
* `lollipop_graph(n,m)` generate the graph with an n-node tail and m-node clique
* `lollipop_graph(n,m,Val{true})` produce and return xy coordinates as well.

Examples
--------

"""
function lollipop_graph end
lollipop_graph(n::Integer) = lollipop_graph(n, n)
lollipop_graph(n::Integer, ::Type{Val{false}}) = lollipop_graph(n, n)
lollipop_graph(n::Integer, ::Type{Val{true}}) = lollipop_graph(n, n, Val{true})
lollipop_graph(n::Integer, m::Integer) = lollipop_graph(n, m, Val{false})
function lollipop_graph(n::Integer, m::Integer, ::Type{Val{false}})
    n >= 0 || throw(ArgumentError("n=$(n) must be larger than 0"))
    m >= 0 || throw(ArgumentError("m=$(m) must be larger than 0"))
    line1 = 1:n
    line2 = 2:n+1
    clique1 = map(x -> x[1], IterTools.subsets(n+1:n+m,2))
    clique2 = map(x -> x[2], IterTools.subsets(n+1:n+m,2))
    ei = [line1; line2; clique1; clique2]
    ej = [line2; line1; clique2; clique1]
    return _matrix_network_direct(sparse(ei,ej,1,n+m,n+m))
end
function lollipop_graph(n::Integer, m::Integer, ::Type{Val{true}})
    A = lollipop_graph(n,m, Val{false})
    xy = [-n:-1 zeros(n);
        (-sqrt(m)*cos.(2*pi*(m:-1:1)/m).+(sqrt(m)+1)) sqrt(m)*sin.(2*pi*(m:-1:1)/m)]
    return A, xy
end

function rewire_graph(A::MatrixNetwork, k::Integer)
    if is_undirected(A)
        ei,ej = undirected_edges(A)
    else
        ei,ej = directed_edges(A)
    end
end

#=
function random_symmetric_edge(A::MatrixNetwork)
    ei,ej,ind1 = random_edge(A)
    ind2 = searchsortedfirst(@view A.ci[A.rp[ej]:A.rp[ej+1]-1], ei)
end
random_undirected_edge(A::MatrixNetwork)

function reverse_edge_index(A::MatrixNetwork, ei::Integer, ej::Integer)

end
=#

# TODO Add chung-lu for general floating point weights
# via union of ER graphs add

"""
`forest_fire_graph`
===================

Create an instance of a forest fire graph. The forest fire
model consists of ...

Usage
-----
`A = forest_fire_graph(c, k, p)` returns a forest-fire sample.
`A = forest_fire_graph(A0, k, p)`
"""


text:
using Plots, Printf, LinearAlgebra

# enable plotting by default
if !@isdefined do_visu; do_visu = true end

# finite-difference support functions
@views av_xi(A) = 0.5*(A[1:end-1,2:end-1].+A[2:end,2:end-1])
@views av_yi(A) = 0.5*(A[2:end-1,1:end-1].+A[2:end-1,2:end])
@views   inn(A) = A[2:end-1,2:end-1]

@views function diffusion_2D_impl(; do_visu=true, save_fig=false)
    # Physics
    lx, ly = 10.0, 10.0   # domain size
    npow   = 3            # power-law exponent
    ttot   = 1.0          # total simulation time
    dt     = 0.2          # physical time step
    # Numerics
    nx, ny = 128, 128     # numerical grid resolution
    tol    = 1e-6         # tolerance
    itMax  = 1e5          # max number of iterations
    # Derived numerics
    dx, dy = lx/nx, ly/ny # grid size
    xc, yc = LinRange(dx/2, lx-dx/2, nx), LinRange(dy/2, ly-dy/2, ny)
    # Array allocation
    qHx    = zeros(nx-1, ny-2) # on staggered grid
    qHy    = zeros(nx-2, ny-1) # on staggered grid
    ResH   = zeros(nx-2, ny-2) # normal grid, without boundary points
    dHdtau = zeros(nx-2, ny-2) # normal grid, without boundary points
    dtau   = zeros(nx-2, ny-2) # normal grid, without boundary points
    # Initial condition
    H      = exp.(.-(xc.-lx/2).^2 .-(yc.-ly/2)'.^2)
    Hold   = copy(H)
    t = 0.0; it = 0; ittot = 0
    # Physical time loop
    while t<ttot
        iter = 0; err = 2*tol
        # Picard-type iteration
        while err>tol && iter<itMax
            qHx    .= -av_xi(H).^npow.*diff(H[:,2:end-1], dims=1)/dx  # flux
            qHy    .= -av_yi(H).^npow.*diff(H[2:end-1,:], dims=2)/dy  # flux
            ResH   .= -(inn(H) - inn(Hold))/dt +
                       (-diff(qHx, dims=1)/dx -diff(qHy, dims=2)/dy)  # residual of the PDE
            dHdtau .= ResH                                            # rate of change
            dtau   .= (1.0./(min(dx, dy)^2 ./inn(H).^npow./4.1) .+ 1.0/dt).^-1  # time step (obeys ~CFL condition)
            H[2:end-1,2:end-1] .= inn(H) .+ dtau.*dHdtau              # update rule, sets the BC as H[1]=H[end]=0
            iter += 1; err = norm(ResH)/length(ResH)
        end
        ittot += iter; it += 1; t += dt
        Hold .= H
    end
    @printf("Total time = %1.2f, time steps = %d, iterations tot = %d \n", round(ttot, sigdigits=2), it, ittot)
    # Visualize
    if do_visu
        fontsize = 12
        opts = (aspect_ratio=1, yaxis=font(fontsize, "Courier"), xaxis=font(fontsize, "Courier"),
                ticks=nothing, framestyle=:box, titlefontsize=fontsize, titlefont="Courier", colorbar_title="",
                xlabel="Lx", ylabel="Ly", xlims=(xc[1],xc[end]), ylims=(yc[1],yc[end]), clims=(0.,1.))
        display(heatmap(xc, yc, H'; c=:davos, title="implicit diffusion (nt=$it, iters=$ittot)", opts...))
        if save_fig savefig("diff2D_impl.png") end
    end
    return xc, yc, H
end

diffusion_2D_impl(; do_visu=do_visu);


text:
ROOT = pwd()
include(joinpath(ROOT, "src_PowSysMod", "PowSysMod_body.jl"))

folder = "Phase_0_IEEE14_1Scenario"
scenario = "scenario_1"
filenames = [folder, scenario]
filenames_path = joinpath("data_GOC", folder, scenario)
instance_path = filenames_path

## Build global pb
OPFpbs = load_OPFproblems(GOCInput, instance_path)
introduce_Sgenvariables!(OPFpbs)
pb_global = build_globalpb!(OPFpbs)

## Build solution
global_point = solution_point(instance_path)
println("minslack of GOC global solution: $(get_minslack(pb_global, global_point))")

## Export problem
pb_global_real = pb_cplx2real(pb_global)
pt_global_real = cplx2real(global_point)

# output_file = joinpath("toto", "real_minlp_instance.dat")
export_to_dat(pb_global_real, "src_ampl", pt_global_real)
# pb_global_real = pb_cplx2real(pb_global)
# pt_global_real = cplx2real(global_point)
#
# output_file = joinpath(amplexportpath, "real_minlp_instance.dat")
# export_to_dat(pb_global_real, output_file, pt_global_real)

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
function bgk!(du, u, p, t)
    dx, e2f, f2e, velo, weights, δ, deg, ll, lr, lpdm, dgl, dgr = p

    ncell = size(u, 1)
    nu = size(u, 2)
    nsp = size(u, 3)

    M = similar(u, ncell, nu, nsp)
    for i = 1:ncell, k = 1:nsp
        w = moments_conserve(u[i, :, k], velo, weights)
        prim = conserve_prim(w, 3.0)
        M[i, :, k] .= maxwellian(velo, prim)
    end
    τ = 1e-2

    f = similar(u)
    for i = 1:ncell, j = 1:nu, k = 1:nsp
        J = 0.5 * dx[i]
        f[i, j, k] = velo[j] * u[i, j, k] / J
    end

    f_face = zeros(eltype(u), ncell, nu, 2)
    #=for i = 1:ncell, j = 1:nu, k = 1:nsp
        # right face of element i
        f_face[i, j, 1] += f[i, j, k] * lr[k]

        # left face of element i
        f_face[i, j, 2] += f[i, j, k] * ll[k]
    end=#

    @views for i = 1:ncell, j = 1:nu
        FluxRC.interp_interface!(f_face[i, j, :], f[i, j, :], ll, lr)
    end

    f_interaction = similar(u, nface, nu)
    for i = 1:nface
        @. f_interaction[i, :] =
            f_face[f2e[i, 1], :, 1] * (1.0 - δ) + f_face[f2e[i, 2], :, 2] * δ
    end

    rhs1 = zeros(eltype(u), ncell, nu, nsp)
    #for i = 1:ncell, j = 1:nu, ppp1 = 1:nsp, k = 1:nsp
    #   rhs1[i, j, ppp1] += f[i, j, k] * lpdm[ppp1, k]
    #end

    @views for i = 1:ncell, j = 1:nu
        FluxRC.poly_derivative!(rhs1[i, j, :], f[i, j, :], lpdm)
    end

    #@views for i = 1:ncell, j = 1:nu, k = 1:nsp
    #    rhs1[i, j, k] = dot(f[i, j, :], lpdm[k, :])
    #end

    for i = 1:ncell, j = 1:nu, ppp1 = 1:nsp
        du[i, j, ppp1] =
            -(
                rhs1[i, j, ppp1] +
                (f_interaction[e2f[i, 2], j] - f_face[i, j, 1]) * dgl[ppp1] +
                (f_interaction[e2f[i, 1], j] - f_face[i, j, 2]) * dgr[ppp1]
            ) + (M[i, j, ppp1] - u[i, j, ppp1]) / τ
    end
end


text:
using Mmap
using DataFrames
using Statistics: mean, std

export Marker,
    Sample,
    PlinkReader,
    samples,
    markers,
    markersDF,
    nsamples,
    nmarkers,
    sample_index,
    marker_index,
    getindex,
    applyGeneticMap!,
    dosageMatrix

"""
    Marker(chrom, id, cm, pos, a1, a2)

PLINK marker info
"""
struct Marker
    chrom::String
    id::String
    cm::Float64
    pos::Int64
    a1::String
    a2::String
end

"""
    load_bim(path)

Load marker info from PLINK .bim file, adding 'chr' prefix if missing.
Return `Vector{Marker}`.
"""
function load_bim(path)

    markers = Marker[]
    bimpath = endswith(path, ".bim") ? path : path * ".bim"
    fh = open(bimpath, "r")
    for line in eachline(fh)
        cols = split(chomp(line))
        push!(markers, Marker(startswith(cols[1], "chr") ? cols[1] : "chr"*cols[1],
                              cols[2],
                              parse(Float64, cols[3]),
                              parse(Int64, cols[4]),
                              cols[5],
                              cols[6]))
    end
    close(fh)

    markers
end

"""
    Sample(fid, iid, father, mother, sex)

PLINK sample info
"""
struct Sample
    fid::String
    iid::String
    father::String
    mother::String
    sex::UInt8
end

"""
    load_fam(path)

Load sample info from PLINK .fam file. Return `Vector{Sample}`.
"""
function load_fam(path)

    samples = Sample[]
    fampath = endswith(path, ".fam") ? path : path * ".fam"
    fh = open(fampath, "r")
    for line in eachline(fh)
        cols = split(chomp(line))
        push!(samples, Sample(cols[1],
                              cols[2],
                              cols[3],
                              cols[4],
                              parse(UInt8, cols[5])))
    end
    close(fh)

    samples
end

"""
PLINK files (.bim, .fam, .bed)
"""
mutable struct PlinkReader
    path::String
    nmarkers::Int64
    nsamples::Int64
    markers::Vector{Marker}
    markerIndex::Dict{String, Int64}
    samples::Vector{Sample}
    sampleIndex::Dict{String, Int64}
    data::Array{UInt8, 2}
end

"""
    PlinkReader(path; markerIndex = false, sampleIndex = false)

Structure to access PLINK .bed, .bim, and .fam files at `path`.

# Arguments

- `markerIndex::Bool` should index marker -> idx be created?
- `sampleIndex::Bool` should index sample -> idx be created?
"""
function PlinkReader(path; markerIndex = false, sampleIndex = false)

    samples = load_fam(path)
    markers = load_bim(path)

    nsamples = length(samples)
    nmarkers = length(markers)

    fh = open(path * ".bed", "r")
    bedheader = read(fh, 3)

    @assert bedheader[1] == 0b01101100 && bedheader[2] == 0b00011011 && bedheader[3] == 0b01 "Only snp major BED v1.0 supported"

    #plinkdata = Mmap.mmap(fh, BitArray, (2, 4*ceil(Int64, 0.25 * nsamples), nmarkers))
    # see source code for Mmap.mmap(io, BitArray, dims)
    nrows = (nsamples + 3) >> 2
    data = Mmap.mmap(fh, Vector{UInt8}, nrows * nmarkers)

    PlinkReader(path,
                nmarkers,
                nsamples,
                markers,
                markerIndex ? Dict([m.id => idx for (idx, m) in enumerate(markers)]) : Dict(),
                samples,
                sampleIndex ? Dict([s.iid => idx for (idx, s) in enumerate(samples)]) : Dict(),
                reshape(data, (nrows, nmarkers)))
end

function Base.show(io::IO, p::PlinkReader)
    print(io, "<PLINK file ($(p.nsamples) samples x $(p.nmarkers) markers) at $(p.path)>")
end

"""
    samples(p::PlinkReader)

Return Array of `Sample` structs for samples in Plink file.
"""
samples(p::PlinkReader) = p.samples

"""
    markers(p::PlinkReader)

Return array of `Marker` structs for markers in Plink file.
"""
markers(p::PlinkReader) = p.markers

"""
    nsamples(p::PlinkReader)

Number of samples in Plink file.
"""
nsamples(p::PlinkReader) = p.nsamples

"""
    nmarkers(p::PlinkReader)

Number of markers in Plink file.
"""
nmarkers(p::PlinkReader) = p.nmarkers

"""
    markersDF(p::PlinkReader)

Get data frame with marker info.

Returns data frame with columns

- Chrom
- Name
- cM
- Pos
- A1
- A2
- Idx
"""
function markersDF(p::PlinkReader)
    m = DataFrame(Chrom = String[], Name = String[], cM = Float64[],
                  Pos = Int64[], A1 = String[], A2 = String[],
                  Idx = Int64[])

    for (idx, mrk) in enumerate(markers(p))
        push!(m, (mrk.chrom, mrk.id, mrk.cm, mrk.pos, mrk.a1, mrk.a2, idx))
    end
    m
end

"""
    sample_index(p::PlinkReader, iid::String)

Return index of sample `iid` in `PlinkReader`.
Uses index if available, otherwise linear search.
"""
function sample_index(p::PlinkReader, iid::String)

    if length(p.sampleIndex) > 0
        return p.sampleIndex[iid]
    end

    s = samples(p)
    for i = 1:length(s)
        if s[i].iid == iid
            return i
        end
    end
    error("No such sample: $iid")
end

"""
    marker_index(p::PlinkReader, id::String)

Return index of marker `id` in `PlinkReader`.
Uses index if available, otherwise linear search.
"""
function marker_index(p::PlinkReader, id::String)

    if length(p.markerIndex) > 0
        return p.markerIndex[id]
    end

    m = markers(p)
    for i = 1:length(m)
        if m[i].id == id
            return i
        end
    end
    error("No such marker: $iid")
end

"""
    getindex(p::PlinkReader, s::Int, m::Int)

Retrieve genotype of sample `s` and marker `m`.

Plink convention for representation of genotypes is as follows:

- `0b00` Hom1
- `0b01` Het
- `0b10` missing
- `0b11` Hom2
"""
@inline function Base.getindex(p::PlinkReader, s::Int, m::Int)
    sp3 = s + 3
    (p.data[sp3 >> 2, m] >> ((sp3 & 0x03) << 1)) & 0x03
end

Base.size(p::PlinkReader) = (p.nsamples, p.nmarkers)

abstract type AbstractGeneticMap end

function applyGeneticMap!(gm::AbstractGeneticMap, p::PlinkReader)

    function applyGeneticMap(gm::GeneticMap, m::Marker)
        Marker(m.chrom,
               m.id,
               gm(m.chrom, m.pos),
               m.pos,
               m.a1,
               m.a2)
    end

    p.markers = [applyGeneticMap(gm, m) for m in p.markers]
    p
end

"""
    dosageMatrix(p::PlinkReader, markerIdx, sampleIdx = nothing;
                   normalize = true)

Create dosage matrix (samples x markers) for `markerIdx` and
`sampleIdx` from `PlinkReader`. Here dosage is expressed as count of
alternative allele.

Replace missing genotypes with mean dosage. If `normalize = true`,
normalize markers to mean μ = 0 and standard deviation σ = 1.
"""
function dosageMatrix(p::PlinkReader, markerIdx, sampleIdx = nothing;
                        normalize = true)

    sampleIdx == nothing && (sampleIdx = 1:nsamples(p))

    nm = length(markerIdx)
    ns = length(sampleIdx)
    g = zeros(Float64, ns, nm)

    # a2counts[gt + 1] = counts of allele2 for gt
    # Plink convention is
    # 0b00 Hom1, 0b01 Het, 0b10 missing, 0b11 Hom2
    a2counts = [0, 1, -1, 2]

    for (midx, m) in enumerate(markerIdx)
        data = view(p.data,:, m)
        counts = 0
        total = 0
        miss = 0
        for (sidx, s) in enumerate(sampleIdx)
            sp3 = s + 3
            plgt = data[sp3 >> 2] >> ((sp3 & 0x03) << 1) & 0x03

            a2count = a2counts[plgt + 1]
            g[sidx, midx] = a2count

            if a2count == -1
                miss += 1
            else
                counts += a2count
                total += 1
            end
        end

        # replace missing genotype with mean
        if miss > 0
            freq = counts/float(total)
            for s = 1:ns
                if g[s, midx] < 0.0
                    g[s, midx] = freq
                end
            end
        end

    end

    if normalize
        sd = std(g; dims = 1)
        sd[sd .== 0.0] .= 1.0 # avoid division by zero
        g ./= sd
        g .-= mean(g; dims = 1)
    end

    g

end


text:
module Moments

export make_mon_expo,
       get_ℝ_block_diag,
       make_xxᵀ_tens_yyᵀ,
       get_ℂ_block_diag,
       make_xx̄ᵀ_tens_yȳᵀ,
       var_kron,
       var_kron_C,
       eᵢ,
       get_ℜℑααᶥββᶥᴿ,
       get_xx̄yȳMM_blocks,
       get_γγᶥζζᶥ_δδᶥηηᶥ

## Utils
"""The standard basis vector eᵢ in dimension n"""
eᵢ(n::Int,i::Int) = [Int(j==i) for j in 1:n]

""" input: A,B (arrays of integer tupples) output: A ⊗ B """
function var_kron(A,B)
   C_temp = [ a + b for  b in B, a ∈ A]
   rcs =  [(i,j) for i in 1:size(A)[1] for j in 1:size(B)[1]]
   return [C_temp[ij[2],kl[2],ij[1],kl[1]] for ij in rcs , kl in rcs ]
end

function var_kron_C(A,B)
   C_temp = [ [a] + b for  b in B, a ∈ A]
   rcs =  [(i,j) for i in 1:size(A)[1] for j in 1:size(B)[1]]
   return [C_temp[ij[2],kl[2],ij[1],kl[1]] for ij in rcs , kl in rcs ]
end

## Moments
"""[x]≦ₜ or [x]₌ₜ"""
function make_mon_expo(n::Int,t::Int; isle::Bool = true)
    @assert typeof(t) == Int64
    t == 0 ? (return [eᵢ(n,0)]) : 0
    tmp = make_mon_expo(n,t-1;isle=isle)
    M_vec = reshape([m + eᵢ(n,i) for i ∈ 1:n, m ∈ tmp],:,1)
    return unique(isle ? vcat(tmp,M_vec) : M_vec)
end

"""[x]≦ₜ[x]ᵀ≦ₜ or [x]₌ₜ[x]ᵀ₌ₜ"""
function make_mon_expo(n::Int,t::Tuple{Int,Int}; isle::Bool = true)
    M_vec1      = make_mon_expo(n,t[1]; isle=isle)
    M_vec2      = make_mon_expo(n,t[2]; isle=isle)
    M_mat_flat  = [ mi+ mj for mi in M_vec1 for mj in M_vec2]
    M_mat       = reshape(M_mat_flat, (length(M_vec1), length(M_vec2)) )
    return M_mat
end

## Complex moments
"""[x,̄x]≦ₜ or [x,̄x]ᵀ₌ₜ"""
make_mon_expo(d::Tuple{Int,Int},t::Int; isle::Bool = true) = make_mon_expo(sum(2 .* d),t; isle=isle)

"""[x,̄x]≦ₜ[x,̄x]*≦ₜ or [x,̄x]₌ₜ[x,̄x]*₌ₜ
but is actually [xᵣₑ,xᵢₘ]≦ₜ[xᵣₑ,xᵢₘ]ᵀ≦ₜ or [xᵣₑ,xᵢₘ]₌ₜ[xᵣₑ,xᵢₘ]ᵀ₌ₜ
"""
function make_mon_expo(d::Tuple{Int,Int},t::Tuple{Int,Int}; isle::Bool = true)
    M_vec1      = make_mon_expo(d,t[1]; isle=isle)
    M_vec2      = make_mon_expo(d,t[1]; isle=isle)
    M_mat      = [ mi+ mj for mi in M_vec1, mj in M_vec2]
    # M_mat       = reshape(M_mat_flat, (length(M_vec1), length(M_vec2)) )
    return M_mat
end
## Real
"""xxᵀ⊗yyᵀ"""
function make_xxᵀ_tens_yyᵀ(d)
    n = sum(d)
    pre_expo = make_mon_expo(n,1;isle = false)
    x, y = pre_expo[1:d[1]], pre_expo[d[1]+1:n]
    return var_kron(x .+ reshape(x,1,d[1]), y .+ reshape(y,1,d[2]))
end

"""Returns diagonal block of the momoment matrix with partition:
    xᵃ⁺ᶜyᵇ⁺ᵈ where |a+c|,|b+d| ∈ 2N """
function get_ℝ_block_diag(d::Tuple{Int,Int},t::Tuple{Int,Int})
    mom_mat = Moments.make_mon_expo(sum(d),t)
    # return Dict("Default" => mom_mat)
    mom_vec = mom_mat[1,:]
    halfsum(arr) = [sum(arr[1:d[1]]),sum(arr[(d[1]+1):end])]
    tf(arr)      = map(x -> (isodd(x[1]),isodd(x[2])), halfsum.(arr))
    tf2(arr)     = map(p -> p[1], findall(arr))

    isevev(pair) = pair == (false,false)
    isodod(pair) = pair == (true,true)
    isevod(pair) = pair == (false,true)
    isodev(pair) = pair == (true,false)

    ee = tf2(isevev.(tf(mom_vec)))
    oo = tf2(isodod.(tf(mom_vec)))
    eo = tf2(isevod.(tf(mom_vec)))
    oe = tf2(isodev.(tf(mom_vec)))

    ℝ_block_diag = Dict("ee" => mom_mat[ee,ee],
                        "oo" => mom_mat[oo,oo],
                        "eo" => mom_mat[eo,eo],
                        "oe" => mom_mat[oe,oe])

    for key in ["ee","oo","eo","oe"]
         isempty(ℝ_block_diag[key]) ? delete!(ℝ_block_diag,key) : continue
    end
    return ℝ_block_diag
end

## Complex

""" ℜe(xx*⊗yy*), ℑm(xx*⊗yy*)
xx*⊗yy* = (xᵣₑ + i xᵢₘ)(xᵣₑ - i xᵢₘ)ᵀ ⊗ (yᵣₑ + i yᵢₘ)(yᵣₑ - i yᵢₘ)ᵀ
= (xᵣₑxᵣₑᵀ + xᵢₘxᵢₘᵀ + i (xᵢₘxᵣₑᵀ - xᵣₑxᵢₘᵀ)) ⊗ (yᵣₑyᵣₑᵀ + yᵢₘyᵢₘᵀ + i (yᵢₘyᵣₑᵀ - yᵣₑyᵢₘᵀ))
⟹ ℜe(xx*⊗yy*) = (xᵣₑxᵣₑᵀ + xᵢₘxᵢₘᵀ ) ⊗ (yᵣₑyᵣₑᵀ + yᵢₘyᵢₘᵀ ) - (xᵢₘxᵣₑᵀ - xᵣₑxᵢₘᵀ) ⊗ (yᵢₘyᵣₑᵀ - yᵣₑyᵢₘᵀ)
⟹ ℑm(xx*⊗yy*) = (xᵣₑxᵣₑᵀ + xᵢₘxᵢₘᵀ) ⊗ (yᵢₘyᵣₑᵀ - yᵣₑyᵢₘᵀ) + (xᵢₘxᵣₑᵀ - xᵣₑxᵢₘᵀ)⊗(yᵣₑyᵣₑᵀ + yᵢₘyᵢₘᵀ) """
U(vec,dₐ) = vec .+ reshape(vec,1,dₐ)
W(vec1,vec2,dₐ) = vec1 .+ reshape(vec2,1,dₐ)
function make_xx̄ᵀ_tens_yȳᵀ(d)
    d₁,d₂ = d
    n = sum(2 .* d)
    pre_expo = make_mon_expo(n,1;isle =false)
    xᵣₑ, xᵢₘ = pre_expo[1:d₁], pre_expo[d₁+1:2*d₁]
    yᵣₑ, yᵢₘ = pre_expo[2*d₁+1:2*d₁+d₂], pre_expo[2*d₁+d₂+1:2*d₁+2*d₂]

    xᵣₑxᵣₑᵀ, xᵢₘxᵢₘᵀ, yᵣₑyᵣₑᵀ, yᵢₘyᵢₘᵀ = U(xᵣₑ,d₁), U(xᵢₘ,d₁), U(yᵣₑ,d₂), U(yᵢₘ,d₂)
    xᵢₘxᵣₑᵀ, xᵣₑxᵢₘᵀ, yᵢₘyᵣₑᵀ, yᵣₑyᵢₘᵀ = W(xᵢₘ,xᵣₑ,d₁), W(xᵣₑ,xᵢₘ,d₁), W(yᵢₘ,yᵣₑ,d₂), W(yᵣₑ,yᵢₘ,d₂)

    Real_dict =  Dict(("+ 1")  => var_kron(xᵣₑxᵣₑᵀ,yᵣₑyᵣₑᵀ),
                      ("+ 2")  => var_kron(xᵣₑxᵣₑᵀ,yᵢₘyᵢₘᵀ),
                      ("+ 3")  => var_kron(xᵢₘxᵢₘᵀ,yᵣₑyᵣₑᵀ),
                      ("+ 4")  => var_kron(xᵢₘxᵢₘᵀ,yᵢₘyᵢₘᵀ),
                      ("- 5")  => var_kron(xᵢₘxᵣₑᵀ,yᵢₘyᵣₑᵀ),
                      ("+ 6")  => var_kron(xᵢₘxᵣₑᵀ,yᵣₑyᵢₘᵀ),
                      ("+ 7")  => var_kron(xᵣₑxᵢₘᵀ,yᵢₘyᵣₑᵀ),
                      ("- 8")  => var_kron(xᵣₑxᵢₘᵀ,yᵣₑyᵢₘᵀ))

    Imag_dict =  Dict(("+ 1")  => var_kron(xᵣₑxᵣₑᵀ,yᵢₘyᵣₑᵀ),
                      ("- 2") => var_kron(xᵣₑxᵣₑᵀ,yᵣₑyᵢₘᵀ),
                      ("+ 3")  => var_kron(xᵢₘxᵢₘᵀ,yᵢₘyᵣₑᵀ),
                      ("- 4") => var_kron(xᵢₘxᵢₘᵀ,yᵣₑyᵢₘᵀ),
                      ("+ 5")  => var_kron(xᵢₘxᵣₑᵀ,yᵣₑyᵣₑᵀ),
                      ("+ 6")  => var_kron(xᵢₘxᵣₑᵀ,yᵢₘyᵢₘᵀ),
                      ("- 7") => var_kron(xᵣₑxᵢₘᵀ,yᵣₑyᵣₑᵀ),
                      ("- 8") => var_kron(xᵣₑxᵢₘᵀ,yᵢₘyᵢₘᵀ))

    return Dict("real" => Real_dict,
                "imag" => Imag_dict)
end

## Block diagonalization
# The complex blocks
"""Splits a vector of exponents into the register components ααᶥββᶥ →  (α,αᶥ,β,βᶥ)"""
split_expo(ααᶥββᶥ,d) =     (ααᶥββᶥ[1:d[1]],
                            ααᶥββᶥ[d[1]+1:2*d[1]],
                            ααᶥββᶥ[1+2*d[1]:2*d[1]+d[2]],
                            ααᶥββᶥ[1+2*d[1]+d[2]:end])

"""Gets the principal block of the moment matrix after block diagonalization"""
function get_xx̄yȳMM_blocks(d,t)
    Iᵗ = Moments.make_mon_expo(d,t[1])
    Iᵗ_temp = map(x ->  sum.(split_expo(x,d)),Iᵗ)
    r_list  = map(v -> v[1]-v[2],Iᵗ_temp) ; s_list  = map(v -> v[3]-v[4],Iᵗ_temp)

    Iᵗ_d = Dict()
    for r ∈ -t[1]:t[1], s ∈ -t[1]:t[1]
        T = Iᵗ[(r_list .== r) .& (s_list .== s) ]
        isempty(T) ? continue : Iᵗ_d[r,s] = T
    end
    return Dict(zip(keys(Iᵗ_d),[Iᵗ_d[k] .+ reshape(Iᵗ_d[k],1,:) for k in keys(Iᵗ_d)]))
end

## The real analogue
"""Returns all  (γ,γᶥ,ζ,ζᶥ),(δ,δᶥ,η,ηᶥ) ∈ (ℕᵈ)⁴ s.t.
(γ,γᶥ,ζ,ζᶥ)+(δ,δᶥ,η,ηᶥ)=(α,αᶥ,β,βᶥ)  ∀  ααᶥββᶥ ∈ Moment matrixₜ
"""
function get_γγᶥζζᶥ_δδᶥηηᶥ(d,t)
    M = Moments.make_mon_expo(d,t)
    M_vec = M[:,1]
    γγᶥζζᶥ_δδᶥηηᶥ_dict = Dict()
    for k in Moments.make_mon_expo(d,2*t[1])# Define variables in the moment matrix.
        fd = hcat(map(x -> [x[1],x[2]],findall(M .== [k]))...)
        γγᶥζζᶥ_δδᶥηηᶥ_dict[k] = [M_vec[fd[1,:]], M_vec[fd[2,:]]]
    end
    return γγᶥζζᶥ_δδᶥηηᶥ_dict
end

"""The coefficients function"""
pf(x) = prod(factorial.(x))
function get_coef(part1,part2,d)
    γ,γᶥ,ζ,ζᶥ = split_expo(part1,d) ; δ,δᶥ,η,ηᶥ = split_expo(part2,d)
    a = ((-1.0)^sum(δᶥ+ηᶥ))*((1.0im)^sum(δ+δᶥ+η+ηᶥ))*prod(pf.([[γ,γᶥ,ζ,ζᶥ]+[δ,δᶥ,η,ηᶥ]...]))
    return a/prod(pf.([γ,γᶥ,ζ,ζᶥ,δ,δᶥ,η,ηᶥ]))
end
get_coef(γ_δ_pair,d) = map((x,y)->get_coef(x,y,d),γ_δ_pair...)

""" get_expo(γγᶥδδᶥ,ζζᶥηηᶥ,d) = [γ+γᶥ,δ+δᶥ,η+ηᶥ,ζ+ζᶥ] """
function get_expo(γ_,δ_,d)
    γ,γᶥ,ζ,ζᶥ = split_expo(γ_,d) ; δ,δᶥ,η,ηᶥ = split_expo(δ_,d)
    return vcat(γ+γᶥ,δ+δᶥ,ζ+ζᶥ,η+ηᶥ)
end
get_expo(γ_δ_,d) = map((x,y)->get_expo(x,y,d),γ_δ_...)

"""f"""
function get_ℜℑααᶥββᶥᴿ(d,ααᶥββᶥ::Array{Int64,1},γdict)
    γ_δ_ = γdict[ααᶥββᶥ]
    return   get_expo(γ_δ_,d),get_coef(γ_δ_,d)
end

function get_ℜℑααᶥββᶥᴿ(d,B,γdict)
    MMex = map(x->get_ℜℑααᶥββᶥᴿ(d,x,γdict)[1],B)
    MMCoef = map(x->get_ℜℑααᶥββᶥᴿ(d,x,γdict)[2],B)

    real_mask = isreal.(MMCoef)

    RMMex = copy(MMex) ; IMMex = copy(MMex)
    RMMex[.!real_mask] .= [[]]
    IMMex[real_mask] .= [[]]

    RMMCoef = copy(MMCoef) ; IMMCoef = copy(MMCoef)
    RMMCoef[.!real_mask] .= [[0.0]] ;RMMCoef = real.(RMMCoef)
    IMMCoef[real_mask] .= [[0.0]] ; IMMCoef = imag.(IMMCoef)

    MMexᴿ = hcat(vcat(RMMex,IMMex),vcat(IMMex,RMMex))
    MMCoefᴿ = hcat(vcat(RMMCoef,IMMCoef),vcat((-1.0)*IMMCoef,RMMCoef))
    return MMexᴿ,MMCoefᴿ
end

""""""
function get_ℂ_block_diag(d,t;noBlock = false)
    if noBlock
        MMexᴿ = Dict("Default" => map(x->[x],make_mon_expo(d,t)))
        MMCoefᴿ = Dict("Default" => fill([1.0],size(MMexᴿ["Default"])...))
        return MMexᴿ,MMCoefᴿ
    end
    xx̄yȳMM_blocks = get_xx̄yȳMM_blocks(d,t)
    γγᶥζζᶥ_δδᶥηηᶥ = get_γγᶥζζᶥ_δδᶥηηᶥ(d,t)

    MMexᴿ= Dict() ; MMCoefᴿ = Dict()
    for B in keys(xx̄yȳMM_blocks)
        MMexᴿ[B], MMCoefᴿ[B] = get_ℜℑααᶥββᶥᴿ(d,xx̄yȳMM_blocks[B],γγᶥζζᶥ_δδᶥηηᶥ)
    end
    return MMexᴿ,MMCoefᴿ
end

end

max_stars_repo_name: JAndriesJ/sep-rank-bounding

max_issues_repo_name: JAndriesJ/sep-rank-bounding

max_forks_repo_name: JAndriesJ/sep-rank-bounding


text:
using Test
using DataStructures

as = rand(10)
bs = rand(10)

dict1 = Dict(as .=> bs)
dict2 = SortedDict( as .=> bs)

f1 = DictFunction(dict1)
f2 = DictFunction(dict2)

@test f1.(as) == f2.(as)
@test f1(rand()) == 0.0
@test f2(rand()) == 0.0

f3 = DictFunction(dict1, 3.0)

@test f3.(as) == f1.(as)
@test f3(rand()) == 3.0


text:
Dict_to_NamedTuple(tmp::Dict) = (; zip(Symbol.(keys(tmp)), values(tmp))...)

"""
    UnitGrid(γ::gcmgrid;option="minimal")

Generate a unit grid, where every grid spacing and area is 1, according to γ.
"""
function UnitGrid(γ::gcmgrid;option="minimal")
    nFaces=γ.nFaces
    ioSize=(γ.ioSize[1],γ.ioSize[2])

    Γ=Dict()

    pc=fill(0.5,2); pg=fill(0.0,2); pu=[0.,0.5]; pv=[0.5,0.];
    if option=="full"
        list_n=("XC","XG","YC","YG","RAC","RAW","RAS","RAZ","DXC","DXG","DYC","DYG","Depth","hFacC","hFacS","hFacW");
        list_u=(u"m",u"m",u"m",u"m",u"m^2",u"m^2",u"m^2",u"m^2",u"m",u"m",u"m",u"m",u"m",1.0,1.0,1.0)
        list_p=(pc,pg,pc,pg,pc,pu,pv,pg,pu,pv,pv,pu,pc,fill(0.5,3),[0.,0.5,0.5],[0.5,0.,0.5])
    else
        list_n=("XC","YC");
        list_u=(u"°",u"°")
        list_p=(pc,pc)
    end

    for ii=1:length(list_n);
        tmp1=fill(1.,(ioSize[:]))
        m=varmeta(list_u[ii],list_p[ii],missing,list_n[ii],list_n[ii]);
        tmp1=γ.read(tmp1,MeshArray(γ,Float64;meta=m));
        Γ[list_n[ii]]=tmp1
    end

    for i in 1:nFaces
        (np,nq)=γ.fSize[i]
        Γ["XC"][i]=vec(0.5:1.0:np-0.5)*ones(1,nq)
        option=="full" ? Γ["XG"][i]=vec(0.0:1.0:np-1.0)*ones(1,nq) : nothing
        Γ["YC"][i]=ones(np,1)*transpose(vec(0.5:1.0:nq-0.5))
        option=="full" ? Γ["YG"][i]=ones(np,1)*transpose(vec(0.0:1.0:nq-1.0)) : nothing
    end

    Γ=Dict_to_NamedTuple(Γ)

    return Γ

end

"""
    UnitGrid(ioSize, tileSize; option="minimal")

Generate a unit grid, where every grid spacing and area is 1, according to `ioSize, tileSize`.

Since `ioSize, tileSize` defines a one to one mapping from global to tiled array, here we use
`read_tiles, write_tiles` instead of the default `read, write`. And we overwrite local `XC,YC`
etc accordingly with global `XC,YC` etc .
"""
function UnitGrid(ioSize::NTuple{2, Int},tileSize::NTuple{2, Int}; option="minimal")

    nF=div(prod(ioSize),prod(tileSize))
    fSize=fill(tileSize,nF)

    γ=gcmgrid("","PeriodicDomain",nF,fSize, ioSize, Float32,
        MeshArrays.read_tiles, MeshArrays.write_tiles)
    Γ=UnitGrid(γ;option=option)

    Γ.XC[:]=γ.read([i-0.5 for i in 1:ioSize[1], j in 1:ioSize[2]],γ)
    Γ.YC[:]=γ.read([j-0.5 for i in 1:ioSize[1], j in 1:ioSize[2]],γ)
    if option=="full"
        Γ.XG[:]=γ.read([i-1.0 for i in 1:ioSize[1], j in 1:ioSize[2]],γ)
        Γ.YG[:]=γ.read([j-1.0 for i in 1:ioSize[1], j in 1:ioSize[2]],γ)
    end

    return Γ,γ
end

"""
    simple_periodic_domain(np::Integer,nq=missing)

Set up a simple periodic domain of size np x nq

```jldoctest; output = false
using MeshArrays
np=16 #domain size is np x np
Γ=simple_periodic_domain(np)
isa(Γ.XC,MeshArray)

# output

true
```

"""
function simple_periodic_domain(np::Integer,nq=missing)
    ismissing(nq) ? nq=np : nothing

    nFaces=1
    ioSize=[np nq]
    facesSize=Array{NTuple{2, Int},1}(undef,nFaces)
    facesSize[:].=[(np,nq)]
    ioPrec=Float32
    γ=gcmgrid("","PeriodicDomain",1,facesSize, ioSize, ioPrec, read, write)

    return UnitGrid(γ)
end

## GridSpec function with default GridName argument:

GridSpec() = GridSpec("PeriodicDomain","./")

## GridSpec function with GridName argument:

"""
    GridSpec(GridName,GridParentDir="./")

Select one of the pre-defined grids (by `GridName`) and return
the corresponding `gmcgrid` -- a global grid specification
which contains the grid files location (`GridParentDir`).

Possible choices for `GridName`:

- `"PeriodicDomain"`
- `"PeriodicChannel"`
- `"CubeSphere"`
- `"LatLonCap"``

```jldoctest; output = false
using MeshArrays
g = GridSpec()
g = GridSpec("PeriodicChannel",MeshArrays.GRID_LL360)
g = GridSpec("CubeSphere",MeshArrays.GRID_CS32)
g = GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
isa(g,gcmgrid)

# output

true
```
"""
function GridSpec(GridName,GridParentDir="./")

grDir=GridParentDir
if GridName=="LatLonCap"
    nFaces=5
    grTopo="LatLonCap"
    ioSize=[90 1170]
    facesSize=[(90, 270), (90, 270), (90, 90), (270, 90), (270, 90)]
    ioPrec=Float64
elseif GridName=="CubeSphere"
    nFaces=6
    grTopo="CubeSphere"
    ioSize=[32 192]
    facesSize=[(32, 32), (32, 32), (32, 32), (32, 32), (32, 32), (32, 32)]
    ioPrec=Float32
elseif GridName=="PeriodicChannel"
    nFaces=1
    grTopo="PeriodicChannel"
    ioSize=[360 160]
    facesSize=[(360, 160)]
    ioPrec=Float32
elseif GridName=="PeriodicDomain"
    nFaces=4
    grTopo="PeriodicDomain"
    ioSize=[80 42]
    facesSize=[(40, 21), (40, 21), (40, 21), (40, 21)]
    ioPrec=Float32
else
    error("unknown GridName case")
end

return gcmgrid(grDir,grTopo,nFaces,facesSize, ioSize, ioPrec, read, write)

end

## GridLoad function

"""
    GridLoad(γ::gcmgrid;option="minimal")

Return a `NamedTuple` of grid variables read from files located in `γ.path` (see `?GridSpec`).

By default, option="minimal" means that only grid cell center positions (XC, YC) are loaded.

The "full" option provides a complete set of grid variables.

Based on the MITgcm naming convention, grid variables are:

- XC, XG, YC, YG, AngleCS, AngleSN, hFacC, hFacS, hFacW, Depth.
- RAC, RAW, RAS, RAZ, DXC, DXG, DYC, DYG.
- DRC, DRF, RC, RF (one-dimensional)

```jldoctest; output = false
using MeshArrays

γ = GridSpec("CubeSphere",MeshArrays.GRID_CS32)
#γ = GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
#γ = GridSpec("PeriodicChannel",MeshArrays.GRID_LL360)

Γ = GridLoad(γ;option="full")

isa(Γ.XC,MeshArray)

# output

true
```
"""
function GridLoad(γ::gcmgrid;option="minimal")

    γ.path==GRID_CS32 ? GRID_CS32_download() : nothing
    γ.path==GRID_LL360 ? GRID_LL360_download() : nothing
    γ.path==GRID_LLC90 ? GRID_LLC90_download() : nothing

    Γ=Dict()

    if option=="full"
        list_n=("XC","XG","YC","YG","RAC","RAW","RAS","RAZ","DXC","DXG","DYC","DYG","Depth");
        if (!isempty(filter(x -> occursin("AngleCS",x), readdir(γ.path))))
            list_n=(list_n...,"AngleCS","AngleSN");
        end
        list_n=(list_n...,"DRC","DRF","RC","RF")
        list_n=(list_n...,"hFacC","hFacS","hFacW");
    else
        list_n=("XC","YC");
    end

    [Γ[ii]=GridLoadVar(ii,γ) for ii in list_n]
    option=="full" ? GridAddWS!(Γ) : nothing
    return Dict_to_NamedTuple(Γ)
end

"""
    GridLoadVar(nam::String,γ::gcmgrid)

Return a grid variable read from files located in `γ.path` (see `?GridSpec`, `?GridLoad`).

Based on the MITgcm naming convention, grid variables are:

- XC, XG, YC, YG, AngleCS, AngleSN, hFacC, hFacS, hFacW, Depth.
- RAC, RAW, RAS, RAZ, DXC, DXG, DYC, DYG.
- DRC, DRF, RC, RF (one-dimensional)

```jldoctest; output = false
using MeshArrays

γ = GridSpec("CubeSphere",MeshArrays.GRID_CS32)
XC = GridLoadVar("XC",γ)

isa(XC,MeshArray)

# output

true
```
"""
function GridLoadVar(nam::String,γ::gcmgrid)
    pc=fill(0.5,2); pg=fill(0.0,2); pu=[0.,0.5]; pv=[0.5,0.];
    list_n=("XC","XG","YC","YG","RAC","RAW","RAS","RAZ","DXC","DXG","DYC","DYG","Depth","AngleCS","AngleSN")
    list_u=(u"°",u"°",u"°",u"°",u"m^2",u"m^2",u"m^2",u"m^2",u"m",u"m",u"m",u"m",u"m",1.0,1.0)
    list_p=(pc,pg,pc,pg,pc,pu,pv,pg,pu,pv,pv,pu,pc,pc,pc)
    #
    list3d_n=("hFacC","hFacS","hFacW");
    list3d_u=(1.0,1.0,1.0)
    list3d_p=(fill(0.5,3),[0.,0.5,0.5],[0.5,0.,0.5])
    #
    list1d_n=("DRC","DRF","RC","RF")

    if sum(nam.==list_n)==1
        ii=findall(nam.==list_n)[1]
        m=varmeta(list_u[ii],list_p[ii],missing,list_n[ii],list_n[ii])
        tmp1=γ.read(joinpath(γ.path,list_n[ii]*".data"),MeshArray(γ,γ.ioPrec;meta=m))
    elseif sum(nam.==list1d_n)==1
        fil=joinpath(γ.path,nam*".data")
        γ.ioPrec==Float64 ? reclen=8 : reclen=4
        n3=Int64(stat(fil).size/reclen)

        fid = open(fil)
        tmp1 = Array{γ.ioPrec,1}(undef,n3)
        read!(fid,tmp1)
        tmp1 = hton.(tmp1)
    elseif sum(nam.==list3d_n)==1
        fil=joinpath(γ.path,"RC.data")
        γ.ioPrec==Float64 ? reclen=8 : reclen=4
        n3=Int64(stat(fil).size/reclen)

        ii=findall(nam.==list3d_n)[1]
        m=varmeta(list3d_u[ii],list3d_p[ii],missing,list3d_n[ii],list3d_n[ii]);
        tmp1=γ.read(joinpath(γ.path,list3d_n[ii]*".data"),MeshArray(γ,γ.ioPrec,n3;meta=m))
    else
        tmp1=missing
    end
end

"""
    GridOfOnes(grTp,nF,nP;option="minimal")

Define all-ones grid variables instead of using `GridSpec` & `GridLoad`. E.g.

```
γ,Γ=GridOfOnes("CubeSphere",6,20);
```
"""
function GridOfOnes(grTp,nF,nP;option="minimal")

    grDir=""
    grTopo=grTp
    nFaces=nF
    if grTopo=="LatLonCap"
        ioSize=[nP nP*nF]
    elseif grTopo=="CubeSphere"
        ioSize=[nP nP*nF]
    elseif grTopo=="PeriodicChannel"
        ioSize=[nP nP]
    elseif grTopo=="PeriodicDomain"
        nFsqrt=Int(sqrt(nF))
        ioSize=[nP*nFsqrt nP*nFsqrt]
    else
        error("unknown configuration (grTopo)")
    end
    facesSize=Array{NTuple{2, Int},1}(undef,nFaces)
    facesSize[:].=[(nP,nP)]
    ioPrec=Float32

    γ=gcmgrid(grDir,grTopo,nFaces,facesSize, ioSize, ioPrec, read, write)

    return γ, UnitGrid(γ;option=option)
end

"""
    Tiles(γ::gcmgrid,ni::Int,nj::Int)

Define sudomain `tiles` of size `ni,nj`. Each tile is defined by a `Dict` where
`tile,face,i,j` correspond to tile ID, face ID, index ranges.

```jldoctest; output = false
using MeshArrays
γ=GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
τ=Tiles(γ,30,30)

isa(τ[1],NamedTuple)

# output

true
```
"""
function Tiles(γ::gcmgrid,ni::Int,nj::Int)
    nt=Int(prod(γ.ioSize)/ni/nj)
    τ=Array{Dict,1}(undef,nt)
    #
    cnt=0
    for iF=1:γ.nFaces
        for jj=Int.(1:γ.fSize[iF][2]/nj)
            for ii=Int.(1:γ.fSize[iF][1]/ni)
                cnt=cnt+1
                i=(1:ni).+ni*(ii-1)
                j=(1:nj).+nj*(jj-1)
                τ[cnt]=Dict("tile" => cnt, "face" => iF, "i" => i, "j" => j)
            end
        end
    end
    #
    return Dict_to_NamedTuple.(τ)
end

"""
    Tiles(τ::Array{Dict},x::MeshArray)

Return an `Array` of tiles which cover `x` according to tile partition `τ`.

```jldoctest; output = false
using MeshArrays
γ=GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
d=γ.read(γ.path*"Depth.data",MeshArray(γ,γ.ioPrec))
τ=Tiles(γ,30,30)
td=Tiles(τ,d)

isa(td[1],Array)

# output

true
```
"""
function Tiles(τ::Array,x::MeshArray)
    nt=length(τ)
    tx=Array{typeof(x[1]),1}(undef,nt)
    dn=size(x[1],1)-x.fSize[1][1]
    for ii=1:nt
        f=τ[ii].face
        i0=minimum(τ[ii].i)
        i1=maximum(τ[ii].i)
        j0=minimum(τ[ii].j)
        j1=maximum(τ[ii].j)
        tx[ii]=view(x[f],i0:i1+dn,j0:j1+dn)
    end
    return tx
end

"""
    GridAddWS!(Γ::Dict)

Compute XW, YW, XS, and YS (vector field locations) from XC, YC (tracer
field locations) and add them to Γ.
"""
function GridAddWS!(Γ::Dict)

    XC=exchange(Γ["XC"])
    YC=exchange(Γ["YC"])
    nFaces=XC.grid.nFaces
    uX=XC.meta.unit
    uY=YC.meta.unit

    XW=MeshArrays.gcmarray(XC.grid,eltype(XC);meta=varmeta(uX,[0.,0.5],missing,"XW","XW"))
    YW=MeshArrays.gcmarray(XC.grid,eltype(XC);meta=varmeta(uY,[0.,0.5],missing,"YW","YW"))
    XS=MeshArrays.gcmarray(XC.grid,eltype(XC);meta=varmeta(uX,[0.5,0.],missing,"XS","XS"))
    YS=MeshArrays.gcmarray(XC.grid,eltype(XC);meta=varmeta(uY,[0.5,0.],missing,"YS","YS"))

    for ff=1:nFaces
        tmp1=XC[ff][1:end-2,2:end-1]
        tmp2=XC[ff][2:end-1,2:end-1]
        tmp2[tmp2.-tmp1.>180]=tmp2[tmp2.-tmp1.>180].-360;
        tmp2[tmp1.-tmp2.>180]=tmp2[tmp1.-tmp2.>180].+360;
        XW[ff]=(tmp1.+tmp2)./2;
       #
        tmp1=XC[ff][2:end-1,1:end-2]
        tmp2=XC[ff][2:end-1,2:end-1]
        tmp2[tmp2.-tmp1.>180]=tmp2[tmp2.-tmp1.>180].-360;
        tmp2[tmp1.-tmp2.>180]=tmp2[tmp1.-tmp2.>180].+360;
        XS[ff]=(tmp1.+tmp2)./2;
       #
        tmp1=YC[ff][1:end-2,2:end-1]
        tmp2=YC[ff][2:end-1,2:end-1]
        YW[ff]=(tmp1.+tmp2)./2;
       #
        tmp1=YC[ff][2:end-1,1:end-2]
        tmp2=YC[ff][2:end-1,2:end-1]
        YS[ff]=(tmp1.+tmp2)./2;
    end;

    Xmax=180; Xmin=-180;
    XS[findall(XS.<Xmin)]=XS[findall(XS.<Xmin)].+360;
    XS[findall(XS.>Xmax)]=XS[findall(XS.>Xmax)].-360;
    XW[findall(XW.<Xmin)]=XW[findall(XW.<Xmin)].+360;
    XW[findall(XW.>Xmax)]=XW[findall(XW.>Xmax)].-360;

    Γ["XW"]=XW
    Γ["XS"]=XS
    Γ["YW"]=YW
    Γ["YS"]=YS
    return Γ
end


text:
Lazy.@forward ModelSwitches.pvec (
    is_calibrated, calibrate!, fix!, param_value, param_default_value, retrieve
    );

has_pvector(switches :: ModelSwitches) = false;

function get_pvector(switches :: ModelSwitches)
    if has_pvector(switches)
        return switches.pvec;
    else
        return nothing;
    end
end

"""
	$(SIGNATURES)

Does object contain a `ParamVector`. Override for user defined types.
By default, ModelObjects are assumed to have `ParamVector`s.
"""
has_pvector(o :: ModelObject) = true;
has_pvector(o) = false;

function ModelObjectsLH.get_object_id(switches :: ModelSwitches)
    if has_pvector(switches)
        return get_object_id(get_pvector(switches))
    else
        return nothing
    end
end

"""
	$(SIGNATURES)

Retrieve the switches that govern the `ModelObject`'s behavior. Nothing if not found.
"""
function get_switches(o :: ModelObject)
    if isdefined(o, :switches)
        return o.switches;
    else
        return nothing;
    end
end
get_switches(o) = nothing;

has_switches(o) = !isnothing(get_switches(o));

"""
	$(SIGNATURES)

Find the ParamVector in a `ModelObject`. Return empty if not found.
Override for user types.
"""
function get_pvector(o)
    # Try the default field
    if has_pvector(o)  &&  isdefined(o, :pvec)
        pvec = o.pvec;
    elseif has_pvector(get_switches(o))
        pvec = get_switches(o).pvec;
    else
        pvec = ParamVector(o.objId);
    end
    return pvec :: ParamVector
end

"""
	$(SIGNATURES)

Retrieve a `Param`.
"""
function retrieve(o, pName :: Symbol)
    pvec = get_pvector(o);
    return retrieve(pvec, pName)
end

# Code permits objects without ParamVector or child objects
function validate(o :: T1) where T1 <: ModelObject
    isValid = true;
    if !isa(get_object_id(o), ObjectId)
        @warn "Invalid ObjectId"
        isValid = false;
    end
    if has_pvector(o)
        if !isa(get_pvector(o), ParamVector)
            @warn "Invalid ParamVector"
            isValid = false;
        end
    end
    return isValid
end

function validate_all_params(o :: ModelObject; silent = true)
    return validate_pvec(collect_pvectors(o); silent);
end

## -----------  Check that params are consistent with `ParamVector`s

function check_param_values(x :: ModelObject, isCalibrated :: Bool)
    isValid = true;
    for o in collect_model_objects(x)
        if !check_own_param_values(o, get_pvector(o), isCalibrated);
            isValid = false;
        end
    end
    return isValid
end

check_calibrated_params(x :: ModelObject) = check_param_values(x, true);
check_fixed_params(x :: ModelObject) = check_param_values(x, false);

"""
    $(SIGNATURES)

Check that param vector values are consistent with object values.
Does not reach into child objects.
"""
function check_own_param_values(x :: ModelObject, pvec, isCalibrated :: Bool)
    # d = make_dict(pvec; isCalibrated = isCalibrated, useValues = isCalibrated);
    pList = calibrated_params(pvec; isCalibrated);
    valid = true;
    # for (pName, pValue) in d
    for p in pList
        pName = name(p);
        if isCalibrated
            pValue = value(p);
        else
            pValue = default_value(p);
        end
        isValid = (getproperty(x, pName) ≈ pValue);
        if !isValid
            valid = false;
            propValue = getproperty(x, pName);
            @warn "Invalid value: $pName: $pValue vs. $propValue";
        end
    end
    return valid
end

"""
	$(SIGNATURES)

Check that calibrated param vector values are consistent with object values.
Does not reach into child objects.
"""
check_own_calibrated_params(x :: ModelObject, pvec) =
    check_own_param_values(x, pvec, true);

"""
    $(SIGNATURES)

Check that all fixed parameters have the correct values
Does not reach into child objects.
"""
check_own_fixed_params(x :: ModelObject, pvec) =
    check_own_param_values(x, pvec, false);

"""
	$(SIGNATURES)

Do two model objects have the same calibrated parameters?
Includes child objects.
This is mainly useful for testing (e.g., loading parameters).
"""
function params_equal(x1 :: ModelObject, x2 :: ModelObject)
    guess1 = make_guess(x1);
    guess2 = make_guess(x2);
    v1 = get_values(x1, guess1);
    v2 = get_values(x2, guess2);
    return isapprox(v1, v2; atol = 1e-6)
end

## -----------  Change values

"""
	$(SIGNATURES)

Set bounds for a parameter. Input must support `retrieve` for a `Param`.
"""
function set_bounds!(pvec, pName :: Symbol; lb = nothing, ub = nothing)
    p = retrieve(pvec, pName);
    set_bounds!(p; lb, ub);
end

"""
	$(SIGNATURES)

Change value of a field in a `ModelObject` and its `ParamVector`.

# Example
```julia
change_value!(mObj, :utilityFunction, :sigma, 2.0);
get_value(mObj, :utilityFunction, :sigma) ≈ 2.0
```
"""
function change_value!(x :: ModelObject, oName :: Symbol, pName :: Symbol,  newValue)
    obj = find_only_object(x, oName);
    pvec = get_pvector(obj);
    @assert length(pvec) > 0  "No ParamVector in $oName / $pName"
    oldValue = change_value!(pvec, pName, newValue);
    # Set value in object as well
    setfield!(obj, pName, deepcopy(newValue));
    return oldValue
end

"""
	$(SIGNATURES)

Retrieve a parameter value in a `ModelObject`.
"""
function get_value(x :: ModelObject, oName :: Symbol, pName :: Symbol)
    obj = find_only_object(x, oName);
    pvec = get_pvector(obj);
    @assert length(pvec) > 0  "No ParamVector in $oName / $pName"
    @assert param_exists(pvec, pName)  "$pName not found";
    p = retrieve(pvec, pName);
    return value(p)
end

## Set fields in struct from param vector (using values, not defaults)
# Does not reach into child objects.
function set_own_values_from_pvec!(x :: ModelObject, isCalibrated :: Bool;
    alwaysUseDefaultValue = false)
    pvec = get_pvector(x);
    # d = make_dict(pvec; isCalibrated = isCalibrated, useValues = isCalibrated);
    # set_own_values_from_dict!(x, d);

    pList = calibrated_params(pvec; isCalibrated);
    for p in pList
        set_own_value_from_param!(x, p; alwaysUseDefaultValue);
    end
    return nothing
end

# Ignores if no corresponding property exists in x.
function set_own_value_from_param!(x :: ModelObject, p :: AbstractParam;
    alwaysUseDefaultValue = false)
    pName = name(p);
    if pName ∈ propertynames(x)
        if is_calibrated(p)  &&  (!alwaysUseDefaultValue)
            val = value(p);
        else
            val = default_value(p);
        end
        setfield!(x, pName, val);
    end
end

"""
	$(SIGNATURES)

Set fields in a struct from a Dict{Symbol, Any}.
Does not change `ParamVector` inside `x` (if any).
Does not change child objects.
"""
function set_own_values_from_dict!(x :: ModelObject,  d :: D1) where D1 <: AbstractDict
    for (k, val) in d
        if k ∈ propertynames(x)
            setfield!(x, k, val);
        else
            @warn "Field $k not found"
        end
    end
    return nothing
end

## Set default values from param vector
# Typically for non-calibrated parameters
function set_own_default_values!(x :: ModelObject, isCalibrated :: Bool)
    set_own_values_from_pvec!(x, isCalibrated; alwaysUseDefaultValue = true);
    # pvec = get_pvector(x);
    # # Last arg: use default values
    # d = make_dict(pvec; isCalibrated = isCalibrated, useValues = false);
    # set_own_values_from_dict!(x, d);
    return nothing
end

"""
	$(SIGNATURES)

Sync all values from all objects' param vectors into objects.
"""
function sync_values!(x :: ModelObject)
    for o in collect_model_objects(x)
        if has_pvector(o)
            sync_own_values!(o);
        end
    end
end

"""
    $(SIGNATURES)

Sync all values from object's param vector into object.
"""
function sync_own_values!(x :: ModelObject)
    # Set calibrated parameters
    set_own_values_from_pvec!(x, true);
    # Set fixed parameters
    set_own_default_values!(x, false);
end

"""
	$(SIGNATURES)

Set all parameters to calibrated or fixed. Useful for experiments that calibrate only select parameters. Recurses into child objects.
"""
function set_calibration_status_all_params!(x :: ModelObject, isCalibrated :: Bool)
    pvecV = collect_pvectors(x);
    for (_, pvec) in pvecV
        set_calibration_status_all_params!(pvec, isCalibrated);
    end
end

"""
	$(SIGNATURES)

Set default values for all parameters to values.
Useful for fixing parameters at previously calibrated values.
"""
function set_default_values_all_params!(x :: ModelObject)
    pvecV = collect_pvectors(x);
    for (_, pvec) in pvecV
        set_default_values_all_params!(pvec);
    end
end

"""
	$(SIGNATURES)

Return all `Param`s in a `ModelObject`. Optionally filtered by calibration status.
Returns `Vector{AbstractParam}`. Note that names of parameters may not be unique.
"""
function all_params(x :: ModelObject; isCalibrated = nothing)
    pvecV = collect_pvectors(x);
    pList = Vector{AbstractParam}();
    for (_, pvec) in pvecV
        pvecList = all_params(pvec; isCalibrated);
        isempty(pvecList)  ||  append!(pList, pvecList);
    end
    return pList
end

"""
	$(SIGNATURES)

Compare parameters across objects. Optionally filtered on calibration status.
Returns list of
- missing on o1
- missing in o2
- different between the objects

Each is a `Dict` that contains another `Dict` for each `ObjectId` with differences.
"""
function compare_params(o1, o2; ignoreCalibrationStatus :: Bool = true)
    pvec1V = collect_pvectors(o1);
    pvec2V = collect_pvectors(o2);

    pMiss1 = Dict{ObjectId, Any}();
    pMiss2 = Dict{ObjectId, Any}();
    pDiff = Dict{ObjectId, Any}();
    for (objId, pvec1) in pvec1V
        if has_pvector(pvec2V, objId)
            pvec2 = find_pvector(pvec2V, objId);
            miss1, miss2, dif = compare_params(pvec1, pvec2; ignoreCalibrationStatus);
            isempty(miss1)  ||  (pMiss1[objId] = miss1);
            isempty(miss2)  ||  (pMiss2[objId] = miss2);
            isempty(dif)    ||  (pDiff[objId] = dif);
        else
            @warn "ParamVector $objId not found."
        end
    end

    return pMiss1, pMiss2, pDiff
end

"""
	$(SIGNATURES)

Dict with parameter info for an entire `ModelObject` and its children.
"""
function make_dict(o :: ModelObject, isCalibrated :: Bool)
    pvecV = collect_pvectors(o);
    return make_dict(pvecV; isCalibrated, valueType = :value);
end

"""
	$(SIGNATURES)

Check that param values match across two model objects (including children).
Allow for differences in the objects listed by own name (Symbol) in `differentIdV`.
"""
function check_params_match(m1, m2,
    allowedIdV :: AbstractVector{ObjectId};
    ignoreCalibrationStatus = true)

    isValid = true;
    miss1, miss2, pDiff = compare_params(m1, m2; ignoreCalibrationStatus);
    if !isempty(miss1)
        isValid = false;
        @warn "Missing params in m1:  $miss1";
    end
    if !isempty(miss2)
        isValid = false;
        @warn "Missing params in m2:  $miss2";
    end

    # Check that differences occur only in permitted objects.
    diffIdV = collect(keys(pDiff));
    if !isempty(diffIdV)
        for oId in diffIdV
            if !any_isequal(oId, allowedIdV)
                isValid = false;
                @warn "Param differences in $oId";
            end
        end
    end
    return isValid
end

# function collect_object_ids(m, differentIdV)
#     return [get_object_id(find_only_object(m, oName))
#         for oName in differentIdV];
#     # for oName in differentIdV
#     #     o = find_object(m, oName);
#     #     if isempty(o)
#     #         isValid = false;
#     #         @warn "Object $oName not found";
#     #     else
#     #         oId = get_object_id(o);
#     #     end
#     # end
# end

# Because broadcasting ObjectId not implemented yet.
function any_isequal(oId, oIdV)
    found = false;
    for oId2 in oIdV
        if isequal(oId, oId2)
            found = true;
            break;
        end
    end
    return found
end

# -----------------