text:
"""
#### Project variables or derived quantities from the **hydro-dataset**:
- projection to a grid related to the maximum level of the loaded data
- overview the list of predefined quantities with: projection()
- select variable(s) and their unit(s)
- limit to a maximum range
- select a coarser grid than the maximum resolution of the loaded data (maps with both resolutions are created)
- give the spatial center (with units) of the data within the box (relevant e.g. for radius dependency)
- relate the coordinates to a direction (x,y,z)
- select between mass (default) and volume weighting
- pass a mask to exclude elements (cells/particles/...) from the calculation
- toggle verbose mode

```julia
projection(   dataobject::HydroDataType, vars::Array{Symbol,1};
                        units::Array{Symbol,1}=[:standard],
                        lmax::Real=dataobject.lmax,
                        mask=[false],
                        direction::Symbol=:z,
                        weighting::Symbol=:mass,
                        xrange::Array{<:Any,1}=[missing, missing],
                        yrange::Array{<:Any,1}=[missing, missing],
                        zrange::Array{<:Any,1}=[missing, missing],
                        center::Array{<:Any,1}=[0., 0., 0.],
                        range_unit::Symbol=:standard,
                        data_center::Array{<:Any,1}=[missing, missing, missing],
                        data_center_unit::Symbol=:standard,
                        verbose::Bool=verbose_mode)

return HydroMapsType

```

#### Arguments
##### Required:
- **`dataobject`:** needs to be of type: "HydroDataType"
- **`var(s)`:** select a variable from the database or a predefined quantity (see field: info, function projection(), dataobject.data)
##### Predefined/Optional Keywords:
- **`unit(s)`:** return the variable in given units
- **`lmax`:** create maps with coarser grid than provided by the maximum level of the loaded data
- **`xrange`:** the range between [xmin, xmax] in units given by argument `range_unit` and relative to the given `center`; zero length for xmin=xmax=0. is converted to maximum possible length
- **`yrange`:** the range between [ymin, ymax] in units given by argument `range_unit` and relative to the given `center`; zero length for ymin=ymax=0. is converted to maximum possible length
- **`zrange`:** the range between [zmin, zmax] in units given by argument `range_unit` and relative to the given `center`; zero length for zmin=zmax=0. is converted to maximum possible length
- **`range_unit`:** the units of the given ranges: :standard (code units), :Mpc, :kpc, :pc, :mpc, :ly, :au , :km, :cm (of typye Symbol) ..etc. ; see for defined length-scales viewfields(info.scale)
- **`center`:** in units given by argument `range_unit`; by default [0., 0., 0.]; the box-center can be selected by e.g. [:bc], [:boxcenter], [value, :bc, :bc], etc..
- **`weighting`:** select between `:mass` weighting (default) and `:volume` weighting
- **`data_center`:** to calculate the data relative to the data_center; in units given by argument `data_center_unit`; by default the argument data_center = center ;
- **`data_center_unit`:** :standard (code units), :Mpc, :kpc, :pc, :mpc, :ly, :au , :km, :cm (of typye Symbol) ..etc. ; see for defined length-scales viewfields(info.scale)
- **`direction`:** select between: :x, :y, :z
- **`mask`:** needs to be of type MaskType which is a supertype of Array{Bool,1} or BitArray{1} with the length of the database (rows)

### Defined Methods - function defined for different arguments

- projection( dataobject::HydroDataType, var::Symbol; ...) # one given variable
- projection( dataobject::HydroDataType, var::Symbol, unit::Symbol; ...) # one given variable with its unit
- projection( dataobject::HydroDataType, vars::Array{Symbol,1}; ...) # several given variables -> array needed
- projection( dataobject::HydroDataType, vars::Array{Symbol,1}, units::Array{Symbol,1}; ...) # several given variables and their corresponding units -> both arrays
- projection( dataobject::HydroDataType, vars::Array{Symbol,1}, unit::Symbol; ...)  # several given variables that have the same unit -> array for the variables and a single Symbol for the unit

#### Examples
...
"""
function projection(   dataobject::HydroDataType, var::Symbol;
                        unit::Symbol=:standard,
                        lmax::Real=dataobject.lmax,
                        mask=[false],
                        direction::Symbol=:z,
                        #plane_orientation::Symbol=:perpendicular,
                        weighting::Symbol=:mass,
                        xrange::Array{<:Any,1}=[missing, missing],
                        yrange::Array{<:Any,1}=[missing, missing],
                        zrange::Array{<:Any,1}=[missing, missing],
                        center::Array{<:Any,1}=[0., 0., 0.],
                        range_unit::Symbol=:standard,
                        data_center::Array{<:Any,1}=[missing, missing, missing],
                        data_center_unit::Symbol=:standard,
                        verbose::Bool=verbose_mode)

    return projection(dataobject, [var], units=[unit],
                            lmax=lmax,
                            mask=mask,
                            direction=direction,
                            #plane_orientation=plane_orientation,
                            weighting=weighting,
                            xrange=xrange,
                            yrange=yrange,
                            zrange=zrange,
                            center=center,
                            range_unit=range_unit,
                            data_center=data_center,
                            data_center_unit=data_center_unit,
                            verbose=verbose)

end

function projection(   dataobject::HydroDataType, var::Symbol, unit::Symbol;
                        lmax::Real=dataobject.lmax,
                        mask=[false],
                        direction::Symbol=:z,
                        #plane_orientation::Symbol=:perpendicular,
                        weighting::Symbol=:mass,
                        xrange::Array{<:Any,1}=[missing, missing],
                        yrange::Array{<:Any,1}=[missing, missing],
                        zrange::Array{<:Any,1}=[missing, missing],
                        center::Array{<:Any,1}=[0., 0., 0.],
                        range_unit::Symbol=:standard,
                        data_center::Array{<:Any,1}=[missing, missing, missing],
                        data_center_unit::Symbol=:standard,
                        verbose::Bool=verbose_mode)

    return projection(dataobject, [var], units=[unit],
                            lmax=lmax,
                            mask=mask,
                            direction=direction,
                            #plane_orientation=plane_orientation,
                            weighting=weighting,
                            xrange=xrange,
                            yrange=yrange,
                            zrange=zrange,
                            center=center,
                            range_unit=range_unit,
                            data_center=data_center,
                            data_center_unit=data_center_unit,
                            verbose=verbose)

end

function projection(   dataobject::HydroDataType, vars::Array{Symbol,1}, units::Array{Symbol,1};
                        lmax::Real=dataobject.lmax,
                        mask=[false],
                        direction::Symbol=:z,
                        #plane_orientation::Symbol=:perpendicular,
                        weighting::Symbol=:mass,
                        xrange::Array{<:Any,1}=[missing, missing],
                        yrange::Array{<:Any,1}=[missing, missing],
                        zrange::Array{<:Any,1}=[missing, missing],
                        center::Array{<:Any,1}=[0., 0., 0.],
                        range_unit::Symbol=:standard,
                        data_center::Array{<:Any,1}=[missing, missing, missing],
                        data_center_unit::Symbol=:standard,
                        verbose::Bool=verbose_mode)

    return projection(dataobject, vars, units=units,
                                                lmax=lmax,
                                                mask=mask,
                                                direction=direction,
                                                #plane_orientation=plane_orientation,
                                                weighting=weighting,
                                                xrange=xrange,
                                                yrange=yrange,
                                                zrange=zrange,
                                                center=center,
                                                range_unit=range_unit,
                                                data_center=data_center,
                                                data_center_unit=data_center_unit,
                                                verbose=verbose)

end

function projection(   dataobject::HydroDataType, vars::Array{Symbol,1}, unit::Symbol;
                        lmax::Real=dataobject.lmax,
                        mask=[false],
                        direction::Symbol=:z,
                        #plane_orientation::Symbol=:perpendicular,
                        weighting::Symbol=:mass,
                        xrange::Array{<:Any,1}=[missing, missing],
                        yrange::Array{<:Any,1}=[missing, missing],
                        zrange::Array{<:Any,1}=[missing, missing],
                        center::Array{<:Any,1}=[0., 0., 0.],
                        range_unit::Symbol=:standard,
                        data_center::Array{<:Any,1}=[missing, missing, missing],
                        data_center_unit::Symbol=:standard,
                        verbose::Bool=verbose_mode)

    return projection(dataobject, vars, units=fill(unit, length(vars)),
                                                lmax=lmax,
                                                mask=mask,
                                                direction=direction,
                                                #plane_orientation=plane_orientation,
                                                weighting=weighting,
                                                xrange=xrange,
                                                yrange=yrange,
                                                zrange=zrange,
                                                center=center,
                                                range_unit=range_unit,
                                                data_center=data_center,
                                                data_center_unit=data_center_unit,
                                                verbose=verbose)

end

#todo: check for uniform grid
function projection(   dataobject::HydroDataType, vars::Array{Symbol,1};
                        units::Array{Symbol,1}=[:standard],
                        lmax::Real=dataobject.lmax,
                        mask=[false],
                        direction::Symbol=:z,
                        #plane_orientation::Symbol=:perpendicular,
                        weighting::Symbol=:mass,
                        xrange::Array{<:Any,1}=[missing, missing],
                        yrange::Array{<:Any,1}=[missing, missing],
                        zrange::Array{<:Any,1}=[missing, missing],
                        center::Array{<:Any,1}=[0., 0., 0.],
                        range_unit::Symbol=:standard,
                        data_center::Array{<:Any,1}=[missing, missing, missing],
                        data_center_unit::Symbol=:standard,
                        verbose::Bool=verbose_mode)

    printtime("", verbose)

    lmin = dataobject.lmin
    #lmax = dataobject.lmax
    simlmax=dataobject.lmax
    #simlmax=lmax
    Nlevel = simlmax-lmin
    boxlen = dataobject.boxlen

    #ranges = [xrange[1],xrange[1],yrange[1],yrange[1],zrange[1],zrange[1]]
    scale = dataobject.scale
    nvarh = dataobject.info.nvarh
    lmax_projected = lmax
    isamr = checkuniformgrid(dataobject, lmax)

    selected_vars = vars #unique(vars)

    sd_names = [:sd, :Σ, :surfacedensity]
    density_names = [:density, :rho, :ρ]

    # checks to use maps instead of projections
    rcheck = [:r_cylinder, :r_sphere]
    anglecheck = [:ϕ]
    ranglecheck = [rcheck..., anglecheck...]
    # for velocity dispersion add necessary velocity components
    # ========================================================
    σcheck = [:σx, :σy, :σz, :σ, :σr_cylinder, :σϕ_cylinder]
    rσanglecheck = [rcheck...,σcheck...,anglecheck...]

    σ_to_v = SortedDict(  :σx => [:vx, :vx2],
                          :σy => [:vy, :vy2],
                          :σz => [:vz, :vz2],
                          :σ  => [:v,  :v2],
                          :σr_cylinder => [:vr_cylinder, :vr_cylinder2],
                          :σϕ_cylinder => [:vϕ_cylinder, :vϕ_cylinder2] )

    for i in σcheck
        idx = findall(x->x==i, selected_vars) #[1]
        if length(idx) >= 1
            selected_v = σ_to_v[i]
            for j in selected_v
                jdx = findall(x->x==j, selected_vars)
                if length(jdx) == 0
                    append!(selected_vars, [j])
                end
            end
        end
    end
    # ========================================================

    notonly_ranglecheck_vars = checkformaps(selected_vars, ranglecheck)

    if weighting == :mass
        # only add :sd if there are also other variables than in ranglecheck
        if !in(:sd, selected_vars) && notonly_ranglecheck_vars
            append!(selected_vars, [:sd])
        end

        if !in(:rho, keys(dataobject.data[1]) )
            error("""[Mera]: For mass weighting variable "rho" is necessary.""")
        end
    end

    # convert given ranges and print overview on screen
    ranges = prepranges(dataobject.info,range_unit, verbose, xrange, yrange, zrange, center, dataranges=dataobject.ranges)

    data_centerm = prepdatacenter(dataobject.info, center, range_unit, data_center, data_center_unit)

    xmin, xmax, ymin, ymax, zmin, zmax = ranges

    if verbose
        println("Selected var(s)=$(tuple(selected_vars...)) ")
        println()
    end

    # rebin data on the maximum loaded grid
    r1 = floor(Int, ranges[1] * (2^simlmax)) + 1
    r2 = ceil(Int, ranges[2] * (2^simlmax))  + 1
    r3 = floor(Int, ranges[3] * (2^simlmax)) + 1
    r4 = ceil(Int, ranges[4] * (2^simlmax))  + 1
    r5 = floor(Int, ranges[5] * (2^simlmax)) + 1
    r6 = ceil(Int, ranges[6] * (2^simlmax))  + 1

    x_coord = :cx
    y_coord = :cy
    rl = data_centerm .* 2^simlmax #.* dataobject.boxlen

    if direction == :z
        # range on maximum used grid

        newrange1 = range(r1, stop=r2, length=(r2-r1))
        newrange2 = range(r3, stop=r4, length=(r4-r3))
        #println(newrange1)
        #println(newrange2)

        x_coord = :cx
        y_coord = :cy
        extent=[r1-1,r2-1,r3-1,r4-1]
        ratio = (extent[2]-extent[1]) / (extent[4]-extent[3])
        extent_center= [extent[1]-rl[1], extent[2]-rl[1], extent[3]-rl[2], extent[4]-rl[2]] .* boxlen ./ 2^simlmax
        extent = extent .* boxlen ./ 2^simlmax
        length1_center = (data_centerm[1] -xmin ) * boxlen
        length2_center = (data_centerm[2] -ymin ) * boxlen

    elseif direction == :y
        # range on maximum used grid
        newrange1 = range(r1, stop=r2, length=(r2-r1))
        newrange2 = range(r5, stop=r6, length=(r6-r5))
        #println(newrange1)
        #println(newrange2)s

        x_coord = :cx
        y_coord = :cz
        extent=[r1-1,r2-1,r5-1,r6-1]
        ratio = (extent[2]-extent[1]) / (extent[4]-extent[3])
        extent_center= [extent[1]-rl[1], extent[2]-rl[1], extent[3]-rl[3], extent[4]-rl[3]] .* boxlen ./ 2^simlmax
        extent = extent .* boxlen ./ 2^simlmax
        length1_center = (data_centerm[1] - xmin  ) * boxlen
        length2_center = (data_centerm[3] - zmin) * boxlen

    elseif direction == :x
        # range on maximum used grid
        newrange1 = range(r3, stop=r4, length=(r4-r3))
        newrange2 = range(r5, stop=r6, length=(r6-r5))
        #println(newrange1)
        #println(newrange2)
        x_coord = :cy
        y_coord = :cz
        extent=[r3-1,r4-1,r5-1,r6-1]
        ratio = (extent[2]-extent[1]) / (extent[4]-extent[3])
        extent_center= [extent[1]-rl[2], extent[2]-rl[2], extent[3]-rl[3], extent[4]-rl[3]] .* boxlen ./ 2^simlmax
        extent = extent .* boxlen ./ 2^simlmax
        length1_center = (data_centerm[2] -ymin ) * boxlen
        length2_center = (data_centerm[3] -zmin) * boxlen
    end

    length1=length( newrange1)
    length2=length( newrange2)
    map = zeros(Float64, length1, length2, length(selected_vars)  )
    map_weight = zeros(Float64, length1 , length2, length(selected_vars) );

    rows = length(dataobject.data)
    mera_mask_inserted = false
    if length(mask) > 1
        if length(mask) !== rows
            error("[Mera] ",now()," : array-mask length: $(length(mask)) does not match with data-table length: $(rows)")
        else
            if in(:mask, colnames(dataobject.data))
                if verbose
                    println(":mask provided by datatable")
                    println()
                end
            else
                Nafter = JuliaDB.ncols(dataobject.data)
                dataobject.data = JuliaDB.insertcolsafter(dataobject.data, Nafter, :mask => mask)
                if verbose
                    println(":mask provided by function")
                    println()
                end
                mera_mask_inserted = true
            end
        end
    end

    alt_shift = 2
    closed=:left
    maps = SortedDict( )
    maps_unit = SortedDict( )
    maps_lmax = SortedDict( )
    maps_mode = SortedDict( )
    if notonly_ranglecheck_vars
        @showprogress 1 "" for level = lmin:simlmax

            first_time_level = fill(1, length(selected_vars) )
            if isamr
                #level_data = filter(row->row.level == level, dataobject.data)
                level_data = filter(p-> p.level == level &&
                                        p.cx >=floor(Int, 2^p.level * xmin) &&
                                        p.cx <=ceil(Int,  2^p.level * xmax) &&
                                        p.cy >=floor(Int, 2^p.level * ymin) &&
                                        p.cy <=ceil(Int,  2^p.level * ymax) &&
                                        p.cz >=floor(Int, 2^p.level * zmin) &&
                                        p.cz <=ceil(Int,  2^p.level * zmax), dataobject.data)
            else # for uniform grid
                #level_data = dataobject.data
                level_data = filter(p-> p.cx >=floor(Int, 2^lmax * xmin) &&
                                        p.cx <=ceil(Int,  2^lmax * xmax) &&
                                        p.cy >=floor(Int, 2^lmax * ymin) &&
                                        p.cy <=ceil(Int,  2^lmax * ymax) &&
                                        p.cz >=floor(Int, 2^lmax * zmin) &&
                                        p.cz <=ceil(Int,  2^lmax * zmax), dataobject.data)
            end

            # rebin data on the used level grid
            rl1 = floor(Int, ranges[1] * (2^level ))  + 1
            rl2 = ceil(Int, ranges[2] * (2^level ))   + 1
            rl3 = floor(Int, ranges[3] * (2^level ))  + 1
            rl4 = ceil(Int, ranges[4] * (2^level ))   + 1
            rl5 = floor(Int, ranges[5] * (2^level)) + 1
            rl6 = ceil(Int, ranges[6] * (2^level))  + 1
            #println("Rebin data on the used maximum used grid")
            #println("xrange: ",rl1, " ", rl2)
            #println("yrange: ",rl3, " ", rl4)
            #println()

            if level == simlmax
                alt_shift = 1
            end

            if direction == :z
                # range on maximum used grid
                new_level_range1 = range(rl1, stop=rl2, length=(rl2-rl1) +alt_shift)
                new_level_range2 = range(rl3, stop=rl4, length=(rl4-rl3) +alt_shift)
                #println(newrange1)
                #println(newrange2)

            elseif direction == :y
                # range on maximum used grid
                new_level_range1 = range(rl1, stop=rl2, length=(rl2-rl1)+alt_shift)
                new_level_range2 = range(rl5, stop=rl6, length=(rl6-rl5)+alt_shift)
                #println(newrange1)
                #println(newrange2)s

            elseif direction == :x
                # range on maximum used grid
                new_level_range1 = range(rl3, stop=rl4, length=(rl4-rl3)+alt_shift)
                new_level_range2 = range(rl5, stop=rl6, length=(rl6-rl5)+alt_shift)
                #println(newrange1)
                #println(newrange2)

            end

            # range on maximum used level grid
            length_level1=length( new_level_range1 )
            length_level2=length( new_level_range2 )

                # needed for mass weighting and/or sd,rho, ekin projections
                if in(:rho, selected_vars) || in(:ρ, selected_vars) || in(:density, selected_vars) ||
                    in(:sd, selected_vars) || in(:Σ, selected_vars) || in(:surfacedensity, selected_vars) ||
                    in(:ekin, selected_vars) ||
                    weighting == :mass

                    if length(mask) == 1
                        global h = fit(Histogram, ( select(level_data, x_coord), select(level_data, y_coord) ),
                                            weights( select(level_data, :rho) ),
                                            closed=closed,
                                            (new_level_range1, new_level_range2)  )
                    else
                        global h = fit(Histogram, ( select(level_data, x_coord), select(level_data, y_coord) ),
                                            weights( select(level_data, :rho) .* select(level_data, :mask) ),
                                            closed=closed,
                                            (new_level_range1, new_level_range2)  )
                    end

                end

                counter = 0
                for ivar in selected_vars
                    counter = counter + 1

                    if !in(ivar, rσanglecheck)  # exclude velocity dispersion symbols and radius/angle maps

                        # non derived variables, density weighted (per level)
                        if !in(ivar, density_names) && !in(ivar, sd_names)
                        #if ivar != :sd  &&  ivar!=:rho
                             #&& ivar != :Σ && ivar != :surfacedensity  #&& ivar != :ρ && ivar != :density
                            if weighting == :mass
                                if length(mask) == 1
                                    #println(ivar, " ", counter)
                                    h_var = fit(Histogram, ( select(level_data, x_coord), select(level_data, y_coord) ),
                                                    weights( getvar(dataobject, ivar, filtered_db=level_data, center=data_centerm, direction=direction) .* select(level_data, :rho)  ),
                                                    closed=closed,
                                                    (new_level_range1, new_level_range2) )
                                else
                                    h_var = fit(Histogram, ( select(level_data, x_coord), select(level_data, y_coord) ),
                                                    weights( getvar(dataobject, ivar, filtered_db=level_data, center=data_centerm, direction=direction) .* select(level_data, :rho) .* select(level_data, :mask)  ),
                                                    closed=closed,
                                                    (new_level_range1, new_level_range2) )
                                end
                            elseif weighting == :volume
                                if length(mask) == 1
                                    h_var = fit(Histogram, ( select(level_data, x_coord), select(level_data, y_coord) ),
                                                    weights( getvar(dataobject, ivar, filtered_db=level_data, center=data_centerm, direction=direction)  ),
                                                    closed=closed,
                                                    (new_level_range1, new_level_range2) )
                                else
                                    h_var = fit(Histogram, ( select(level_data, x_coord), select(level_data, y_coord) ),
                                                    weights( getvar(dataobject, ivar, filtered_db=level_data, center=data_centerm, direction=direction)  .* select(level_data, :mask) ),
                                                    closed=closed,
                                                    (new_level_range1, new_level_range2) )
                                end
                            end

                        end

                        # scale to current levels
                        if in(ivar, sd_names) #ivar == :sd #|| ivar == :Σ || ivar == :surfacedensity
                            map_buffer = h.weights .* (boxlen / 2^level)

                        elseif in(ivar, density_names) #ivar == :rho #|| ivar == :ρ || ivar == :density
                            map_buffer = h.weights .* (boxlen / 2^level)

                        elseif !in(ivar, density_names) && !in(ivar, sd_names) && weighting == :volume
                        #elseif ivar !== :sd && ivar != :rho && weighting == false
                            map_buffer = h_var.weights .* (boxlen / 2^level)

                        else # for any kind of mass weighted projection
                            map_buffer = h_var.weights .* (boxlen / 2^level)^3
                            #if first_time == 1
                            map_buffer_weight = h.weights .* (boxlen / 2^level)^3
                            #end

                        end

                        s = size(map_buffer)
                        #println( s )
                        # remap on selected gridsize
                        if level != simlmax
                            ratio_level = 2^simlmax / 2^level
                            scaled_length1 = Int(ratio_level * (length_level1-alt_shift) )
                            scaled_length2 = Int(ratio_level * (length_level2-alt_shift) )

                            r1_diff = f_min(new_level_range1[1], ratio_level)-newrange1[1]
                            r2_diff = f_max(new_level_range1[end]-1., ratio_level)-newrange1[end]+1.

                            r3_diff = f_min(new_level_range2[1], ratio_level)-newrange2[1]
                            r4_diff = f_max(new_level_range2[end]-1., ratio_level)-newrange2[end]+1.

                            # map on lmax grid
                            remapped = [map_buffer[floor(Int,x),floor(Int,y)]  for x in range(1, stop=s[1], length=scaled_length1) , y in range(1, stop=s[2], length=scaled_length2)]

                            if !in(ivar, density_names) && !in(ivar, sd_names) && weighting == :mass
                            #if ivar != :sd && ivar != :rho && weighting == true

                                if first_time_level[counter] == 1
                                    map[:, :, counter] .+= remapped[Int(1-r1_diff):Int(end-r2_diff), Int(1-r3_diff):Int(end-r4_diff)] ./ (2^simlmax / 2^level  )^2

                                    remap_weight = [map_buffer_weight[floor(Int,x),floor(Int,y)]  for x in range(1, stop=s[1], length=scaled_length1), y in range(1, stop=s[2], length=scaled_length2)]
                                    map_weight[:, :, counter] .+= remap_weight[Int(1-r1_diff):Int(end-r2_diff), Int(1-r3_diff):Int(end-r4_diff)] ./ (2^simlmax / 2^level  )^2  # and correct for remapping
                                    #println(counter, " ", level, " ", sum(map_weight[:, :, counter]))
                                    first_time_level[counter] = 0
                                end
                            else
                                map[:, :, counter] .+= remapped[Int(1-r1_diff):Int(end-r2_diff), Int(1-r3_diff):Int(end-r4_diff)]
                            end
                        elseif level == simlmax
                            map[:,:,counter] .+= map_buffer

                            if !in(ivar, density_names) && !in(ivar, sd_names) && weighting == :mass && first_time_level[counter] == 1
                            #if ivar != :sd && ivar != :rho && weighting == true && first_time_level[counter] == 1
                                map_weight[:, :, counter] .+= map_buffer_weight
                                #println(counter, " ", level, " ", sum(map_weight[:, :, counter]))
                                first_time_level[counter] = 0
                            end
                        end

                    end

                end

        end
    end # notonly_ranglecheck_vars

    # calc density maps & reverse weighting
    counter = 0
    for ivar in selected_vars
        counter = counter + 1
        if !in(ivar, rσanglecheck) # exclude velocity dispersion symbols and radius/angle maps
            if !in(ivar, density_names) && !in(ivar, sd_names) && weighting == :mass

            #if ivar != :sd && ivar != :Σ && ivar != :surfacedensity && ivar != :rho && ivar != :ρ && ivar != :density && weighting == true
                #map[:,:,counter] = map[:,:,counter] ./ map_weight
                selected_unit, unit_name= getunit(dataobject, ivar, selected_vars, units, uname=true)
                #println(counter, " ", size(map), " ", ivar, " ", selected_unit, " ", unit_name)
                #println(ivar, " ", counter, " ", unit_name)
                maps[Symbol(ivar)] = map[:,:, counter] ./ map_weight[:, :, counter] .* selected_unit
                #println(sum(map[:,:, counter]), " ", sum(map_weight[:, :, counter]), " ", sum(map[:,:, counter] ./ map_weight[:, :, counter]))
                #println()
                maps_unit[Symbol( ivar )] = unit_name

            else
                selected_unit, unit_name= getunit(dataobject, ivar, selected_vars, units, uname=true)
                #println(counter, " ", size(map), " ", ivar, " ", selected_unit, " ", unit_name)
                maps[Symbol(ivar)] = map[:,:, counter]  .* selected_unit
                #println(sum(map[:,:, counter]), " ", sum(map_weight[:, :, counter]), " ", sum(map[:,:, counter] ./ map_weight[:, :, counter]))
                #println()
                maps_unit[Symbol( string(ivar)  )] = unit_name
            end

        end
    end

    # create radius map
    for ivar in selected_vars
        if in(ivar, rcheck)
            selected_unit, unit_name= getunit(dataobject, ivar, selected_vars, units, uname=true)
            map_R = zeros(Float64, length1, length2 );
            for i = 1:(length1)
                for j = 1:(length2)
                    x = i * dataobject.boxlen / 2^simlmax
                    y = j * dataobject.boxlen / 2^simlmax
                    radius = sqrt((x-length1_center)^2 + (y-length2_center)^2)
                    map_R[i,j] = radius * selected_unit
                end
            end

            maps[Symbol(ivar)] = map_R
            maps_unit[Symbol( string(ivar)  )] = unit_name
        end
    end

    # create ϕ-angle map
    for ivar in selected_vars
        if in(ivar, anglecheck)
            map_ϕ = zeros(Float64, length1, length2 );
            for i = 1:(length1)
                for j = 1:(length2)
                    x = i * dataobject.boxlen / 2^simlmax - length1_center
                    y = j * dataobject.boxlen / 2^simlmax - length2_center
                    if x > 0. && y >= 0.
                        map_ϕ[i,j] = atan(y / x)
                    elseif x > 0. && y < 0.
                        map_ϕ[i,j] = atan(y / x) + 2. * pi
                    elseif x < 0.
                        map_ϕ[i,j] = atan(y / x) + pi
                    elseif x==0 && y > 0
                        map_ϕ[i,j] = pi/2.
                    elseif x==0 && y < 0
                        map_ϕ[i,j] = 3. * pi/2.
                    end
                end
            end

            maps[Symbol(ivar)] = map_ϕ
            maps_unit[Symbol( string(ivar)  )] = :radian
        end
    end

    # create velocity dispersion maps, after all other maps are created
    counter = 0
    for ivar in selected_vars
        counter = counter + 1

        if in(ivar, σcheck)
            #for iσ in σcheck
                selected_unit, unit_name= getunit(dataobject, ivar, selected_vars, units, uname=true)

                    selected_v = σ_to_v[ivar]
                    iv  = maps[selected_v[1]]
                    iv_unit = maps_unit[Symbol( string(selected_v[1])  )]
                    iv2 = maps[selected_v[2]]
                    iv2_unit = maps_unit[Symbol( string(selected_v[2])  )]
                    if iv_unit == iv2_unit
                        if iv_unit == unit_name
                            maps[Symbol(ivar)] = sqrt.( iv2 .- iv .^2 )
                        elseif iv_unit == :standard
                            maps[Symbol(ivar)] = sqrt.( iv2 .- iv .^2 )  .* selected_unit
                        elseif iv_unit == :km_s
                            maps[Symbol(ivar)] = sqrt.( iv2 .- iv .^2 )  ./ dataobject.info.scale.km_s
                        end
                    elseif iv_unit != iv2_unit
                        if iv_unit == :km_s && unit_name == :standard
                            iv = iv ./ dataobject.info.scale.km_s
                        elseif iv_unit == :standard && unit_name == :km_s
                            iv = iv .* dataobject.info.scale.km_s
                        end
                        if iv2_unit == :km_s && unit_name == :standard
                            iv2 = iv2 ./ dataobject.info.scale.km_s.^2
                        elseif iv2_unit == :standard && unit_name == :km_s
                            iv2 = iv2 .* dataobject.info.scale.km_s.^2
                        end

                        # overwrite NaN due to radius = 0
                        #iv2 = iv2[isnan.(iv2)] .= 0
                        #iv  = iv[isnan.(iv)] .= 0

                        maps[Symbol(ivar)] = sqrt.( iv2 .- iv .^2 )
                    end

                    maps_unit[Symbol( string(ivar)  )] = unit_name
                #end
            #end
        end
    end

    if mera_mask_inserted # delete column :mask
        dataobject.data = select(dataobject.data, Not(:mask))
    end

    finalmaps = HydroMapsType(maps, maps_unit, maps_lmax, maps_mode, lmax_projected, lmin, simlmax, ranges, extent, extent_center, ratio, boxlen, dataobject.smallr, dataobject.smallc, dataobject.scale, dataobject.info)

    if simlmax > lmax
            return remap(finalmaps, lmax, weighting=weighting)
    else
        return finalmaps
    end

end


text:
#=
    InterpolatedRejectionSampling
    Copyright © 2019 Mark Wells <mwellsa@gmail.com>

    Distributed under terms of the AGPL-3.0 license.
=#

module InterpolatedRejectionSampling

using Base.Iterators
using Interpolations
using Interpolations: Extrapolation
using StatsBase: sample, Weights

export irsample, irsample!

@inline function midpoints(x::AbstractVector{Float64})
    length(x) == 1 && return x
    retval = Vector{Float64}(undef, length(x)-1)
    @fastmath @inbounds @simd for i in eachindex(retval)
        retval[i] = (x[i  ] +
                     x[i+1])/2
    end
    return retval
end

@inline function midpoints(x::AbstractRange{Float64})
    length(x) == 1 && return x
    Δx = 0.5*step(x)
    return range(first(x)+Δx, stop=last(x)-Δx, length=length(x)-1)
end

@inline get_interp(interp::AbstractExtrapolation{Float64,N,ITPT,IT}, val::NTuple{N,Float64}) where {N,ITPT,IT} = interp(val...)

#@inline get_knots(interp::Extrapolation{Float64,1,ITPT,BSpline{Linear},ET}) where {ITPT,ET} = first(interp.itp.ranges)
#@inline get_knots(interp::Extrapolation{Float64,1,ITPT,Gridded{Linear},ET}) where {ITPT,ET} = first(interp.itp.knots)

@inline get_knots(interp::Extrapolation{Float64,N,ITPT,BSpline{Linear},ET}) where {N,ITPT,ET} = interp.itp.ranges
@inline get_knots(interp::Extrapolation{Float64,N,ITPT,Gridded{Linear},ET}) where {N,ITPT,ET} = interp.itp.knots

@inline get_coefs(interp::Extrapolation{Float64,N,ITPT,BSpline{Linear},ET}) where {N,ITPT,ET} = interp.itp.itp.coefs
@inline get_coefs(interp::Extrapolation{Float64,N,ITPT,Gridded{Linear},ET}) where {N,ITPT,ET} = interp.itp.coefs

@inline get_Δ(x::AbstractVector) = length(x) > 1 ? diff(x) : 1

@inline integrate(interp::AbstractExtrapolation{Float64,N,ITPT,IT}, (x, Δx)::NTuple{2,NTuple{N,Float64}}
                  ) where {N,ITPT,IT} = prod(Δx)*get_interp(interp, x)

function integrate(knots::NTuple{N,AbstractVector{Float64}},
                   interp::AbstractExtrapolation{Float64,N,ITPT,IT}
                  ) where {N,ITPT,IT}
    midknots = map(midpoints, knots)
    Δknots = map(get_Δ, knots)
    return sum(x -> integrate(interp, x), zip(product(midknots...), product(Δknots...)))
end

@inline integrate(interp::AbstractExtrapolation{Float64,N,ITPT,IT}
                 ) where {N,ITPT,IT} = integrate(get_knots(interp), interp)

@inline is_normalized(interp::AbstractExtrapolation{Float64,N,ITPT,IT}
                     ) where {N,ITPT,IT} = isapprox(integrate(interp), one(Float64))

function normalize_interp(interp::AbstractExtrapolation{Float64,N,ITPT,IT}) where {N,ITPT,IT}
    knots = get_knots(interp)
    coefs = get_coefs(interp)
    A = integrate(interp)
    coefs ./= A
    return LinearInterpolation(knots, coefs)
end

@inline sliced_knots(k::AbstractVector{Float64}, s::T
                    ) where T<:Union{Missing,Float64} = ismissing(s) ? k : [s]

@inline sliced_knots(knots::NTuple{N,AbstractVector}, slice::AbstractVector{Union{Missing,Float64}}
                    ) where N = ntuple(i -> sliced_knots(knots[i], slice[i]), Val(N))

struct Cells{Float64,N,ITPT,IT,ET}
    knots::NTuple{N,Vector}
    pmass::Array{Float64,N}
    cinds::CartesianIndices{N,NTuple{N,Base.OneTo{Int}}}
    interp::Extrapolation{Float64,N,ITPT,IT,ET}

    function Cells(interp::Extrapolation{Float64,N,ITPT,IT,ET},
                   knots::NTuple{N,AbstractVector} = get_knots(interp)
                  ) where {T<:Union{Missing,Float64},N,ITPT,IT,ET}

        ksz = map(length, knots)
        midpnt = map(midpoints, knots)
        msz = map(length, midpnt)

        pmass = Array{Float64,N}(undef,msz)
        for (i,k) in enumerate(product(midpnt...))
            pmass[i] = get_interp(interp, k)
        end
        pmass ./= sum(pmass)

        new{Float64,N,ITPT,IT,ET}(knots, pmass, CartesianIndices(msz), interp)
    end

    function Cells(interp::Extrapolation{Float64,N,ITPT,IT,ET},
                   slice::AbstractVector{Union{Missing,Float64}}
                  ) where {N,ITPT,IT,ET}
        knots = get_knots(interp)
        sknots = sliced_knots(knots, slice)
        return Cells(interp, sknots)
    end
end

@inline get_interp(C::Cells{Float64,N,ITPT,IT,ET}, val::NTuple{N,Float64}
                  ) where {N,ITPT,IT,ET} = get_interp(C.interp, val)

import StatsBase.sample
@inline sample(C::Cells) = sample(CartesianIndices(C.cinds), Weights(vec(C.pmass)))

@inline _get_xmin(x::AbstractVector{Float64}, i::Int)::Float64 = x[i]
@inline _get_span(x::AbstractVector{Float64}, i::Int)::Float64 = length(x) == 1 ? zero(Float64) : x[i+1] - x[i]

struct Support{Float64,N}
    xmin::NTuple{N,Float64}
    span::NTuple{N,Float64}

    function Support(C::Cells{Float64,N}, cind::CartesianIndex{N}) where {N}
        xmin = ntuple(i -> _get_xmin(C.knots[i], cind.I[i])::Float64, Val(N))
        span = ntuple(i -> _get_span(C.knots[i], cind.I[i])::Float64, Val(N))
        new{Float64,N}(xmin,span)
    end
end

@inline Base.getindex(S::Support{Float64,N}, i::Int) where N = (S.xmin[i], S.span[i])

Base.iterate(S::Support{Float64,N}, state::Int=1) where N = state > N ? nothing : (S[state], state+1)

@inline propose_sample(S::Support{Float64,N}
                      ) where N = ntuple(i -> iszero(last(S[i])) ?
                                         first(S[i]) :
                                         first(S[i]) + rand()*last(S[i]),
                                         Val(N)
                                        )

@inline get_extrema(S::Support{Float64,N}
                   ) where N = product(ntuple(i -> (first(S[i]), first(S[i]) + last(S[i])),
                                              Val(N)
                                             )...
                                      )

@inline maxmapreduce(f,a) = mapreduce(f, max, a)

struct Envelope{Float64,N,ITPT,IT,ET}
    support::Support{Float64,N}
    maxvalue::Float64
    interp::Extrapolation{Float64,N,ITPT,IT,ET}
    function Envelope(C::Cells{Float64,N,ITPT,IT,ET},
                      cind::CartesianIndex{N}
                     ) where {N,ITPT,IT,ET}
        support = Support(C, cind)
        spnts = get_extrema(support)
        maxvalue = maxmapreduce(x -> get_interp(C,x), spnts)
        new{Float64,N,ITPT,IT,ET}(support, maxvalue, C.interp)
    end
end

@inline get_interp(E::Envelope{Float64,N,ITPT,IT,ET},
                   val::NTuple{N,Float64}
                  ) where {N,ITPT,IT,ET} = get_interp(E.interp, val)

@inline propose_sample(E::Envelope) = propose_sample(E.support)

function rsample(E::Envelope{Float64,N,ITPT,IT,ET}) where {N,ITPT,IT,ET}
    while true
        samp = propose_sample(E)
        if rand()*E.maxvalue ≤ get_interp(E, samp)
            return samp
        end
    end
    error("unable to draw a sample after $maxruns runs")
end

function rsample(C::Cells)
    cind = sample(C)
    E = Envelope(C,cind)
    return rsample(E)
end

function rsample(interp::Extrapolation)
    C = Cells(interp)
    return rsample(C)
end

function irsample(knots::NTuple{N,AbstractVector{Float64}},
                  probs::AbstractArray{Float64,N}
                 ) where N
    interp = LinearInterpolation(knots, probs)
    return rsample(interp)
end

function irsample(knots::NTuple{N,AbstractVector{Float64}},
                  probs::AbstractArray{Float64,N},
                  n::Int
                 ) where N
    interp = LinearInterpolation(knots, probs)
    retval = Matrix{Float64}(undef, N, n)
    for s in eachcol(retval)
        s .= rsample(interp)
    end
    return retval
end

function rsample(interp::Extrapolation, slice::AbstractVector{Union{Missing,Float64}})
    C = Cells(interp, slice)
    return rsample(C)
end

function irsample!(slices::AbstractMatrix{Union{Missing,Float64}},
                   knots::NTuple{N,AbstractVector{Float64}},
                   probs::AbstractArray{Float64,N}
                  ) where N
    interp = LinearInterpolation(knots, probs)
    for s in eachcol(slices)
        s .= rsample(interp, s)
    end
end

end

ext: jl


text:
# Model structure
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#

const IntTuple = Union{Tuple{Integer,Integer}, Tuple{Integer,Integer,Integer},Array{Int64,1},Array{Int32,1}}
const RealTuple = Union{Tuple{Real,Real}, Tuple{Real,Real,Real},Array{Float64,1},Array{Float32,1}}

export Model

# Object for velocity/slowness models
mutable struct Model
    n::IntTuple
    d::RealTuple
    o::RealTuple
    nb::Integer # number of absorbing boundaries points on each side
    m   # slowness squared
    rho # density in g / m^3 (1 for water)
end

"""
    Model
        n::IntTuple
        d::RealTuple
        o::RealTuple
        nb::Integer
        m::Array
        rho::Array

Model structure for seismic velocity models.

`n`: number of gridpoints in (x,y,z) for 3D or (x,z) for 2D

`d`: grid spacing in (x,y,z) or (x,z) (in meters)

`o`: origin of coordinate system in (x,y,z) or (x,z) (in meters)

`nb`: number of absorbing boundary points in each direction

`m`: velocity model in slowness squared (s^2/km^2)

`rho`: density (g / m^3)

Constructor
===========

The parameters `n`, `d`, `o` and `m` are mandatory, whith `nb` and `rho` being optional input arguments.

    Model(n, d, o, m; nb=40, rho=ones(n))

"""
function Model(n::IntTuple, d::RealTuple, o::RealTuple, m; rho=[], nb=40)
    isempty(rho) && (rho = 1)
    return Model(n,d,o,nb,m,rho)
end

function Model(n::IntTuple, d::RealTuple, o::RealTuple, m, rho; nb=40)
    isempty(rho) && (rho = 1)
    return Model(n,d,o,nb,m,rho)
end


text:
# Test Fishing Problem - Switched Nonlinear System
function nonlinsystest(solver=IpoptSolver(tol=1e-04, max_iter = 25); objtol = 1e-3, primaltol = 1e-3)
  println("Testing Switched Nonlinear Systems Optimization ", string(typeof(solver)))

  # Define Time Interval
  t0 = 0.0; tf = 12.0

  uvec = [repmat([0.0; 1.0], 4, 1); 0.0]'  # Input Vector

  # Cost Funcction Matrix
  C = [1.0 0.0 -1.0 0.0;
       0.0 1.0 0.0  -1.0]
  Q = C'*C

  # Define Initial State
  x0 = [0.5; 0.7; 1; 1]

  ### Define System Dynamics
  function nldyn(x::Array{Float64,1}, u::Array{Float64,1})
    n = length(x)
    f = zeros(n)
    f[1] = x[1] - x[1]*x[2] - 0.4*x[1]*u[1]
    f[2] = -x[2] + x[1]*x[2] - 0.2*x[2]*u[1]
    return f

  end

  function nldyn_deriv(x::Array{Float64,1}, u::Array{Float64,1})
    df = [1.0-x[2]-0.4*u[1]       -x[1]                   0    0;
          x[2]                     -1+x[1]-0.2*u[1]       0    0;
          0                        0                      0    0;
          0                        0                      0    0]
  end

  m = stoproblem(x0, nldyn, nldyn_deriv, uvec, ngrid = 300, t0=t0, tf=tf, Q=Q, solver=solver)

  solve!(m)
  # @test getstat(m) == :Optimal

  # Test Optimal Solution
  tauopt = gettau(m)
  @test_approx_eq_eps tauopt[1] 2.4434718158206916 primaltol
  @test_approx_eq_eps tauopt[2] 4.122362522221089 primaltol
  @test_approx_eq_eps tauopt[3] 4.433072844726073 primaltol
  @test_approx_eq_eps tauopt[4] 4.68252002334405  primaltol
  @test_approx_eq_eps tauopt[5] 5.201667827750571 primaltol
  @test_approx_eq_eps tauopt[6] 5.369908894975762 primaltol
  @test_approx_eq_eps tauopt[7] 6.376845190917198 primaltol
  @test_approx_eq_eps tauopt[8] 6.47160340206027 primaltol

  # Test Optimal Value
  @test_approx_eq_eps getobjval(m) 1.3454355602054182 objtol

  println("Passed")
end


text:
module TestIM

using IterationManagers
using Test
using LinearAlgebra: norm

include("newton.jl")

end  # module


text:
using Test, LinearAlgebra, Random
using DynamicPolynomials, HomotopyContinuation, StaticArrays
import TreeViews, ProjectiveVectors, PolynomialTestSystems

import PolynomialTestSystems: cyclic, katsura, equations, ipp2, heart, griewank_osborne
const HC = HomotopyContinuation

function test_treeviews(x)
    @test TreeViews.hastreeview(x)
    @test_nowarn TreeViews.treelabel(devnull, x, MIME"application/prs.juno.inline"())
    for i=1:TreeViews.numberofnodes(x)
        @test_nowarn TreeViews.nodelabel(devnull, x, i, MIME"application/prs.juno.inline"())
        @test_nowarn TreeViews.treenode(x, i)
    end
end

function test_show_juno(x)
    @test show(stdout, MIME("application/prs.juno.inline"), x) === x
end

# We order the tests such that isolated things are tested first
@testset "HomotopyContinuation" begin
    include("utilities_test.jl")
    include("affine_patches_test.jl")
    include("problem_test.jl")
    include("systems_test.jl")
    include("homotopies_test.jl")
    include("predictors_test.jl")
    include("correctors_test.jl")
    include("core_tracker_test.jl")
    include("path_tracker_test.jl")
    include("solve_test.jl")
    include("polyhedral_test.jl")
    include("result_test.jl")
    include("root_count_test.jl")
    include("monodromy_test.jl")
    include("nextjournal.jl")
end


text:
#######
# Demo setup

cool(x) = x + 1
cool(x, y) = x + y + 1

# a rule we define so we can test rules
dummy_identity(x) = x
@scalar_rule(dummy_identity(x), One())

#######

_second(t) = Base.tuple_type_head(Base.tuple_type_tail(t))

@testset "frule and rrule" begin
    @test frule(cool, 1) === nothing
    @test frule(cool, 1; iscool=true) === nothing
    @test rrule(cool, 1) === nothing
    @test rrule(cool, 1; iscool=true) === nothing

    # add some methods:
    ChainRulesCore.@scalar_rule(Main.cool(x), one(x))
    @test hasmethod(rrule, Tuple{typeof(cool),Number})
    ChainRulesCore.@scalar_rule(Main.cool(x::String), "wow such dfdx")
    @test hasmethod(rrule, Tuple{typeof(cool),String})
    # Ensure those are the *only* methods that have been defined
    cool_methods = Set(m.sig for m in methods(rrule) if _second(m.sig) == typeof(cool))
    only_methods = Set([Tuple{typeof(rrule),typeof(cool),Number},
                        Tuple{typeof(rrule),typeof(cool),String}])
    @test cool_methods == only_methods

    frx, cool_pushforward = frule(cool, 1)
    @test frx == 2
    @test cool_pushforward(NamedTuple(), 1) == 1
    rrx, cool_pullback = rrule(cool, 1)
    self, rr1 = cool_pullback(1)
    @test self == NO_FIELDS
    @test rrx == 2
    @test rr1 == 1
end

@testset "Basic Wirtinger scalar_rule" begin
    myabs2(x) = abs2(x)
    @scalar_rule(myabs2(x), Wirtinger(x', x))

    @testset "real input" begin
        # even though our rule was define in terms of Wirtinger,
        # pushforward result will be real as real (even if seed is Compex)

        x = rand(Float64)
        f, myabs2_pushforward = frule(myabs2, x)
        @test f === x^2

        Δ = One()
        df = @inferred myabs2_pushforward(NamedTuple(), Δ)
        @test df === x + x

        Δ = rand(Complex{Int64})
        df = @inferred myabs2_pushforward(NamedTuple(), Δ)
        @test df === Δ * (x + x)
    end

    @testset "complex input" begin
        z = rand(Complex{Float64})
        f, myabs2_pushforward = frule(myabs2, z)
        @test f === abs2(z)

        df = @inferred myabs2_pushforward(NamedTuple(), One())
        @test df === Wirtinger(z', z)

        Δ = rand(Complex{Int64})
        df = @inferred myabs2_pushforward(NamedTuple(), Δ)
        @test df === Wirtinger(Δ * z', Δ * z)
    end
end

@testset "Advanced Wirtinger @scalar_rule: abs_to_pow" begin
    # This is based on SimeonSchaub excellent example:
    # https://gist.github.com/simeonschaub/a6dfcd71336d863b3777093b3b8d9c97

    # This is much more complex than the previous case
    # as it has many different types
    # depending on input, and the output types do not always agree

    abs_to_pow(x, p) = abs(x)^p
    @scalar_rule(
        abs_to_pow(x::Real, p),
        (
            p == 0 ? Zero() : p * abs_to_pow(x, p-1) * sign(x),
            Ω * log(abs(x))
        )
    )

    @scalar_rule(
        abs_to_pow(x::Complex, p),
        @setup(u = abs(x)),
        (
            p == 0 ? Zero() : p * u^(p-1) * Wirtinger(x' / 2u, x / 2u),
            Ω * log(abs(x))
        )
    )

    f = abs_to_pow
    @testset "f($x, $p)" for (x, p) in Iterators.product(
        (2, 3.4, -2.1, -10+0im, 2.3-2im),
        (0, 1, 2, 4.3, -2.1, 1+.2im)
    )
        expected_type_df_dx =
            if iszero(p)
                Zero
            elseif typeof(x) <: Complex
                Wirtinger
            elseif typeof(p) <: Complex
                Complex
            else
                Real
            end

        expected_type_df_dp =
            if typeof(p) <: Real
                Real
            else
                Complex
            end

        res = frule(f, x, p)
        @test res !== nothing  # Check the rule was defined
        fx, f_pushforward = res
        df(Δx, Δp) = f_pushforward(NamedTuple(), Δx, Δp)

        df_dx::Thunk = df(One(), Zero())
        df_dp::Thunk = df(Zero(), One())
        @test fx == f(x, p)  # Check we still get the normal value, right
        @test df_dx() isa expected_type_df_dx
        @test df_dp() isa expected_type_df_dp

        res = rrule(f, x, p)
        @test res !== nothing  # Check the rule was defined
        fx, f_pullback = res
        dself, df_dx, df_dp = f_pullback(One())
        @test fx == f(x, p)  # Check we still get the normal value, right
        @test dself == NO_FIELDS
        @test df_dx() isa expected_type_df_dx
        @test df_dp() isa expected_type_df_dp
    end
end


text:
"""
    function get_Pg_thermal_nonzeroslope(args...; kargs...)

Identify generators which cost function has a non zero slope, locate their buses and returns
a matrix with their Name, Power output and Bus number (sorted by Bus). Skips Slack
generators, generators which Active Power is binding and Generators with constant price.
Note: For a quadratic cost function of the form: αPg²+ βPg + γ, a non zero slope generator
will have the parameter α ≠ 0.

# Arguments
- `sys::System`:                            Power system in p.u. (from PowerSystems.jl)
- `res::PowerSimulations.
    OperationsProblemResults`:              Results of the solved OPF for the system
                                             (from PowerSimulations.jl)

# Keywords
- `dual_gen_tol::Float64 = 1e-1`:           Tolerance to identify any binding generators
"""
function get_Pg_thermal_nonzeroslope(
    sys::System,
    res::PowerSimulations.OperationsProblemResults;
    dual_gen_tol::Float64
    )
    # Define elements
    gen_loc = 0
    gen_nonzero = 0 #number of generators with non zero slope (counter starter)
    gens_thermal_nonzeroslope_Pg = Array{Float64}(undef, 0)
    gens_thermal_nonzeroslope_names = Array{String}(undef, 0)
    gens_thermal_nonzeroslope_busnumbers = Array{Int64}(undef, 0)

    # Get all thermal Generators
    gens_thermal = get_components(ThermalStandard, sys)
    all_PGenThermal = get_variables(res)[:P__ThermalStandard] #Optimised PGen

    # Get dual variables from the OPF results
    dualPgub = get_duals(res)[:P_ub__ThermalStandard__RangeConstraint]
    dualPglb = get_duals(res)[:P_lb__ThermalStandard__RangeConstraint]

    # Identify Generators and save their information
    for gen_thermal in gens_thermal
        gen_loc = gen_loc +1
        if get_bustype(gen_thermal.bus) ≠ BusTypes.REF  #skips slack
            bindub = !isapprox(dualPgub[1,gen_loc], 0.0; atol = dual_gen_tol) #is binding?
            bindlb = !isapprox(dualPglb[1,gen_loc], 0.0; atol = dual_gen_tol) #is binding?
            if !bindub && !bindlb #skips binding Pg gens Dual = 0
                α, β = get_cost(get_variable(get_operation_cost(gen_thermal)))
                if α ≠ 0 # Only gens with α(Pg^2)
                    gen_nonzero = gen_nonzero +1
                    gen_thermal_nonzeroslope_name = gen_thermal.name
                    resize!(gens_thermal_nonzeroslope_Pg, gen_nonzero)
                    resize!(gens_thermal_nonzeroslope_names, gen_nonzero)
                    resize!(gens_thermal_nonzeroslope_busnumbers, gen_nonzero)
                    gens_thermal_nonzeroslope_Pg[gen_nonzero] =
                        all_PGenThermal[1, gen_thermal_nonzeroslope_name]
                    gens_thermal_nonzeroslope_names[gen_nonzero] =
                        gen_thermal_nonzeroslope_name
                    gens_thermal_nonzeroslope_busnumbers[gen_nonzero] =
                        gen_thermal.bus.number
                end
            end
        end
    end

    # Matrix of generators with their information sorted
    gens_thermal_nonzeroslope_Matrix = Array{Any}(undef,(gen_nonzero, 3))
    gens_thermal_nonzeroslope_Matrix[:, 1] = gens_thermal_nonzeroslope_names
    gens_thermal_nonzeroslope_Matrix[:, 2] = gens_thermal_nonzeroslope_Pg
    gens_thermal_nonzeroslope_Matrix[:, 3] = gens_thermal_nonzeroslope_busnumbers
    #sort by bus to match with PTDF structure
    gens_thermal_nonzeroslope_Matrix =
       gens_thermal_nonzeroslope_Matrix[sortperm(gens_thermal_nonzeroslope_Matrix[:, 3]), :]

    return gens_thermal_nonzeroslope_Matrix
end

max_stars_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_issues_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_forks_repo_name: invenia/MaximumGeneratorProfitTCRDD


text:
"""
    DGMesh{T} <: Mesh{T}

A structure representing a discretisation scheme to be used for a DiscretisedFluidQueue.

# Arguments:
- `nodes::AbstractArray{Float64, 1}`: The edges of the cells.
- `n_bases::Int`: The number of basis functions used to represent the solution on each cell
"""
struct DGMesh{T} <: Mesh{T}
    nodes::T
    n_bases::Int
    function DGMesh{T}(nodes::T,n_bases::Int) where T
        (n_bases<=0)&&throw(DomainError("n_bases must be positive"))
        (nodes[1]!=0.0)&&throw(DomainError("first node must be 0"))
        return new{T}(nodes,n_bases)
    end
end
# Convenience constructors
DGMesh(nodes::T,n_bases::Int) where T = DGMesh{T}(nodes::T,n_bases::Int)
DGMesh() = DGMesh([0.0],0)

"""

    n_bases_per_cell(mesh::DGMesh)

Number of bases in a cell
"""
n_bases_per_cell(mesh::DGMesh) = mesh.n_bases

"""

    cell_nodes(mesh::DGMesh)

The positions of the GLJ nodes within each cell of a mesh.
"""
function cell_nodes(mesh::DGMesh)
    cellnodes = zeros(Float64, n_bases_per_cell(mesh), n_intervals(mesh))
    if n_bases_per_cell(mesh) > 1
        z = Jacobi.zglj(n_bases_per_cell(mesh), 0, 0) # the LGL nodes
    elseif n_bases_per_cell(mesh) == 1
        z = 0.0
    end
    for i = 1:n_intervals(mesh)
        # Map the LGL nodes on [-1,1] to each cell
        cellnodes[:, i] .= (mesh.nodes[i+1] + mesh.nodes[i]) / 2 .+ (mesh.nodes[i+1] - mesh.nodes[i]) / 2 * z
    end
    cellnodes[1,:] .+= sqrt(eps()) # move the cell edges because funny things can happen at boundaries
    if n_bases_per_cell(mesh)>1
        cellnodes[end,:] .-= sqrt(eps())
    end
    return cellnodes
end

"""
    local_dg_operators(mesh::DGMesh; v::Bool = false)

Construct the block matrices and vectors which define the discretised
generator (i.e. the `blocks`, and `boundary_flux` in `LazyGenerator`).
"""
function local_dg_operators(
    mesh::DGMesh;
    v::Bool = false,
)
    ## Construct local blocks
    V = vandermonde(n_bases_per_cell(mesh))

    Dw = (
        DwInv = LinearAlgebra.diagm(0 => 1.0 ./ V.w),
        Dw = LinearAlgebra.diagm(0 => V.w),
    )# function weights so that we can work in probability land as
    # opposed to density land

    MLocal = Dw.DwInv * V.inv' * V.inv * Dw.Dw
    GLocal = Dw.DwInv * V.inv' * V.inv * (V.D * V.inv) * Dw.Dw
    MInvLocal = Dw.DwInv * V.V * V.V' * Dw.Dw
    Phi = (V.inv*V.V)[[1; end], :]

    PosDiagBlock = -Dw.DwInv * Phi[end, :] * Phi[end, :]' * Dw.Dw
    NegDiagBlock = Dw.DwInv * Phi[1, :] * Phi[1, :]' * Dw.Dw
    UpDiagBlock = Dw.DwInv * Phi[end, :] * Phi[1, :]' * Dw.Dw
    LowDiagBlock = Dw.DwInv * Phi[1, :] * Phi[end, :]' * Dw.Dw

    out = (
        G = GLocal, M = MLocal, MInv = MInvLocal, V = V, Phi = Phi, Dw = Dw,
        PosDiagBlock = PosDiagBlock, NegDiagBlock = NegDiagBlock,
        UpDiagBlock = UpDiagBlock, LowDiagBlock = LowDiagBlock,
        )
    v && println("UPDATE: Local operators created; ", keys(out))
    return out
end

function build_lazy_generator(dq::DiscretisedFluidQueue{DGMesh{T}}; v::Bool = false) where T

    m = local_dg_operators(dq.mesh; v=v)
    blocks = (m.LowDiagBlock*m.MInv*2, (m.G+m.PosDiagBlock)*m.MInv*2,
        -(m.G+m.NegDiagBlock)*m.MInv*2, m.UpDiagBlock*m.MInv*2)

    boundary_flux = BoundaryFlux(
        OneBoundaryFlux((m.Dw.DwInv * m.Phi[end, :]*2)[:],
            (m.Phi[end, :]' * m.Dw.Dw * m.MInv)[:]),
        OneBoundaryFlux((m.Dw.DwInv * m.Phi[1, :]*2)[:],
            (m.Phi[1, :]' * m.Dw.Dw * m.MInv)[:])
    )

    D = Matrix{Float64}(LinearAlgebra.I(n_bases_per_cell(dq)))

    out = LazyGenerator(dq,blocks,boundary_flux,D)
    v && println("UPDATE: LazyGenerator object created with keys ", keys(out))
    return out
end


text:
using Pkg

Pkg.develop(PackageSpec(path=joinpath(dirname(@__FILE__), "..")))
# when first running instantiate
Pkg.instantiate()

using Documenter, ProxAL

makedocs(
    sitename = "ProxAL.jl",
    format = Documenter.HTML(
        prettyurls = Base.get(ENV, "CI", nothing) == "true",
        mathengine = Documenter.KaTeX()
    ),
    strict = true,
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "Formulation" => "man/formulation.md",
            "Algorithm" => "man/algorithm.md",
            "Usage" => "man/usage.md",
        ],
        "Library" => [
            "Model parameters" => "lib/modelparams.md",
            "Algorithm parameters" => "lib/algparams.md",
            "Main functions" => "lib/algorithm.md",
            "Developer reference" => "lib/developer.md",
        ]
    ]
)