text:
    # David Schaffer. Multiple Objective Optimization with Vector Evaluated Genetic Algorithms.
    # In Proceedings of the 1st International Conference on Genetic Algorithms,
    # L. Erlbaum Associates Inc. pp. 93–100, 1985.
    # http://dl.acm.org/citation.cfm?id=645511.657079

    using NSGAIII

    # using PyPlot
    # function plot_pop(P)
    #     clf()
    #     p = plot(map(x -> x.y[1], P), map(x -> x.y[2], P), "bo", markersize=1)
    #     !isinteractive() && show()
    #     sleep(0.2)
    # end

    using UnicodePlots
    function plot_pop(P)
        println()
        display(scatterplot(map(x -> x.y[1], P), map(x -> x.y[2], P)))
        sleep(0.4)
    end

    const d = RealCoding(6, [-10], [10])
    z1(x) = x^2
    z2(x) = (x-2)^2
    z(bits) = begin
        x = decode(bits, d)[1]
        z1(x), z2(x)
    end

    pop = nsga(100, 20, ()->rand(Bool, d.nbbitstotal), z, 2, fplot = plot_pop)
    pop = nsga(100, 20, ()->rand(Bool, d.nbbitstotal), z, [[1.,0.],[0.5,0.5],[0.,1.]], fplot = plot_pop)

    pop = nsga(100, 20, ()->rand(Bool, d.nbbitstotal), z, 30, fplot = plot_pop)

    [(decode(ind.x, d)[1], ind.y) for ind in pop]


text:
module {{{PKG}}}

# Define a handler function that is called by the Lambda runtime
function handle_event(event_data, headers)
    @info "Handling request" event_data headers
    return "Hello World!"
end

end

ext: jl


text:
using ModernGL, GLWindow, GLAbstraction, GLFW, GeometryTypes

const window = GLWindow.create_glcontext("Example", resolution=(512, 512), debugging=true)

const vsh = vert"""
{{GLSL_VERSION}}
in vec2 position;

void main(){
	gl_Position = vec4(position, 0, 1.0);
}
"""

const fsh = frag"""
{{GLSL_VERSION}}
out vec4 outColor;

void main() {
	outColor = vec4(1.0, 1.0, 1.0, 1.0);
}
"""

const triangle = std_renderobject(
	Dict{Symbol, Any}(
        :position => GLBuffer(Point2f0[(0.0, 0.5), (0.5, -0.5), (-0.5,-0.5)]),
    ),
	LazyShader(vsh, fsh)
)

glClearColor(0, 0, 0, 1)

while isopen(window)
  	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)

    render(triangle)

  	swapbuffers(window)
  	pollevents()
end


text:
using CCBlade
using Statistics: mean

# Comparison using geometry and data from)
# https://rotorcraft.arc.nasa.gov/Publications/files/RamasamyGB_ERF10_836.pdf

# add an "overloaded" function to handle the tip loss used just for hover)
function solvehover(rotor, section, op)

	 if section.r/rotor.Rtip > 0.94
        phi = 0.0
        theta = section.theta + op.pitch
        alpha = theta - phi
        _, cd = afeval(section.af, alpha, 1.0, 1.0)
        cl = 0.0

        # resolve into normal and tangential forces
        cn = 0.0
        ct = cd
        F = 1.0
        G = 1.0

        u = 0.0
        v = 0.0
        a = 0.0
        ap = 0.0
        W = op.Vy
        Np = 0.0
        Tp = ct*0.5*op.rho*W^2*section.chord

        return Outputs(Np, Tp, a, ap, u, v, phi, alpha, W, cl, cd, cn, ct, F, G)
    else
    	return solve(rotor, section, op)
    end
end

chord = 0.060
theta = 0.0
Rtip = 0.656
Rhub = 0.19*Rtip

# turbine = false
B = 3

rotor = Rotor(Rhub, Rtip, B, turbine=false)

r = range(Rhub, Rtip, length=30)

af = AlphaAF("naca0012v2.txt", radians=false)
function af2(alpha, Re, M)
	cl, cd = afeval(af, alpha, Re, M)
	return cl, cd+0.014
end
sections = Section.(r, chord, theta, af2)

rho = 1.225
Omega = 800*pi/30
Vinf = 0.0

Vinf2 = 1e-4

nP = 40
pitch = range(1e-4, 20*pi/180, length=nP)

CT = zeros(nP)
CQ = zeros(nP)
FM = zeros(nP)
CT2 = zeros(nP)
CQ2 = zeros(nP)
FM2 = zeros(nP)

for i = 1:nP
	local op = simple_op.(Vinf, Omega, r, rho, pitch=pitch[i])
	local op2 = simple_op.(Vinf2, Omega, r, rho, pitch=pitch[i])
	local outputs = solvehover.(Ref(rotor), sections, op)
	local outputs2 = solvehover.(Ref(rotor), sections, op2)
	local T, Q = thrusttorque(rotor, sections, outputs)
	T2, Q2 = thrusttorque(rotor, sections, outputs2)

	FM[i], CT[i], CQ[i] = nondim(T, Q, Vinf, Omega, rho, rotor, "helicopter")
	FM2[i], CT2[i], CQ2[i] = nondim(T2, Q2, Vinf, Omega, rho, rotor, "helicopter")

end

sigma = B * chord / (pi * Rtip)

using PyPlot
close("all")

rc("axes.spines", right=false, top=false)
rc("font", size=14.0)
rc("legend", frameon=false)
rc("text", usetex="True")

figure()
plot(CT/sigma, CQ/sigma, color="#348ABD")
# plot(CT2/sigma, CQ2/sigma, "k--")
xlim([0, 0.15])
ylim([0, 0.025])
xlabel(L"C_T/\sigma")
ylabel(L"C_Q/\sigma", rotation=0)

data = [
0.0006486486486486538 0.001826767151767151
0.00021621621621621678 0.0019711018711018706
0.0009729729729729721 0.0024036122661122636
0.0021621621621621644 0.002283134095634097
0.010594594594594595 0.0027378378378378337
0.011351351351351346 0.0026174636174636154
0.012540540540540546 0.0024008316008316023
0.011675675675675679 0.003098154885654885
0.012000000000000004 0.0034346153846153867
0.017081081081081077 0.0035535862785862765
0.01935135135135136 0.0036491943866943863
0.017621621621621626 0.0029765332640332605
0.017081081081081077 0.0028324324324324295
0.019135135135135137 0.0027357848232848203
0.05016216216216216 0.004891787941787941
0.05059459459459458 0.005228222453222447
0.054054054054054064 0.005684121621621619
0.05502702702702702 0.006501195426195424
0.05427027027027025 0.006597531185031182
0.06940540540540538 0.007194854469854465
0.07189189189189188 0.007410602910602908
0.07243243243243241 0.007771049896049893
0.07178378378378378 0.00810774428274428
0.07675675675675675 0.008322895010395012
0.0829189189189189 0.008657952182952181
0.08572972972972973 0.009017853430353426
0.08637837837837836 0.0090176975051975
0.088 0.00923365384615384
0.08637837837837836 0.009570582120582115
0.0937297297297297 0.010602468814968808
0.09448648648648647 0.010265748440748436
0.09848648648648647 0.010841709979209978
0.0963243243243243 0.011178768191268186
0.09448648648648647 0.01264555613305613
0.11589189189189186 0.014948102910602907
0.11794594594594593 0.015139916839916836
0.11816216216216216 0.015284095634095628
0.11805405405405403 0.015596621621621619
0.11502702702702702 0.015621387733887729
0.11956756756756756 0.014586642411642408
0.12054054054054059 0.014850831600831596
0.1211891891891892 0.014850675675675672
]
plot(data[:, 1], data[:, 2], "o", color="#A60628")
text(0.133, 0.014, "BEM", color="#348ABD")
text(0.08, 0.017, "experimental", color="#A60628")
gca().yaxis.set_label_coords(-0.09,1.05)
# savefig("hoverCT.pdf")

figure()
plot(CT/sigma, FM, color="#348ABD")
xlim([0, 0.15])
ylim([0, 0.7])
xlabel(L"C_T/\sigma")
ylabel("figure of merit", rotation=0)

data = [
0.0007194244604316561 0.0031982942430702765
0.011305241521068862 0.09786780383795313
0.017163412127440908 0.16439232409381654
0.05056526207605345 0.45223880597014926
0.07194244604316546 0.5462686567164178
0.08581706063720454 0.5827292110874199
0.09362795477903393 0.5641791044776119
0.12127440904419323 0.5936034115138591
0.11582733812949644 0.5533049040511726
]
plot(data[:, 1], data[:, 2], "o", color="#A60628")
text(0.12, 0.67, "BEM", color="#348ABD")
text(0.12, 0.5, "experimental", color="#A60628")
gca().yaxis.set_label_coords(0,1.05)
# savefig("hoverFM.pdf")

# figure()
# plot(CT/sigma, CQ/sigma, color="#348ABD")
# plot(CT2/sigma, CQ2/sigma, "k--")
# xlim([0, 0.15])
# ylim([0, 0.025])
# xlabel(L"C_T/\sigma")
# ylabel(L"C_Q/\sigma", rotation=0)
# gca().yaxis.set_label_coords(-0.09,1.05)
# legend(["modified", "original"])
# # savefig("hoverCTcompare.pdf")

# figure()
# plot(CT/sigma, FM, color="#348ABD")
# plot(CT2/sigma, FM2, "k--")
# xlim([0, 0.15])
# ylim([0, 0.7])
# xlabel(L"C_T/\sigma")
# ylabel("figure of merit", rotation=0)
# gca().yaxis.set_label_coords(0,1.05)
# legend(["modified", "original"])
# # savefig("hoverFMcompare.pdf")

# figure()
# plot(pitch*180/pi, CT/sigma, color="#348ABD")
# plot(pitch*180/pi, CT2/sigma, "k--")
# xlabel("pitch (deg)")
# ylabel(L"C_T/\sigma", rotation=0)
# gca().yaxis.set_label_coords(-0.09,1.05)

figure()
plot(pitch*180/pi, CQ/sigma, color="#348ABD")
plot(pitch*180/pi, CQ2/sigma, "k--")
xlabel("pitch (deg)")
ylabel(L"C_Q/\sigma", rotation=0)
gca().yaxis.set_label_coords(-0.09,1.05)
text(13.7, 0.025, "reformulated", color="#348ABD")
text(18.5, 0.02, "original")
# savefig("CQcompare.pdf")

# ---------- distributed loads --------

# adjust pitch to match thrust
# CT = 0.0013
op = simple_op.(Vinf, Omega, r, rho, pitch=3.5*pi/180)
outputs = solvehover.(Ref(rotor), sections, op)

T, Q = thrusttorque(rotor, sections, outputs)
A = pi*Rtip^2
CT = T / (rho * A * (Omega * Rtip)^2)
println(CT)

figure()
plot(r/Rtip, outputs.u/(Omega*Rtip), color="#348ABD")
xlabel(L"r/R")
ylabel(L"u/(\Omega R)", rotation=0)

data = [
0.20619554695062925 -0.011297248925009579
0.21055179090029041 -0.009261520970755763
0.21926427879961274 -0.00896913482969
0.22434656340755083 -0.013037917332672955
0.24394966118102612 -0.011289932236205238
0.2323330106485963 -0.014780555618091357
0.23814133591481124 -0.013325941601566901
0.25338818973862537 -0.013759033296562279
0.2642787996127783 -0.015065062248137043
0.27516940948693125 -0.016225742362502543
0.28751210067763794 -0.015641955019248532
0.297676669893514 -0.018256264211260947
0.3136495643756051 -0.0188345640379117
0.3238141335914811 -0.02057678020666831
0.3354307841239109 -0.019121040545712426
0.3463213939980639 -0.020572418334496495
0.36592449177153913 -0.02245815416826133
0.37536302032913843 -0.02129353429838582
0.38625363020329134 -0.022744912087169916
0.39932236205227484 -0.023759821247664292
0.41384317521781216 -0.022884914113329907
0.4283639883833494 -0.023754193025507117
0.4428848015488867 -0.02564091379814945
0.46030977734753137 -0.02578288570206444
0.47773475314617614 -0.027087648303653858
0.49661181026137463 -0.028246780656926096
0.5169409486931267 -0.028388189738625377
0.537270087124879 -0.028965645331952575
0.5546950629235237 -0.030125059096332674
0.573572120038722 -0.030993493775186315
0.5793804453049369 -0.029684228595871137
0.5939012584704744 -0.030553507508048375
0.6084220716360116 -0.029823949210923253
0.62729912875121 -0.02996563970373041
0.6454501452081316 -0.028363284855579846
0.6606969990319458 -0.029377771899412444
0.6803000968054211 -0.029664670523874948
0.70135527589545 -0.02980593890002027
0.7187802516940947 -0.030093259641144593
0.739109390125847 -0.03023466872284386
0.7550822846079381 -0.03153971273554111
0.7681510164569214 -0.03269997073324481
0.776863504356244 -0.03168084040613252
0.7899322362052275 -0.031969005380580395
0.7971926427879961 -0.03298504018550623
0.813165537270087 -0.032400549314482566
0.8349467570183928 -0.03268702582228328
0.8480154888673765 -0.03370193498277767
0.8654404646660212 -0.03384390688669267
0.8741529525653435 -0.033260823071208284
0.8785091965150047 -0.03151579279137311
0.8828654404646661 -0.029334715999909985
0.8850435624394966 -0.02729941016231796
0.887947725072604 -0.0251186147819627
0.8952081316553726 -0.021919533307818748
0.9010164569215875 -0.017993989058736143
0.9060987415295256 -0.014213934352416788
0.9126331074540173 -0.010433598234989544
0.9184414327202324 -0.005781309799860454
0.9242497579864473 -0.000983672527522031
0.9300580832526621 0.0032325693959791638
0.9373184898354305 0.007449092730588226
0.9431268151016456 0.011229288142461526
0.9496611810261375 0.01457357774826084
0.952565343659245 0.01588228010536029
0.9612778315585673 0.01588396857200744
0.9699903194578896 0.016321703550282514
0.9787028073572119 0.015015252482045938
0.9903194578896417 0.012982620049978577
1.0019361084220715 0.01327556901326006
0.9837850919651501 0.013998795560458334
]
plot(data[:, 1], -data[:, 2], "o", color="#A60628")
gca().invert_yaxis()
text(0.5, 0.015, "BEM", color="#348ABD")
text(0.5, 0.042, "experimental", color="#A60628")
xlim([0, 1])
ylim([0.1, -0.05])
gca().yaxis.set_label_coords(-0.08,1.05)
# savefig("hoveru1.pdf")

# CT = 0.0076
op = simple_op.(Vinf, Omega, r, rho, pitch=11.8*pi/180)
outputs = solvehover.(Ref(rotor), sections, op)

T, Q = thrusttorque(rotor, sections, outputs)
A = pi*Rtip^2
CT = T / (rho * A * (Omega * Rtip)^2)
println(CT)

figure()
plot(r/Rtip, outputs.u/(Omega*Rtip), color="#348ABD")
xlabel(L"r/R")
ylabel(L"u/(\Omega R)", rotation=0)

data = [
0.19326598413451757 -0.008918299545731909
0.20112550003390056 -0.011346396365855313
0.20542850362736453 -0.01491436707573393
0.20972642213031384 -0.017769069089429795
0.22470811580446126 -0.01920130178317174
0.2147321852328971 -0.01991077361177028
0.22757610685470198 -0.021484846430266435
0.23614143331751297 -0.022914638280561395
0.23758356498745675 -0.025197640518001233
0.2432809004000271 -0.02434388772120144
0.2525733948064275 -0.027771103125635632
0.25826056003796866 -0.025490812936470264
0.26537968675842427 -0.024066987592379144
0.27110956675028813 -0.02777815445114923
0.28323547359142986 -0.02863868736863516
0.29180385110854967 -0.030496440436639766
0.306067530002034 -0.031215133229371472
0.3146318394467421 -0.032502271340429864
0.3274737270323411 -0.03379103668045291
0.34173130381720784 -0.03365380703776527
0.35244457251339073 -0.03636829615567157
0.3659963387348294 -0.037229371482812385
0.37598955861414324 -0.03894501322123532
0.39026340768865686 -0.04109024340633262
0.4045230185097294 -0.0412383212421181
0.41735270187809337 -0.04081524171130245
0.43020170859041273 -0.043102583225981414
0.44304563021221766 -0.04467665604447758
0.4630178317174045 -0.04611078717201167
0.48156213980608853 -0.04725906841141772
0.4986846565868872 -0.048977422198115114
0.5257800528849412 -0.04955834293850429
0.5514516238389042 -0.05042402874771171
0.577113024611838 -0.04986317716455352
0.5956512305919045 -0.050155535968540235
0.6198918570750558 -0.05030741067190994
0.6441436707573395 -0.05202847650688182
0.6826459420977692 -0.05261373652451014
0.70404298596515 -0.05390575632246254
0.7211563495830225 -0.05434022645603091
0.7439751847582886 -0.055062173706691975
0.7610936334666756 -0.056209912536443155
0.7767852057766627 -0.05721445521730287
0.7981822496440435 -0.05850647501525526
0.8081795375957691 -0.06079273171062444
0.816750966167198 -0.06307844599633873
0.8267574072818495 -0.06664858634483693
0.841019052139128 -0.06708197165909553
0.8524096548918569 -0.06480385110854973
0.8609424367753749 -0.06166872330327479
0.8637738151739102 -0.05881673333785341
0.8680188487355073 -0.054253440911248216
0.8708329378262931 -0.04897633737880533
0.8729239270459013 -0.042272425249169424
0.879988134788799 -0.033145297986304154
0.8849298257508984 -0.02629981693674148
0.8912983931113971 -0.01959753203607023
0.8969519967455417 -0.012609668452098446
0.9033215811241437 -0.006050037290663765
0.9089792528307004 0.0003672113363617918
0.9146409926096684 0.006213845006441114
0.9210166790968879 0.01191755373245644
0.926677401857753 0.01790684114177233
0.9330551223811783 0.023325242389314535
0.9394348769408092 0.028458336158383626
0.9443948742287611 0.032736049901688255
0.9486602481524171 0.03444626754356228
0.9543535154925755 0.035870635297308294
0.9629005356295339 0.037008610753271415
0.9728856193640245 0.03643419892874094
0.9786032951386534 0.03443487694080956
0.9885944809817613 0.03300454268085973
0.9957339480642755 0.03157529324021968
0.998601939114516 0.029291748593124965
]
plot(data[:, 1], -data[:, 2], "o", color="#A60628")
gca().invert_yaxis()
text(0.5, 0.068, "BEM", color="#348ABD")
text(0.5, 0.04, "experimental", color="#A60628")
xlim([0, 1])
ylim([0.1, -0.05])
gca().yaxis.set_label_coords(-0.08,1.05)
# savefig("hoveru2.pdf")


text:
using LinearAlgebra, Images, FileIO, Colors

function sinkhorn( a, b, precision, costMatrix, n_iter )
	n, m = length( a ), length( b )
	C = costMatrix
	costMatrix = costMatrix ./ max( costMatrix... ) # renormalizing to avoid underflows in division
	K = exp.( -costMatrix/precision )
	u, v = ones( n ), ones( m )
	for i in 1:n_iter
		# u = a ./ K*v
		u = exp.( log.( a ) - log.( K*v ) )
		v = exp.( log.( b ) - log.( K'*u ) )
	end
	# P = [ u[i]*K[i,j]*v[j] for i in 1:n, j in 1:m ]
	P = u.*K*diagm( v )
	entropy = -sum( P.*log.( P ) )
	return P, sum( P.*C ) - precision*entropy
end

function getImage( path::String, quantization::Int )
	execstring =`convert $path -dither None -colors $quantization $path.temp.png`
	run(execstring)
	img = load( "$path.temp.png" )
	rm( "$path.temp.png" )
	return img
end

function swapPalette( paletteSourceImage::String, paletteTargetImage::String, maxColors=10_000, mode="mean" )
	source = RGB.( getImage( paletteSourceImage, maxColors ) )
	target = RGB.( getImage( paletteTargetImage, maxColors ) )

	scols = unique( source )
	tcols = unique( target )

	# we don't account for frequency of colors, so the distributions are uniform
	scdist = ones( length( scols ) )/ length( scols )
	tcdist = ones( length( tcols ) )/ length( tcols )

	costMatrix = [ colordiff( t, s ) for t in tcols, s in scols ]

	P, d = sinkhorn( tcdist, scdist, 0.01, costMatrix, 100 )

	if mode=="matching"
		perm = [ argmax( P[i,:] ) for i in 1:length( tcols ) ]
		mp = [ tcols[i] => scols[perm][i] for i in 1:length( tcols ) ]
	elseif mode=="mean"
		scols = LCHab.( scols )
		scols = [ [a.l, a.c, a.h] for a in scols ]
		mp = [tcols[i] => LCHab( sum( scols.*P[i,:]/sum( P[i,:] ) )... ) for i in 1:length( tcols )]
	else
		@error "Unknown mode:\t$mode."
	end
	for i in 1:length( tcols )  replace!( target, mp[i] )  end
	return target
end

nargs = length( ARGS )
if 3 <= nargs<= 5
	maxColors = nargs>=4 ? ARGS[4] : 10_000
	mode      = nargs==5 ? ARGS[5] : "mean"
	@time img = swapPalette( ARGS[1], ARGS[2], maxColors, mode )
	save( ARGS[3], img )
else
	println( "Usage: julia paletteSwap.jl paletteSourceImage paletteTargetImage outputFile [maxColors=10_000 mode='mean']" )
end


text:
using JuMP, EAGO, Gurobi, CSV, DataFrames, Ipopt
#m = Model(with_optimizer(Ipopt.Optimizer))

m = Model(with_optimizer(EAGO.Optimizer, relaxed_optimizer = Gurobi.Optimizer(OutputFlag=0)))

x_Idx = Any[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24]
@variable(m, x[x_Idx])
JuMP.set_lower_bound(x[21], 0.0)
JuMP.set_lower_bound(x[1], -0.5)
JuMP.set_upper_bound(x[1], 0.5)
JuMP.set_lower_bound(x[2], 5.4)
JuMP.set_upper_bound(x[2], 6.4)
JuMP.set_lower_bound(x[3], 0.4)
JuMP.set_upper_bound(x[3], 1.4)
JuMP.set_lower_bound(x[4], 4.9)
JuMP.set_upper_bound(x[4], 5.9)
JuMP.set_lower_bound(x[5], 1.3)
JuMP.set_upper_bound(x[5], 2.3)
JuMP.set_lower_bound(x[6], 3.9)
JuMP.set_upper_bound(x[6], 4.9)
JuMP.set_lower_bound(x[7], 2.1)
JuMP.set_upper_bound(x[7], 3.1)
JuMP.set_lower_bound(x[8], 4.1)
JuMP.set_upper_bound(x[8], 5.1)
JuMP.set_lower_bound(x[9], 2.8)
JuMP.set_upper_bound(x[9], 3.8)
JuMP.set_lower_bound(x[10], 3.0)
JuMP.set_upper_bound(x[10], 4.0)
JuMP.set_lower_bound(x[11], 3.9)
JuMP.set_upper_bound(x[11], 4.9)
JuMP.set_lower_bound(x[12], 3.2)
JuMP.set_upper_bound(x[12], 4.2)
JuMP.set_lower_bound(x[13], 4.7)
JuMP.set_upper_bound(x[13], 5.7)
JuMP.set_lower_bound(x[14], 2.3)
JuMP.set_upper_bound(x[14], 3.3)
JuMP.set_lower_bound(x[15], 5.6)
JuMP.set_upper_bound(x[15], 6.6)
JuMP.set_lower_bound(x[16], 2.3)
JuMP.set_upper_bound(x[16], 3.3)
JuMP.set_lower_bound(x[17], 6.0)
JuMP.set_upper_bound(x[17], 7.0)
JuMP.set_lower_bound(x[18], 1.9)
JuMP.set_upper_bound(x[18], 2.9)
JuMP.set_lower_bound(x[19], 6.9)
JuMP.set_upper_bound(x[19], 7.9)
JuMP.set_lower_bound(x[20], 1.0)
JuMP.set_upper_bound(x[20], 2.0)
JuMP.set_upper_bound(x[21], 10.0)
JuMP.set_lower_bound(x[22], -2.0)
JuMP.set_upper_bound(x[22], 2.0)
JuMP.set_lower_bound(x[23], -2.0)
JuMP.set_upper_bound(x[23], 2.0)
JuMP.set_lower_bound(x[24], -2.0)
JuMP.set_upper_bound(x[24], 2.0)

# ----- Constraints ----- #
@NLconstraint(m, e2, x[22]*x[1]+ (x[1])^2*x[23]+ (x[1])^3*x[24]-x[2]+x[21] == 0.0)
@NLconstraint(m, e3, x[22]*x[3]+ (x[3])^2*x[23]+ (x[3])^3*x[24]-x[4]+x[21] == 0.0)
@NLconstraint(m, e4, x[22]*x[5]+ (x[5])^2*x[23]+ (x[5])^3*x[24]-x[6]+x[21] == 0.0)
@NLconstraint(m, e5, x[22]*x[7]+ (x[7])^2*x[23]+ (x[7])^3*x[24]-x[8]+x[21] == 0.0)
@NLconstraint(m, e6, x[22]*x[9]+ (x[9])^2*x[23]+ (x[9])^3*x[24]-x[10]+x[21] == 0.0)
@NLconstraint(m, e7, x[22]*x[11]+ (x[11])^2*x[23]+ (x[11])^3*x[24]-x[12]+x[21] == 0.0)
@NLconstraint(m, e8, x[22]*x[13]+ (x[13])^2*x[23]+ (x[13])^3*x[24]-x[14]+x[21] == 0.0)
@NLconstraint(m, e9, x[22]*x[15]+ (x[15])^2*x[23]+ (x[15])^3*x[24]-x[16]+x[21] == 0.0)
@NLconstraint(m, e10, x[22]*x[17]+ (x[17])^2*x[23]+ (x[17])^3*x[24]-x[18]+x[21] == 0.0)
@NLconstraint(m, e11, x[22]*x[19]+ (x[19])^2*x[23]+ (x[19])^3*x[24]-x[20]+x[21] == 0.0)

# ----- Objective ----- #
@NLobjective(m, Min, ( (x[1])^2+ ((-5.9)+x[2])^2+ ((-0.9)+x[3])^2+ ((-5.4)+x[4])^2+ ((-1.8)+x[5])^2+ ((-4.4)+x[6])^2+ ((-2.6)+x[7])^2+ ((-4.6)+x[8])^2+ ((-3.3)+x[9])^2+ ((-3.5)+x[10])^2+ ((-4.4)+x[11])^2+ ((-3.7)+x[12])^2+ ((-5.2)+x[13])^2+ ((-2.8)+x[14])^2+ ((-6.1)+x[15])^2+ ((-2.8)+x[16])^2+ ((-6.5)+x[17])^2+ ((-2.4)+x[18])^2+ ((-7.4)+x[19])^2+ ((-1.5)+x[20])^2))

JuMP.optimize!(m)

backend_opt = backend(m).optimizer.model.optimizer
last_relaxed_opt = backend_opt.relaxed_optimizer
run_time = backend_opt._run_time
println("run time: $run_time")


text:
using CatViews
using Test
using Random: randn!

@testset "CatView tests" begin

    @testset "Array properties" begin
        x = CatView([1,2,3,4],[5,6,7,8])
        @test length(x) == 8
        @test ndims(x) == 1
        [ @test x[i] == i for i in 1:length(x) ]
        @test x[end] == 8
        @test_throws BoundsError x[0]
        @test_throws BoundsError x[9]
    end

    @testset "mutations" begin
        a = randn(10,10)
        b = randn(9,9)

        x = CatView(view(a,:,2),view(b,:,3))
        @test length(x) == 19
        [ @test x[i] == a[i,2] for i in 1:size(a,1) ]
        [ @test x[size(a,1)+i] == b[i,3] for i in 1:size(b,1) ]

        ## Test mutates original object
        x[1] = 55
        @test a[1,2] == 55.0
        x[1:3] = [55, 56, 57]
        @test a[1,2] == 55.0
        @test a[2,2] == 56.0
        @test a[3,2] == 57.0

        randn!(x)
        @test isapprox(a[:,2],x[1:10])
        @test isapprox(b[:,3],x[11:end])
    end

    @testset "Matrix multiply" begin
        a = randn(10,10)
        b = randn(9,9)
        x = CatView(view(a,:,2),view(b,:,3))

        c = vcat(a[:,2],b[:,3])
        A = randn(19,19)
        @test isapprox(A*c,A*x)
    end

    @testset "Other constructors" begin
        x = CatView([1,2,3,4],[5,6,7,8])
        x = CatView(x,[9,10,11,12])
        @test length(x) == 12
        [ @test x[i] == i for i in 1:length(x) ]

        y = CatView([1,2,3,4],[5,6,7,8],[9,10,11,12])
        @test length(y) == 12
        [ @test y[i] == i for i in 1:length(y) ]

        z = ([1,2,3], @view([0.5, 0.3, 0.2][1:2]))
        @test vcat(z...) == collect(CatView(z))
    end

    @testset "Iteration" begin
        x = CatView([1,2],[3,4])
        @test all(x .== [1,2,3,4])
    end
end

@testset "splitview tests" begin
    @testset "without preallocation" begin

        x,(A,B),s,e = splitview(Int64,(2,3),(4,5))

        # test sizes
        @test length(x) == 26
        @test size(A) == (2,3)
        @test size(B) == (4,5)
        @test eltype(A) == Int64
        @test eltype(B) == Int64

        # test partititions between vectors
        @test s == (1,7)
        @test e == (6,26)

        # set A to 1:6 by mutating x
        [ x[i] = i for i = s[1]:e[1] ]
        @test all(A == [1 3 5; 2 4 6])

        # set B to 1:20 by mutating it directly
        [ B[i] = i for i in 1:length(B) ]
        [ @test x[i+6] == i for i in 1:20 ]

        # test on a higher-order array
        sz1 = (2,2,2)
        sz2 = (2,3,2)
        sz3 = (3,2,2)
        x,(A,B,C),s,e = splitview(Int64,sz1,sz2,sz3)

        @test size(A) == sz1
        @test size(B) == sz2
        @test size(C) == sz3
        @test s == (1, length(A)+1, length(A)+length(B)+1)
        @test e == (length(A), length(A)+length(B), length(A)+length(B)+length(C))
    end

    @testset "with preallocation" begin

        # preallocate x
        x = collect(1:26)
        (A,B),s,e = splitview(x,(2,3),(4,5))

        #test sizes
        @test size(A) == (2,3)
        @test size(B) == (4,5)
        @test s == (1,7)
        @test e == (6,26)

        # test that A and B faithfully view into x
        for i = 1:6
            @test A[i] == i
            @test B[i] == 6+i
        end

        # test on higher-order array
        sz1 = (2,2,2)
        sz2 = (2,3,2)
        sz3 = (3,2,2)
        N = prod(sz1)+prod(sz2)+prod(sz3)
        x = collect(1:N)
        (A,B,C),s,e = splitview(x,sz1,sz2,sz3)

        @test eltype(A) == Int
        @test size(A) == sz1
        @test size(B) == sz2
        @test size(C) == sz3
        @test s == (1, length(A)+1, length(A)+length(B)+1)
        @test e == (length(A), length(A)+length(B), length(A)+length(B)+length(C))

        # test that arrays correctly hold 1:N
        for (i,j,X) in zip(s,e,(A,B,C))
            for (xi,Xi) in zip(i:j,eachindex(X))
                @test x[xi] == X[Xi] == xi
            end
        end

    end
end

@testset "indexing" begin

    x,children, = splitview(Int64,(3,3),(2,2,2,2,2,2),(4,5,2),(10,),(1,10),(3,3))
    copyto!(x,1:length(x))

    for child in children
        for child_idx = eachindex(child)
            parent_idx = vecidx(child, child_idx)
            @test x[parent_idx] == child[child_idx]
        end
    end

    for child in children
        d = size(child)
        for i = 1:prod(d)
            idx = Tuple(CartesianIndices(d)[i])
            pi1 = vecidx(child, idx)
            pi2 = vecidx(child, idx...)
            @test x[pi1] == x[pi2] == child[idx...] == child[i]
        end
    end

end

@testset "mapreduce" begin
    orig = collect(1:100)
    xs = CatView(orig[21:80], orig[1:20], orig[81:end])

    @test mapreduce(x->2x, +, xs) ==  sum(2*orig)
    @test sum(xs) == sum(orig)
end


text:
struct UranusBarycenter <: abstractBarycenter end

naif_id(::Type{UranusBarycenter}) = 7
parent(::Type{UranusBarycenter}) = SSB


text:
"""
$(SIGNATURES)

Construct the matrix of a random 2D Ising lattice ``∑ᵢⱼJᵢⱼZᵢZⱼ``, where ``Jᵢⱼ``s are uniformly distributed between [-1, 1). Generate sparse matrix when `sp` is set to true.
"""
function random_ising(num_qubits::Integer; sp=false)
    J = Vector{Float64}()
    idx = Vector{Vector{Int64}}()
    for i ∈ 1:num_qubits
        for j ∈ i+1:num_qubits
            push!(J, 2*rand()-1)
            push!(idx, [i, j])
        end
    end
    two_local_term(J, idx, num_qubits, sp=sp)
end

"""
$(SIGNATURES)

Construct the matrix of the alternating sectors chain model described in (https://www.nature.com/articles/s41467-018-05239-9#Sec1). Generate sparse matrix when `sp` is set to true.
"""
function alt_sec_chain(w1, w2, n, num_qubits; sp=false)
    J = Vector{Float64}()
    idx = Vector{Vector{Int64}}()
    for i in 1:num_qubits-1
        push!(J, isodd(ceil(Int, i/n)) ? w1 : w2)
        push!(idx, [i, i+1])
    end
    two_local_term(J, idx, num_qubits, sp=sp)
end

function build_example_hamiltonian(num_qubits; sp=false)
    A = (s) -> (1 - s)
    B = (s) -> s
    if num_qubits == 1
        sp ? SparseHamiltonian([A, B], [spσx, spσz]) : DenseHamiltonian([A, B], [σx, σz])
    elseif num_qubits == 4
        Hd = -standard_driver(4, sp = sp)
        Hp = -q_translate("ZIII+IZII+IIZI+IIIZ", sp = sp)
        sp ? SparseHamiltonian([A, B], [Hd, Hp]) : DenseHamiltonian([A, B], [Hd, Hp])
    end
end


text:
using Distributions, StatsBase
using JLD, PyPlot

for set in ["set1", "set2", "set3", "set4"]
    file = jldopen("res/power_$(set).jld", "r")
    power = read(file, "power")
    close(file)

    fig = figure(figsize=(3.14961, 3.14961), dpi=1000)

    scatter(-.75:.05:.75, power[1, :], s=5, marker="o")
    scatter(-.75:.05:.75, power[2, :], s=5, marker="v")

    minorticks_on()
    grid()
    grid(which="minor", ls="dotted", lw=".5")

    xlabel("δ", size="xx-small")
    ylabel("proportion of rejections", size="xx-small")

    ylim(0, 1.05)

    tick_params("both", labelsize="xx-small")

    tight_layout()

    savefig("res/exper2_$(set).pdf")

    close(fig)
end


text:
## Test grenander and related methods ##
module Test_grenander

using MultipleTesting
using Test

@testset "Grenander estimators" begin

    @testset "Grenander estimator" begin

        pv = [0.2, 0.4, 0.6, 0.8, 1.0]
        p, f, F = MultipleTesting.grenander(pv)
        @test p ≈ pv
        @test f ≈ fill(1.0, size(p))
        @test F ≈ pv

        pv = [0.1, 0.3, 0.5, 0.7, 0.9]
        p, f, F = MultipleTesting.grenander(pv)
        @test p ≈ pv
        @test f ≈ fill(1.0, size(p))
        @test F ≈ pv .+ 0.1

        pv = [0.1, 0.1, 0.3, 0.3, 0.5, 0.5, 0.7, 0.7, 0.9, 0.9]
        p, f, F = MultipleTesting.grenander(pv)
        @test p ≈ unique(pv)
        @test f ≈ fill(1.0, size(p))
        @test F ≈ p .+ 0.1

        pv = [0.1, 0.2, 0.5, 0.6, 0.9]
        p, f, F = MultipleTesting.grenander(pv)
        @test p ≈ pv
        @test f ≈ [2.0, 1.0, 1.0, 2 / 3, 2 / 3]
        @test F ≈ [0.2, 0.4, 0.7, 0.8, 1.0]

    end

    @testset "Isotonic regression" begin

        function isotonic_regression_reference(y::AbstractVector{T}, w::AbstractVector{T}) where T <: AbstractFloat
            y = copy(y)
            w = copy(w)
            cnts = fill(1, length(y))
            i = 2
            while i <= length(y)
                if y[i] < y[i - 1]
                    y[i - 1] = (w[i] * y[i] + w[i - 1] * y[i - 1]) / (w[i] + w[i - 1])
                    w[i - 1] = w[i] + w[i - 1]
                    cnts[i - 1] += cnts[i]
                    deleteat!(y, i)
                    deleteat!(w, i)
                    deleteat!(cnts, i)
                    i = max(i - 2, 1)
                end
                i += 1
            end
            yisotonic = vcat([y[idx] .* fill(1.0, cnt) for (idx, cnt) in enumerate(cnts)]...)
            return yisotonic
        end

        function isotonic_regression_reference(y::AbstractVector{T}) where T <: AbstractFloat
            isotonic_regression_reference(y, fill(1.0, size(y)))
        end

        # pooled adjacent violators example: Robertson, page 10
        @test isotonic_regression_reference([22.5; 23.333; 20.833; 24.25], [3.0;3.0;3.0;2.0]) ≈ [22.222; 22.222; 22.222; 24.25]

        # if input already ordered, then output should be the same
        @test isotonic_regression_reference([1.0; 2.0; 3.0]) ≈ [1.0; 2.0; 3.0]
        @test isotonic_regression_reference([1., 41., 51., 1., 2., 5., 24.], [1., 2., 3., 4., 5., 6., 7.]) ≈ [1.0, 13.95, 13.95, 13.95, 13.95, 13.95, 24.]

        # single value or empty vector remains unchanged
        r = rand(1)
        @test MultipleTesting.isotonic_regression(r) == r
        r = Vector{Float64}()
        @test MultipleTesting.isotonic_regression(r) == r

        r = rand(10)
        @test MultipleTesting.isotonic_regression(r, fill(1.0, size(r))) == MultipleTesting.isotonic_regression(r)
        @test_throws DimensionMismatch MultipleTesting.isotonic_regression(rand(10), ones(5))

    end

end

end