text:
"""
`macbeth(colorchart::Symbol,color::Symbol)`

returns `RSpec` of one of the 24 MacBeth color patches of a given name.

available color charts:

`:Macbeth_chart1`, λ=380:10:730 nm, average of multiple measured charts from BabelColor

`:Macbeth_chart2`, λ=380:4:780 nm

`:Macbeth_chart3`, λ=380:5:780 nm

available color names:

:darkskin, :lightskin, :bluesky, :foliage, :blueflower, :bluishgreen,

:orange, :purplishblue, :moderatered, :purple, :yellowgreen, :orangeyellow,

:blue, :green, :red, :yellow, :magenta, :cyan,

:white95, :neutral8, :neutral65, :neutral5, :neutral35, :black2
"""
function macbeth(colorchart::Symbol,color::Symbol)
    chart = eval(colorchart)
    color == :darkskin ? reflectance_spec(chart[:,1],chart[:,2]) :
    color == :lightskin ? reflectance_spec(chart[:,1],chart[:,3]) :
    color == :bluesky ? reflectance_spec(chart[:,1],chart[:,4]) :
    color == :foliage ? reflectance_spec(chart[:,1],chart[:,5]) :
    color == :blueflower ? reflectance_spec(chart[:,1],chart[:,6]) :
    color == :bluishgreen ? reflectance_spec(chart[:,1],chart[:,7]) :
    color == :orange ? reflectance_spec(chart[:,1],chart[:,8]) :
    color == :purplishblue ? reflectance_spec(chart[:,1],chart[:,9]) :
    color == :moderatered ? reflectance_spec(chart[:,1],chart[:,10]) :
    color == :purple ? reflectance_spec(chart[:,1],chart[:,11]) :
    color == :yellowgreen ? reflectance_spec(chart[:,1],chart[:,12]) :
    color == :orangeyellow ? reflectance_spec(chart[:,1],chart[:,13]) :
    color == :blue ? reflectance_spec(chart[:,1],chart[:,14]) :
    color == :green ? reflectance_spec(chart[:,1],chart[:,15]) :
    color == :red ? reflectance_spec(chart[:,1],chart[:,16]) :
    color == :yellow ? reflectance_spec(chart[:,1],chart[:,17]) :
    color == :magenta ? reflectance_spec(chart[:,1],chart[:,18]) :
    color == :cyan ? reflectance_spec(chart[:,1],chart[:,19]) :
    color == :white95 ? reflectance_spec(chart[:,1],chart[:,20]) :
    color == :neutral8 ? reflectance_spec(chart[:,1],chart[:,21]) :
    color == :neutral65 ? reflectance_spec(chart[:,1],chart[:,22]) :
    color == :neutral5 ? reflectance_spec(chart[:,1],chart[:,23]) :
    color == :neutral35 ? reflectance_spec(chart[:,1],chart[:,24]) :
    color == :black2 ? reflectance_spec(chart[:,1],chart[:,25]) : error("Color does not exist.")
end

"""
`macbeth(colorchart::Symbol,num::Int)`

returns `RSpec` of one of the 24 MacBeth color patches of a given number.

available color charts:

`:Macbeth_chart1`, λ=380:10:730 nm, average of multiple measured charts from BabelColor

`:Macbeth_chart2`, λ=380:4:780 nm

`:Macbeth_chart3`, λ=380:5:780 nm

available color patch numbers: 1 to 24
"""
function macbeth(colorchart::Symbol,num::Int)
    colorchart in (:Macbeth_chart1, :Macbeth_chart2, :Macbeth_chart3) ? nothing : DomainError("$colorchart does not exist.")
    if num < 1 || num > 24
        DomainError("Color index $num does not exist. Index must be in the range 1 to 24")
    else
        i=[:darkskin, :lightskin, :bluesky, :foliage, :blueflower, :bluishgreen, :orange, :purplishblue, :moderatered, :purple, :yellowgreen, :orangeyellow, :blue, :green, :red, :yellow, :magenta, :cyan, :white95, :neutral8, :neutral65, :neutral5, :neutral35, :black2]
        macbeth(colorchart,i[num])
    end
end

"""
`mbspec(colorchart::Symbol,name::Symbol)`

returns `RSpec` adapted to SPECENV limits for λ.
"""
function mbspec(colorchart::Symbol,name::Symbol)
    macbethspec=macbeth(colorchart,name)
    specspline=spline3(macbethspec.λ,macbethspec.r,:natural)
    extrapspline=extrap(specspline,SPECENV.λmin:SPECENV.λmax,:linear)
    lam,s=interp(extrapspline,SPECENV.λmin:SPECENV.Δλ:SPECENV.λmax)
    reflectance_spec(lam,s)
end

"""
`mbspec(colorchart::Symbol, num::Int)`

returns `RSpec` adapted to th 390:830 nm range of the standard CMF
"""
function mbspec(colorchart::Symbol, num::Int)
    macbethspec=macbeth(colorchart,num)
    specspline=spline3(macbethspec.λ,macbethspec.r,:natural)
    extrapspline=extrap(specspline,SPECENV.λmin:SPECENV.Δλ:SPECENV.λmax,:linear)
    lam,s=interp(extrapspline,SPECENV.λmin:SPECENV.Δλ:SPECENV.λmax)
    reflectance_spec(lam,s)
end

"""
`mbcolor(spec::RSpec,illum::LSpec)`

returns XYZ color according to illuminant
"""
function mbcolor(spec::RSpec,illum::LSpec)
    illum*spec*cmf(CMF2012_10)
end

# dark skin
# light skin
# blue sky
# foliage
# blue flower
# bluish green
# orange
# purplish blue
# moderate red
# purple
# yellow green
# orange yellow
# blue
# green
# red
# yellow
# magenta
# cyan
# white 9.5 (.05 D)
# neutral 8 (.23 D)
# neutral 6.5 (.44 D)
# neutral 5 (.70 D)
# neutral 3.5 (1.05 D)
# black 2 (1.5 D)

"""
BabelColor Color Checker spectra, average of 20 color checker measurements
"""
const Macbeth_chart1=
   [380  0.055  0.121  0.141  0.051  0.158  0.145  0.053  0.132  0.098  0.095  0.060  0.062  0.064  0.051  0.049  0.056  0.155  0.114  0.199  0.182  0.152  0.109  0.068  0.031;
    390  0.058  0.148  0.184  0.054  0.209  0.185  0.053  0.171  0.116  0.119  0.061  0.062  0.074  0.053  0.048  0.053  0.202  0.145  0.259  0.240  0.197  0.133  0.076  0.032;
    400  0.061  0.180  0.254  0.055  0.300  0.250  0.053  0.233  0.131  0.148  0.062  0.063  0.094  0.054  0.047  0.052  0.284  0.192  0.421  0.367  0.272  0.163  0.083  0.032;
    410  0.062  0.197  0.307  0.056  0.380  0.299  0.053  0.290  0.136  0.172  0.063  0.064  0.138  0.055  0.047  0.052  0.346  0.235  0.660  0.506  0.330  0.181  0.086  0.032;
    420  0.062  0.201  0.325  0.057  0.412  0.323  0.053  0.329  0.134  0.182  0.064  0.064  0.192  0.057  0.047  0.052  0.362  0.259  0.811  0.566  0.349  0.187  0.088  0.033;
    430  0.062  0.204  0.331  0.059  0.425  0.340  0.054  0.362  0.132  0.176  0.066  0.065  0.239  0.059  0.047  0.054  0.355  0.284  0.863  0.581  0.356  0.191  0.090  0.033;
    440  0.061  0.208  0.334  0.060  0.429  0.358  0.054  0.387  0.131  0.160  0.070  0.066  0.280  0.062  0.047  0.056  0.334  0.316  0.877  0.586  0.360  0.194  0.091  0.033;
    450  0.061  0.216  0.333  0.062  0.429  0.383  0.054  0.399  0.129  0.139  0.075  0.066  0.311  0.066  0.046  0.059  0.305  0.352  0.884  0.587  0.361  0.195  0.091  0.032;
    460  0.061  0.229  0.327  0.062  0.422  0.420  0.055  0.390  0.126  0.118  0.086  0.067  0.312  0.074  0.045  0.066  0.275  0.390  0.890  0.588  0.361  0.194  0.090  0.032;
    470  0.061  0.250  0.314  0.063  0.405  0.466  0.056  0.359  0.121  0.100  0.104  0.069  0.281  0.092  0.045  0.081  0.246  0.426  0.894  0.587  0.359  0.193  0.090  0.032;
    480  0.061  0.277  0.300  0.065  0.381  0.509  0.058  0.313  0.116  0.085  0.137  0.072  0.231  0.123  0.044  0.108  0.217  0.446  0.897  0.585  0.357  0.192  0.089  0.032;
    490  0.062  0.304  0.287  0.068  0.348  0.545  0.060  0.257  0.111  0.074  0.190  0.077  0.175  0.176  0.044  0.154  0.189  0.444  0.901  0.585  0.356  0.192  0.089  0.032;
    500  0.065  0.325  0.271  0.076  0.313  0.567  0.067  0.207  0.106  0.065  0.269  0.090  0.126  0.244  0.044  0.228  0.167  0.423  0.905  0.586  0.356  0.192  0.089  0.032;
    510  0.070  0.330  0.251  0.103  0.282  0.575  0.088  0.167  0.101  0.059  0.377  0.129  0.090  0.306  0.045  0.339  0.148  0.384  0.906  0.586  0.357  0.192  0.089  0.032;
    520  0.076  0.314  0.231  0.146  0.254  0.571  0.123  0.137  0.096  0.056  0.478  0.208  0.066  0.338  0.045  0.464  0.126  0.335  0.908  0.587  0.358  0.192  0.089  0.032;
    530  0.079  0.289  0.215  0.177  0.229  0.553  0.151  0.117  0.093  0.053  0.533  0.303  0.052  0.334  0.046  0.557  0.107  0.280  0.907  0.586  0.358  0.192  0.089  0.032;
    540  0.080  0.277  0.200  0.183  0.214  0.524  0.172  0.105  0.093  0.051  0.551  0.378  0.045  0.316  0.047  0.613  0.099  0.229  0.907  0.586  0.358  0.192  0.089  0.032;
    550  0.084  0.279  0.184  0.170  0.208  0.488  0.200  0.097  0.094  0.051  0.547  0.426  0.041  0.293  0.049  0.647  0.101  0.183  0.910  0.586  0.358  0.192  0.089  0.032;
    560  0.090  0.280  0.167  0.149  0.202  0.443  0.250  0.090  0.097  0.052  0.529  0.468  0.039  0.261  0.052  0.670  0.103  0.144  0.910  0.585  0.358  0.192  0.089  0.032;
    570  0.102  0.294  0.156  0.132  0.195  0.398  0.338  0.086  0.109  0.052  0.505  0.521  0.038  0.228  0.058  0.692  0.109  0.117  0.912  0.587  0.359  0.192  0.089  0.032;
    580  0.119  0.344  0.150  0.121  0.194  0.349  0.445  0.084  0.157  0.051  0.471  0.575  0.038  0.196  0.071  0.709  0.136  0.100  0.912  0.587  0.359  0.192  0.089  0.031;
    590  0.134  0.424  0.147  0.114  0.202  0.299  0.536  0.084  0.268  0.053  0.427  0.613  0.038  0.164  0.103  0.722  0.199  0.089  0.912  0.586  0.359  0.191  0.088  0.031;
    600  0.143  0.489  0.144  0.109  0.216  0.253  0.583  0.084  0.401  0.059  0.381  0.631  0.038  0.134  0.177  0.731  0.290  0.081  0.910  0.584  0.357  0.190  0.088  0.031;
    610  0.147  0.523  0.141  0.104  0.231  0.223  0.590  0.084  0.501  0.073  0.346  0.637  0.039  0.114  0.313  0.739  0.400  0.076  0.912  0.582  0.356  0.189  0.087  0.032;
    620  0.151  0.542  0.140  0.103  0.241  0.206  0.586  0.083  0.557  0.095  0.327  0.637  0.040  0.102  0.471  0.747  0.514  0.074  0.915  0.580  0.354  0.188  0.087  0.032;
    630  0.158  0.558  0.140  0.104  0.253  0.198  0.582  0.084  0.579  0.117  0.317  0.638  0.040  0.096  0.586  0.753  0.611  0.073  0.916  0.577  0.351  0.186  0.086  0.032;
    640  0.168  0.576  0.141  0.107  0.276  0.192  0.579  0.088  0.587  0.139  0.312  0.641  0.040  0.092  0.651  0.759  0.682  0.073  0.918  0.575  0.348  0.184  0.085  0.032;
    650  0.179  0.594  0.145  0.109  0.310  0.190  0.579  0.095  0.589  0.162  0.309  0.644  0.041  0.090  0.682  0.764  0.726  0.073  0.920  0.573  0.346  0.182  0.084  0.032;
    660  0.188  0.611  0.150  0.111  0.345  0.192  0.584  0.105  0.592  0.189  0.314  0.646  0.042  0.090  0.698  0.769  0.754  0.075  0.921  0.572  0.344  0.181  0.084  0.032;
    670  0.190  0.623  0.151  0.111  0.365  0.200  0.595  0.117  0.595  0.221  0.327  0.648  0.043  0.093  0.707  0.772  0.769  0.076  0.920  0.570  0.341  0.179  0.083  0.032;
    680  0.188  0.634  0.147  0.111  0.367  0.212  0.611  0.133  0.601  0.256  0.345  0.653  0.043  0.099  0.715  0.777  0.779  0.076  0.921  0.569  0.339  0.178  0.083  0.032;
    690  0.185  0.651  0.141  0.112  0.363  0.223  0.628  0.155  0.607  0.295  0.362  0.661  0.043  0.104  0.725  0.784  0.789  0.075  0.923  0.568  0.337  0.177  0.082  0.032;
    700  0.186  0.672  0.134  0.117  0.362  0.231  0.644  0.186  0.614  0.336  0.376  0.671  0.044  0.109  0.734  0.792  0.800  0.072  0.926  0.567  0.334  0.175  0.081  0.032;
    710  0.192  0.693  0.131  0.123  0.368  0.233  0.653  0.218  0.617  0.370  0.380  0.679  0.047  0.111  0.740  0.798  0.807  0.072  0.928  0.567  0.333  0.174  0.081  0.032;
    720  0.200  0.710  0.133  0.129  0.377  0.229  0.654  0.255  0.617  0.404  0.378  0.684  0.050  0.110  0.744  0.801  0.813  0.073  0.929  0.565  0.331  0.173  0.081  0.032;
    730  0.214  0.728  0.144  0.135  0.394  0.229  0.659  0.296  0.618  0.445  0.380  0.689  0.055  0.110  0.748  0.805  0.821  0.079  0.932  0.565  0.330  0.172  0.080  0.032]

# Macbeth color checker spectra from 380 to 780 nm in 4 nm steps.
# Measurement info:
# t=r n=101 o=380 s=4
# Reflectance (IE, signal DIVIDED by reference)
# Mon Oct 24 19:32:41 1994 |  Reflectance | Macbeth - Dark Skin
# Mon Oct 24 19:32:41 1994 |  Reflectance | Macbeth - Light Skin
# Mon Oct 24 19:32:41 1994 |  Reflectance | Macbeth - Blue Sky
# Mon Oct 24 19:32:41 1994 |  Reflectance | Macbeth - Foilage
# Mon Oct 24 19:32:41 1994 |  Reflectance | Macbeth - Blue Flower
# Mon Oct 24 19:32:41 1994 |  Reflectance | Macbeth - Bluish Green
# Mon Oct 24 19:32:41 1994 |  Reflectance | Macbeth - Orange
# Mon Oct 24 19:32:41 1994 |  Reflectance | Macbeth - Purplish Blue
# Mon Oct 24 19:32:41 1994 |  Reflectance | Macbeth - Moderate Red
# Mon Oct 24 19:32:41 1994 |  Reflectance | Macbeth - Purple
# Mon Oct 24 19:32:41 1994 |  Reflectance | Macbeth - Yellow Green
# Mon Oct 24 19:32:41 1994 |  Reflectance | Macbeth - Orange Yellow
# Mon Oct 24 19:32:41 1994 |  Reflectance | Macbeth - Blue
# Mon Oct 24 19:32:41 1994 |  Reflectance | Macbeth - Green
# Mon Oct 24 19:32:41 1994 |  Reflectance | Macbeth - Red
# Mon Oct 24 19:32:41 1994 |  Reflectance | Macbeth - Yellow
# Mon Oct 24 19:32:41 1994 |  Reflectance | Macbeth - Magenta
# Mon Oct 24 19:32:41 1994 |  Reflectance | Macbeth - Cyan
# Mon Oct 24 19:32:41 1994 |  Reflectance | Macbeth - White
# Mon Oct 24 19:32:41 1994 |  Reflectance | Macbeth - Neutral-8
# Mon Oct 24 19:32:41 1994 |  Reflectance | Macbeth - Neutral-6.5
# Mon Oct 24 19:32:41 1994 |  Reflectance | Macbeth - Neutral-5
# Mon Oct 24 19:32:41 1994 |  Reflectance | Macbeth - Neutral-3.5
# Mon Oct 24 19:32:41 1994 |  Reflectance | Macbeth - Black
const Macbeth_chart2=
 [380  6.4071e-02  1.2475e-01  1.1468e-01  5.0734e-02  1.2553e-01  1.2553e-01  8.9655e-02  9.6830e-02  1.2497e-01  9.4605e-02  8.1090e-02  9.1491e-02  6.2681e-02  6.6018e-02  7.6196e-02  9.1101e-02  1.2608e-01  1.0545e-01  2.0117e-01  1.7297e-01  1.5973e-01  1.2814e-01  8.4316e-02  2.9522e-02;
  384  6.4733e-02  1.3434e-01  1.2972e-01  5.0986e-02  1.4153e-01  1.3814e-01  8.5062e-02  1.0816e-01  1.2623e-01  1.0308e-01  7.8747e-02  8.6499e-02  6.7197e-02  6.8378e-02  7.6232e-02  7.6283e-02  1.3778e-01  1.1407e-01  2.0472e-01  1.8778e-01  1.7059e-01  1.3696e-01  8.6242e-02  2.9476e-02;
  388  6.5438e-02  1.4751e-01  1.5111e-01  5.1014e-02  1.6553e-01  1.5424e-01  8.0092e-02  1.2567e-01  1.2963e-01  1.2005e-01  7.9493e-02  8.1751e-02  7.3594e-02  6.9355e-02  7.3502e-02  7.6866e-02  1.5355e-01  1.3069e-01  2.3700e-01  2.1728e-01  1.9147e-01  1.4756e-01  8.9124e-02  2.9949e-02;
  392  6.8866e-02  1.6097e-01  1.8383e-01  5.1617e-02  2.0643e-01  1.8176e-01  7.7533e-02  1.5045e-01  1.4135e-01  1.3903e-01  7.8439e-02  7.9301e-02  8.4605e-02  7.2747e-02  6.9168e-02  7.5119e-02  1.8900e-01  1.5179e-01  2.8685e-01  2.6904e-01  2.2531e-01  1.6218e-01  9.2367e-02  3.0768e-02;
  396  7.3362e-02  1.7731e-01  2.2234e-01  5.1973e-02  2.5702e-01  2.1669e-01  7.7190e-02  1.8441e-01  1.5110e-01  1.5856e-01  7.8611e-02  7.6401e-02  9.7869e-02  7.0442e-02  6.8666e-02  7.1744e-02  2.2727e-01  1.7419e-01  3.5533e-01  3.2928e-01  2.6484e-01  1.8003e-01  9.5659e-02  3.1405e-02;
  400  7.6425e-02  1.9042e-01  2.5508e-01  5.1462e-02  3.0526e-01  2.4664e-01  7.7339e-02  2.1495e-01  1.5599e-01  1.7496e-01  7.8838e-02  7.6499e-02  1.1115e-01  7.1491e-02  6.3779e-02  6.9591e-02  2.6572e-01  1.9631e-01  4.4189e-01  3.9547e-01  2.9536e-01  1.8995e-01  9.6455e-02  3.1542e-02;
  404  7.8657e-02  2.0003e-01  2.8255e-01  5.2083e-02  3.4699e-01  2.7036e-01  7.4286e-02  2.4286e-01  1.5639e-01  1.8960e-01  7.7590e-02  7.5525e-02  1.2558e-01  7.2427e-02  6.4337e-02  6.9019e-02  2.9546e-01  2.1477e-01  5.3735e-01  4.6024e-01  3.2138e-01  1.9831e-01  9.7315e-02  3.1666e-02;
  408  7.8971e-02  2.0450e-01  3.0071e-01  5.2572e-02  3.7787e-01  2.8667e-01  7.2468e-02  2.6610e-01  1.5438e-01  2.0058e-01  7.8324e-02  7.5639e-02  1.4031e-01  7.3277e-02  6.5157e-02  6.5189e-02  3.1770e-01  2.3048e-01  6.3410e-01  5.1246e-01  3.3873e-01  2.0281e-01  9.8318e-02  3.1773e-02;
  412  7.7607e-02  2.0428e-01  3.1043e-01  5.2140e-02  3.9662e-01  2.9632e-01  7.1278e-02  2.8345e-01  1.5145e-01  2.0702e-01  7.6552e-02  7.4231e-02  1.5401e-01  7.2664e-02  6.3472e-02  6.2839e-02  3.3122e-01  2.4063e-01  7.1368e-01  5.4340e-01  3.4599e-01  2.0491e-01  9.8493e-02  3.1826e-02;
  416  7.5724e-02  2.0475e-01  3.1526e-01  5.2234e-02  4.0629e-01  3.0121e-01  6.8699e-02  2.9559e-01  1.4819e-01  2.1056e-01  7.6033e-02  7.3869e-02  1.6609e-01  7.2829e-02  6.2686e-02  6.4007e-02  3.3802e-01  2.4633e-01  7.6904e-01  5.5774e-01  3.4897e-01  2.0531e-01  9.8455e-02  3.1750e-02;
  420  7.4192e-02  2.0503e-01  3.1744e-01  5.3259e-02  4.1123e-01  3.0631e-01  6.8919e-02  3.0525e-01  1.4682e-01  2.1139e-01  7.7186e-02  7.3794e-02  1.7684e-01  7.3715e-02  6.1765e-02  6.1659e-02  3.4314e-01  2.5246e-01  8.0101e-01  5.6783e-01  3.5109e-01  2.0480e-01  9.8569e-02  3.1558e-02;
  424  7.2052e-02  2.0479e-01  3.1887e-01  5.3628e-02  4.1386e-01  3.1189e-01  6.8337e-02  3.1688e-01  1.4373e-01  2.0723e-01  7.7312e-02  7.2924e-02  1.9090e-01  7.4196e-02  6.1331e-02  6.1007e-02  3.4106e-01  2.6328e-01  8.1825e-01  5.6993e-01  3.5303e-01  2.0529e-01  9.8604e-02  3.1463e-02;
  428  7.0115e-02  2.0549e-01  3.1911e-01  5.3901e-02  4.1362e-01  3.1676e-01  6.8442e-02  3.2760e-01  1.4183e-01  2.0000e-01  7.7634e-02  7.2614e-02  2.0643e-01  7.4711e-02  6.1042e-02  5.9439e-02  3.3679e-01  2.7103e-01  8.2489e-01  5.6964e-01  3.5376e-01  2.0497e-01  9.8680e-02  3.1369e-02;
  432  6.8481e-02  2.0694e-01  3.1854e-01  5.4762e-02  4.1312e-01  3.2121e-01  6.7812e-02  3.3817e-01  1.4249e-01  1.9101e-01  7.8363e-02  7.2251e-02  2.2253e-01  7.5173e-02  5.9826e-02  5.9135e-02  3.3081e-01  2.8062e-01  8.2936e-01  5.7082e-01  3.5445e-01  2.0497e-01  9.8952e-02  3.1229e-02;
  436  6.6574e-02  2.0762e-01  3.1756e-01  5.5717e-02  4.1113e-01  3.2719e-01  6.6510e-02  3.4839e-01  1.4017e-01  1.8148e-01  7.8844e-02  7.1542e-02  2.4047e-01  7.6424e-02  5.8672e-02  6.0000e-02  3.2248e-01  2.9336e-01  8.3512e-01  5.7173e-01  3.5546e-01  2.0505e-01  9.9336e-02  3.1156e-02;
  440  6.5338e-02  2.0963e-01  3.1627e-01  5.6677e-02  4.0753e-01  3.3306e-01  6.6431e-02  3.5876e-01  1.3911e-01  1.7072e-01  8.0130e-02  7.1388e-02  2.6204e-01  7.7438e-02  5.9227e-02  6.0198e-02  3.1303e-01  3.0716e-01  8.3772e-01  5.7163e-01  3.5613e-01  2.0538e-01  9.9575e-02  3.1081e-02;
  444  6.4376e-02  2.1343e-01  3.1486e-01  5.8045e-02  4.0358e-01  3.3930e-01  6.6551e-02  3.6798e-01  1.3826e-01  1.5945e-01  8.1563e-02  7.1459e-02  2.8483e-01  7.9138e-02  5.9045e-02  6.0527e-02  3.0254e-01  3.2121e-01  8.3988e-01  5.7313e-01  3.5662e-01  2.0554e-01  9.9788e-02  3.1120e-02;
  448  6.3572e-02  2.1660e-01  3.1274e-01  5.9623e-02  3.9938e-01  3.4820e-01  6.6149e-02  3.7342e-01  1.3705e-01  1.4789e-01  8.3440e-02  7.1632e-02  3.0433e-01  8.0808e-02  5.9276e-02  6.1762e-02  2.9007e-01  3.3522e-01  8.4281e-01  5.7375e-01  3.5716e-01  2.0508e-01  9.9689e-02  3.1018e-02;
  452  6.2950e-02  2.2073e-01  3.0962e-01  6.0985e-02  3.9485e-01  3.5937e-01  6.6394e-02  3.7311e-01  1.3550e-01  1.3670e-01  8.6224e-02  7.2482e-02  3.1553e-01  8.2762e-02  5.9002e-02  6.1802e-02  2.7779e-01  3.4963e-01  8.4397e-01  5.7401e-01  3.5676e-01  2.0491e-01  9.9391e-02  3.0840e-02;
  456  6.2683e-02  2.2633e-01  3.0617e-01  6.2383e-02  3.8983e-01  3.7233e-01  6.5800e-02  3.6867e-01  1.3420e-01  1.2680e-01  8.9350e-02  7.3250e-02  3.1800e-01  8.5733e-02  5.8250e-02  6.4700e-02  2.6650e-01  3.6450e-01  8.4633e-01  5.7333e-01  3.5633e-01  2.0467e-01  9.9017e-02  3.0750e-02;
  460  6.2584e-02  2.3194e-01  3.0242e-01  6.3586e-02  3.8498e-01  3.8912e-01  6.5956e-02  3.6080e-01  1.3291e-01  1.1742e-01  9.3843e-02  7.4531e-02  3.1339e-01  8.9373e-02  5.7079e-02  6.8008e-02  2.5549e-01  3.8180e-01  8.4808e-01  5.7413e-01  3.5651e-01  2.0426e-01  9.8791e-02  3.0719e-02;
  464  6.2576e-02  2.3796e-01  2.9817e-01  6.4543e-02  3.7927e-01  4.0945e-01  6.6067e-02  3.5015e-01  1.3261e-01  1.0860e-01  9.9726e-02  7.6570e-02  3.0213e-01  9.3735e-02  5.6845e-02  7.3308e-02  2.4436e-01  3.9954e-01  8.5000e-01  5.7363e-01  3.5595e-01  2.0351e-01  9.8445e-02  3.0625e-02;
  468  6.2513e-02  2.4470e-01  2.9350e-01  6.5508e-02  3.7297e-01  4.3272e-01  6.6180e-02  3.3791e-01  1.3271e-01  1.0082e-01  1.0668e-01  7.9489e-02  2.8590e-01  9.8451e-02  5.6669e-02  8.2805e-02  2.3346e-01  4.1592e-01  8.5186e-01  5.7341e-01  3.5573e-01  2.0321e-01  9.8101e-02  3.0533e-02;
  472  6.2661e-02  2.5056e-01  2.8828e-01  6.6265e-02  3.6540e-01  4.5780e-01  6.7120e-02  3.2389e-01  1.3166e-01  9.4588e-02  1.1420e-01  8.3090e-02  2.6570e-01  1.0423e-01  5.6099e-02  9.8906e-02  2.2252e-01  4.2947e-01  8.5264e-01  5.7305e-01  3.5558e-01  2.0275e-01  9.7967e-02  3.0496e-02;
  476  6.3237e-02  2.5771e-01  2.8177e-01  6.7347e-02  3.5668e-01  4.8594e-01  6.7766e-02  3.0530e-01  1.2810e-01  8.8345e-02  1.2604e-01  8.8440e-02  2.3864e-01  1.1232e-01  5.5043e-02  1.2432e-01  2.1052e-01  4.4010e-01  8.5438e-01  5.7261e-01  3.5506e-01  2.0254e-01  9.7756e-02  3.0422e-02;
  480  6.3656e-02  2.6197e-01  2.7546e-01  6.8146e-02  3.4787e-01  5.0863e-01  6.8042e-02  2.8597e-01  1.2458e-01  8.3729e-02  1.3909e-01  9.3681e-02  2.1253e-01  1.2057e-01  5.5261e-02  1.5479e-01  1.9969e-01  4.4531e-01  8.5520e-01  5.7182e-01  3.5500e-01  2.0215e-01  9.7639e-02  3.0362e-02;
  484  6.4841e-02  2.6557e-01  2.6894e-01  6.9042e-02  3.4059e-01  5.2701e-01  6.8617e-02  2.6507e-01  1.2054e-01  7.9732e-02  1.5454e-01  9.8687e-02  1.8642e-01  1.2953e-01  5.4726e-02  1.9067e-01  1.8880e-01  4.4661e-01  8.5634e-01  5.7214e-01  3.5459e-01  2.0193e-01  9.7549e-02  3.0320e-02;
  488  6.6060e-02  2.6773e-01  2.6219e-01  6.9975e-02  3.3361e-01  5.4065e-01  6.9228e-02  2.4389e-01  1.1674e-01  7.5960e-02  1.7319e-01  1.0343e-01  1.6163e-01  1.3947e-01  5.4402e-02  2.2968e-01  1.7969e-01  4.4430e-01  8.5764e-01  5.7244e-01  3.5433e-01  2.0185e-01  9.7507e-02  3.0266e-02;
  492  6.6829e-02  2.6968e-01  2.5556e-01  7.1470e-02  3.2616e-01  5.4896e-01  7.0069e-02  2.2407e-01  1.1360e-01  7.2535e-02  1.9606e-01  1.0795e-01  1.3912e-01  1.5069e-01  5.3704e-02  2.6667e-01  1.7130e-01  4.3831e-01  8.5845e-01  5.7222e-01  3.5417e-01  2.0150e-01  9.7477e-02  3.0220e-02;
  496  6.7875e-02  2.7235e-01  2.4908e-01  7.3861e-02  3.1879e-01  5.5332e-01  7.1422e-02  2.0656e-01  1.1130e-01  6.9576e-02  2.2401e-01  1.1301e-01  1.1973e-01  1.6359e-01  5.3385e-02  2.9943e-01  1.6404e-01  4.2911e-01  8.5957e-01  5.7189e-01  3.5437e-01  2.0130e-01  9.7527e-02  3.0211e-02;
  500  6.9173e-02  2.7749e-01  2.4217e-01  7.8808e-02  3.1054e-01  5.5487e-01  7.2942e-02  1.8989e-01  1.0985e-01  6.6721e-02  2.6248e-01  1.1903e-01  1.0174e-01  1.8157e-01  5.3230e-02  3.2955e-01  1.5705e-01  4.1521e-01  8.5990e-01  5.7194e-01  3.5450e-01  2.0091e-01  9.7634e-02  3.0179e-02;
  504  7.0555e-02  2.8515e-01  2.3623e-01  8.6976e-02  3.0145e-01  5.5420e-01  7.4673e-02  1.7654e-01  1.0839e-01  6.4858e-02  3.0302e-01  1.2574e-01  8.8544e-02  2.0278e-01  5.2712e-02  3.5288e-01  1.5094e-01  4.0023e-01  8.6119e-01  5.7186e-01  3.5445e-01  2.0100e-01  9.7753e-02  3.0124e-02;
  508  7.2031e-02  2.9439e-01  2.3010e-01  9.9209e-02  2.9084e-01  5.5207e-01  7.7720e-02  1.6469e-01  1.0760e-01  6.3716e-02  3.4642e-01  1.3447e-01  7.7974e-02  2.2939e-01  5.2216e-02  3.7532e-01  1.4512e-01  3.8333e-01  8.6310e-01  5.7236e-01  3.5493e-01  2.0110e-01  9.7840e-02  3.0109e-02;
  512  7.3919e-02  3.0157e-01  2.2348e-01  1.1483e-01  2.7741e-01  5.4843e-01  8.2053e-02  1.5314e-01  1.0599e-01  6.2829e-02  3.9008e-01  1.4646e-01  6.9450e-02  2.6031e-01  5.1866e-02  3.9764e-01  1.3831e-01  3.6464e-01  8.6346e-01  5.7240e-01  3.5521e-01  2.0098e-01  9.7878e-02  3.0088e-02;
  516  7.6004e-02  3.0371e-01  2.1675e-01  1.3159e-01  2.6194e-01  5.4377e-01  8.7821e-02  1.4196e-01  1.0333e-01  6.1979e-02  4.3007e-01  1.6185e-01  6.2693e-02  2.9220e-01  5.1998e-02  4.2160e-01  1.3073e-01  3.4443e-01  8.6403e-01  5.7212e-01  3.5528e-01  2.0076e-01  9.7907e-02  3.0067e-02;
  520  7.8044e-02  2.9991e-01  2.0973e-01  1.4717e-01  2.4569e-01  5.3735e-01  9.5644e-02  1.3142e-01  1.0046e-01  6.1047e-02  4.6348e-01  1.8174e-01  5.7507e-02  3.2048e-01  5.2076e-02  4.4801e-01  1.2280e-01  3.2326e-01  8.6460e-01  5.7210e-01  3.5505e-01  2.0056e-01  9.7961e-02  3.0037e-02;
  524  8.0135e-02  2.9100e-01  2.0252e-01  1.5968e-01  2.2889e-01  5.2925e-01  1.0666e-01  1.2178e-01  9.8470e-02  6.0045e-02  4.9172e-01  2.1044e-01  5.3267e-02  3.4302e-01  5.2169e-02  4.7813e-01  1.1548e-01  2.9946e-01  8.6553e-01  5.7201e-01  3.5545e-01  2.0045e-01  9.8110e-02  3.0036e-02;
  528  8.1538e-02  2.8269e-01  1.9694e-01  1.6582e-01  2.1600e-01  5.2054e-01  1.1923e-01  1.1521e-01  9.7290e-02  5.9073e-02  5.0892e-01  2.4240e-01  5.0437e-02  3.5385e-01  5.2142e-02  5.0568e-01  1.1031e-01  2.7841e-01  8.6661e-01  5.7203e-01  3.5542e-01  2.0044e-01  9.8252e-02  3.0009e-02;
  532  8.2706e-02  2.7677e-01  1.9226e-01  1.6783e-01  2.0570e-01  5.1021e-01  1.3464e-01  1.1021e-01  9.7447e-02  5.8298e-02  5.1991e-01  2.7864e-01  4.8298e-02  3.5719e-01  5.2434e-02  5.3157e-01  1.0681e-01  2.5796e-01  8.6723e-01  5.7157e-01  3.5540e-01  2.0034e-01  9.8298e-02  2.9983e-02;
  536  8.3679e-02  2.7440e-01  1.8865e-01  1.6653e-01  1.9826e-01  4.9843e-01  1.5220e-01  1.0638e-01  9.8343e-02  5.7705e-02  5.2568e-01  3.1715e-01  4.6628e-02  3.5427e-01  5.2684e-02  5.5592e-01  1.0539e-01  2.3819e-01  8.6827e-01  5.7175e-01  3.5534e-01  2.0025e-01  9.8343e-02  2.9992e-02;
  540  8.4435e-02  2.7516e-01  1.8548e-01  1.6234e-01  1.9323e-01  4.8492e-01  1.7145e-01  1.0339e-01  1.0000e-01  5.7234e-02  5.2694e-01  3.5540e-01  4.5210e-02  3.4613e-01  5.2516e-02  5.7605e-01  1.0565e-01  2.1903e-01  8.6855e-01  5.7153e-01  3.5524e-01  1.9992e-01  9.8306e-02  2.9968e-02;
  544  8.5287e-02  2.7828e-01  1.8269e-01  1.5586e-01  1.9103e-01  4.6971e-01  1.9386e-01  1.0063e-01  1.0299e-01  5.6932e-02  5.2526e-01  3.9371e-01  4.4020e-02  3.3399e-01  5.2801e-02  5.9504e-01  1.0716e-01  1.9976e-01  8.6939e-01  5.7207e-01  3.5578e-01  1.9992e-01  9.8348e-02  2.9953e-02;
  548  8.6416e-02  2.8097e-01  1.7982e-01  1.4858e-01  1.9025e-01  4.5288e-01  2.1750e-01  9.7698e-02  1.0652e-01  5.7252e-02  5.2049e-01  4.2655e-01  4.3116e-02  3.1911e-01  5.3492e-02  6.1113e-01  1.0883e-01  1.8189e-01  8.7030e-01  5.7245e-01  3.5595e-01  2.0000e-01  9.8465e-02  2.9962e-02;
  552  8.7958e-02  2.8055e-01  1.7674e-01  1.4144e-01  1.8923e-01  4.3620e-01  2.4136e-01  9.4615e-02  1.0912e-01  5.7532e-02  5.1361e-01  4.5228e-01  4.2341e-02  3.0337e-01  5.4174e-02  6.2393e-01  1.0957e-01  1.6604e-01  8.7061e-01  5.7225e-01  3.5610e-01  2.0007e-01  9.8579e-02  2.9963e-02;
  556  9.0110e-02  2.7832e-01  1.7319e-01  1.3414e-01  1.8623e-01  4.1905e-01  2.6696e-01  9.1136e-02  1.1062e-01  5.7890e-02  5.0520e-01  4.7355e-01  4.1685e-02  2.8703e-01  5.5033e-02  6.3436e-01  1.0938e-01  1.5165e-01  8.7033e-01  5.7231e-01  3.5619e-01  2.0007e-01  9.8681e-02  2.9963e-02;
  560  9.2770e-02  2.7552e-01  1.6943e-01  1.2720e-01  1.8182e-01  4.0200e-01  2.9477e-01  8.7688e-02  1.1124e-01  5.8246e-02  4.9613e-01  4.9184e-01  4.1102e-02  2.7065e-01  5.5855e-02  6.4388e-01  1.0909e-01  1.3923e-01  8.7044e-01  5.7258e-01  3.5684e-01  2.0007e-01  9.8712e-02  2.9957e-02;
  564  9.6496e-02  2.7526e-01  1.6538e-01  1.2165e-01  1.7800e-01  3.8514e-01  3.2551e-01  8.4933e-02  1.1156e-01  5.8255e-02  4.8584e-01  5.0827e-01  4.0708e-02  2.5410e-01  5.7281e-02  6.5312e-01  1.0967e-01  1.2817e-01  8.7036e-01  5.7344e-01  3.5767e-01  2.0035e-01  9.8879e-02  2.9986e-02;
  568  1.0158e-01  2.8118e-01  1.6120e-01  1.1882e-01  1.7630e-01  3.6751e-01  3.5941e-01  8.3036e-02  1.1367e-01  5.8084e-02  4.7411e-01  5.2370e-01  4.0412e-02  2.3668e-01  6.0041e-02  6.6305e-01  1.1284e-01  1.1854e-01  8.7088e-01  5.7438e-01  3.5886e-01  2.0096e-01  9.9176e-02  3.0076e-02;
  572  1.0666e-01  2.9300e-01  1.5774e-01  1.1824e-01  1.7725e-01  3.5195e-01  3.8903e-01  8.2369e-02  1.2026e-01  5.7739e-02  4.6272e-01  5.3580e-01  4.0209e-02  2.2127e-01  6.3264e-02  6.7073e-01  1.1857e-01  1.1144e-01  8.7147e-01  5.7497e-01  3.5949e-01  2.0135e-01  9.9260e-02  3.0121e-02;
  576  1.1207e-01  3.1088e-01  1.5462e-01  1.1735e-01  1.8028e-01  3.3661e-01  4.1623e-01  8.2388e-02  1.3371e-01  5.7203e-02  4.5066e-01  5.4631e-01  4.0073e-02  2.0633e-01  6.8206e-02  6.7744e-01  1.2757e-01  1.0574e-01  8.7203e-01  5.7546e-01  3.6022e-01  2.0158e-01  9.9472e-02  3.0172e-02;
  580  1.1741e-01  3.3353e-01  1.5172e-01  1.1560e-01  1.8492e-01  3.2142e-01  4.3968e-01  8.2977e-02  1.5631e-01  5.6874e-02  4.3767e-01  5.5540e-01  4.0000e-02  1.9178e-01  7.5016e-02  6.8350e-01  1.3968e-01  1.0117e-01  8.7184e-01  5.7592e-01  3.6078e-01  2.0168e-01  9.9612e-02  3.0220e-02;
  584  1.2256e-01  3.5951e-01  1.4911e-01  1.1337e-01  1.9062e-01  3.0627e-01  4.6020e-01  8.4094e-02  1.8986e-01  5.6787e-02  4.2421e-01  5.6401e-01  3.9930e-02  1.7782e-01  8.4664e-02  6.9011e-01  1.5494e-01  9.7529e-02  8.7262e-01  5.7662e-01  3.6128e-01  2.0196e-01  9.9747e-02  3.0272e-02;
  588  1.2750e-01  3.8858e-01  1.4618e-01  1.0987e-01  1.9696e-01  2.8908e-01  4.7827e-01  8.5351e-02  2.3811e-01  5.7008e-02  4.0764e-01  5.7176e-01  3.9832e-02  1.6288e-01  1.0056e-01  6.9460e-01  1.7480e-01  9.4103e-02  8.7213e-01  5.7641e-01  3.6114e-01  2.0199e-01  9.9690e-02  3.0286e-02;
  592  1.3138e-01  4.1252e-01  1.4389e-01  1.0641e-01  2.0232e-01  2.7405e-01  4.9164e-01  8.6386e-02  2.8736e-01  5.7790e-02  3.9286e-01  5.7821e-01  3.9853e-02  1.5085e-01  1.2082e-01  6.9963e-01  1.9579e-01  9.1392e-02  8.7302e-01  5.7668e-01  3.6148e-01  2.0183e-01  9.9695e-02  3.0336e-02;
  596  1.3469e-01  4.3337e-01  1.4166e-01  1.0270e-01  2.0684e-01  2.5978e-01  5.0258e-01  8.6855e-02  3.3818e-01  5.9292e-02  3.7797e-01  5.8367e-01  3.9832e-02  1.3968e-01  1.4778e-01  7.0348e-01  2.2071e-01  8.8836e-02  8.7275e-01  5.7629e-01  3.6128e-01  2.0156e-01  9.9520e-02  3.0342e-02;
  600  1.3774e-01  4.5103e-01  1.3957e-01  9.9350e-02  2.1116e-01  2.4719e-01  5.1170e-01  8.6946e-02  3.8677e-01  6.1725e-02  3.6403e-01  5.8836e-01  3.9829e-02  1.2995e-01  1.8287e-01  7.0703e-01  2.5062e-01  8.6533e-02  8.7301e-01  5.7626e-01  3.6113e-01  2.0148e-01  9.9409e-02  3.0390e-02;
  604  1.4056e-01  4.6560e-01  1.3742e-01  9.6514e-02  2.1580e-01  2.3649e-01  5.1952e-01  8.6578e-02  4.3173e-01  6.5253e-02  3.5195e-01  5.9326e-01  3.9831e-02  1.2156e-01  2.2841e-01  7.1005e-01  2.8786e-01  8.4718e-02  8.7333e-01  5.7623e-01  3.6101e-01  2.0122e-01  9.9303e-02  3.0424e-02;
  608  1.4340e-01  4.7767e-01  1.3472e-01  9.4289e-02  2.1982e-01  2.2713e-01  5.2633e-01  8.5894e-02  4.7219e-01  7.0645e-02  3.4129e-01  5.9737e-01  3.9914e-02  1.1428e-01  2.8892e-01  7.1274e-01  3.3575e-01  8.3095e-02  8.7379e-01  5.7567e-01  3.6037e-01  2.0074e-01  9.9143e-02  3.0457e-02;
  612  1.4611e-01  4.8658e-01  1.3236e-01  9.3517e-02  2.2108e-01  2.2103e-01  5.3269e-01  8.5344e-02  5.0068e-01  7.7227e-02  3.3416e-01  6.0147e-01  4.0034e-02  1.0936e-01  3.4814e-01  7.1533e-01  3.8382e-01  8.2187e-02  8.7373e-01  5.7554e-01  3.6009e-01  2.0051e-01  9.8985e-02  3.0490e-02;
  616  1.4861e-01  4.9363e-01  1.2982e-01  9.3681e-02  2.2012e-01  2.1641e-01  5.3780e-01  8.4867e-02  5.2190e-01  8.5255e-02  3.2871e-01  6.0477e-01  4.0161e-02  1.0565e-01  4.0732e-01  7.1729e-01  4.3376e-01  8.1596e-02  8.7417e-01  5.7539e-01  3.5959e-01  2.0022e-01  9.8836e-02  3.0516e-02;
  620  1.5142e-01  4.9995e-01  1.2719e-01  9.4584e-02  2.1876e-01  2.1280e-01  5.4245e-01  8.4737e-02  5.3802e-01  9.4311e-02  3.2473e-01  6.0832e-01  4.0312e-02  1.0284e-01  4.6204e-01  7.1937e-01  4.8293e-01  8.1291e-02  8.7418e-01  5.7495e-01  3.5892e-01  1.9989e-01  9.8632e-02  3.0553e-02;
  624  1.5434e-01  5.0529e-01  1.2434e-01  9.5683e-02  2.1829e-01  2.1025e-01  5.4722e-01  8.5321e-02  5.5046e-01  1.0437e-01  3.2164e-01  6.1198e-01  4.0491e-02  1.0059e-01  5.1176e-01  7.2153e-01  5.3162e-01  8.1166e-02  8.7534e-01  5.7474e-01  3.5855e-01  1.9951e-01  9.8381e-02  3.0588e-02;
  628  1.5757e-01  5.1023e-01  1.2154e-01  9.6109e-02  2.2036e-01  2.0805e-01  5.5117e-01  8.6834e-02  5.5864e-01  1.1461e-01  3.1914e-01  6.1567e-01  4.0725e-02  9.8881e-02  5.5117e-01  7.2388e-01  5.7463e-01  8.1130e-02  8.7527e-01  5.7409e-01  3.5805e-01  1.9904e-01  9.8134e-02  3.0597e-02;
  632  1.6083e-01  5.1462e-01  1.1909e-01  9.5794e-02  2.2508e-01  2.0636e-01  5.5521e-01  8.9169e-02  5.6362e-01  1.2361e-01  3.1740e-01  6.1830e-01  4.0931e-02  9.7634e-02  5.7834e-01  7.2608e-01  6.0726e-01  8.1178e-02  8.7539e-01  5.7361e-01  3.5747e-01  1.9863e-01  9.8002e-02  3.0668e-02;
  636  1.6442e-01  5.1901e-01  1.1662e-01  9.4909e-02  2.3236e-01  2.0462e-01  5.5844e-01  9.2156e-02  5.6727e-01  1.3205e-01  3.1584e-01  6.2078e-01  4.1122e-02  9.6468e-02  5.9792e-01  7.2727e-01  6.3481e-01  8.1195e-02  8.7532e-01  5.7299e-01  3.5673e-01  1.9818e-01  9.7818e-02  3.0696e-02;
  640  1.6845e-01  5.2335e-01  1.1432e-01  9.3843e-02  2.4305e-01  2.0308e-01  5.6131e-01  9.5793e-02  5.6952e-01  1.4007e-01  3.1411e-01  6.2340e-01  4.1288e-02  9.5382e-02  6.1262e-01  7.2909e-01  6.5675e-01  8.1529e-02  8.7532e-01  5.7260e-01  3.5618e-01  1.9759e-01  9.7640e-02  3.0723e-02;
  644  1.7380e-01  5.2864e-01  1.1201e-01  9.2701e-02  2.6067e-01  2.0182e-01  5.6462e-01  1.0071e-01  5.7121e-01  1.4911e-01  3.1262e-01  6.2646e-01  4.1526e-02  9.4222e-02  6.2392e-01  7.3036e-01  6.7714e-01  8.2108e-02  8.7532e-01  5.7223e-01  3.5560e-01  1.9716e-01  9.7466e-02  3.0816e-02;
  648  1.7939e-01  5.3443e-01  1.1013e-01  9.1816e-02  2.8124e-01  2.0120e-01  5.6786e-01  1.0564e-01  5.7236e-01  1.5749e-01  3.1203e-01  6.2874e-01  4.1737e-02  9.3563e-02  6.3124e-01  7.3204e-01  6.9162e-01  8.2735e-02  8.7625e-01  5.7186e-01  3.5509e-01  1.9681e-01  9.7305e-02  3.0913e-02;
  652  1.8531e-01  5.4028e-01  1.0835e-01  9.1356e-02  3.0462e-01  2.0133e-01  5.7073e-01  1.1036e-01  5.7318e-01  1.6577e-01  3.1257e-01  6.3114e-01  4.1920e-02  9.2829e-02  6.3654e-01  7.3330e-01  7.0383e-01  8.3448e-02  8.7672e-01  5.7122e-01  3.5452e-01  1.9627e-01  9.7102e-02  3.0968e-02;
  656  1.9176e-01  5.4700e-01  1.0669e-01  9.1231e-02  3.3009e-01  2.0237e-01  5.7364e-01  1.1483e-01  5.7364e-01  1.7442e-01  3.1434e-01  6.3372e-01  4.2020e-02  9.2490e-02  6.4050e-01  7.3450e-01  7.1366e-01  8.4205e-02  8.7645e-01  5.7122e-01  3.5388e-01  1.9583e-01  9.6948e-02  3.1003e-02;
  660  1.9933e-01  5.5444e-01  1.0518e-01  9.1799e-02  3.5842e-01  2.0456e-01  5.7698e-01  1.1890e-01  5.7410e-01  1.8451e-01  3.1789e-01  6.3597e-01  4.2211e-02  9.2518e-02  6.4365e-01  7.3669e-01  7.2182e-01  8.4988e-02  8.7674e-01  5.7074e-01  3.5329e-01  1.9564e-01  9.6835e-02  3.1127e-02;
  664  2.0726e-01  5.6262e-01  1.0375e-01  9.3454e-02  3.8729e-01  2.0778e-01  5.7922e-01  1.2158e-01  5.7448e-01  1.9573e-01  3.2258e-01  6.3852e-01  4.2324e-02  9.2837e-02  6.4564e-01  7.3719e-01  7.2865e-01  8.5531e-02  8.7666e-01  5.7021e-01  3.5280e-01  1.9516e-01  9.6584e-02  3.1157e-02;
  668  2.1460e-01  5.7109e-01  1.0254e-01  9.6422e-02  4.1163e-01  2.1121e-01  5.8098e-01  1.2246e-01  5.7486e-01  2.0654e-01  3.2745e-01  6.3983e-01  4.2406e-02  9.3362e-02  6.4736e-01  7.3776e-01  7.3305e-01  8.5640e-02  8.7665e-01  5.6921e-01  3.5212e-01  1.9468e-01  9.6375e-02  3.1224e-02;
  672  2.2209e-01  5.8104e-01  1.0159e-01  1.0112e-01  4.3321e-01  2.1518e-01  5.8384e-01  1.2209e-01  5.7543e-01  2.1840e-01  3.3321e-01  6.4222e-01  4.2550e-02  9.4162e-02  6.4970e-01  7.3844e-01  7.3704e-01  8.5661e-02  8.7716e-01  5.6936e-01  3.5180e-01  1.9449e-01  9.6357e-02  3.1336e-02;
  676  2.2968e-01  5.9131e-01  1.0055e-01  1.0818e-01  4.5113e-01  2.1923e-01  5.8576e-01  1.2043e-01  5.7466e-01  2.3135e-01  3.3893e-01  6.4355e-01  4.2557e-02  9.5377e-02  6.5141e-01  7.3833e-01  7.4018e-01  8.5391e-02  8.7656e-01  5.6819e-01  3.5104e-01  1.9417e-01  9.6255e-02  3.1415e-02;
  680  2.3841e-01  6.0357e-01  9.9497e-02  1.1980e-01  4.6731e-01  2.2359e-01  5.8802e-01  1.1797e-01  5.7476e-01  2.4719e-01  3.4545e-01  6.4563e-01  4.2721e-02  9.7074e-02  6.5341e-01  7.3708e-01  7.4348e-01  8.4865e-02  8.7654e-01  5.6790e-01  3.5053e-01  1.9364e-01  9.6159e-02  3.1509e-02;
  684  2.4606e-01  6.1493e-01  9.8778e-02  1.3353e-01  4.7783e-01  2.2710e-01  5.9050e-01  1.1575e-01  5.7511e-01  2.6222e-01  3.5100e-01  6.4796e-01  4.2878e-02  9.8688e-02  6.5520e-01  7.3620e-01  7.4570e-01  8.4299e-02  8.7647e-01  5.6742e-01  3.5014e-01  1.9308e-01  9.5928e-02  3.1611e-02;
  688  2.5405e-01  6.2741e-01  9.8209e-02  1.5052e-01  4.8634e-01  2.3059e-01  5.9292e-01  1.1411e-01  5.7456e-01  2.7810e-01  3.5611e-01  6.5025e-01  4.3153e-02  1.0036e-01  6.5652e-01  7.3668e-01  7.4832e-01  8.3923e-02  8.7685e-01  5.6695e-01  3.4984e-01  1.9279e-01  9.5880e-02  3.1729e-02;
  692  2.6230e-01  6.3890e-01  9.7607e-02  1.6996e-01  4.9180e-01  2.3372e-01  5.9548e-01  1.1303e-01  5.7466e-01  2.9433e-01  3.6048e-01  6.5175e-01  4.3505e-02  1.0182e-01  6.5751e-01  7.3815e-01  7.5011e-01  8.3119e-02  8.7683e-01  5.6624e-01  3.4896e-01  1.9247e-01  9.5658e-02  3.1772e-02;
  696  2.7222e-01  6.5336e-01  9.7104e-02  1.9438e-01  4.9671e-01  2.3681e-01  5.9807e-01  1.1286e-01  5.7481e-01  3.1373e-01  3.6477e-01  6.5467e-01  4.4054e-02  1.0373e-01  6.5862e-01  7.4111e-01  7.5296e-01  8.2361e-02  8.7670e-01  5.6560e-01  3.4835e-01  1.9206e-01  9.5437e-02  3.1856e-02;
  700  2.8121e-01  6.6551e-01  9.6746e-02  2.1505e-01  4.9935e-01  2.3861e-01  6.0043e-01  1.1354e-01  5.7440e-01  3.3046e-01  3.6759e-01  6.5640e-01  4.4685e-02  1.0503e-01  6.5944e-01  7.4403e-01  7.5401e-01  8.1779e-02  8.7679e-01  5.6573e-01  3.4781e-01  1.9171e-01  9.5358e-02  3.1991e-02;
  704  2.9028e-01  6.7703e-01  9.6373e-02  2.3385e-01  5.0130e-01  2.3951e-01  6.0233e-01  1.1477e-01  5.7427e-01  3.4702e-01  3.6913e-01  6.5803e-01  4.5423e-02  1.0613e-01  6.5976e-01  7.4611e-01  7.5475e-01  8.1477e-02  8.7651e-01  5.6520e-01  3.4741e-01  1.9136e-01  9.5207e-02  3.2081e-02;
  708  2.9987e-01  6.8798e-01  9.6137e-02  2.4991e-01  5.0258e-01  2.3918e-01  6.0429e-01  1.1652e-01  5.7468e-01  3.6412e-01  3.6948e-01  6.5966e-01  4.6223e-02  1.0687e-01  6.6009e-01  7.4807e-01  7.5579e-01  8.1545e-02  8.7639e-01  5.6438e-01  3.4700e-01  1.9099e-01  9.5064e-02  3.2197e-02;
  712  3.1220e-01  7.0021e-01  9.5778e-02  2.6397e-01  5.0362e-01  2.3770e-01  6.0640e-01  1.1906e-01  5.7441e-01  3.8443e-01  3.6823e-01  6.6141e-01  4.7207e-02  1.0691e-01  6.6013e-01  7.4925e-01  7.5778e-01  8.2260e-02  8.7633e-01  5.6375e-01  3.4623e-01  1.9049e-01  9.4925e-02  3.2307e-02;
  716  3.2450e-01  7.1011e-01  9.5599e-02  2.7283e-01  5.0487e-01  2.3555e-01  6.0855e-01  1.2167e-01  5.7427e-01  4.0271e-01  3.6644e-01  6.6314e-01  4.8201e-02  1.0656e-01  6.6060e-01  7.5074e-01  7.5963e-01  8.3538e-02  8.7601e-01  5.6327e-01  3.4566e-01  1.9006e-01  9.4710e-02  3.2366e-02;
  720  3.3865e-01  7.1916e-01  9.5705e-02  2.7941e-01  5.0526e-01  2.3406e-01  6.1053e-01  1.2467e-01  5.7432e-01  4.2189e-01  3.6543e-01  6.6526e-01  4.9347e-02  1.0648e-01  6.6105e-01  7.5200e-01  7.6253e-01  8.5053e-02  8.7663e-01  5.6337e-01  3.4531e-01  1.8973e-01  9.4611e-02  3.2476e-02;
  724  3.5390e-01  7.2651e-01  9.5699e-02  2.8388e-01  5.0564e-01  2.3319e-01  6.1211e-01  1.2806e-01  5.7370e-01  4.4092e-01  3.6518e-01  6.6639e-01  5.0605e-02  1.0664e-01  6.6138e-01  7.5240e-01  7.6409e-01  8.7182e-02  8.7641e-01  5.6284e-01  3.4434e-01  1.8956e-01  9.4405e-02  3.2526e-02;
  728  3.7210e-01  7.3433e-01  9.5600e-02  2.8755e-01  5.0643e-01  2.3441e-01  6.1436e-01  1.3304e-01  5.7410e-01  4.6368e-01  3.6700e-01  6.6833e-01  5.2470e-02  1.0739e-01  6.6252e-01  7.5384e-01  7.6754e-01  9.0328e-02  8.7588e-01  5.6206e-01  3.4355e-01  1.8921e-01  9.4230e-02  3.2661e-02;
  732  3.8801e-01  7.3960e-01  9.5632e-02  2.8970e-01  5.0639e-01  2.3704e-01  6.1599e-01  1.3861e-01  5.7355e-01  4.8331e-01  3.7075e-01  6.6996e-01  5.4553e-02  1.0902e-01  6.6296e-01  7.5443e-01  7.6967e-01  9.3737e-02  8.7680e-01  5.6119e-01  3.4330e-01  1.8871e-01  9.4232e-02  3.2769e-02;
  736  4.0416e-01  7.4531e-01  9.5714e-02  2.9143e-01  5.0735e-01  2.4151e-01  6.1878e-01  1.4608e-01  5.7347e-01  5.0367e-01  3.7714e-01  6.7143e-01  5.7429e-02  1.1159e-01  6.6449e-01  7.5551e-01  7.7306e-01  9.8408e-02  8.7673e-01  5.6041e-01  3.4265e-01  1.8869e-01  9.4122e-02  3.2845e-02;
  740  4.2033e-01  7.5000e-01  9.6057e-02  2.9305e-01  5.0772e-01  2.4789e-01  6.2073e-01  1.5630e-01  5.7358e-01  5.2439e-01  3.8549e-01  6.7317e-01  6.1301e-02  1.1480e-01  6.6504e-01  7.5610e-01  7.7642e-01  1.0561e-01  8.7520e-01  5.5976e-01  3.4252e-01  1.8797e-01  9.4065e-02  3.2927e-02;
  744  4.3742e-01  7.5334e-01  9.6193e-02  2.9441e-01  5.0749e-01  2.5642e-01  6.2252e-01  1.7027e-01  5.7311e-01  5.4638e-01  3.9603e-01  6.7436e-01  6.7031e-02  1.1871e-01  6.6626e-01  7.5577e-01  7.7926e-01  1.1596e-01  8.7485e-01  5.5893e-01  3.4176e-01  1.8696e-01  9.3884e-02  3.3030e-02;
  748  4.5107e-01  7.5554e-01  9.6416e-02  2.9525e-01  5.0745e-01  2.6432e-01  6.2344e-01  1.8498e-01  5.7269e-01  5.6424e-01  4.0596e-01  6.7539e-01  7.3822e-02  1.2267e-01  6.6734e-01  7.5554e-01  7.8131e-01  1.2815e-01  8.7515e-01  5.5820e-01  3.4068e-01  1.8675e-01  9.3919e-02  3.3117e-02;
  752  4.6391e-01  7.5822e-01  9.6792e-02  2.9615e-01  5.0762e-01  2.7193e-01  6.2510e-01  2.0160e-01  5.7217e-01  5.8019e-01  4.1500e-01  6.7602e-01  8.1877e-02  1.2642e-01  6.6840e-01  7.5581e-01  7.8308e-01  1.4246e-01  8.7410e-01  5.5774e-01  3.3974e-01  1.8669e-01  9.3905e-02  3.3216e-02;
  756  4.7486e-01  7.6097e-01  9.7167e-02  2.9705e-01  5.0838e-01  2.7953e-01  6.2610e-01  2.2079e-01  5.7263e-01  5.9417e-01  4.2378e-01  6.7678e-01  9.1780e-02  1.3017e-01  6.7079e-01  7.5579e-01  7.8492e-01  1.6018e-01  8.7470e-01  5.5666e-01  3.3966e-01  1.8579e-01  9.3615e-02  3.3228e-02;
  760  4.8310e-01  7.6143e-01  9.7734e-02  2.9805e-01  5.0775e-01  2.8565e-01  6.2704e-01  2.4087e-01  5.7296e-01  6.0437e-01  4.3221e-01  6.7793e-01  1.0282e-01  1.3348e-01  6.7276e-01  7.5586e-01  7.8529e-01  1.7797e-01  8.7356e-01  5.5547e-01  3.3909e-01  1.8553e-01  9.3519e-02  3.3237e-02;
  764  4.8773e-01  7.6247e-01  9.8060e-02  2.9893e-01  5.0792e-01  2.8994e-01  6.2827e-01  2.6021e-01  5.7205e-01  6.1203e-01  4.3705e-01  6.7933e-01  1.1401e-01  1.3646e-01  6.7458e-01  7.5534e-01  7.8543e-01  1.9188e-01  8.7173e-01  5.5542e-01  3.3789e-01  1.8547e-01  9.3547e-02  3.3365e-02;
  768  4.9093e-01  7.6222e-01  9.8541e-02  2.9976e-01  5.0749e-01  2.9243e-01  6.2815e-01  2.8028e-01  5.7177e-01  6.1790e-01  4.3967e-01  6.7981e-01  1.2603e-01  1.3876e-01  6.7547e-01  7.5473e-01  7.8667e-01  2.0343e-01  8.7027e-01  5.5521e-01  3.3699e-01  1.8529e-01  9.3375e-02  3.3462e-02;
  772  4.9159e-01  7.6105e-01  9.9257e-02  3.0043e-01  5.0684e-01  2.9441e-01  6.2886e-01  3.0493e-01  5.7059e-01  6.2221e-01  4.4192e-01  6.7853e-01  1.4075e-01  1.4005e-01  6.7697e-01  7.5401e-01  7.8569e-01  2.1326e-01  8.6977e-01  5.5377e-01  3.3551e-01  1.8447e-01  9.3195e-02  3.3465e-02;
  776  4.8987e-01  7.6072e-01  9.9571e-02  3.0090e-01  5.0585e-01  2.9521e-01  6.2822e-01  3.2767e-01  5.6976e-01  6.2471e-01  4.4232e-01  6.7771e-01  1.5460e-01  1.4069e-01  6.7966e-01  7.5370e-01  7.8293e-01  2.1929e-01  8.6906e-01  5.5144e-01  3.3449e-01  1.8383e-01  9.2907e-02  3.3484e-02;
  780  4.8817e-01  7.5921e-01  1.0047e-01  3.0163e-01  5.0485e-01  2.9616e-01  6.2737e-01  3.5211e-01  5.6883e-01  6.2699e-01  4.4281e-01  6.7701e-01  1.6937e-01  1.4114e-01  6.8050e-01  7.5145e-01  7.8131e-01  2.2404e-01  8.6661e-01  5.4983e-01  3.3466e-01  1.8379e-01  9.3021e-02  3.3501e-02]

#=
    Data from from Noboru Ohta, Alan R. Robertson, “Colorimetry: Fundamentals and Applications”, (May 2006), DOI:10.1002/0470094745
=#
const Macbeth_chart3 =
   [380.0  0.048  0.103  0.113  0.048  0.123  0.110  0.053  0.099  0.096  0.101  0.056  0.060  0.069  0.055  0.052  0.054  0.118  0.093  0.153  0.150  0.138  0.113  0.074  0.032;
    385.0  0.051  0.120  0.138  0.049  0.152  0.133  0.054  0.120  0.108  0.115  0.058  0.061  0.081  0.056  0.052  0.053  0.142  0.110  0.189  0.184  0.167  0.131  0.079  0.033;
    390.0  0.055  0.141  0.174  0.049  0.197  0.167  0.054  0.150  0.123  0.135  0.059  0.063  0.096  0.057  0.052  0.054  0.179  0.134  0.245  0.235  0.206  0.150  0.084  0.033;
    395.0  0.060  0.163  0.219  0.049  0.258  0.208  0.054  0.189  0.135  0.157  0.059  0.064  0.114  0.058  0.052  0.053  0.228  0.164  0.319  0.299  0.249  0.169  0.088  0.034;
    400.0  0.065  0.182  0.266  0.050  0.328  0.252  0.054  0.231  0.144  0.177  0.060  0.065  0.136  0.058  0.051  0.053  0.283  0.195  0.409  0.372  0.289  0.183  0.091  0.035;
    405.0  0.068  0.192  0.300  0.049  0.385  0.284  0.054  0.268  0.145  0.191  0.061  0.065  0.156  0.058  0.051  0.053  0.322  0.220  0.536  0.459  0.324  0.193  0.093  0.035;
    410.0  0.068  0.197  0.320  0.049  0.418  0.303  0.053  0.293  0.144  0.199  0.061  0.064  0.175  0.059  0.050  0.053  0.343  0.238  0.671  0.529  0.346  0.199  0.094  0.036;
    415.0  0.067  0.199  0.330  0.050  0.437  0.314  0.053  0.311  0.141  0.203  0.061  0.064  0.193  0.059  0.050  0.052  0.354  0.249  0.772  0.564  0.354  0.201  0.094  0.036;
    420.0  0.064  0.201  0.336  0.050  0.446  0.322  0.052  0.324  0.138  0.206  0.062  0.064  0.208  0.059  0.049  0.052  0.359  0.258  0.840  0.580  0.357  0.202  0.094  0.036;
    425.0  0.062  0.203  0.337  0.051  0.448  0.329  0.052  0.335  0.134  0.198  0.063  0.064  0.224  0.060  0.049  0.052  0.357  0.270  0.868  0.584  0.358  0.203  0.094  0.036;
    430.0  0.059  0.205  0.337  0.052  0.448  0.336  0.052  0.348  0.132  0.190  0.064  0.064  0.244  0.062  0.049  0.053  0.350  0.281  0.878  0.585  0.359  0.203  0.094  0.036;
    435.0  0.057  0.208  0.337  0.053  0.447  0.344  0.052  0.361  0.132  0.179  0.066  0.065  0.265  0.063  0.049  0.053  0.339  0.296  0.882  0.587  0.360  0.204  0.095  0.036;
    440.0  0.055  0.212  0.335  0.054  0.444  0.353  0.052  0.373  0.131  0.168  0.068  0.065  0.290  0.065  0.049  0.053  0.327  0.315  0.883  0.587  0.361  0.205  0.095  0.035;
    445.0  0.054  0.217  0.334  0.056  0.440  0.363  0.052  0.383  0.131  0.156  0.071  0.066  0.316  0.067  0.049  0.054  0.313  0.334  0.885  0.588  0.362  0.205  0.095  0.035;
    450.0  0.053  0.224  0.331  0.058  0.434  0.375  0.052  0.387  0.129  0.144  0.075  0.067  0.335  0.070  0.049  0.055  0.298  0.352  0.886  0.588  0.362  0.205  0.095  0.035;
    455.0  0.053  0.231  0.327  0.060  0.428  0.390  0.052  0.383  0.128  0.132  0.079  0.068  0.342  0.074  0.048  0.056  0.282  0.370  0.886  0.587  0.361  0.205  0.094  0.035;
    460.0  0.052  0.240  0.322  0.061  0.421  0.408  0.052  0.374  0.126  0.120  0.085  0.069  0.338  0.078  0.048  0.059  0.267  0.391  0.887  0.586  0.361  0.204  0.094  0.035;
    465.0  0.052  0.251  0.316  0.063  0.413  0.433  0.052  0.361  0.126  0.110  0.093  0.073  0.324  0.084  0.047  0.065  0.253  0.414  0.888  0.585  0.359  0.204  0.094  0.035;
    470.0  0.052  0.262  0.310  0.064  0.405  0.460  0.053  0.345  0.125  0.101  0.104  0.077  0.302  0.091  0.047  0.075  0.239  0.434  0.888  0.583  0.358  0.203  0.094  0.035;
    475.0  0.053  0.273  0.302  0.065  0.394  0.492  0.054  0.325  0.123  0.093  0.118  0.084  0.273  0.101  0.046  0.093  0.225  0.449  0.888  0.582  0.358  0.203  0.093  0.035;
    480.0  0.054  0.282  0.293  0.067  0.381  0.523  0.055  0.301  0.119  0.086  0.135  0.092  0.239  0.113  0.045  0.121  0.209  0.458  0.888  0.581  0.357  0.202  0.093  0.034;
    485.0  0.055  0.289  0.285  0.068  0.372  0.548  0.056  0.275  0.114  0.080  0.157  0.100  0.205  0.125  0.045  0.157  0.195  0.461  0.888  0.580  0.356  0.202  0.093  0.034;
    490.0  0.057  0.293  0.276  0.070  0.362  0.566  0.057  0.247  0.109  0.075  0.185  0.107  0.172  0.140  0.044  0.202  0.182  0.457  0.888  0.580  0.356  0.202  0.093  0.034;
    495.0  0.059  0.296  0.268  0.072  0.352  0.577  0.059  0.223  0.105  0.070  0.221  0.115  0.144  0.157  0.044  0.252  0.172  0.447  0.888  0.580  0.356  0.202  0.092  0.034;
    500.0  0.061  0.301  0.260  0.078  0.342  0.582  0.061  0.202  0.103  0.067  0.269  0.123  0.120  0.180  0.044  0.303  0.163  0.433  0.887  0.580  0.356  0.202  0.092  0.034;
    505.0  0.062  0.310  0.251  0.088  0.330  0.583  0.064  0.184  0.102  0.063  0.326  0.133  0.101  0.208  0.044  0.351  0.155  0.414  0.887  0.580  0.356  0.202  0.093  0.034;
    510.0  0.065  0.321  0.243  0.106  0.314  0.580  0.068  0.167  0.100  0.061  0.384  0.146  0.086  0.244  0.044  0.394  0.146  0.392  0.887  0.580  0.356  0.202  0.093  0.034;
    515.0  0.067  0.326  0.234  0.130  0.294  0.576  0.076  0.152  0.097  0.059  0.440  0.166  0.074  0.286  0.044  0.436  0.135  0.366  0.887  0.581  0.356  0.202  0.093  0.034;
    520.0  0.070  0.322  0.225  0.155  0.271  0.569  0.086  0.137  0.094  0.058  0.484  0.193  0.066  0.324  0.044  0.475  0.124  0.339  0.887  0.581  0.357  0.202  0.093  0.034;
    525.0  0.072  0.310  0.215  0.173  0.249  0.560  0.101  0.125  0.091  0.056  0.516  0.229  0.059  0.351  0.044  0.512  0.113  0.310  0.887  0.582  0.357  0.202  0.093  0.034;
    530.0  0.074  0.298  0.208  0.181  0.231  0.549  0.120  0.116  0.089  0.054  0.534  0.273  0.054  0.363  0.044  0.544  0.106  0.282  0.887  0.582  0.357  0.203  0.093  0.034;
    535.0  0.075  0.291  0.203  0.182  0.219  0.535  0.143  0.110  0.090  0.053  0.542  0.323  0.051  0.363  0.044  0.572  0.102  0.255  0.887  0.582  0.358  0.203  0.093  0.034;
    540.0  0.076  0.292  0.198  0.177  0.211  0.519  0.170  0.106  0.092  0.052  0.545  0.374  0.048  0.355  0.045  0.597  0.102  0.228  0.887  0.583  0.358  0.203  0.093  0.034;
    545.0  0.078  0.297  0.195  0.168  0.209  0.501  0.198  0.103  0.096  0.052  0.541  0.418  0.046  0.342  0.046  0.615  0.105  0.204  0.886  0.583  0.358  0.203  0.093  0.034;
    550.0  0.079  0.300  0.191  0.157  0.209  0.480  0.228  0.099  0.102  0.053  0.533  0.456  0.045  0.323  0.047  0.630  0.107  0.180  0.886  0.583  0.358  0.203  0.093  0.034;
    555.0  0.082  0.298  0.188  0.147  0.207  0.458  0.260  0.094  0.106  0.054  0.524  0.487  0.044  0.303  0.048  0.645  0.107  0.159  0.887  0.584  0.358  0.203  0.092  0.034;
    560.0  0.087  0.295  0.183  0.137  0.201  0.436  0.297  0.090  0.108  0.055  0.513  0.512  0.043  0.281  0.050  0.660  0.106  0.141  0.887  0.584  0.359  0.203  0.093  0.033;
    565.0  0.092  0.295  0.177  0.129  0.196  0.414  0.338  0.086  0.109  0.055  0.501  0.534  0.042  0.260  0.053  0.673  0.107  0.126  0.887  0.585  0.359  0.203  0.093  0.033;
    570.0  0.100  0.305  0.172  0.126  0.196  0.392  0.380  0.083  0.112  0.054  0.487  0.554  0.041  0.238  0.057  0.686  0.112  0.114  0.888  0.586  0.360  0.204  0.093  0.033;
    575.0  0.107  0.326  0.167  0.125  0.199  0.369  0.418  0.083  0.126  0.053  0.472  0.570  0.041  0.217  0.063  0.698  0.123  0.104  0.888  0.587  0.361  0.204  0.093  0.033;
    580.0  0.115  0.358  0.163  0.122  0.206  0.346  0.452  0.083  0.157  0.052  0.454  0.584  0.040  0.196  0.072  0.708  0.141  0.097  0.887  0.588  0.361  0.205  0.093  0.033;
    585.0  0.122  0.397  0.160  0.119  0.215  0.324  0.481  0.085  0.208  0.052  0.436  0.598  0.040  0.177  0.086  0.718  0.166  0.092  0.886  0.588  0.361  0.205  0.093  0.033;
    590.0  0.129  0.435  0.157  0.115  0.223  0.302  0.503  0.086  0.274  0.053  0.416  0.609  0.040  0.158  0.109  0.726  0.198  0.088  0.886  0.588  0.361  0.205  0.093  0.033;
    595.0  0.134  0.468  0.153  0.109  0.229  0.279  0.520  0.087  0.346  0.055  0.394  0.617  0.040  0.140  0.143  0.732  0.235  0.083  0.886  0.588  0.361  0.205  0.092  0.033;
    600.0  0.138  0.494  0.150  0.104  0.235  0.260  0.532  0.087  0.415  0.059  0.374  0.624  0.039  0.124  0.192  0.737  0.279  0.080  0.887  0.588  0.360  0.204  0.092  0.033;
    605.0  0.142  0.514  0.147  0.100  0.241  0.245  0.543  0.086  0.473  0.065  0.358  0.630  0.039  0.111  0.256  0.742  0.333  0.077  0.888  0.587  0.360  0.204  0.092  0.033;
    610.0  0.146  0.530  0.144  0.098  0.245  0.234  0.552  0.085  0.517  0.074  0.346  0.635  0.040  0.101  0.332  0.746  0.394  0.075  0.889  0.586  0.359  0.204  0.092  0.033;
    615.0  0.150  0.541  0.141  0.097  0.245  0.226  0.560  0.084  0.547  0.086  0.337  0.640  0.040  0.094  0.413  0.749  0.460  0.074  0.890  0.586  0.358  0.203  0.091  0.033;
    620.0  0.154  0.550  0.137  0.098  0.243  0.221  0.566  0.084  0.567  0.099  0.331  0.645  0.040  0.089  0.486  0.753  0.522  0.073  0.891  0.585  0.357  0.203  0.091  0.033;
    625.0  0.158  0.557  0.133  0.100  0.243  0.217  0.572  0.085  0.582  0.113  0.328  0.650  0.040  0.086  0.550  0.757  0.580  0.073  0.891  0.584  0.356  0.202  0.091  0.033;
    630.0  0.163  0.564  0.130  0.100  0.247  0.215  0.578  0.088  0.591  0.126  0.325  0.654  0.041  0.084  0.598  0.761  0.628  0.073  0.891  0.583  0.355  0.201  0.090  0.033;
    635.0  0.167  0.569  0.126  0.099  0.254  0.212  0.583  0.092  0.597  0.138  0.322  0.658  0.041  0.082  0.631  0.765  0.666  0.073  0.891  0.581  0.354  0.201  0.090  0.033;
    640.0  0.173  0.574  0.123  0.097  0.269  0.210  0.587  0.098  0.601  0.149  0.320  0.662  0.042  0.080  0.654  0.768  0.696  0.073  0.890  0.580  0.353  0.200  0.090  0.033;
    645.0  0.180  0.582  0.120  0.096  0.291  0.209  0.593  0.105  0.604  0.161  0.319  0.667  0.042  0.078  0.672  0.772  0.722  0.073  0.889  0.579  0.352  0.199  0.090  0.033;
    650.0  0.188  0.590  0.118  0.095  0.318  0.208  0.599  0.111  0.607  0.172  0.319  0.672  0.042  0.077  0.686  0.777  0.742  0.074  0.889  0.578  0.351  0.198  0.089  0.033;
    655.0  0.196  0.597  0.115  0.095  0.351  0.209  0.602  0.118  0.608  0.182  0.320  0.675  0.043  0.076  0.694  0.779  0.756  0.075  0.889  0.577  0.350  0.198  0.089  0.033;
    660.0  0.204  0.605  0.112  0.095  0.384  0.211  0.604  0.123  0.607  0.193  0.324  0.676  0.043  0.075  0.700  0.780  0.766  0.076  0.889  0.576  0.349  0.197  0.089  0.033;
    665.0  0.213  0.614  0.110  0.097  0.417  0.215  0.606  0.126  0.606  0.205  0.330  0.677  0.043  0.075  0.704  0.780  0.774  0.076  0.889  0.575  0.348  0.197  0.088  0.033;
    670.0  0.222  0.624  0.108  0.101  0.446  0.220  0.608  0.126  0.605  0.217  0.337  0.678  0.044  0.075  0.707  0.781  0.780  0.077  0.888  0.574  0.346  0.196  0.088  0.033;
    675.0  0.231  0.637  0.106  0.110  0.470  0.227  0.611  0.124  0.605  0.232  0.345  0.681  0.044  0.077  0.712  0.782  0.785  0.076  0.888  0.573  0.346  0.195  0.088  0.033;
    680.0  0.242  0.652  0.105  0.125  0.490  0.233  0.615  0.120  0.605  0.248  0.354  0.685  0.044  0.078  0.718  0.785  0.791  0.075  0.888  0.572  0.345  0.195  0.087  0.033;
    685.0  0.251  0.668  0.104  0.147  0.504  0.239  0.619  0.117  0.604  0.266  0.362  0.688  0.044  0.080  0.721  0.785  0.794  0.074  0.888  0.571  0.344  0.194  0.087  0.033;
    690.0  0.261  0.682  0.104  0.174  0.511  0.244  0.622  0.115  0.605  0.282  0.368  0.690  0.045  0.082  0.724  0.787  0.798  0.074  0.888  0.570  0.343  0.194  0.087  0.032;
    695.0  0.271  0.697  0.103  0.210  0.517  0.249  0.625  0.115  0.606  0.301  0.375  0.693  0.046  0.085  0.727  0.789  0.801  0.073  0.888  0.569  0.342  0.193  0.087  0.032;
    700.0  0.282  0.713  0.103  0.247  0.520  0.252  0.628  0.116  0.606  0.319  0.379  0.696  0.048  0.088  0.729  0.792  0.804  0.072  0.888  0.568  0.341  0.192  0.086  0.032;
    705.0  0.294  0.728  0.102  0.283  0.522  0.252  0.630  0.118  0.604  0.338  0.381  0.698  0.050  0.089  0.730  0.792  0.806  0.072  0.887  0.567  0.340  0.192  0.086  0.032;
    710.0  0.305  0.745  0.102  0.311  0.523  0.250  0.633  0.120  0.602  0.355  0.379  0.698  0.051  0.089  0.730  0.793  0.807  0.071  0.886  0.566  0.339  0.191  0.086  0.032;
    715.0  0.318  0.753  0.102  0.329  0.522  0.248  0.633  0.124  0.601  0.371  0.376  0.698  0.053  0.090  0.729  0.792  0.807  0.073  0.886  0.565  0.338  0.191  0.086  0.032;
    720.0  0.334  0.762  0.102  0.343  0.521  0.244  0.633  0.128  0.599  0.388  0.373  0.698  0.056  0.090  0.727  0.790  0.807  0.075  0.886  0.564  0.337  0.190  0.085  0.032;
    725.0  0.354  0.774  0.102  0.353  0.521  0.245  0.636  0.133  0.598  0.406  0.372  0.700  0.060  0.090  0.728  0.792  0.810  0.078  0.885  0.562  0.336  0.189  0.085  0.032;
    730.0  0.372  0.783  0.102  0.358  0.522  0.245  0.637  0.139  0.596  0.422  0.375  0.701  0.064  0.089  0.729  0.792  0.813  0.082  0.885  0.562  0.335  0.189  0.085  0.032;
    735.0  0.392  0.788  0.104  0.362  0.521  0.251  0.639  0.149  0.595  0.436  0.382  0.701  0.070  0.092  0.729  0.790  0.814  0.090  0.885  0.560  0.334  0.188  0.085  0.032;
    740.0  0.409  0.791  0.104  0.364  0.521  0.260  0.638  0.162  0.593  0.451  0.392  0.701  0.079  0.094  0.727  0.787  0.813  0.100  0.884  0.560  0.333  0.188  0.085  0.032;
    745.0  0.420  0.787  0.104  0.360  0.516  0.269  0.633  0.178  0.587  0.460  0.401  0.695  0.091  0.097  0.723  0.782  0.810  0.116  0.884  0.558  0.332  0.187  0.084  0.032;
    750.0  0.436  0.789  0.104  0.362  0.514  0.278  0.633  0.197  0.584  0.471  0.412  0.694  0.104  0.102  0.721  0.778  0.808  0.133  0.883  0.557  0.331  0.187  0.084  0.032;
    755.0  0.450  0.794  0.106  0.364  0.514  0.288  0.636  0.219  0.584  0.481  0.422  0.696  0.120  0.106  0.724  0.780  0.811  0.154  0.882  0.556  0.330  0.186  0.084  0.032;
    760.0  0.462  0.801  0.106  0.368  0.517  0.297  0.641  0.242  0.586  0.492  0.433  0.700  0.138  0.110  0.728  0.782  0.814  0.176  0.882  0.555  0.329  0.185  0.084  0.032;
    765.0  0.465  0.799  0.107  0.368  0.515  0.301  0.639  0.259  0.584  0.495  0.436  0.698  0.154  0.111  0.727  0.781  0.813  0.191  0.881  0.554  0.328  0.185  0.084  0.032;
    770.0  0.448  0.771  0.110  0.355  0.500  0.297  0.616  0.275  0.566  0.482  0.426  0.673  0.168  0.112  0.702  0.752  0.785  0.200  0.880  0.553  0.327  0.184  0.083  0.032;
    775.0  0.432  0.747  0.115  0.346  0.491  0.296  0.598  0.294  0.551  0.471  0.413  0.653  0.186  0.112  0.680  0.728  0.765  0.208  0.880  0.551  0.326  0.184  0.083  0.032;
    780.0  0.421  0.734  0.120  0.341  0.487  0.296  0.582  0.316  0.540  0.467  0.404  0.639  0.204  0.112  0.664  0.710  0.752  0.214  0.879  0.550  0.325  0.183  0.083  0.032]


text:
# IndexCounter

eval(Expr(
    :struct, true, :IndexCounter, Expr(:block,
        [:($(t.fieldname)::Int) for t in plexostables if isnothing(t.identifier)]...
    )
))

IndexCounter() = IndexCounter(zeros(Int, length(fieldnames(IndexCounter)))...)

function increment!(x::IndexCounter, fieldname::Symbol)
    idx = getfield(x, fieldname) + 1
    setfield!(x, fieldname, idx)
    return idx
end

# PLEXOSSolutionDataset

eval(Expr(
    :struct, false, :(PLEXOSSolutionDataset <: AbstractDataset), Expr(:block,
        [:($(t.fieldname)::Vector{$(t.fieldtype)}) for t in plexostables]...
    )
))

function PLEXOSSolutionDataset(
    zippath::String, xmlname::String=defaultxml(zippath)
)

    resultsarchive = _open_plexoszip(zippath)
    xml = parsexml(resultsarchive[xmlname])
    return PLEXOSSolutionDataset(xml)
end

function PLEXOSSolutionDataset(xml::Document)

    summary = PLEXOSSolutionDatasetSummary(xml)
    result = PLEXOSSolutionDataset(summary, consolidated=false)
    idxcounter = IndexCounter()

    for loadorder in 1:7
        for element in eachelement(xml.root)

            # Ignore the band table
            element.name == "t_band" && continue

            table = plexostables_lookup[element.name]
            table.loadorder == loadorder || continue

            idx = if isnothing(table.identifier)
                      increment!(idxcounter, table.fieldname)
                  else
                      getchildint(table.identifier, element) + table.indexoffset
                  end

            getfield(result, table.fieldname)[idx] =
                eval(table.fieldtype)(element, result)

        end
    end

    return consolidate(result, summary)

end

function PLEXOSSolutionDataset(
    summary::PLEXOSSolutionDatasetSummary;
    consolidated::Bool=false)

    selector = consolidated ? first : last

    return PLEXOSSolutionDataset((
        Vector{eval(t.fieldtype)}(undef, selector(getfield(summary, t.fieldname)))
        for t in plexostables)...)

end

function consolidate(
    unconsolidated::PLEXOSSolutionDataset,
    summary::PLEXOSSolutionDatasetSummary)

    result = PLEXOSSolutionDataset(summary, consolidated=true)

    for name in fieldnames(PLEXOSSolutionDataset)
        idx = 0
        vec = getfield(unconsolidated, name)
        for i in 1:length(vec)
            if isassigned(vec, i)
                idx += 1
                getfield(result, name)[idx] = vec[i]
            end
        end
    end

    return result

end

ext: jl


text:
"""
    classification_report(i, yhat, y)

    Returns multiple classification metrics based on predictions yhat and
    ground truth values y, for run i.
 """
function classification_report(i, yhat, y)
    P′ = [pdf(yhat[j], 1.0) for j ∈ 1:size(yhat)[1]]
    yₚ = [p >= 0.5 for p ∈ P′]
    out = DataFrame(
        Run = i,
        TN  = true_negative(yₚ, y),
        FP  = false_positive(yₚ, y),
        FN  = false_negative(yₚ, y),
        TP  = true_positive(yₚ, y),
        Sensitivity = sensitivity(yₚ, y),
        Specificity = specificity(yₚ, y),
        PPV = ppv(yₚ, y),
        NPV = npv(yₚ, y),
        Accuracy = accuracy(yₚ, y),
        AUC = auc(yhat, y),
        MCC = mcc(yₚ, y),
        F1  = f1score(yₚ, y),
        Brier = mean((yₚ .- float(y)).^2)
    )
    return out
end

"""
    classification_predictions(fold, train_sites, test_sites, ytrue, ypred)

    Returns a table with individual level predictions
"""
function classification_predictions(fold, train_site, test_sites, ytrue, ypred)
    return DataFrame(
        fold = fold,
        train_site = train_site,
        test_site = vec(Matrix(test_sites)),
        ytrue = int(ytrue, type=Int64) .- 1,
        ypred = map(y->pdf(y, 1.0), ypred))
end

"""
    compute_reliability_curve(predictions, n_bins)

"""
function compute_reliability_curve(
    predictions::DataFrame,
    n_bins::Int64=10)

    sort!(predictions, :ypred)
    N = size(predictions)[1]
    increment = convert(Int, ceil(N/n_bins))
    indices = collect(1:increment:N)
    push!(indices, N)
    fₚ = []; fₒ = []
    for i ∈ 2:length(indices)
        pred_i = predictions[indices[i-1]:indices[i],:]
        push!(fₚ, pred_i[:,:ypred] |> mean)
        push!(fₒ, pred_i[:,:ytrue] |> mean)
    end

    return fₚ, fₒ
end

"""
    plot_reliability_curves(predictions; color, alpha, title, xlabel, ylabel, figsize)
"""
function plot_reliability_curves(
    predictions::DataFrame;
    color::Symbol=:black,
    alpha::Float64=1.0,
    title::Union{String, LaTeXString}="Reliability",
    xlabel::Union{String, LaTeXString}=L"f_p",
    ylabel::Union{String, LaTeXString}=L"f_o",
    figsize::Tuple=(250,250))
    fig = plot([0,1], [0,1], c=:black,
        ls=:dash, label=nothing, size=figsize,
        title=title, xlabel=xlabel, ylabel=ylabel)
    for i ∈ unique(predictions.fold)
        fₚ, fₒ = compute_reliability_curve(predictions[predictions.fold .== i,:])
        fig = plot!(fₚ, fₒ, c=color, alpha=alpha, label=nothing)
    end
    return fig
end

"""
    plot_roc_curves(roc_curves; color, alpha, title, xlabel, ylabel, figsize)

    Plots ROC curves.
"""
function plot_roc_curves(
    roc_curves::Vector;
    color::Symbol=:black,
    alpha::Float64=0.5,
    title::Union{String, LaTeXString}="ROC Curve",
    xlabel::Union{String, LaTeXString}=L"FPR",
    ylabel::Union{String, LaTeXString}=L"TPR",
    figsize::Tuple{Int64,Int64}=(250,250))

    fig = plot([0,1], [0,1], c=color, ls=:dash,
                label=nothing, title=title,
                ylabel=ylabel, xlabel=xlabel, size=figsize)
    K = size(roc_curves)[1]
    for i ∈ 1:K
        plot!(roc_curves[i][1], roc_curves[i][2],
        c=color, alpha=alpha, label=nothing)
    end
    return fig
end

"""
    store_classification_predictions(df, fold, train_sites, test_sites, ytrue, ypred)

    Returns a table with individual level predictions
"""
function store_classification_predictions(df, fold, train_site, test_sites, ytrue, ypred)
    return [df; classification_predictions(fold, train_site, test_sites, ytrue, ypred)]
end

"""
    store_classification_results(df, fold, yhat, ytrue)

    This function takes the results dataframe and adds to it a row for the
    classification results in the current fold.
"""
function store_classification_results(df, fold, yhat, ytrue)
    return [df; classification_report(fold, yhat, coerce(ytrue, OrderedFactor))]
end

"""
    summarize_classification_results(results, sigdigits)

    Returns a formatted version of the results with a column showing mean and 95%ci
    rounded to `sigdigits`.
"""
function summarize_classification_results(results::DataFrame, sigdigits::Int64=4)
    summary_res = describe(results, :mean, :std)
    summary_res[:, :mean] = round.(summary_res[:, :mean]; sigdigits=sigdigits)
    summary_res[:, :se] = round.(summary_res[:, :std] ./ sqrt(K); sigdigits=sigdigits)
    summary_res[:, :lci] = round.(summary_res[:, :mean] .+ quantile(Normal(), 0.025).*summary_res[:, :se]; sigdigits=sigdigits)
    summary_res[:, :uci] = round.(summary_res[:, :mean] .+ quantile(Normal(), 1-0.025).*summary_res[:, :se]; sigdigits=sigdigits)
    summary_res[:, :Result] = map((x,y,z) -> "$x ($y, $z)", summary_res.mean, summary_res.lci, summary_res.uci)
    return summary_res
end


text:
##
## Transformations are
## * givens rotations  Find [c s; -conj(s) conj(c)][a,b] = [x 0];
## * fuse U,V -> UV or UV,D (when complex real)
## * turnover: [   [ -->  [
##               [      [   [
## to this we add abstract "passthrough!" functions

##################################################

"""
    c,s,r  = givensrot(a,b)

where `[c s; -s conj(c)] * [a,b] = [r, 0]
"""
@inline function givensrot(a::T, b::T) where {T <: Real}
    G, r = givens(a,b,1,2)
    s = r >= 0 ? one(T) : -one(T)   #sign(r)
    s*G.c, s*G.s, s*r
end

## This leaves s real
## a, b= rand(Complex{Float64},2)
## G,r = givens(b,a, 1,2)
## c, s = G.s, real(G.c)
## [c s; -s conj(c)] * [a,b]
function givensrot(a::Complex{T},b::Complex{T}) where {T <: Real}
    G, r = givens(b, a, 1, 2)
    G.s, real(G.c), r
end

givensrot(a::Complex{T},b::T) where {T <: Real} =  givensrot(a, complex(b, zero(T)))

## This is suggested in the book. A bit slower, but only one function needed...
## function givensrotX(a::S, b::T) where {T <: Real, S <: Union{T, Complex{T}}}
##     m = max(norm(a), norm(b))
##     x,y = a/m, b/m
##     r = sqrt(norm(x)^2 + y^2)
##     c = x/r
##     s = y/r
##     r = m*r
##     return conj(c),real(s),r
## end

# avoid square root when norm is near 1
#
# 5.1. Backward Stability of the Turnover. Suggests
# this modification to c,s to ensure unitarity
# Assuming that c and s have been produced by valid formulas, the
# magnitude of the correction will be extremely tiny, on the order of
# the unit roundoff. This seemingly trivial correction is crucial. If it
# is neglected, the core transformations will gradually drift away from
# unitarity over the course of many operations, and the algorithm will
# fail. If it is not neglected, the programmer has a good chance of
# producing an accurate turnover.
function approx_givensrot(a::S,b::T) where {S,T}
    delta::T = abs(a*conj(a)) + b^2 - one(T)
    r::T = one(T) - delta/(one(T)+one(T))  # taylor approx for 1/sqrt(1 + delta) when delta ~ 0
    c = a*r
    s = b*r
    return conj(c),s

end

## extra polish; not used
function polish_givens(c::S, s::T) where {S, T}

    delta = real(conj(c)*c) + s*s - one(T)
    r = one(T) -  delta/(one(T) + one(T))
    r * c, r * s

end

##################################################
## Fuse
## fuse combines two rotations, a and b, into one,
## For general rotation, we have no phase consideration, so return an identity rotator
@inline function fuse(a::AbstractRotator{T,S}, b::AbstractRotator{T,S}) where {T,S <: Real}
    #    idx(a) == idx(b) || error("can't fuse")
    i = idx(a)

    c1, s1 = vals(a)
    c2, s2 = vals(b)

    u = c1 * c2 - s1 * conj(s2)
    v = c1 * s2 + s1 * conj(c2)

    u,v = polish_givens(u,v)
    return Rotator(u,v,i), IdentityRotator{T, S}(i)

end

## U * V -> UV [ alpha 0; 0 conj(alpha)]

## For ComplexRealRotator, the result of a*b will not have a real sign
## we output by rotating by alpha.
## We have U*V = (UV) * D
## To get D * (UV) a flip is needed
@inline function fuse(a::Rotator{T,S}, b::Rotator{T,S}) where {T, S <: Complex{T}}

    i = idx(a)
    #    idx(a) == idx(b) || error("can't fuse")
    c1, s1 = vals(a)
    c2, s2 = vals(b)

    u = c1 * c2 - s1 * s2
    v = c1*s2 + s1 * conj(c2)

    s = norm(v)
    alpha =  iszero(v) ? one(Complex{T}) : v/s
    c = u * alpha

    c, s = polish_givens(c, s)

    Rotator(c, s, i), DiagonalRotator(conj(alpha), i)
end

## Fuse for diagonal matrices
fuse(D::IdentityRotator, U) = U, D
function fuse(D::DiagonalRotator, U)
    alpha,_ = vals(D)
    i = idx(D)
    @assert i == idx(U)
    c, s = vals(U)
    Rotator(c*alpha, s, i)
end

##################################################
# Turnover: Q1    Q3   | x x x |      U1
#              Q2    = | x x x | = U3    U2
#                      | x x x |
#
# This is the key computation once matrices are written as rotators
# This function returns
# c4,s4,c5,s5,c6,s6

# following pp15-17 in book
# and the paper
# M = U    W
#       V
# we find U1 by computing a givens rotation to clear [3,1]
#
#     * M =| a x x |
# U1'      | b x x |
#          | 0 x x |
# Then use fact that a^2 + b^2 ~ 1 to use approx givens to find V1 to clear [2,1]
# V1' *     * M  = | 1  0 0   |
#       U1'        | 0  c s   |
#                  | 0 -s cbar|
# The 1 in [1,1] must be 1 and not -1;  so the  lower part is a rotator
# from here [2:3,2:3] is unitary, we pick [3,3] for c and use
# the invariant s1*s2 = s5*s6 to find s6 unless s5=0
#
# question: polishing slows things down, and doesn't seem necessary,
# as approx_givensrot does this; but does make Wilkinson case a bit better
@inline function _turnover(c1::S, s1::T, c2::S, s2::T, c3::S, s3::T) where {S,T}

    ## UVW has first row [UVW11, UVW21, UVW31]
    UVW21 = -c2 * s3 * conj(c1) - c3*s1
    UVW31 = s2 * s3

    if iszero(UVW21) && iszero(UVW31)
        if iszero(s2*conj(c3)) ## this is the case of diagonal matrix Q1*Q2*Q3
            c4 = one(S)
            c5 = c1*c3 - c2*s1*s3
            c6 = c2
            return c4, zero(T), c5, zero(T), c6, zero(T)
        end
    end

    c4, s4, r4 = givensrot(UVW21, UVW31)

    c4, s4 = conj(c4), -s4  # conjugate, as [c4 s4; -s4 conj(s4)]*[a,b] = [r,0]
#    c4, s4 =  polish_givens(c4, s4)

    UVW11 = c1*c3 - c2*s1*s3
    c5, s5::T = approx_givensrot(UVW11, real(r4))
    c5, s5 = conj(c5), -s5
#    c5, s5 =  polish_givens(c5, s5)

    ##  use s1*s2 = s5 * s6 if we can
    ## get a from  M  = V1' * U1' * U *  V  *  W; a = M[3,3]
    if !iszero(s5)
        a = conj(c1) * s2 * s4 + conj(c2) * c4
        b = s1*s2/s5
        c6, s6 = approx_givensrot(a, b)
    else
        a = c2*conj(c1)
        b = -s2
        c6, s6 = givensrot(a, b)
    end
#    c6, s6 = approx_givensrot(a, b)
#    c6, s6 =  polish_givens(c6, s6)

    return (c4, s4, c5, s5, c6, s6)

end

##
##  Turnover interface for rotators
##
function turnover(Q1::R1t,
                  Q2::R2t,
                  Q3::R3t) where {R1t, R2t, R3t}

    c1, s1 = vals(Q1); c2, s2 = vals(Q2); c3,s3 = vals(Q3)
    i,j,k = idx(Q1), idx(Q2), idx(Q3)

    # @assert i == k && (abs(j-i) == 1)

    c4,s4,c5,s5,c6,s6 = _turnover(c1,s1, c2,s2, c3,s3)

    R1 = R3t(c4, s4, j)
    R2 = R2t(c5, s5, i)
    R3 = R1t(c6, s6, j)

    R1, R2, R3

end


text:
# A typing force to keep @eval out of your code.

struct New{T}
    invoke::T
end

unbox(x) = x
unbox(x::New{T}) where T = x.invoke
unbox(args::Tuple) = Tuple(unbox(i) for i in args)

import Base.getindex
getindex(::Type{New},x::Type) = Union{x,New{x}}

f(x::Int,y::Int) = 3

@generated function f(args...)
    @show args

    if findfirst(x->x<:New,args)==nothing
        return quote
            error("A method error. Please define what to do for f$args")
        end
    else
        return quote
            return f(unbox(args)...)
        end
    end
end


text:
##################################################################
# Muestra en pantalla los diagramas de bigotes
##################################################################
using PlotlyJS
using CSV
using DataFrames
limit = 16
path = "Memoria/img/7-algoritmo-inicializar-pesos/experimento/"

# Imprimimos  tabla bigotes tiempo
df = DataFrame(CSV.File("Experimentos/inicializacion-pesos-red-neuronal/resultados/2_air_self_noise/tiempos.csv"))
table_df = Tables.matrix(df)
traces = [
    box(
        y = table_df[i, 2:16],
        name = table_df[i, 1]
    )
    for i in 1:2
]
ref = plot(traces, Layout(yaxis_title="Tiempos en segundos", boxmode="group"))
savefig(ref, path*"grafico-bigotes-tiempo.png")

# Error cuadrático medio
file_path = "Experimentos/inicializacion-pesos-red-neuronal/resultados/2_air_self_noise/error_entrenamiento.csv"
df = DataFrame(CSV.File(file_path))
table_df = Tables.matrix(df)
traces = [
    box(
        y = table_df[i, 2:limit],
        name = table_df[i, 1]
    )
    for i in 1:2
]
ref = plot(traces, Layout(yaxis_title="Error cuadrático medio", boxmode="group"))
savefig(ref, path*"grafico-bigotes-error_entrenamiento.png")

# Error cuadrático medio
file_path = "Experimentos/inicializacion-pesos-red-neuronal/resultados/2_air_self_noise/error_test.csv"
df = DataFrame(CSV.File(file_path))
table_df = Tables.matrix(df)
traces = [
    box(
        y = table_df[i, 2:limit],
        name = table_df[i, 1]
    )
    for i in 1:2
]
ref = plot(traces, Layout(yaxis_title="Error cuadrático medio", boxmode="group"))
path = "Memoria/img/7-algoritmo-inicializar-pesos/experimento/"
savefig(ref, path*"grafico-bigotes-error_test.png")


text:
@testset "EpsBoxArchive" begin
    # check that frontier elements are mutually nondominated and
    # ther are no epsilon-index duplicates
    function check_frontier(a::EpsBoxArchive)
        for (i, el_i) in enumerate(pareto_frontier(a))
            for (j, el_j) in enumerate(pareto_frontier(a))
                if j > i
                    comp, index_match = hat_compare(fitness(el_i), fitness(el_j), fitness_scheme(a))
                    if comp != 0 || index_match
                        @warn "Dominated pair: fitness($i)=$(fitness(el_i)) fitness($j)=$(fitness(el_j))"
                        return comp, index_match # return if violating pair is found
                    end
                end
            end
        end
        return (0, false) # all non-dominated, no index matches
    end

    scheme = EpsBoxDominanceFitnessScheme{2}(0.1, is_minimizing=true)

    @testset "EpsBoxFrontierIndividual" begin
        arch_indiv = BlackBoxOptim.EpsBoxFrontierIndividual(
                        convert(IndexedTupleFitness, (2.21, 1.12), scheme), [3.0, 2.0], 3, 5, 2)
        @test arch_indiv.fitness.orig == (2.21, 1.12)
        @test arch_indiv.fitness.index == (22, 11)
        @test arch_indiv.timestamp > 0
        @test arch_indiv.tag == 3
        @test arch_indiv.num_fevals == 5
        @test arch_indiv.n_restarts == 2
    end

    @testset "Constructing a small archive and adding to it" begin
        a = EpsBoxArchive(scheme, max_size=100)

        @test capacity(a) == 100
        @test length(a)   == 0
        @test isempty(a)
        @test best_candidate(a) === nothing
        @test isnafitness(best_fitness(a), fitness_scheme(a))
        @test BlackBoxOptim.noprogress_streak(a) == 0
        @test BlackBoxOptim.tagcounts(a) == Dict{Int,Int}()

        BlackBoxOptim.add_candidate!(a, convert(IndexedTupleFitness, (2.21, 1.12), scheme), [0.0, 5.0], 1, 8)
        @test capacity(a)         == 100
        @test length(a)           == 1
        @test !isempty(a)
        @test best_fitness(a) == fitness(first(pareto_frontier(a)))
        @test best_candidate(a) === params(first(pareto_frontier(a)))
        @test BlackBoxOptim.best_front_elem(a) == first(pareto_frontier(a))
        @test first(pareto_frontier(a)).num_fevals == 8
        @test BlackBoxOptim.noprogress_streak(a) == 0
        @test BlackBoxOptim.tagcounts(a) == Dict(1=>1)

        # dominated candidate is ignored and best candidate unchanged
        BlackBoxOptim.add_candidate!(a, (3.21, 1.12), [1.0, 5.0], 2)
        @test capacity(a)         == 100
        @test length(a)           == 1
        @test BlackBoxOptim.best_front_elem(a) == first(pareto_frontier(a))
        @test first(pareto_frontier(a)).num_fevals == 8
        @test BlackBoxOptim.noprogress_streak(a) == 1
        @test BlackBoxOptim.tagcounts(a) == Dict(1=>1)

        # non-dominated candidate is appended and best candidate updated
        BlackBoxOptim.add_candidate!(a, (1.21, 2.11), [1.2, 3.0], 3)
        @test capacity(a)         == 100
        @test length(a)           == 2
        @test findfirst(a.frontier, SI.Point((12,21))) !== nothing
        candi_leaf, candi_ix = findfirst(a.frontier, SI.Point((12,21)))
        @test candi_leaf[candi_ix].num_fevals == 3
        @test best_fitness(a).orig == (1.21, 2.11)
        @test best_candidate(a) == [1.2, 3.0]
        @test BlackBoxOptim.best_front_elem(a) == candi_leaf[candi_ix]
        @test BlackBoxOptim.noprogress_streak(a) == 0
        @test BlackBoxOptim.tagcounts(a) == Dict(1=>1, 3=>1)

        # another non-dominated candidate (no best update)
        BlackBoxOptim.add_candidate!(a, (0.52, 3.15), [1.5, 3.0], 3)
        @test capacity(a)         == 100
        @test length(a)           == 3
        @test best_fitness(a).orig == (1.21, 2.11)
        @test BlackBoxOptim.noprogress_streak(a) == 0
        @test BlackBoxOptim.tagcounts(a) == Dict(1=>1, 3=>2)

        # the candidate that dominates all current frontier
        BlackBoxOptim.add_candidate!(a, (0.22, 0.22), [1.5, 3.0], 5)
        @test capacity(a)         == 100
        @test length(a)           == 1
        @test best_fitness(a).orig == (0.22, 0.22)
        @test BlackBoxOptim.noprogress_streak(a) == 0
        @test BlackBoxOptim.tagcounts(a) == Dict(5=>1)

        BlackBoxOptim.add_candidate!(a, (0.21, 0.21), [1.5, 3.0], 6)
        @test capacity(a)         == 100
        @test length(a)           == 1
        @test best_fitness(a).orig == (0.21, 0.21)
        @test BlackBoxOptim.noprogress_streak(a) == 1
        @test BlackBoxOptim.tagcounts(a) == Dict(6=>1)

        @testset "noprogress_streak(since_restart=0) is reset after restart" begin
            @test BlackBoxOptim.noprogress_streak(a, since_restart=false) == 1
            @test BlackBoxOptim.noprogress_streak(a, since_restart=true) == 1
            BlackBoxOptim.notify!(a, :restart)
            @test BlackBoxOptim.noprogress_streak(a, since_restart=false) == 1
            @test BlackBoxOptim.noprogress_streak(a, since_restart=true) == 0
            BlackBoxOptim.add_candidate!(a, (0.2, 0.2), [1.5, 3.0], 6)
            @test BlackBoxOptim.noprogress_streak(a, since_restart=false) == 2
            @test BlackBoxOptim.noprogress_streak(a, since_restart=true) == 1
        end

        @testset "updating best candidate index if old and new non-dominated" begin
            a = EpsBoxArchive(scheme, max_size=100)

            BlackBoxOptim.add_candidate!(a, convert(IndexedTupleFitness, (1.0, 0.0), scheme), [0.0, 1.0], 1)
            @test best_fitness(a).orig == (1.0, 0.0)
            @test length(a)           == 1
            BlackBoxOptim.add_candidate!(a, convert(IndexedTupleFitness, (0.6, 0.6), scheme), [0.0, 2.0], 2)
            @test best_fitness(a).orig == (1.0, 0.0)
            @test best_candidate(a) == [0.0, 1.0]
            @test length(a)           == 2
            BlackBoxOptim.add_candidate!(a, convert(IndexedTupleFitness, (0.0, 1.1), scheme), [0.0, 3.0], 3)
            @test best_fitness(a).orig == (1.0, 0.0)
            @test length(a)           == 3
            BlackBoxOptim.add_candidate!(a, convert(IndexedTupleFitness, (0.0, 0.9), scheme), [0.0, 4.0], 4)
            @test best_fitness(a).orig == (0.0, 0.9)
            @test best_candidate(a) == [0.0, 4.0]
            @test length(a)           == 3
        end

        @testset "handling dulicate elements" begin
            a = EpsBoxArchive(scheme, max_size=100)

            BlackBoxOptim.add_candidate!(a, convert(IndexedTupleFitness, (1.25, 0.0), scheme), [0.0, 1.0], 1)
            @test BlackBoxOptim.best_front_elem(a) == first(pareto_frontier(a))
            @test length(a)         == 1
            BlackBoxOptim.add_candidate!(a, convert(IndexedTupleFitness, (0.6, 0.6), scheme), [0.0, 2.0], 2)
            @test best_fitness(a).orig == (0.6, 0.6)
            @test length(a)         == 2
            BlackBoxOptim.add_candidate!(a, convert(IndexedTupleFitness, (0.0, 1.3), scheme), [0.0, 3.0], 3)
            @test best_fitness(a).orig == (0.6, 0.6)
            @test best_candidate(a) == [0.0, 2.0] # no change
            @test length(a)         == 3
            BlackBoxOptim.add_candidate!(a, convert(IndexedTupleFitness, (0.6, 0.6), scheme), [0.0, 4.0], 4)
            @test best_fitness(a).orig == (0.6, 0.6)
            @test best_candidate(a) == [0.0, 2.0] # same fitness, so no change
            @test length(a)         == 3
        end
    end

    @testset "shaffer1() Pareto frontier" begin
        schaffer1(x) = (sum(abs2, x), sum(xx -> abs2(xx - 2.0), x))
        scheme = EpsBoxDominanceFitnessScheme{2}(0.1, is_minimizing=true)
        a = EpsBoxArchive(scheme)

        for t in range(0, stop=2, length=1000)
            params = [t, t]
            BlackBoxOptim.add_candidate!(a, schaffer1(params), params)
        end

        @test length(a) == 40
    end

    @testset "archive copies the individuals" begin
        a = EpsBoxArchive(scheme, max_size=100)
        @test length(a) == 0

        indiv = [1.5, 2.0]
        BlackBoxOptim.add_candidate!(a, (1.5, 2.0), indiv)
        # equal but not identical
        @test first(a.frontier).params == indiv
        @test first(a.frontier).params !== indiv

        # modify the vector
        indiv[1] = 5.0
        # test that archived version is still the same
        @test first(a.frontier).params == [1.5, 2.0]
    end

    @testset "simulate optimization" begin
        # generate random fitness from (0..2, 0..2, 0..2)
        # and lying inside the (1.0-r1)^1.5 + (1.0-r2)^1.5 + (1.0-r3)^1.5 <= 8.0
        function rfitness()
            while true
                r = ntuple(_ -> rand(), 3)
                if sum(x -> (1 - x)^1.5, r) <= 1.0
                    return 2.0 .* r
                end
            end
        end

        scheme3 = EpsBoxDominanceFitnessScheme{3}(0.1, is_minimizing=true)
        a = EpsBoxArchive(scheme3, max_size=1000)
        indiv = [0.0, 0.0]
        for i in 1:10000
            BlackBoxOptim.add_candidate!(a, rfitness(), indiv)
            # periodically check that pareto_frontier(a) is indeed
            # epsilon-index frontier
            if mod(i, 100) == 0
                @test length(a) > 0
                @test check_frontier(a) == (0, false)
            end
        end
    end
end


text:
#Possible activation functions
function sigmoid(x::Number)
    return exp(x)/(exp(x)+1)
end

function relu(x::Number)
    if x > 0
        return x
    else
        return 0
    end
end

#Not sure if any of the below is in
#the right direction
struct Layer
    nodes :: Number
end

mutable struct buildabob
    activation :: Function
    hiddenlayers :: Number
    firstlayer :: Layer
    lastlayer :: Layer
end

max_stars_repo_name: ian-double-u/NNfromScratch

max_issues_repo_name: ian-double-u/NNfromScratch

max_forks_repo_name: ian-double-u/NNfromScratch


text:
using NeuralNetworkAnalysis

@taylorize function sys!(dx, x, p, t)
    dx[1] = x[2]
    dx[2] = zero(x[2])
    return dx
end;

function oscillator(X::LazySet)
    if high(X)[1] > 1
        if low(X)[1] < 1
            res = Interval(-0.1, 0.1)
        else
            res = Interval(-0.1, -0.1)
        end
    else
        res = Interval(0.1, 0.1)
    end
    return res
end;

function oscillator(x::Vector)
    if x[1] > 1
        res = [-0.1]
    else
        res = [0.1]
    end
    return res
end;
controller = BlackBoxController(oscillator);

X₀ = Interval(0.95, 0.99);
U = ZeroSet(1);
vars_idx = Dict(:state_vars=>[1], :control_vars=>2);
ivp = @ivp(x' = sys!(x), dim: 2, x(0) ∈ X₀ × U);
period = 0.1;
prob = ControlledPlant(ivp, controller, vars_idx, period);
alg = TMJets21b(abstol=1e-10, orderT=8, orderQ=2, adaptive=true);
alg_nn = BlackBoxSolver();
N = 10;
T = N*period;

sol_raw = solve(prob, T=T, alg_nn=alg_nn, alg=alg);

import DifferentialEquations
sim = simulate(prob, T=T; trajectories=10, include_vertices=true);

# using Plots
# fig = plot()
# vars = (0, 1)
# plot!(fig, sol_raw, vars=vars, c=:yellow, lw=2.0, ls=:dash)
# plot_simulation!(fig, sim; vars=vars, color=:red, lab="");
# savefig("oscillator_$(N)_steps")


text:
using Distributions

include("create_maze.jl")

struct CollectPointsCpuCfg
    maze_columns::Int
    maze_rows::Int
    maze_cell_size::Int
    agent_radius::Int
    point_radius::Int
    agent_movement_range::Float32
    number_of_sensors::Int
    reward_per_collected_positive_point::Float32
    reward_per_collected_negative_point::Float32
    number_time_steps::Int

    function CollectPointsCpuCfg(config::OrderedDict)

        new(
            config["maze_columns"],
            config["maze_rows"],
            config["maze_cell_size"],
            config["agent_radius"],
            config["point_radius"],
            config["agent_movement_range"],
            4,
            config["reward_per_collected_positive_point"],
            config["reward_per_collected_negative_point"],
            config["number_time_steps"],
        )

    end
end

mutable struct CollectPointsCpu
    config::CollectPointsCpuCfg
    screen_width::Int
    screen_height::Int
    agent_position_x::Int
    agent_position_y::Int
    point_x::Int
    point_y::Int
    negative_point_x::Int
    negative_point_y::Int
    maze::Any
    sensor_top::Float32
    sensor_bottom::Float32
    sensor_left::Float32
    sensor_right::Float32
    t::Float32

    function CollectPointsCpu(config::CollectPointsCpuCfg)

        screen_width = config.maze_cell_size * config.maze_columns
        screen_height = config.maze_cell_size * config.maze_rows

        # Agent coordinates
        agent_position_x, agent_position_y = place_randomly_in_maze(config, config.agent_radius)

        # Point coordinates for positive point
        point_x, point_y = place_randomly_in_maze(config, config.point_radius)

        # Point coordinates for negative point
        negative_point_x, negative_point_y = place_randomly_in_maze(config, config.point_radius)

        # Create Maze
        maze = make_maze(config.maze_columns, config.maze_rows)

        sensor_top = 0.0
        sensor_bottom = 0.0
        sensor_left = 0.0
        sensor_right = 0.0

        t = 0

        new(
            configuration,
            screen_width,
            screen_height,
            agent_position_x,
            agent_position_y,
            point_x,
            point_y,
            negative_point_x,
            negative_point_y,
            maze,
            sensor_top,
            sensor_bottom,
            sensor_left,
            sensor_right,
            t,
        )

    end
end

function step(env::CollectPointsCpu)

    action = rand(Uniform(-1.0,1.0), 2)

    # Movement range for agent
    mr = env.config.agent_movement_range

    # Move agent
    env.agent_position_x += convert(Int, round(clamp(action[1] * mr, -mr, mr)))
    env.agent_position_y += convert(Int, round(clamp(action[2] * mr, -mr, mr)))

    # Check agent collisions with outer walls
    env.agent_position_y = max(env.agent_position_y, env.config.agent_radius)  # Upper border
    env.agent_position_y = min(env.agent_position_y, env.screen_height - env.config.agent_radius)  # Lower bord.
    env.agent_position_x = min(env.agent_position_x, env.screen_width - env.config.agent_radius)  # Right border
    env.agent_position_x = max(env.agent_position_x, env.config.agent_radius)  # Left border

    # Get cell indizes of agents current position
    cell_x = convert(Int, ceil(env.agent_position_x / env.config.maze_cell_size))
    cell_y = convert(Int, ceil(env.agent_position_y / env.config.maze_cell_size))

    # Get current cell
    cell = env.maze[cell_x, cell_y]

    # Get coordinates of current cell
    x_left, x_right, y_top, y_bottom = get_coordinates_maze_cell(cell_x, cell_y, env.config.maze_cell_size)

    # Check agent collisions with maze walls
    if cell.walls[North] == true
        env.agent_position_y = max(env.agent_position_y, y_top + env.config.agent_radius)
    end
    if cell.walls[South] == true
        env.agent_position_y = min(env.agent_position_y, y_bottom - env.config.agent_radius)
    end
    if cell.walls[East] == true
        env.agent_position_x = min(env.agent_position_x, x_right - env.config.agent_radius)
    end
    if cell.walls[West] == true
        env.agent_position_x = max(env.agent_position_x, x_left + env.config.agent_radius)
    end

    reward = 0.0

    # Collect positive point in reach
    distance = sqrt((env.point_x - env.agent_position_x)^2 + (env.point_y - env.agent_position_y)^2)
    if distance <= env.config.point_radius + env.config.agent_radius
        env.point_x, env.point_y = place_randomly_in_maze(env.config, env.config.point_radius)
        reward += env.config.reward_per_collected_positive_point
    end

    # Collect negative point in reach
    distance = sqrt((env.negative_point_x - env.agent_position_x)^2 + (env.negative_point_y - env.agent_position_y)^2)
    if distance <= env.config.point_radius + env.config.agent_radius
        env.negative_point_x, env.negative_point_y = place_randomly_in_maze(env.config, env.config.point_radius)
        reward += env.config.reward_per_collected_negative_point
    end

    println(reward)

    return

end

function place_randomly_in_maze(config::CollectPointsCpuCfg, radius::Int)

    x = rand(radius:(config.maze_cell_size - radius)) + rand(0:config.maze_columns-1) * config.maze_cell_size
    y = rand(radius:(config.maze_cell_size - radius)) + rand(0:config.maze_rows-1) * config.maze_cell_size

    return x, y
end


text:
#InspectDR: Basic math operations
#-------------------------------------------------------------------------------
#=NOTE:
These tools should eventually be moved to a separate unit.
=#

#==Useful tests
===============================================================================#

#Verifies that v is strictly increasing (no repeating values):
#TODO: support repeating values (non-strict)
function isincreasing(v::Vector)
	if length(v) < 1; return true; end
	prev = v[1]

	for x in v[2:end]
		if !(x > prev) #Make sure works for NaN
			return false
		end
	end
	return true
end

isincreasing(r::AbstractRange) = (step(r) > 0)

#==Basic operations
===============================================================================#

#Safe version of extrema (returns DNaN on error):
#Also ignores NaNs.  TODO: optimize?
function extrema_nan(v::Vector)
	vmin = DInf; vmax = -DInf
	for elem in v
		if !isnan(elem)
			if elem > vmax
				vmax = elem
			end
			if elem < vmin
				vmin = elem
			end
		end
	end
	if vmin>vmax
		vmin = vmax = DNaN
	end
	return (vmin, vmax)
end
#Last line