text:
include("ndgrid.jl")

####################################################################
## Structs for forward KLT tracking

mutable struct F_ImageTrackerVariables

    orgI::Array{Float64,2}                      #the first frame in the image sequence, does not need to be updated after every tracking itteration
    I_nextFrame::Array{Float64,2}               #the next frame in the image sequence, needs to be updated after every tracking itteration
    p_reference::Array{CartesianIndex{2}}       #feature pixel coordinates. p_reference[1,:] = coordinate of the reference feature on orgI. p_reference[2,:] = coordinate of the new feature on I_nextFrame

end

struct F_ImageTrackerConstants

    TrackingType::String
    windowSize::Int8
    number_features::Int8
    x::Array{Float64,2}
    y::Array{Float64,2}
    xder::Array{Int64,2}
    yder::Array{Int64,2}
    DGaussx::Array{Float64,2}
    DGaussy::Array{Float64,2}
    W_Jacobian::Array{Float64,3}

end

## Structs for forward KLT tracking
####################################################################
## Structs for forward KLT tracking with pyramid images

mutable struct FP_ImageTrackerVariables

    orgI::Array{Float64,2}                      #the first frame in the image sequence, does not need to be updated after every tracking itteration
    I_nextFrame::Array{Float64,2}               #the next frame in the image sequence, needs to be updated after every tracking itteration
    orgI_downsample::Array{Float64,2}           #when usning pyramid images orgI_downsample = orgI downsampled by a factor of downsampleFactor
    I_nextFrame_downsample::Array{Float64,2}    #when usning pyramid images I_nextFrame_downsample = I_nextFrame downsampled by a factor of downsampleFactor
    p_reference::Array{CartesianIndex{2}}       #feature pixel coordinates. p_reference[1,:] = coordinate of the reference feature on orgI. p_reference[2,:] = coordinate of the new feature on I_nextFrame

end

struct FP_ImageTrackerConstants

    TrackingType::String
    windowSize::Int8
    number_features::Int8
    downsampleFactor::Int8
    x::Array{Float64,2}
    y::Array{Float64,2}
    xder::Array{Int64,2}
    yder::Array{Int64,2}
    DGaussx::Array{Float64,2}
    DGaussy::Array{Float64,2}
    W_Jacobian::Array{Float64,3}

end

## Structs for forward KLT tracking with pyramid images
####################################################################
## Structs for inversed KLT tracking

mutable struct I_ImageTemplate

    T::Array{Float64,2}
    I_steepest::Array{Float64,2}
    H_inv::Array{Float64,2}

end

mutable struct I_ImageTrackerVariables

    orgI::Array{Float64,2}                      #the first frame in the image sequence, does not need to be updated after every tracking itteration
    I_nextFrame::Array{Float64,2}               #the next frame in the image sequence, needs to be updated after every tracking itteration
    Template::Array{I_ImageTemplate,1}            #template information for the tracked feature
    p_reference::Array{CartesianIndex{2}}       #feature pixel coordinates. p_reference[1,:] = coordinate of the reference feature on orgI. p_reference[2,:] = coordinate of the new feature on I_nextFrame

end

struct I_ImageTrackerConstants

    TrackingType::String
    windowSize::Int8
    number_features::Int8
    x::Array{Float64,2}
    y::Array{Float64,2}
    xder::Array{Int64,2}
    yder::Array{Int64,2}
    DGaussx::Array{Float64,2}
    DGaussy::Array{Float64,2}
    W_Jacobian::Array{Float64,3}

end

## Structs for inversed KLT tracking
####################################################################
## Structs for inversed KLT tracking with pyramid images

mutable struct IP_ImageTemplate

    T::Array{Float64,2}
    T_downsample::Array{Float64,2}
    I_steepest::Array{Float64,2}
    I_steepest_downsample::Array{Float64,2}
    H_inv::Array{Float64,2}
    H_inv_downsample::Array{Float64,2}

end

mutable struct IP_ImageTrackerVariables

    orgI::Array{Float64,2}                      #the first frame in the image sequence, does not need to be updated after every tracking itteration
    I_nextFrame::Array{Float64,2}               #the next frame in the image sequence, needs to be updated after every tracking itteration
    orgI_downsample::Array{Float64,2}           #when usning pyramid images orgI_downsample = orgI downsampled by a factor of downsampleFactor
    I_nextFrame_downsample::Array{Float64,2}    #when usning pyramid images I_nextFrame_downsample = I_nextFrame downsampled by a factor of downsampleFactor
    Template::Array{IP_ImageTemplate,1}            #template information for the tracked feature
    p_reference::Array{CartesianIndex{2}}       #feature pixel coordinates. p_reference[1,:] = coordinate of the reference feature on orgI. p_reference[2,:] = coordinate of the new feature on I_nextFrame

end

struct IP_ImageTrackerConstants

    TrackingType::String
    windowSize::Int8
    number_features::Int8
    downsampleFactor::Int8
    x::Array{Float64,2}
    y::Array{Float64,2}
    xder::Array{Int64,2}
    yder::Array{Int64,2}
    DGaussx::Array{Float64,2}
    DGaussy::Array{Float64,2}
    W_Jacobian::Array{Float64,3}

end

## Structs for inversed KLT tracking with pyramid images
####################################################################

##
"""
    ImageTrackerSetup(orgI_setup, corners,[windowSize = 20, downsampleFactor_setup = 2])

orgI_setup:             Image to do tracking on
corners:                List of corners, preferably a returned list from getApproxBestHarrisInWindow()
windowSize              Half the size of the window to look for a matching feature. Actual window size is windowSize * 2 + 1, windowSize = 20 will result in a 41x41 window
downsampleFactor_setup  downsampleFactor_setup = 2: use pyramid images (first track on a downsampled image, then on the regular image with a smaller windowsize), increases speed. Set downsampleFactor_setup = 1 to use only regular image

Returns two structures, one of type ImageTrackerVariables and one of type ImageTrackerConstants
"""

function ImageTrackerSetup(orgI_setup, corners; windowSize = 20, TrackingType_setup = "Inversed")

    if (TrackingType_setup == "Pyramid" || TrackingType_setup == "InversedPyramid")
        downsampleFactor_setup = 2
    else
        downsampleFactor_setup = 1
    end

    ref_point = CartesianIndex{2}
    ref_point_downsampled = CartesianIndex{2}

    windowSize_setup = round(Int64, windowSize/downsampleFactor_setup)
    orgI_downsample_setup = imresize(orgI_setup, (Int64(length(orgI_setup[:,1])/downsampleFactor_setup),Int64(length(orgI_setup[1,:])/downsampleFactor_setup)));
    number_features_setup = length(corners[:,1])

    # Pixel coordinate matrix
    p_reference_setup = Array{CartesianIndex{2}}(2,number_features_setup)
    for count = 1:number_features_setup
        p_reference_setup[1, count] = corners[count].keypoint
        p_reference_setup[2, count] = corners[count].keypoint
    end

    # Make derivatives kernels
    x_setup, y_setup = ndgrid(0:2*windowSize_setup,0:2*windowSize_setup)

    x_setup=x_setup.-(windowSize_setup+1)
    y_setup=y_setup.-(windowSize_setup+1)
    # x_setup=x_setup.-(windowSize_setup)
    # y_setup=y_setup.-(windowSize_setup)

    sigma = 2;
    xder_setup,yder_setup =ndgrid(floor(-3*sigma):ceil(3*sigma),floor(-3*sigma):ceil(3*sigma));
    DGaussx_setup =-(xder_setup./(2*pi*sigma^4)).*exp(-(xder_setup.^2+yder_setup.^2)/(2*sigma^2));
    DGaussy_setup =-(yder_setup./(2*pi*sigma^4)).*exp(-(xder_setup.^2+yder_setup.^2)/(2*sigma^2));

    DGaussx_reflected_setup = parent(reflect(DGaussx_setup))
    DGaussy_reflected_setup = parent(reflect(DGaussy_setup))

    # Evaluate the Jacobian
    W_Jacobian_x_setup=[x_setup[:] zeros(size(x_setup[:])) y_setup[:] zeros(size(x_setup[:])) ones(size(x_setup[:])) zeros(size(x_setup[:]))]
    W_Jacobian_y_setup=[zeros(size(x_setup[:])) x_setup[:] zeros(size(x_setup[:])) y_setup[:] zeros(size(x_setup[:])) ones(size(x_setup[:]))]

    W_Jacobian_setup = Array{Float64}(2,6,length(x_setup))
    I_steepest_setup = Array{Float64}(length(x_setup),6)

    for j1=1:length(x_setup)
        W_Jacobian_setup[:,:,j1]=[W_Jacobian_x_setup[j1,:] W_Jacobian_y_setup[j1,:]]'
    end

    T = zeros(windowSize_setup * 2 + 1, windowSize_setup * 2 + 1 )
    Ix = zeros(windowSize_setup * 2 + 1, windowSize_setup * 2 + 1 )
    Iy = zeros(windowSize_setup * 2 + 1, windowSize_setup * 2 + 1 )
    T_downsample = zeros(windowSize_setup * 2 + 1, windowSize_setup * 2 + 1 )
    Ix_downsample = zeros(windowSize_setup * 2 + 1, windowSize_setup * 2 + 1 )
    Iy_downsample = zeros(windowSize_setup * 2 + 1, windowSize_setup * 2 + 1 )

    if (TrackingType_setup == "Inversed")
        single_template_setup = I_ImageTemplate(zeros(windowSize_setup*2+1, windowSize_setup*2+1), zeros(length(x_setup), 6), zeros(6, 6))
        Template_setup = Array{I_ImageTemplate,1}(number_features_setup)
    else
        single_template_setup = IP_ImageTemplate(zeros(windowSize_setup*2+1, windowSize_setup*2+1), zeros(windowSize_setup*2+1, windowSize_setup*2+1), zeros(length(x_setup), 6), zeros(length(x_setup), 6), zeros(6, 6), zeros(6, 6))
        Template_setup = Array{IP_ImageTemplate,1}(number_features_setup)
    end
    # Template_setup = Array{ImageTemplate,1}(number_features_setup)

    Ix_grad = imfilter(orgI_setup, centered(DGaussx_reflected_setup));
    Iy_grad = imfilter(orgI_setup, centered(DGaussy_reflected_setup));

    Ix_grad_downsample = imfilter(orgI_downsample_setup, centered(DGaussx_reflected_setup));
    Iy_grad_downsample = imfilter(orgI_downsample_setup, centered(DGaussy_reflected_setup));

    for count = 1:number_features_setup

        x_min = round(Int64,p_reference_setup[1,count][1]-windowSize_setup)
        x_max = round(Int64,p_reference_setup[1,count][1]+windowSize_setup)
        y_min = round(Int64,p_reference_setup[1,count][2]-windowSize_setup)
        y_max = round(Int64,p_reference_setup[1,count][2]+windowSize_setup)

        x_min_downsample = round(Int64,p_reference_setup[1,count][1]/downsampleFactor_setup-windowSize_setup)
        x_max_downsample = round(Int64,p_reference_setup[1,count][1]/downsampleFactor_setup+windowSize_setup)
        y_min_downsample = round(Int64,p_reference_setup[1,count][2]/downsampleFactor_setup-windowSize_setup)
        y_max_downsample = round(Int64,p_reference_setup[1,count][2]/downsampleFactor_setup+windowSize_setup)

        T = orgI_setup[x_min : x_max , y_min : y_max]
        T_downsample = orgI_downsample_setup[x_min_downsample : x_max_downsample , y_min_downsample : y_max_downsample]

        single_template_setup.T = T
        if (TrackingType_setup == "InversedPyramid")
            single_template_setup.T_downsample = T_downsample
        end

        Ix = Ix_grad[x_min : x_max , y_min : y_max]
        Iy = Iy_grad[x_min : x_max , y_min : y_max]

        Ix_downsample = Ix_grad_downsample[x_min_downsample : x_max_downsample , y_min_downsample : y_max_downsample]
        Iy_downsample = Iy_grad_downsample[x_min_downsample : x_max_downsample , y_min_downsample : y_max_downsample]

        # W_p = [ 1 0 p_reference_setup[1,count][1]; 0 1 p_reference_setup[1,count][2]];

        Gradient1 = 0
        for j1=1:length(x_setup)
            Gradient1=[Ix[j1] Iy[j1]];
            Gradient1_downsample=[Ix_downsample[j1] Iy_downsample[j1]];

            single_template_setup.I_steepest[j1,1:6] = Gradient1 * W_Jacobian_setup[:,:,j1];
            if (TrackingType_setup == "InversedPyramid")
                single_template_setup.I_steepest_downsample[j1,1:6] = Gradient1_downsample * W_Jacobian_setup[:,:,j1];
            end
        end

        H=zeros(6,6);
        H_downsample=zeros(6,6);
        for j2=1:length(x_setup)
            H = H + single_template_setup.I_steepest[j2,:] * single_template_setup.I_steepest[j2,:]'
            if (TrackingType_setup == "InversedPyramid")
                H_downsample = H_downsample + single_template_setup.I_steepest_downsample[j2,:] * single_template_setup.I_steepest_downsample[j2,:]'
            end
        end

        single_template_setup.H_inv = inv(H)
        if (TrackingType_setup == "InversedPyramid")
            single_template_setup.H_inv_downsample = inv(H_downsample)
        end

        Template_setup[count] = deepcopy(single_template_setup)

    end

    # ImageTrackerVariables
    orgI_constructor                    = orgI_setup
    I_nextFrame_constructor             = orgI_setup
    orgI_downsample_constructor         = orgI_downsample_setup
    I_nextFrame_downsample_constructor  = orgI_downsample_setup
    Template_constructor                = Template_setup
    p_reference_constructor             = p_reference_setup

    # ImageTrackerConstants
    TrackingType_constructor            = TrackingType_setup
    windowSize_constructor              = windowSize_setup
    number_features_constructor         = number_features_setup
    downsampleFactor_constructor        = downsampleFactor_setup
    x_constructor                       = x_setup
    y_constructor                       = y_setup
    xder_constructor                    = xder_setup
    yder_constructor                    = yder_setup
    DGaussx_constructor                 = DGaussx_reflected_setup
    DGaussy_constructor                 = DGaussy_reflected_setup
    W_Jacobian_constructor              = W_Jacobian_setup

    if (TrackingType_setup == "Forward")

        ITVar = F_ImageTrackerVariables(    orgI_constructor,
                                            I_nextFrame_constructor,
                                            p_reference_constructor)

        #
        ITConst = F_ImageTrackerConstants(  TrackingType_constructor,
                                            windowSize_constructor,
                                            number_features_constructor,
                                            x_constructor,
                                            y_constructor,
                                            xder_constructor,
                                            yder_constructor,
                                            DGaussx_constructor,
                                            DGaussy_constructor,
                                            W_Jacobian_constructor)

    elseif (TrackingType_setup == "Pyramid")

        ITVar = FP_ImageTrackerVariables(   orgI_constructor,
                                            I_nextFrame_constructor,
                                            orgI_downsample_constructor,
                                            I_nextFrame_downsample_constructor,
                                            p_reference_constructor)

        #
        ITConst = FP_ImageTrackerConstants( TrackingType_constructor,
                                            windowSize_constructor,
                                            number_features_constructor,
                                            downsampleFactor_constructor,
                                            x_constructor,
                                            y_constructor,
                                            xder_constructor,
                                            yder_constructor,
                                            DGaussx_constructor,
                                            DGaussy_constructor,
                                            W_Jacobian_constructor)

    elseif (TrackingType_setup == "Inversed")

        ITVar = I_ImageTrackerVariables(    orgI_constructor,
                                            I_nextFrame_constructor,
                                            Template_constructor ,
                                            p_reference_constructor)

        #
        ITConst = I_ImageTrackerConstants(  TrackingType_constructor,
                                            windowSize_constructor,
                                            number_features_constructor,
                                            x_constructor,
                                            y_constructor,
                                            xder_constructor,
                                            yder_constructor,
                                            DGaussx_constructor,
                                            DGaussy_constructor,
                                            W_Jacobian_constructor)

    else

        ITVar = IP_ImageTrackerVariables(   orgI_constructor,
                                            I_nextFrame_constructor,
                                            orgI_downsample_constructor,
                                            I_nextFrame_downsample_constructor,
                                            Template_constructor ,
                                            p_reference_constructor)

        #
        ITConst = IP_ImageTrackerConstants( TrackingType_constructor,
                                            windowSize_constructor,
                                            number_features_constructor,
                                            downsampleFactor_constructor,
                                            x_constructor,
                                            y_constructor,
                                            xder_constructor,
                                            yder_constructor,
                                            DGaussx_constructor,
                                            DGaussy_constructor,
                                            W_Jacobian_constructor)

    end

    return ITVar, ITConst

end
#####################
function fillNewImageTemplates!(ITVar, ITConst)

    for count = 1:ITConst.number_features
        testedges(ITVar, ITConst, count)
    end
    #update reference image and reference feature coordinates
    ITVar.orgI = ITVar.I_nextFrame
    if (ITConst.TrackingType == "Pyramid" || ITConst.TrackingType == "InversedPyramid")
        ITVar.orgI_downsample = ITVar.I_nextFrame_downsample
    end
    ITVar.p_reference[1, :] = ITVar.p_reference[2, :]

    Ix = zeros(ITConst.windowSize * 2 + 1, ITConst.windowSize * 2 + 1 )
    Iy = zeros(ITConst.windowSize * 2 + 1, ITConst.windowSize * 2 + 1 )

    # Evaluate the gradient of the template
    Ix_grad = imfilter(ITVar.orgI, centered(ITConst.DGaussx));
    Iy_grad = imfilter(ITVar.orgI, centered(ITConst.DGaussy));

    if (ITConst.TrackingType == "Pyramid" || ITConst.TrackingType == "InversedPyramid")
        Ix_grad_downsample = imfilter(ITVar.orgI_downsample, centered(ITConst.DGaussx));
        Iy_grad_downsample = imfilter(ITVar.orgI_downsample, centered(ITConst.DGaussy));
    end

    if (ITConst.TrackingType == "Inversed" || ITConst.TrackingType == "InversedPyramid")

        for count = 1:ITConst.number_features

            if (ITVar.p_reference[2, count][1] == 0)

            else

                x_min = round(Int64, ITVar.p_reference[2, count][1])-ITConst.windowSize
                x_max = round(Int64, ITVar.p_reference[2, count][1])+ITConst.windowSize
                y_min = round(Int64, ITVar.p_reference[2, count][2])-ITConst.windowSize
                y_max = round(Int64, ITVar.p_reference[2, count][2])+ITConst.windowSize

                if (ITConst.TrackingType == "InversedPyramid")
                    x_min_downsample = round(Int64, ITVar.p_reference[2, count][1]/ITConst.downsampleFactor)-ITConst.windowSize
                    x_max_downsample = round(Int64, ITVar.p_reference[2, count][1]/ITConst.downsampleFactor)+ITConst.windowSize
                    y_min_downsample = round(Int64, ITVar.p_reference[2, count][2]/ITConst.downsampleFactor)-ITConst.windowSize
                    y_max_downsample = round(Int64, ITVar.p_reference[2, count][2]/ITConst.downsampleFactor)+ITConst.windowSize
                end

                ITVar.Template[count].T = ITVar.orgI[x_min:x_max,y_min:y_max]
                if (ITConst.TrackingType == "InversedPyramid")
                    ITVar.Template[count].T_downsample = ITVar.orgI_downsample[x_min_downsample:x_max_downsample,y_min_downsample:y_max_downsample]
                end

                Ix = Ix_grad[x_min:x_max,y_min:y_max]
                Iy = Iy_grad[x_min:x_max,y_min:y_max]

                if (ITConst.TrackingType == "InversedPyramid")
                    Ix_downsample = Ix_grad_downsample[x_min_downsample:x_max_downsample,y_min_downsample:y_max_downsample]
                    Iy_downsample = Iy_grad_downsample[x_min_downsample:x_max_downsample,y_min_downsample:y_max_downsample]
                end

                Gradient1 = 0
                for j1=1:length(ITConst.x)
                    Gradient1=[Ix[j1] Iy[j1]];
                    if (ITConst.TrackingType == "InversedPyramid")
                        Gradient1_downsample=[Ix_downsample[j1] Iy_downsample[j1]];
                    end

                    ITVar.Template[count].I_steepest[j1,1:6] = Gradient1 * ITConst.W_Jacobian[:,:,j1];
                    if (ITConst.TrackingType == "InversedPyramid")
                        ITVar.Template[count].I_steepest_downsample[j1,1:6] = Gradient1_downsample * ITConst.W_Jacobian[:,:,j1];
                    end
                end

                H=zeros(6,6);
                if (ITConst.TrackingType == "InversedPyramid")
                    H_downsample=zeros(6,6);
                end
                for j2=1:length(ITConst.x)
                    H = H + ITVar.Template[count].I_steepest[j2,:] * ITVar.Template[count].I_steepest[j2,:]'
                    if (ITConst.TrackingType == "InversedPyramid")
                        H_downsample = H_downsample + ITVar.Template[count].I_steepest_downsample[j2,:] * ITVar.Template[count].I_steepest_downsample[j2,:]'
                    end
                end

                ITVar.Template[count].H_inv = inv(H)
                if (ITConst.TrackingType == "InversedPyramid")
                    ITVar.Template[count].H_inv_downsample = inv(H_downsample)
                end
            end
        end
    end

end
####################
"""
    testedges(ITVar, ITConst, corner_i)

ITVar:   construct of variables of type ImageTrackerVariables
ITConst: construct of constants of type ImageTrackerConstants
corner_i: Feature to test

Returns true if the frame is within the bounds of the image, false otherwise
"""
function testedges(ITVar, ITConst, corner_i)

    if (ITConst.TrackingType == "Forward" || ITConst.TrackingType == "Inversed")
        if (ITVar.p_reference[2, corner_i][1]-ITConst.windowSize > 0 && ITVar.p_reference[2, corner_i][1]+ITConst.windowSize <= size(ITVar.orgI, 1) && ITVar.p_reference[2, corner_i][2]-ITConst.windowSize > 0 && ITVar.p_reference[2, corner_i][2]+ITConst.windowSize < size(ITVar.orgI, 2))
            return true
        else
            ITVar.p_reference[2, corner_i] = CartesianIndex(0, 0)
            ITVar.p_reference[1, corner_i] = CartesianIndex(0, 0)
            return false
        end
    elseif (ITConst.TrackingType == "Pyramid" || ITConst.TrackingType == "InversedPyramid")
        if (ITVar.p_reference[2, corner_i][1]-(ITConst.windowSize*2+1) > 0 && ITVar.p_reference[2, corner_i][1]+(ITConst.windowSize*2+1) < size(ITVar.orgI, 1) && ITVar.p_reference[2, corner_i][2]-(ITConst.windowSize*2+1) > 0 && ITVar.p_reference[2, corner_i][2]+(ITConst.windowSize*2+1) < size(ITVar.orgI, 2))
            return true
        else
            ITVar.p_reference[2, corner_i] = CartesianIndex(0, 0)
            ITVar.p_reference[1, corner_i] = CartesianIndex(0, 0)
            return false
        end
    end
end

#####################
"""
    KTL_Tracker!(ITVar, ITConst)

ITVar:   construct of variables of type ImageTrackerVariables
ITConst: construct of constants of type ImageTrackerConstants

Returns nothing, updates p_reference[2,:] with the new coordinates of the tracked features
"""

function KTL_Tracker!(ITVar, ITConst)

    for corner_i = 1:ITConst.number_features
        if (testedges(ITVar, ITConst, corner_i))
            if (ITConst.TrackingType == "Forward")
                # slower code, does not use pyramid images============================================
                trackOneFeature(ITVar, ITConst, corner_i)
                # slower code, does not use pyramid images============================================
            elseif (ITConst.TrackingType == "Pyramid")
                # faster code, use pyramid images=====================================================
                trackOneFeaturePyramid(ITVar, ITConst, corner_i)
                # faster code, use pyramid images=====================================================
            elseif (ITConst.TrackingType == "Inversed")
                # slower code, does not use pyramid images============================================
                trackOneFeatureInverse(ITVar, ITConst, corner_i)
                # slower code, does not use pyramid images============================================
            elseif (ITConst.TrackingType == "InversedPyramid")
                # slower code, does not use pyramid images============================================
                trackOneFeatureInversePyramid(ITVar, ITConst, corner_i)
                # slower code, does not use pyramid images============================================
            end

        else
            ITVar.p_reference[2, corner_i] = CartesianIndex(0, 0)
        end
    end
    return nothing
end

####################
"""
    trackOneFeatureInversePyramid(ITVar, ITConst, corner_i)

ITVar:   construct of variables of type ImageTrackerVariables
ITConst: construct of constants of type ImageTrackerConstants
corner_i:index of selected corner in ITVar.p_reference to do tracking on

Returns nothing, updates p_reference[2,corner_i] with the new coordinates of the tracked feature
"""

function trackOneFeatureInversePyramid(ITVar, ITConst, corner_i)
    Threshold = 0.2
    converged = true

    I_warped = zeros(ITConst.windowSize * 2 + 1, ITConst.windowSize * 2 + 1 )
    Ix = zeros(ITConst.windowSize * 2 + 1, ITConst.windowSize * 2 + 1 )
    Iy = zeros(ITConst.windowSize * 2 + 1, ITConst.windowSize * 2 + 1 )

    p = [0 0 0 0 round(Int64, ITVar.p_reference[2, corner_i][1]/2) round(Int64, ITVar.p_reference[2, corner_i][2]/2)]

    delta_p = [0 0 0 0 100 100]

    for count_pyramid = 1:2

        delta_p = [0 0 0 0 100 100]

        counter = 0;

        while ( norm(delta_p) > Threshold)
            counter= counter + 1;
            converged = true
            if(counter > 100)
                converged = false
                break;
            end

            #The affine matrix for template rotation and translation
            W_p = [ 1+p[1] p[3] p[5]; p[2] 1+p[4] p[6]];

            #1 Warp I with w
            if count_pyramid == 1
                warping!(I_warped, ITVar.I_nextFrame_downsample, ITConst.x, ITConst.y, W_p)
                I_error= I_warped - ITVar.Template[corner_i].T_downsample
            else
                warping!(I_warped, ITVar.I_nextFrame, ITConst.x, ITConst.y, W_p)
                I_error= I_warped - ITVar.Template[corner_i].T
            end

            # Break if outside image
            if (!testedges(ITVar, ITConst, corner_i))
                converged = false
                break;
            end

            #2 Multiply steepest descend with error

            if count_pyramid == 1
                total=zeros(6,1);
                for j3=1:length(ITConst.x)
                    total = total + (ITVar.Template[corner_i].I_steepest_downsample[j3,:]' * I_error[j3])'
                end
            else
                total=zeros(6,1);
                for j3=1:length(ITConst.x)
                    total = total + (ITVar.Template[corner_i].I_steepest[j3,:]' * I_error[j3])'
                end
            end

            #3 Computer delta_p

            if count_pyramid == 1
                delta_p = ITVar.Template[corner_i].H_inv_downsample*total;

            else
                delta_p = ITVar.Template[corner_i].H_inv*total;
            end

            #4 Update the parameters p <- p + delta_p
            p = p - 0.1 * delta_p';
        end
        if count_pyramid ==1
            p = [0 0 0 0 round(Int64, p[5]*2) round(Int64, p[6]*2)]
        end
    end
    if (converged == true)
        if (testedges(ITVar, ITConst, corner_i) == true)
            ITVar.p_reference[2, corner_i] = CartesianIndex(round(Int64,p[5]), round(Int64,p[6]))
        else
            ITVar.p_reference[2, corner_i] = CartesianIndex(0, 0)
        end
    else
        ITVar.p_reference[2, corner_i] = CartesianIndex(0, 0)
    end
end

#############
####################
"""
    trackOneFeatureInverse(ITVar, ITConst, corner_i)

ITVar:   construct of variables of type ImageTrackerVariables
ITConst: construct of constants of type ImageTrackerConstants
corner_i:index of selected corner in ITVar.p_reference to do tracking on

Returns nothing, updates p_reference[2,corner_i] with the new coordinates of the tracked feature
Does not use orgI_downsample or I_nextFrame_downsample at all, downsampleFactor must be 1
"""

function trackOneFeatureInverse(ITVar, ITConst, corner_i)
    Threshold = 0.2
    converged = true

    I_warped = zeros(ITConst.windowSize * 2 + 1, ITConst.windowSize * 2 + 1 )
    Ix = zeros(ITConst.windowSize * 2 + 1, ITConst.windowSize * 2 + 1 )
    Iy = zeros(ITConst.windowSize * 2 + 1, ITConst.windowSize * 2 + 1 )

    p = [0 0 0 0 ITVar.p_reference[2, corner_i][1] ITVar.p_reference[2, corner_i][2]]

    delta_p = [0 0 0 0 100 100]

    counter = 0;

    while ( norm(delta_p) > Threshold)
        counter= counter + 1;
        converged = true
        if(counter > 100)
            converged = false
            break;
        end

        #The affine matrix for template rotation and translation
        W_p = [ 1+p[1] p[3] p[5]; p[2] 1+p[4] p[6]];

        #1 Warp I with w
        warping!(I_warped, ITVar.I_nextFrame, ITConst.x, ITConst.y, W_p)

        #2 Subtract T from I
        I_error= I_warped - ITVar.Template[corner_i].T

        # Break if outside image
        if((p[5]>(size(ITVar.I_nextFrame,1))-1)||(p[6]>(size(ITVar.I_nextFrame,2)-1))||(p[5]<0)||(p[6]<0))
            break;#put back
        end

        #3 Multiply steepest descend with error
        total=zeros(6,1);
        for j3=1:length(ITConst.x)
            total = total + (ITVar.Template[corner_i].I_steepest[j3,:]' * I_error[j3])'
        end

        #4 Computer delta_p
        delta_p = ITVar.Template[corner_i].H_inv*total;

        #5 Update the parameters p <- p - delta_p
        p = p - 0.1 * delta_p';

    end
    if (converged == true)
        # @show converged
        if (testedges(ITVar, ITConst, corner_i))
            ITVar.p_reference[2, corner_i] = CartesianIndex(round(Int,p[5]), round(Int,p[6]))
        else
            ITVar.p_reference[2, corner_i] = CartesianIndex(0, 0)
        end
    else
        # @show converged
        ITVar.p_reference[2, corner_i] = CartesianIndex(0, 0)
    end
end

#############
####################
"""
    trackOneFeaturePyramid(ITVar, ITConst, corner_i)

ITVar:   construct of variables of type ImageTrackerVariables
ITConst: construct of constants of type ImageTrackerConstants
corner_i:index of selected corner in ITVar.p_reference to do tracking on

Returns nothing, updates p_reference[2,corner_i] with the new coordinates of the tracked feature
Use orgI_downsample and I_nextFrame_downsample to track features, downsampleFactor must be 2 for now, can be increased in future work
"""

function trackOneFeaturePyramid(ITVar, ITConst, corner_i)
    Threshold = 0.2
    converged = true

    I_warped = zeros(ITConst.windowSize * 2 + 1, ITConst.windowSize * 2 + 1 )
    Ix = zeros(ITConst.windowSize * 2 + 1, ITConst.windowSize * 2 + 1 )
    Iy = zeros(ITConst.windowSize * 2 + 1, ITConst.windowSize * 2 + 1 )

    p = [0 0 0 0 round(Int, ITVar.p_reference[2, corner_i][1]/2) round(Int, ITVar.p_reference[2, corner_i][2]/2)]

    for count_pyramid = 1:2

        if count_pyramid == 1
            if (testedges(ITVar, ITConst, corner_i))
                T = ITVar.orgI_downsample[round(Int64, ITVar.p_reference[1, corner_i][1]/2)-ITConst.windowSize:round(Int64, ITVar.p_reference[1, corner_i][1]/2)+ITConst.windowSize,round(Int64, ITVar.p_reference[1, corner_i][2]/2)-ITConst.windowSize:round(Int64, ITVar.p_reference[1, corner_i][2]/2)+ITConst.windowSize]
                T = convert(Array{Float64}, T)
            else
                break
            end
        else
            T = ITVar.orgI[round(Int64, ITVar.p_reference[1, corner_i][1])-ITConst.windowSize:round(Int64, ITVar.p_reference[1, corner_i][1])+ITConst.windowSize,round(Int64, ITVar.p_reference[1, corner_i][2])-ITConst.windowSize:round(Int64, ITVar.p_reference[1, corner_i][2])+ITConst.windowSize]
            T = convert(Array{Float64}, T)
        end

        delta_p = [0 0 0 0 100 100]

        # Filter the images to get the derivatives
        if count_pyramid == 1
            Ix_grad = imfilter(ITVar.I_nextFrame_downsample, centered(ITConst.DGaussx));
            Iy_grad = imfilter(ITVar.I_nextFrame_downsample, centered(ITConst.DGaussy));
        else
            Ix_grad = imfilter(ITVar.I_nextFrame, centered(ITConst.DGaussx));
            Iy_grad = imfilter(ITVar.I_nextFrame, centered(ITConst.DGaussy));
        end
        counter = 0;

        while ( norm(delta_p) > Threshold)
            counter= counter + 1;
            converged = true
            if(counter > 100)
                converged = false
                break;
            end

            #The affine matrix for template rotation and translation
            W_p = [ 1+p[1] p[3] p[5]; p[2] 1+p[4] p[6]];

            #1 Warp I with w
            if count_pyramid == 1
                warping!(I_warped, ITVar.I_nextFrame_downsample, ITConst.x, ITConst.y, W_p)
            else
                warping!(I_warped, ITVar.I_nextFrame, ITConst.x, ITConst.y, W_p)
            end

            #2 Subtract I from T
            I_error= T - I_warped

            # Break if outside image
            if((p[5]>(size(ITVar.I_nextFrame,1))-1)||(p[6]>(size(ITVar.I_nextFrame,2)-1))||(p[5]<0)||(p[6]<0))
                break;
            end

            #3 Warp the gradient
            warping!(Ix, Ix_grad, ITConst.x, ITConst.y, W_p);
            warping!(Iy, Iy_grad, ITConst.x, ITConst.y, W_p);

            #4 Compute steepest descent
            I_steepest=zeros(length(ITConst.x),6);
            Gradient1 = 0
            W_Jacobian = 0
            for j1=1:length(ITConst.x)
                # W_Jacobian=[ITConst.W_Jacobian_x[j1,:] ITConst.W_Jacobian_y[j1,:]]';
                Gradient1=[Ix[j1] Iy[j1]];
                I_steepest[j1,1:6] = Gradient1 * ITConst.W_Jacobian[:,:,j1];
            end

            #5 Compute Hessian
            H=zeros(6,6);
            for j2=1:length(ITConst.x)
                H = H + I_steepest[j2,:] * I_steepest[j2,:]'
            end

            #6 Multiply steepest descend with error
            total=zeros(6,1);
            for j3=1:length(ITConst.x)
                total = total + (I_steepest[j3,:]' * I_error[j3])'
            end

            #7 Computer delta_p
            delta_p=H\total;

            #8 Update the parameters p <- p + delta_p
            p = p + 0.1 * delta_p';
        end
        if count_pyramid ==1
            p = [0 0 0 0 round(Int64, p[5]*2) round(Int64, p[6]*2)]
        end
    end
    if (converged == true)
        if (testedges(ITVar, ITConst, corner_i))
            ITVar.p_reference[2, corner_i] = CartesianIndex(round(Int64,p[5]), round(Int64,p[6]))
        else
            ITVar.p_reference[2, corner_i] = CartesianIndex(0, 0)
        end
    else
        ITVar.p_reference[2, corner_i] = CartesianIndex(0, 0)
    end
end

####################
####################
"""
    trackOneFeature(ITVar, ITConst, corner_i)

ITVar:   construct of variables of type ImageTrackerVariables
ITConst: construct of constants of type ImageTrackerConstants
corner_i:index of selected corner in ITVar.p_reference to do tracking on

Returns nothing, updates p_reference[2,corner_i] with the new coordinates of the tracked feature
Does not use orgI_downsample or I_nextFrame_downsample at all, downsampleFactor must be 1
"""

function trackOneFeature(ITVar, ITConst, corner_i)
    Threshold = 0.2
    converged = true

    I_warped = zeros(ITConst.windowSize * 2 + 1, ITConst.windowSize * 2 + 1 )
    Ix = zeros(ITConst.windowSize * 2 + 1, ITConst.windowSize * 2 + 1 )
    Iy = zeros(ITConst.windowSize * 2 + 1, ITConst.windowSize * 2 + 1 )

    p = [0 0 0 0 ITVar.p_reference[2, corner_i][1] ITVar.p_reference[2, corner_i][2]]

    T = ITVar.orgI[round(Int64, ITVar.p_reference[1, corner_i][1])-ITConst.windowSize:round(Int64, ITVar.p_reference[1, corner_i][1])+ITConst.windowSize,round(Int64, ITVar.p_reference[1, corner_i][2])-ITConst.windowSize:round(Int64, ITVar.p_reference[1, corner_i][2])+ITConst.windowSize]
    T = convert(Array{Float64}, T)

    delta_p = [0 0 0 0 100 100]

    # Filter the images to get the derivatives
    Ix_grad = imfilter(ITVar.I_nextFrame, centered(ITConst.DGaussx));
    Iy_grad = imfilter(ITVar.I_nextFrame, centered(ITConst.DGaussy));

    counter = 0;

    while ( norm(delta_p) > Threshold)
        counter= counter + 1;
        converged = true
        if(counter > 100)
            converged = false
            break;
        end

        #The affine matrix for template rotation and translation
        W_p = [ 1+p[1] p[3] p[5]; p[2] 1+p[4] p[6]];

        #1 Warp I with w
        warping!(I_warped, ITVar.I_nextFrame, ITConst.x, ITConst.y, W_p)

        #2 Subtract I from T
        I_error= T - I_warped

        # Break if outside image
        if((p[5]>(size(ITVar.I_nextFrame,1))-1)||(p[6]>(size(ITVar.I_nextFrame,2)-1))||(p[5]<0)||(p[6]<0))
            # @show converged = false
            break;#put back
        end

        #3 Warp the gradient
        warping!(Ix, Ix_grad, ITConst.x, ITConst.y, W_p);
        warping!(Iy, Iy_grad, ITConst.x, ITConst.y, W_p);

        #4 Compute steepest descent
        I_steepest=zeros(length(ITConst.x),6);
        Gradient1 = 0
        W_Jacobian = 0
        for j1=1:length(ITConst.x)
            # W_Jacobian=[ITConst.W_Jacobian_x[j1,:] ITConst.W_Jacobian_y[j1,:]]';
            Gradient1=[Ix[j1] Iy[j1]];
            I_steepest[j1,1:6] = Gradient1 * ITConst.W_Jacobian[:,:,j1]#W_Jacobian;
        end

        #5 Compute Hessian
        H=zeros(6,6);
        for j2=1:length(ITConst.x)
            H = H + I_steepest[j2,:] * I_steepest[j2,:]'
        end

        #6 Multiply steepest descend with error
        total=zeros(6,1);
        for j3=1:length(ITConst.x)
            total = total + (I_steepest[j3,:]' * I_error[j3])'
        end

        #7 Computer delta_p
        delta_p=H\total;

        #8 Update the parameters p <- p + delta_p
        p = p + 0.1 * delta_p';
    end
    if (converged == true)
        if (testedges(ITVar, ITConst, corner_i))
            ITVar.p_reference[2, corner_i] = CartesianIndex(round(Int64,p[5]), round(Int64,p[6]))
        else
            ITVar.p_reference[2, corner_i] = CartesianIndex(0, 0)
        end
    else
        ITVar.p_reference[2, corner_i] = CartesianIndex(0, 0)
    end
end

#############
#############
"""
    warping!(Iout,Iin,x,y,M)

Iout:   returned warped image
Iin:    image to be warped
x:      derivatives kernel on x axis
y:      derivatives kernel on y axis
M:      affine matrix for template rotation and translation

Return a warped image
"""
function warping!(Iout,Iin,x,y,M)
    # Affine transformation function (Rotation, Translation, Resize)
    # Calculate the Transformed coordinates

    Tlocalx =  M[1,1] .* x .+ M[1,2] .*y .+ M[1,3]
    Tlocaly =  M[2,1] .* x .+ M[2,2] .*y .+ M[2,3]

    #Iout  = interp2(Iin, Tlocalx, Tlocaly,'*linear');
    # All the neighborh pixels involved in linear interpolation.
    xBas0=floor.(Int64, Tlocalx)
    yBas0=floor.(Int64, Tlocaly)
    xBas1=xBas0+1
    yBas1=yBas0+1

    # Linear interpolation constants (percentages)
    xCom=Tlocalx-xBas0
    yCom=Tlocaly-yBas0
    perc0=(1-xCom).*(1-yCom)
    perc1=(1-xCom).*yCom
    perc2=xCom.*(1-yCom)
    perc3=xCom.*yCom

    # limit indexes to boundaries
    check_xBas0=(xBas0.<0).|(xBas0.>(size(Iin,1)-1))
    check_yBas0=(yBas0.<0).|(yBas0.>(size(Iin,2)-1))
    xBas0[check_xBas0]=0
    yBas0[check_yBas0]=0
    check_xBas1=(xBas1.<0).|(xBas1.>(size(Iin,1).-1))
    check_yBas1=(yBas1.<0).|(yBas1.>(size(Iin,2).-1))
    xBas1[check_xBas1]=0
    yBas1[check_yBas1]=0

    # Iout=zeros(Float64, size(x, 1), size(x, 2))

    Iin_one = Iin
    # Get the intensities
    intensity_xyz0=Iin_one[1+xBas0+yBas0*size(Iin,1)];
    intensity_xyz1=Iin_one[1+xBas0+yBas1*size(Iin,1)];
    intensity_xyz2=Iin_one[1+xBas1+yBas0*size(Iin,1)];
    intensity_xyz3=Iin_one[1+xBas1+yBas1*size(Iin,1)];

    Iout[:,:] .= intensity_xyz0 .* perc0 .+ intensity_xyz1 .* perc1 .+ intensity_xyz2 .* perc2 .+ intensity_xyz3 .* perc3;
    # Iout[:,:,i]=reshape(Iout_one, (size(x,1) size(x,2)));
    # Iout[:,:]=reshape(Iout_one, size(x,1), size(x,2))

    # end
    return nothing

end

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

max_stars_repo_name: dehann/SensorFeatureTracking

max_issues_repo_name: dehann/SensorFeatureTracking

max_forks_repo_name: dehann/SensorFeatureTracking


text:
# When scaling of the amplitude is allowed, we call it's shape.
# TODO Technically, it should be supertype or subtype of amplitude constrain? Just now I do it for easyness as subtype.

"""
ConstrainedByShape{T,N}

Set of abstract arryas with elemet typ `T` and dimensions `N defined by the shape constraint `|x| = s amp` for some `s`.
Field n contains the sum of square of all elements.
"""
struct ConstrainedByShape{T,N} <: AmplitudeConstrainedSet where {T <: Real, N}
    amp::Array{T,N}  #
    n::T
end

"""
ConstrainedByShape(a::AbstractArray{T,N})

Construct set defined by the shape constraint `|x| = s a` for some scaling `s`. Type and dimension of the set are inhereited from the array.
"""
function ConstrainedByShape(a::AbstractArray{T,N}) where {T,N}
    ConstrainedByShape{T,N}(a, sum(a .^2))
end

# function ConstrainedByShape(a::AbstractArray{Union{T,Nothing},N}) where {T,N}
#     ConstrainedByShape{Union{T,Nothing},N}(a)
# end

# function ConstrainedByAmplitude(a::AbstractArray{T,N})
#     ConstrainedByAmplitude{T,N}(a)
# end
export ConstrainedByShape

"""
ConstrainedByAmplitudeMasked(a, mask::Vector)
ConstrainedByAmplitudeMasked(a, AbstractArray{Bool})

The amplitude constrained set only in the indexes given by mask:  `|xᵢ| = aᵢ` for `i ∈ mask`.
"""
struct ConstrainedByShapeMasked{T,N} <: AmplitudeConstrainedSet where {T <: Real, N}
    amp::Array{T,N}  #
    mask::Union{Vector{CartesianIndex{N}}, Vector{Int}} #because 1D arrays are indexed as Vector
    n::T
end

ConstrainedByShapeMasked(amp, mask::Vector) = ConstrainedByShapeMasked(amp, mask, sum(abs2,amp[mask]))

ConstrainedByShapeMasked(a, m::AbstractArray{Bool}) = ConstrainedByShapeMasked(a, findall(m))

export ConstrainedByShapeMasked

"""
ConstrainedByShapeSaturated(a, mask::Vector)
ConstrainedByShapeSaturated(a, AbstractArray{Bool})

The amplitude constrained set only in the indexes given by mask:  `|xᵢ| = s aᵢ` for `i ∈ mask` for some `s`.
Outside the mask function should be larger than the satruation level `|xᵢ| > s ` for `i ∉ mask`.

For this set, the amplitude should be provided in the range from 0 to 1 (1 corresponding to the saturated values).
"""
struct ConstrainedByShapeSaturated{T,N} <: AmplitudeConstrainedSet where {T <: Real, N}
    amp::Array{T,N}  #
    mask::Union{Vector{CartesianIndex{N}}, Vector{Int}} #because 1D arrays are indexed as Vector
    sat::Union{Vector{CartesianIndex{N}}, Vector{Int}} # complementary set to mask
    n::T
end

ConstrainedByShapeSaturated(amp, mask::Vector, sat::Vector) = ConstrainedByShapeSaturated(amp, mask,sat, sum(abs2,amp[mask]))

ConstrainedByShapeSaturated(a, m::AbstractArray{Bool}) = ConstrainedByShapeSaturated(a, findall(m), findall(.!(m)))

export ConstrainedByShapeSaturated

function project!(xp, x, feasset::ConstrainedByShape)
    s = abs.(x)[:]' * feasset.amp[:] / feasset.n
    @inbounds for i in eachindex(xp)
        xp[i] = update_amplitude(s * feasset.amp[i], x[i])
        end
    return xp
end

function project!(x, feasset::ConstrainedByShape)
    s = abs.(x)[:]' * feasset.amp[:] / feasset.n
    for i in eachindex(x)
        x[i] = update_amplitude(s * feasset.amp[i], x[i])
    end

    return x
end

function project!(xp, x, feasset::ConstrainedByShapeMasked)
    s = abs.(x)[feasset.mask]' * feasset.amp[feasset.mask] / feasset.n
    # s = x[feasset.mask]' * feasset.amp[feasset.mask] / feasset.n
    @inbounds for i in feasset.mask
        xp[i] = update_amplitude(s * feasset.amp[i], x[i])
    end
    return xp
end

function project!(x, feasset::ConstrainedByShapeMasked)
    s = abs.(x)[feasset.mask]' * feasset.amp[feasset.mask] / feasset.n
    @inbounds for i in feasset.mask
        x[i] = update_amplitude(s * feasset.amp[i], x[i])
    end
    return x
end

using Roots

function project!(xp, x, feasset::ConstrainedByShapeSaturated)
    s0 = abs.(x)[feasset.mask]' * feasset.amp[feasset.mask] / feasset.n

    xs = sort(abs.(x[feasset.sat]))
    function rp(s)
        j = searchsortedlast(xs,s)
        return (s - s0) * feasset.n + j*s - sum(xs[1:j])
    end
    # println(s0, extrema(xs)) # debug
    smin, smax = extrema(xs)
    sopt = find_zero(rp,(min(smin,s0), max(smax,s0)))

    @inbounds for i in feasset.mask
        xp[i] = update_amplitude(sopt * feasset.amp[i], x[i])
    end
    @inbounds for i in feasset.sat
        xp[i] = update_amplitude(max.(abs(x[i]),sopt), x[i])
    end

    return xp
end

project!(x, feasset::ConstrainedByShapeSaturated) = project!(x, x, feasset)

ext: jl


text:
using Pkg
using TOML
Pkg.add("PyPlot")
Pkg.add(url = "https://github.com/russelljjarvis/SpikeSynchrony.jl")
Pkg.add(url="https://github.com/paulmthompson/SpikeSorting.jl.git")
Pkg.add(url="https://github.com/paulmthompson/Spikes.jl.git")
Pkg.add("Conda")
Pkg.add("PyCall")
run(`$(which pip) install -networkunit`)
valid_test_folder = joinpath(@__DIR__, "testfiles", "valid")

tml = TOML.tryparsefile("../Project.toml")


text:
module TinyTimeModels

using Optim

export fit

include("kalman.jl")
include("parameterestimation.jl")

end # module


text:
using Plots

f(x) = 100*x - 60
g(x) = 30 - 60*x
h(x) = 0
v(x) = max.(f(x),g(x),h(x))
x = 0:0.1:1

plot(x, f.(x),
    xaxis = ("prob. (mu)", [0.0, 0.5, 0.6, 1.0]),
    yaxis = ("expected value", (-10,40)),
    label = "outdoor")
plot!(x, g.(x),
    label = "indoor")
plot!(x, h.(x),
    label = "none")
plot!(x, v.(x),
    label = "max{o,i,n}")

max_stars_repo_name: pedropessoa/git-intro

max_issues_repo_name: pedropessoa/git-intro

max_forks_repo_name: pedropessoa/git-intro


text:
using Statistics: median

test = """
    [({(<(())[]>[[{[]{<()<>>
    [(()[<>])]({[<{<<[]>>(
    {([(<{}[<>[]}>{[]{[(<()>
    (((({<>}<{<{<>}{[]{[]{}
    [[<[([]))<([[{}[[()]]]
    [{[{({}]{}}([{[{{{}}([]
    {<[[]]>}<{[{[{[]{()[[[]
    [<(<(<(<{}))><([]([]()
    <{([([[(<>()){}]>(<<{{
    <{([{{}}[<[[[<>{}]]]>[]]
    """

const PAIRS = Dict('(' => ')', '[' => ']', '{' => '}', '<' => '>')
const CORRUPT_SCORES = Dict(')' => 3, ']' => 57, '}' => 1197, '>' => 25137)
const MISSING_SCORES = Dict(')' => 1, ']' => 2, '}' => 3, '>' => 4)

function parse_input(input::AbstractString)::Vector{String}
    return split(input, r"\s+"; keepempty = false)
end

function score(line::AbstractString)
    open = []
    for c in line
        if c in keys(PAIRS)
            push!(open, PAIRS[c])
        else
            pop!(open) == c || return (CORRUPT_SCORES[c], 0)
        end
    end
    s = foldr((c, s) -> 5s + MISSING_SCORES[c], open; init = 0)
    return (0, s)
end

q1(lines::AbstractVector)::Int = lines .|> score .|> first |> sum

function q2(lines::AbstractVector)::Int
    return filter(!iszero, last.(score.(lines))) |> median |> Int
end

let lines = parse_input(test)
    @assert q1(lines) == 26397
    @assert q2(lines) == 288957
end

lines = parse_input(read("day10.in", String))
println("Q1: ", q1(lines))
println("Q2: ", q2(lines))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
module TestSynchronous

#println("\nSynchronousExamples: Demonstrating the ability to simulate models with synchronous semantics")

using TinyModia
@usingModiaPlot
using DifferentialEquations

MassWithSpringDamper = Model(
    m = 1.0,
    k = 1.0,           # spring constant
    d = 0.1,           # damping coefficient
    x = Var(init=0.0), # Position
    v = Var(init=5.0), # Velocity
    equations1 = :[
        der(x) = v
        m*der(v) = f - k*x - d*v
    ]
)

SpeedControl = MassWithSpringDamper | Model(
    k    = 0.0,
    K    = 5.0,   # Gain of speed P controller
    vref = 100.0, # Speed ref.
    vd = Var(start=0.0),
    u  = Var(start=0.0),
    equations = :[
        vd = sample(v, Clock(0.1))    # speed sensor
        u  = K*(vref-vd)              # P controller for speed
        f  = hold(u, Clock(0.1))      # force actuator
    ]
)

speedControl = @instantiateModel(SpeedControl)
simulate!(speedControl, Tsit5(), stopTime=1.5, logEvents=false, requiredFinalStates = [133.3845202465569, 98.03694943140056] )
plot(speedControl, [("v", "vd"), ("u","f")], heading="SpeedControl", figure=1)

SpeedControlPI = MassWithSpringDamper | Model(
    k    = 0.0,
    K    = 5.0,   # Gain of speed P controller
	Ti = 2, 	  # Integral time
    vref = 100.0, # Speed ref.
	dt = 0.1,     # sampling interval
    vd = Var(start=0.0),
    u  = Var(start=0.0),
	e = Var(start=0.0),
	intE = Var(start=10.0),  # To check that previous(intE) is corrected initialized
    equations = :[
		vd = sample(v, Clock(0.1))     # speed sensor

		# PI controller for speed
		e = vref-vd
        previous_intE = previous(intE, Clock(0.1))
		intE = previous_intE + e
		u = K*(e + intE/Ti)

		# force actuator
		f = hold(u, Clock(0.1))
    ]
)

speedControlPI = @instantiateModel(SpeedControlPI)
simulate!(speedControlPI, Tsit5(), stopTime=1.5, logEvents=false, requiredFinalStates = [145.828878491462, 99.95723678839984])
plot(speedControlPI, [("v", "vd"), ("u","f"), ("previous_intE", "intE")], heading="SpeedControlPI", figure=2)

end


text:
param = DriftScan(
m = m_He3,
theta = 45.75,
D = 5*mel,
h = 0.1,
xi = 0.38/A,
a = 4.3*A,
energies = range(4*mel,5*mel,length=10),
z = (-1*A,6*A)
)

@test_nowarn driftscan(param)


text:
using OrderedCollections
using StaticArrays
using Test

using CLIMAParameters
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

using ClimateMachine
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.Grids
using ClimateMachine.MPIStateArrays
using ClimateMachine.VariableTemplates
using ClimateMachine.SingleStackUtils
using ClimateMachine.DGMethods: LocalGeometry

using ClimateMachine.BalanceLaws:
    BalanceLaw, Auxiliary, Prognostic, Gradient, GradientFlux
import ClimateMachine.BalanceLaws:
    vars_state, nodal_init_state_auxiliary!, init_state_prognostic!

struct EmptyBalLaw{FT, PS} <: BalanceLaw
    "Parameters"
    param_set::PS
    "Domain height"
    zmax::FT

end
EmptyBalLaw(param_set, zmax) =
    EmptyBalLaw{typeof(zmax), typeof(param_set)}(param_set, zmax)

vars_state(::EmptyBalLaw, ::Auxiliary, FT) = @vars(x::FT, y::FT, z::FT)
vars_state(::EmptyBalLaw, ::Prognostic, FT) = @vars(ρ::FT)
vars_state(::EmptyBalLaw, ::Gradient, FT) = @vars()
vars_state(::EmptyBalLaw, ::GradientFlux, FT) = @vars()

function nodal_init_state_auxiliary!(
    m::EmptyBalLaw,
    aux::Vars,
    tmp::Vars,
    geom::LocalGeometry,
)
    aux.x = geom.coord[1]
    aux.y = geom.coord[2]
    aux.z = geom.coord[3]
end

function init_state_prognostic!(
    m::EmptyBalLaw,
    state::Vars,
    aux::Vars,
    coords,
    t::Real,
)
    z = aux.z
    x = aux.x
    y = aux.y
    state.ρ = (1 - 4 * (z - m.zmax / 2)^2) * (2 - x - y)
end

function test_hmean(
    grid::DiscontinuousSpectralElementGrid{T, dim, N},
    Q::MPIStateArray,
    vars,
) where {T, dim, N}
    state_vars_avg = get_horizontal_mean(grid, Q, vars)
    target = target_meanprof(grid)
    @test state_vars_avg["ρ"] ≈ target
end

function test_hvar(
    grid::DiscontinuousSpectralElementGrid{T, dim, N},
    Q::MPIStateArray,
    vars,
) where {T, dim, N}
    state_vars_var = get_horizontal_variance(grid, Q, vars)
    target = target_varprof(grid)
    @test state_vars_var["ρ"] ≈ target
end

function test_horizontally_ave(
    grid::DiscontinuousSpectralElementGrid{T, dim, N},
    Q_in::MPIStateArray,
    vars,
) where {T, dim, N}
    Q = deepcopy(Q_in)
    state_vars_var = get_horizontal_variance(grid, Q, vars)
    i_vars = varsindex(vars, :ρ)
    horizontally_average!(grid, Q, i_vars)
    FT = eltype(Q)
    state_vars_var = get_horizontal_variance(grid, Q, vars)
    @test all(isapprox.(state_vars_var["ρ"], 0, atol = 10 * eps(FT)))
end

function target_meanprof(
    grid::DiscontinuousSpectralElementGrid{T, dim, N},
) where {T, dim, N}
    Nq = N + 1
    Ntot = Nq * grid.topology.stacksize
    z = Array(get_z(grid))
    target =
        SVector{Ntot, T}([1.0 - 4.0 * (z_i - z[Ntot] / 2.0)^2 for z_i in z])
    return target
end

function target_varprof(
    grid::DiscontinuousSpectralElementGrid{T, dim, N},
) where {T, dim, N}
    Nq = N + 1
    nvertelem = grid.topology.stacksize
    Ntot = Nq * nvertelem
    z = Array(get_z(grid))
    x = z[1:Nq] * nvertelem
    scaled_var = 0.0
    for i in 1:Nq
        for j in 1:Nq
            scaled_var = scaled_var + (2 - x[i] - x[j]) * (2 - x[i] - x[j])
        end
    end
    target = SVector{Ntot, Float64}([
        (1.0 - 4.0 * (z_i - z[Ntot] / 2.0)^2) *
        (1.0 - 4.0 * (z_i - z[Ntot] / 2.0)^2) *
        (scaled_var / Nq / Nq - 1) for z_i in z
    ])
    return target
end

function main()
    FT = Float64
    ClimateMachine.init()

    m = EmptyBalLaw(param_set, FT(1))

    # Prescribe polynomial order of basis functions in finite elements
    N_poly = 5
    # Specify the number of vertical elements
    nelem_vert = 20
    # Specify the domain height
    zmax = m.zmax
    # Initial and final times
    t0 = 0.0
    timeend = 1.0
    dt = 0.1
    # Establish a `ClimateMachine` single stack configuration
    driver_config = ClimateMachine.SingleStackConfiguration(
        "SingleStackUtilsTest",
        N_poly,
        nelem_vert,
        zmax,
        param_set,
        m,
    )
    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        ode_dt = dt,
    )

    # tests
    test_hmean(
        driver_config.grid,
        solver_config.Q,
        vars_state(m, Prognostic(), FT),
    )
    test_hvar(
        driver_config.grid,
        solver_config.Q,
        vars_state(m, Prognostic(), FT),
    )

    r1, z1 = reduce_nodal_stack(
        max,
        solver_config.dg.grid,
        solver_config.Q,
        vars_state(m, Prognostic(), FT),
        "ρ",
        i = 6,
        j = 6,
    )
    @test r1 ≈ 8.880558532968455e-16 && z1 == 10
    r2, z2 = reduce_nodal_stack(
        +,
        solver_config.dg.grid,
        solver_config.Q,
        vars_state(m, Prognostic(), FT),
        "ρ",
        i = 3,
        j = 3,
    )
    @test r2 ≈ 102.73283921735293 && z2 == 20
    ns = reduce_element_stack(
        +,
        solver_config.dg.grid,
        solver_config.Q,
        vars_state(m, Prognostic(), FT),
        "ρ",
    )
    (r3, z3) = let
        f(a, b) = (a[1] + b[1], b[2])
        reduce(f, ns)
    end
    @test r3 ≈ FT(2877.6) && z3 == 20

    test_horizontally_ave(
        driver_config.grid,
        solver_config.Q,
        vars_state(m, Prognostic(), FT),
    )

end
main()

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
##  Jacobi Operator

recA(::Type{T},S::Ultraspherical,k) where {T} = (2*(k+order(S)))/(k+one(T))   # one(T) ensures we get correct type
recB(::Type{T},::Ultraspherical,_) where {T} = zero(T)
recC(::Type{T},S::Ultraspherical,k) where {T} = (k-one(T)+2order(S))/(k+one(T))   # one(T) ensures we get correct type

# x p_k
recα(::Type{T},::Ultraspherical,_) where {T} = zero(T)
recβ(::Type{T},S::Ultraspherical,k) where {T} = k/(2*(k-one(T)+order(S)))   # one(T) ensures we get correct type
recγ(::Type{T},S::Ultraspherical,k) where {T} = (k-2+2order(S))/(2*(k-one(T)+order(S)))   # one(T) ensures we get correct type

## Multiplication
# these are special cases

Base.stride(M::ConcreteMultiplication{U,V}) where {U<:Chebyshev,V<:Ultraspherical} =
    stride(M.f)
Base.stride(M::ConcreteMultiplication{U,V}) where {U<:Ultraspherical,V<:Chebyshev} =
    stride(M.f)
Base.stride(M::ConcreteMultiplication{U,V}) where {U<:Ultraspherical,V<:Ultraspherical} =
    stride(M.f)

function Multiplication(f::Fun{C},sp::Ultraspherical{Int}) where C<:Chebyshev
    if order(sp) == 1
        cfs = f.coefficients
        MultiplicationWrapper(f,
            SpaceOperator(SymToeplitzOperator(cfs/2) +
                                HankelOperator(view(cfs,3:length(cfs))/(-2)),
                          sp,sp))

    else
        ConcreteMultiplication(f,sp)
    end
end

## Derivative

#Derivative(k::Integer,d::IntervalOrSegment)=Derivative(k-1:k,d)
#Derivative(d::IntervalOrSegment)=Derivative(1,d)

Derivative(sp::Ultraspherical{LT,DD},m::Integer) where {LT,DD<:IntervalOrSegment} =
    ConcreteDerivative(sp,m)
function Integral(sp::Ultraspherical{LT,DD},m::Integer) where {LT,DD<:IntervalOrSegment}
    λ = order(sp)
    if m ≤ λ
        ConcreteIntegral(sp,m)
    else # Convert up
        nsp = Ultraspherical(λ+1,domain(sp))
        IntegralWrapper(Integral(nsp,m)*Conversion(sp,nsp),m)
    end
end

rangespace(D::ConcreteDerivative{Ultraspherical{LT,DD,RR}}) where {LT,DD<:IntervalOrSegment,RR} =
    Ultraspherical(order(domainspace(D))+D.order,domain(D))

bandwidths(D::ConcreteDerivative{Ultraspherical{LT,DD,RR}}) where {LT,DD<:IntervalOrSegment,RR} = -D.order,D.order
bandwidths(D::ConcreteIntegral{Ultraspherical{LT,DD,RR}}) where {LT,DD<:IntervalOrSegment,RR} = D.order,-D.order
Base.stride(D::ConcreteDerivative{Ultraspherical{LT,DD,RR}}) where {LT,DD<:IntervalOrSegment,RR} = D.order

function getindex(D::ConcreteDerivative{Ultraspherical{TT,DD,RR},K,T},
               k::Integer,j::Integer) where {TT,DD<:IntervalOrSegment,RR,K,T}
    m=D.order
    d=domain(D)
    λ=order(domainspace(D))

    if j==k+m
        convert(T,(pochhammer(one(T)*λ,m)*(4/complexlength(d)).^m))
    else
        zero(T)
    end
end

## Integral

linesum(f::Fun{Ultraspherical{LT,DD,RR}}) where {LT,DD<:IntervalOrSegment,RR} =
    sum(setcanonicaldomain(f))*arclength(d)/2

rangespace(D::ConcreteIntegral{Ultraspherical{LT,DD,RR}}) where {LT,DD<:IntervalOrSegment,RR} =
    order(domainspace(D)) == 1 ? Chebyshev(domain(D)) : Ultraspherical(order(domainspace(D))-D.order,domain(D))

function getindex(Q::ConcreteIntegral{Ultraspherical{LT,DD,RR}},k::Integer,j::Integer) where {LT,DD<:IntervalOrSegment,RR}
    T=eltype(Q)
    m=Q.order
    d=domain(Q)
    λ=order(domainspace(Q))
    @assert m<=λ

    if λ == 1 && k==j+1
        C = complexlength(d)/2
        convert(T,C./(k-1))
    elseif λ > 1 && k==j+m
        convert(T,pochhammer(one(T)*λ,-m)*(complexlength(d)/4)^m)
    else
        zero(T)
    end
end

## Conversion Operator

function Conversion(A::Chebyshev,B::Ultraspherical)
    if order(B) ≤ 1
        ConcreteConversion(A,B)
    else
        d=domain(A)
        US=Ultraspherical(order(B)-1,d)
        ConversionWrapper(TimesOperator(Conversion(US,B),
                                        Conversion(Chebyshev(d),US)))
    end
end

function Conversion(A::Ultraspherical,B::Chebyshev)
    if order(A) == 1//2
        ConcreteConversion(A,B)
    else
        error("Not implemented")
    end
end

maxspace_rule(A::Ultraspherical,B::Chebyshev) = A

function Conversion(A::Ultraspherical,B::Ultraspherical)
    a=order(A); b=order(B)
    if b==a
        ConversionWrapper(Operator(I,A))
    elseif a<b≤a+1  || b<a≤b+1
        ConcreteConversion(A,B)
    elseif b ≠ 1
        d=domain(A)
        ConversionWrapper(TimesOperator(Conversion(Ultraspherical(b-1,d),B),
                                        Conversion(A,Ultraspherical(b-1,d))))
    else
        error("Cannot convert from $A to $B")
    end
end

maxspace_rule(A::Ultraspherical,B::Ultraspherical) = order(A) > order(B) ? A : B

function getindex(M::ConcreteConversion{C,Ultraspherical{Int,DD,RR},T},
               k::Integer,j::Integer) where {DD,RR,C<:Chebyshev,T}
   # order must be 1
    if k==j==1
        one(T)
    elseif k==j
        one(T)/2
    elseif j==k+2
        -one(T)/2
    else
        zero(T)
    end
end

function getindex(M::ConcreteConversion{Ultraspherical{Int,DD,RR},Ultraspherical{Int,DD,RR},T},
                   k::Integer,j::Integer) where {DD,RR,T}
    #  we can assume that λ==m+1
    λ=order(rangespace(M))
    c=λ-one(T)  # this supports big types
    if k==j
        c/(k - 2 + λ)
    elseif j==k+2
        -c/(k + λ)
    else
        zero(T)
    end
end

function getindex(M::ConcreteConversion{Ultraspherical{LT,DD,RR},Ultraspherical{LT,DD,RR},T},
                k::Integer,j::Integer) where {LT,DD,RR,T}
    λ=order(rangespace(M))
    if order(domainspace(M))+1==λ
        c=λ-one(T)  # this supports big types
        if k==j
            c/(k - 2 + λ)
        elseif j==k+2
            -c/(k + λ)
        else
            zero(T)
        end
    else
        error("Not implemented")
    end
end

bandwidths(C::ConcreteConversion{<:Chebyshev,<:Ultraspherical{Int}}) = 0,2  # order == 1
bandwidths(C::ConcreteConversion{<:Ultraspherical{Int},<:Ultraspherical{Int}}) = 0,2

bandwidths(C::ConcreteConversion{<:Chebyshev,<:Ultraspherical}) =
    0,order(rangespace(C))==1 ? 2 : ∞
bandwidths(C::ConcreteConversion{<:Ultraspherical,<:Chebyshev}) =
    0,order(domainspace(C))==1 ? 2 : ∞

bandwidths(C::ConcreteConversion{<:Ultraspherical,<:Ultraspherical}) =
    0,order(domainspace(C))+1==order(rangespace(C)) ? 2 : ∞

Base.stride(C::ConcreteConversion{<:Chebyshev,<:Ultraspherical{Int}}) = 2
Base.stride(C::ConcreteConversion{<:Ultraspherical,<:Ultraspherical}) = 2

## coefficients

# return the space that has banded Conversion to the other
conversion_rule(a::Chebyshev,b::Ultraspherical{Int}) =
    if domainscompatible(a,b)
        a
    else
        NoSpace()
    end

conversion_rule(a::Ultraspherical{LT},b::Ultraspherical{LT}) where {LT} =
    if domainscompatible(a,b) && isapproxinteger(order(a)-order(b))
        order(a) < order(b) ? a : b
    else
        NoSpace()
    end

function coefficients(g::AbstractVector,sp::Ultraspherical{Int},C::Chebyshev)
    if order(sp) == 1
        ultraiconversion(g)
    else
        # do one at a time
        coefficients(g,sp,Ultraspherical(1,domain(sp)),C)
    end
end
function coefficients(g::AbstractVector,C::Chebyshev,sp::Ultraspherical)
    if order(sp) == 1
        ultraconversion(g)
    else
        # do one at a time
        coefficients(g,C,Ultraspherical(1,domain(sp)),sp)
    end
end

# TODO: include in getindex to speed up
Integral(sp::Chebyshev{DD},m::Integer) where {DD<:IntervalOrSegment} =
    IntegralWrapper(TimesOperator([Integral(Ultraspherical(m,domain(sp)),m),
                                   Conversion(sp,Ultraspherical(m,domain(sp)))]),m)

## Non-banded conversions

function getindex(M::ConcreteConversion{C,Ultraspherical{LT,DD,RR},T},
            k::Integer,j::Integer) where {DD,RR,LT,C<:Chebyshev,T}
    λ = order(rangespace(M))
    if λ == 1
        if k==j==1
            one(T)
        elseif k==j
            one(T)/2
        elseif j==k+2
            -one(T)/2
        else
            zero(T)
        end
    elseif λ == 0.5
        # Cheb-to-Leg
        if j==k==1
            one(T)
        elseif j==k
            convert(T,sqrt(π)/(2FastTransforms.Λ(k-1)))
        elseif k < j && iseven(k-j)
            convert(T,-(j-1)*(k-0.5)*(FastTransforms.Λ((j-k-2)/2)/(j-k))*
                            (FastTransforms.Λ((j+k-3)/2)/(j+k-1)))
        else
            zero(T)
        end
    else
        error("Not implemented")
    end
end

function getindex(M::ConcreteConversion{Ultraspherical{LT,DD,RR},C,T},
            k::Integer,j::Integer) where {DD,RR,LT,C<:Chebyshev,T}
    λ = order(domainspace(M))
    if λ == 1
        # order must be 1
        if k==j==1
            one(T)
        elseif k==j
            one(T)/2
        elseif j==k+2
            -one(T)/2
        else
            zero(T)
        end
    elseif λ == 0.5
        if k==1 && isodd(j)
            convert(T,FastTransforms.Λ((j-1)/2)^2/π)
        elseif k ≤ j && iseven(k-j)
            convert(T,FastTransforms.Λ((j-k)/2)*FastTransforms.Λ((k+j-2)/2)*2/π)
        else
            zero(T)
        end
    else
        error("Not implemented")
    end
end

function getindex(M::ConcreteConversion{Ultraspherical{LT,DD,RR},
                                     Ultraspherical{LT2,DD,RR},T},
                     k::Integer,j::Integer) where {DD,RR,LT,LT2,T}
    λ1 = order(domainspace(M))
    λ2 = order(rangespace(M))
    if abs(λ1-λ2) < 1
        if j ≥ k && iseven(k-j)
            convert(T,(λ1 < λ2 && k ≠ j ? -1 : 1) *  # fix sign for lgamma
                exp(lgamma(λ2)+log(k-1+λ2)-lgamma(λ1)-lgamma(λ1-λ2) + lgamma((j-k)/2+λ1-λ2)-
                lgamma((j-k)/2+1)+lgamma((k+j-2)/2+λ1)-lgamma((k+j-2)/2+λ2+1)))
        else
            zero(T)
        end
    else
        error("Not implemented")
    end
end

ReverseOrientation(S::Ultraspherical) = ReverseOrientationWrapper(NegateEven(S,reverseorientation(S)))
Reverse(S::Ultraspherical) = ReverseWrapper(NegateEven(S,S))


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	must_run(EP::Model, inputs::Dict, CapacityReserveMargin::Int)

This function defines the constraints for operation of `must-run' or non-dispatchable resources, such as rooftop solar systems that do not receive dispatch signals, run-of-river hydroelectric facilities without the ability to spill water, or cogeneration systems that must produce a fixed quantity of heat in each time step. This resource type can also be used to model baseloaded or self-committed thermal generators that do not respond to economic dispatch.

For must-run resources ($y\in \mathcal{MR}$) output in each time period $t$ must exactly equal the available capacity factor times the installed capacity, not allowing for curtailment. These resources are also not eligible for contributing to frequency regulation or operating reserve requirements.

```math
\begin{aligned}
\Theta_{y,z,t} = \rho^{max}_{y,z,t}\times \Delta^{total}_{y,z}
\hspace{4 cm}  \forall y \in \mathcal{MR}, z \in \mathcal{Z},t \in \mathcal{T}
\end{aligned}
```
"""
function must_run(EP::Model, inputs::Dict, CapacityReserveMargin::Int)

	println("Must-Run Resources Module")

	dfGen = inputs["dfGen"]

	T = inputs["T"]     # Number of time steps (hours)
	Z = inputs["Z"]     # Number of zones
	G = inputs["G"] 	# Number of generators

	MUST_RUN = inputs["MUST_RUN"]

	### Expressions ###

	## Power Balance Expressions ##

	@expression(EP, ePowerBalanceNdisp[t=1:T, z=1:Z],
		sum(EP[:vP][y,t] for y in intersect(MUST_RUN, dfGen[dfGen[!,:Zone].==z,:][!,:R_ID])))

	EP[:ePowerBalance] += ePowerBalanceNdisp

	# Capacity Reserves Margin policy
	if CapacityReserveMargin > 0
		@expression(EP, eCapResMarBalanceMustRun[res=1:inputs["NCapacityReserveMargin"], t=1:T], sum(dfGen[y,Symbol("CapRes_$res")] * EP[:eTotalCap][y] * inputs["pP_Max"][y,t]  for y in MUST_RUN))
		EP[:eCapResMarBalance] += eCapResMarBalanceMustRun
	end

	### Constratints ###

	@constraint(EP, [y in MUST_RUN, t=1:T], EP[:vP][y,t] == inputs["pP_Max"][y,t]*EP[:eTotalCap][y])
	## CO2 Policy Module Must Run Generation by zone
	@expression(EP, eGenerationByMustRun[z=1:Z, t=1:T], # the unit is GW
		sum(EP[:vP][y,t] for y in intersect(inputs["MUST_RUN"], dfGen[dfGen[!,:Zone].==z,:R_ID]))
	)
	EP[:eGenerationByZone] += eGenerationByMustRun

	return EP
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
#########################################################################
# PyDict: no-copy wrapping of a Julia object around a Python dictionary

# we check for "items" attr since PyMapping_Check doesn't do this (it only
# checks for __getitem__) and PyMapping_Check returns true for some
# scipy scalar array members, grrr.
function is_mapping_object(o::PyObject)
    pyisinstance(o, @pyglobalobj :PyDict_Type) ||
    (pyquery((@pyglobal :PyMapping_Check), o) &&
      ccall((@pysym :PyObject_HasAttrString), Cint, (PyPtr,Ptr{UInt8}), o, "items") == 1)
end

"""
    PyDict(o::PyObject)
    PyDict(d::Dict{K,V})

This returns a PyDict, which is a no-copy wrapper around a Python dictionary.

Alternatively, you can specify the return type of a `pycall` as PyDict.
"""
mutable struct PyDict{K,V,isdict} <: AbstractDict{K,V}
    o::PyObject
    # isdict = true for python dict, otherwise is a generic Mapping object

    function PyDict{K,V,isdict}(o::PyObject) where {K,V,isdict}
        if !isdict && !ispynull(o) && !is_mapping_object(o)
            throw(ArgumentError("only Dict and Mapping objects can be converted to PyDict"))
        end
        return new{K,V,isdict}(o)
    end
end

PyDict{K,V}(o::PyObject) where {K,V} = PyDict{K,V,pyisinstance(o, @pyglobalobj :PyDict_Type)}(o)
PyDict{K,V}() where {K,V} = PyDict{K,V,true}(PyObject(@pycheckn ccall((@pysym :PyDict_New), PyPtr, ())))

PyDict(o::PyObject) = PyDict{PyAny,PyAny}(o)
PyObject(d::PyDict) = d.o
PyDict() = PyDict{PyAny,PyAny}()
PyDict(d::AbstractDict{K,V}) where {K,V} = PyDict{K,V}(PyObject(d))
PyDict(d::AbstractDict{Any,Any}) = PyDict{PyAny,PyAny}(PyObject(d))
PyDict(d::AbstractDict{Any,V}) where {V} = PyDict{PyAny,V}(PyObject(d))
PyDict(d::AbstractDict{K,Any}) where {K} = PyDict{K,PyAny}(PyObject(d))
convert(::Type{PyDict}, o::PyObject) = PyDict(o)
convert(::Type{PyDict{K,V}}, o::PyObject) where {K,V} = PyDict{K,V}(o)
unsafe_convert(::Type{PyPtr}, d::PyDict) = PyPtr(d.o)

haskey(d::PyDict{K,V,true}, key) where {K,V} = 1 == ccall(@pysym(:PyDict_Contains), Cint, (PyPtr, PyPtr), d, PyObject(key))
keys(::Type{T}, d::PyDict{K,V,true}) where {T,K,V} = convert(Vector{T}, PyObject(@pycheckn ccall((@pysym :PyDict_Keys), PyPtr, (PyPtr,), d)))
values(::Type{T}, d::PyDict{K,V,true}) where {T,K,V} = convert(Vector{T}, PyObject(@pycheckn ccall((@pysym :PyDict_Values), PyPtr, (PyPtr,), d)))

keys(::Type{T}, d::PyDict{K,V,false}) where {T,K,V} = convert(Vector{T}, pycall(d.o["keys"], PyObject))
values(::Type{T}, d::PyDict{K,V,false}) where {T,K,V} = convert(Vector{T}, pycall(d.o["values"], PyObject))
haskey(d::PyDict{K,V,false}, key) where {K,V} = 1 == ccall(@pysym(:PyMapping_HasKey), Cint, (PyPtr, PyPtr), d, PyObject(key))

similar(d::PyDict{K,V}) where {K,V} = Dict{pyany_toany(K),pyany_toany(V)}()
eltype(::Type{PyDict{K,V}}) where {K,V} = Pair{pyany_toany(K),pyany_toany(V)}
Base.keytype(::PyDict{K,V}) where {K,V} = pyany_toany(K)
Base.valtype(::PyDict{K,V}) where {K,V} = pyany_toany(V)
Base.keytype(::Type{PyDict{K,V}}) where {K,V} = pyany_toany(K)
Base.valtype(::Type{PyDict{K,V}}) where {K,V} = pyany_toany(V)

function setindex!(d::PyDict, v, k)
    @pycheckz ccall((@pysym :PyObject_SetItem), Cint, (PyPtr, PyPtr, PyPtr),
                     d, PyObject(k), PyObject(v))
    v
end

get(d::PyDict{K,V}, k, default) where {K,V} = get(d.o, V, k, default)

function pop!(d::PyDict{K,V,true}, k) where {K,V}
    v = d[k]
    @pycheckz ccall(@pysym(:PyDict_DelItem), Cint, (PyPtr, PyPtr), d, PyObject(k))
    return v
end
function pop!(d::PyDict{K,V,false}, k) where {K,V}
    v = d[k]
    @pycheckz ccall(@pysym(:PyObject_DelItem), Cint, (PyPtr, PyPtr), d, PyObject(k))
    return v
end

function pop!(d::PyDict, k, default)
    try
        return pop!(d, k)
    catch
        return default
    end
end

function delete!(d::PyDict{K,V,true}, k) where {K,V}
    e = ccall(@pysym(:PyDict_DelItem), Cint, (PyPtr, PyPtr), d, PyObject(k))
    e == -1 && pyerr_clear() # delete! ignores errors in Julia
    return d
end
function delete!(d::PyDict{K,V,false}, k) where {K,V}
    e = ccall(@pysym(:PyObject_DelItem), Cint, (PyPtr, PyPtr), d, PyObject(k))
    e == -1 && pyerr_clear() # delete! ignores errors in Julia
    return d
end

function empty!(d::PyDict{K,V,true}) where {K,V}
    @pycheck ccall((@pysym :PyDict_Clear), Cvoid, (PyPtr,), d)
    return d
end
function empty!(d::PyDict{K,V,false}) where {K,V}
    # for generic Mapping items we must delete keys one by one
    for k in keys(d)
        delete!(d, k)
    end
    return d
end

length(d::PyDict{K,V,true}) where {K,V} = @pycheckz ccall(@pysym(:PyDict_Size), Int, (PyPtr,), d)
length(d::PyDict{K,V,false}) where {K,V} = @pycheckz ccall(@pysym(:PyObject_Size), Int, (PyPtr,), d)
isempty(d::PyDict) = length(d) == 0

struct PyDict_Iterator
    # arrays to pass key, value, and pos pointers to PyDict_Next
    ka::Ref{PyPtr}
    va::Ref{PyPtr}
    pa::Ref{Int}
    i::Int # current position in items list (0-based)
    len::Int # length of items list
end

function Base.iterate(d::PyDict{K,V,true}, itr=PyDict_Iterator(Ref{PyPtr}(), Ref{PyPtr}(), Ref(0), 0, length(d))) where {K,V}
    itr.i >= itr.len && return nothing
    if 0 == ccall((@pysym :PyDict_Next), Cint,
                    (PyPtr, Ref{Int}, Ref{PyPtr}, Ref{PyPtr}),
                    d, itr.pa, itr.ka, itr.va)
        error("unexpected end of PyDict_Next")
    end
    ko = pyincref(itr.ka[]) # PyDict_Next returns
    vo = pyincref(itr.va[]) #   borrowed ref, so incref
    (Pair(convert(K,ko), convert(V,vo)),
        PyDict_Iterator(itr.ka, itr.va, itr.pa, itr.i+1, itr.len))
end

# Iterator for generic mapping, using Python items iterator.
# Our approach is to wrap an iterator over d.o["items"]
# which necessitates including d.o["items"] in the state.
function _start(d::PyDict{K,V,false}) where {K,V}
    d_items = pycall(d.o."items", PyObject)
    (d_items, iterate(d_items))
end
function Base.iterate(d::PyDict{K,V,false}, itr=_start(d)) where {K,V}
    d_items, iter_result = itr
    iter_result === nothing && return nothing
    item, state = iter_result
    iter_result = iterate(d_items, state)
    (item[1] => item[2], (d_items, iter_result))
end


text:
using Test
using LogicCircuits
using ProbabilisticCircuits
using DataFrames
using CUDA

@testset "HCLT construction" begin
    num_vars = 3
    num_cats = 2

    data = DataFrame([true true false; false true false; false false false], :auto)

    pc = hclt(num_vars, num_cats; data, num_hidden_cats = 4)

    @test pc isa ProbCircuit
    @test num_children(pc) == 4
    @test pc.children[1].children[2].children[1].literal == Lit(1)
    @test pc.children[1].children[2].children[2].literal == Lit(-1)
end


text:
import Chipmunk
const Cp = Chipmunk

abstract type Object{T} <: AbstractVector{T} end
Base.IndexStyle(::Type{<:Object}) = IndexLinear()
Base.size(p::Object) = (2,)
Base.length(p::Object) = 2
function Base.getindex(p::Object, i::Int)
    if i == 1
        return Cp.get_position(p.body).x # TODO: needs too long # TODO: tooks too many memory
    elseif i == 2
        return Cp.get_position(p.body).y # TODO: needs too long # TODO: tooks too many memory
    end
    return 0.0
end
# function Base.setindex!(p::Object, v, i::Int)
#     if i == 1
#         Cp.set_position(p.body, Cp.Vect(v, p[2]))
#     elseif i == 2
#         Cp.set_position(p.body, Cp.Vect(p[1], v))
#     end
# end
get_position(p::Object) = [Cp.get_position(p.body).x, Cp.get_position(p.body).y]
set_position!(p::Object, new_position::AbstractVector{<:Number}) = Cp.set_position(p.body, Cp.Vect(new_position...))
get_velocity(p::Object) = [Cp.get_velocity(p.body).x, Cp.get_velocity(p.body).y]
set_velocity!(p::Object, new_velocity::AbstractVector{<:Number}) = Cp.set_velocity(p.body, Cp.Vect(new_velocity...)) # TODO: tooks too many memory
set_velocity!(p::Object, new_velocity_x::Number, new_velocity_y::Number) = Cp.set_velocity(p.body, Cp.Vect(new_velocity_x, new_velocity_y)) # TODO: tooks too many memory
get_angular_velocity(p::Object) = Cp.get_angular_velocity(p.body)
get_size(p::Object) = Cp.get_radius(p.shape)
get_mass(p::Object) = Cp.get_mass(p.shape)
get_friction(p::Object) = Cp.get_friction(p.shape)

"""
    Robot

The current state of an robot.

# Arguments

position::AbstractVector - 2d Vector with the position in the plane
body_orientation::T      - angle (in radiens) relative to x axis with [-π, π], determined from `atan2`
"""
mutable struct Robot{T<:Number} <: Object{T}
    # position::AbstractVector{T}
    body_orientation::T
    head_orientation::T
    # velocity::T
    # angularvelocity::T
    # size::T
    # mass::T
    hfov::T
    max_velocity::T
    max_angularvelocity::T
    body::Cp.Body
    shape::Cp.Shape
    order::Order
end
function Robot(position::AbstractVector{T}, body_orientation::T, head_orientation::T, velocity::AbstractVector{T}, angularvelocity::T, size::T, mass::T, hfov::T, max_velocity::T, max_angularvelocity::T, order::Order) where T
    body = Cp.Body(0, 0) # TODO: KinematicBody
    shape = Cp.CircleShape(body, size, Cp.Vect(0, 0))
    Cp.set_position(body, Cp.Vect(position...))
    Cp.set_friction(shape, 0.8) #TODO: add command line argument
    Cp.set_elasticity(shape, 0) #TODO: add command line argument
    Cp.set_mass(shape, mass)
    Cp.set_collision_type(shape, 0)
    Cp.set_velocity(body, Cp.Vect(velocity...))
    Cp.set_angular_velocity(body, angularvelocity)
    Robot(body_orientation, head_orientation, hfov, max_velocity, max_angularvelocity, body, shape, order)
end
# http://doc.aldebaran.com/2-1/family/robots/index_robots.html
Robot(;position::AbstractVector{T}=[0.0, 0.0],
       body_orientation::T=0.0,
       head_orientation::T=0.0,
       velocity::AbstractVector{T}=[0.0, 0.0],
       angularvelocity::T=0.0,
       size::T=0.15,
       mass::T=5.305350006,
       hfov::T=deg2rad(60.97),
       max_velocity::T=1.0,
       max_angularvelocity::T=Float64(π/2),
       order::Order=NoOrder()) where T<:Number =
    Robot(position, body_orientation, head_orientation, velocity, angularvelocity, size, mass, hfov, max_velocity, max_angularvelocity, order)
Robot(r::Robot) = Robot(get_position(r), r.body_orientation, r.head_orientation, get_velocity(r), get_angular_velocity(r), get_size(r), get_mass(r), r.hfov, r.max_velocity, r.max_angularvelocity, r.order)

Robot(position::AbstractVector{T}, body_orientation::T) where T<:Number =
    Robot(position=position, body_orientation=body_orientation)
NaoV5(position::AbstractVector{T}, body_orientation::T) where T<:Number = Robot(position, body_orientation)
Nao(position::AbstractVector{T}, body_orientation::T) where T<:Number = NaoV5(position, body_orientation)

get_angularvelocity(r::Robot) = Cp.get_angular_velocity(r.body)

apply_order!(robot::Robot) = apply_order!(robot, robot.order)

"""
    Obstacle

The current state of an Obstacle.
"""
mutable struct Obstacle{T<:Number} <: Object{T}
    # position::AbstractVector{T}
    # velocity::T
    # size::T
    # mass::T
    max_velocity::T
    deceleration::T
    friction::T
    fixed::Bool
    body::Cp.Body
    shape::Cp.Shape
    color::AbstractVector{<:Real} # TODO: switch to Color Object
end
function Obstacle(position::AbstractVector{T}, velocity::AbstractVector{T}, size::T, mass::T, max_velocity::T, friction::T, gravity::T, fixed::Bool, color::AbstractVector{<:Real}) where T
    if fixed
        body = Cp.StaticBody()
    else
        body = Cp.Body(0, 0)
        Cp.set_velocity(body, Cp.Vect(velocity...))
    end
    Cp.set_position(body, Cp.Vect(position...))
    Cp.set_velocity(body, Cp.Vect(velocity...))
    shape = Cp.CircleShape(body, size, Cp.Vect(0, 0))
    Cp.set_friction(shape, friction) #TODO: add command line argument
    Cp.set_elasticity(shape, 1.0) #TODO: add command line argument
    Cp.set_collision_type(shape, 1)
    Cp.set_mass(shape, mass)
    Obstacle(max_velocity, friction*gravity, friction, fixed, body, shape, color)
end
Obstacle(;position::AbstractVector{T}=[0.0, 0.0],
          velocity::AbstractVector{T}=[0.0, 0.0],
          size::T=0.05,
          mass::T=0.048,
          max_velocity::T=0.25,
          friction::T=0.1,
          gravity::T=9.81,
          fixed::Bool=false,
          color::AbstractVector{<:Real}=[1.0, 0.1, 0.1, 1.0]) where T<:Number =
    Obstacle(position, velocity, size, mass, max_velocity, friction, gravity, fixed, color)
Obstacle(o::Obstacle) = Obstacle(get_position(o), get_velocity(o), get_size(o), get_mass(o), o.max_velocity, o.friction, o.deceleration/o.friction, o.fixed, o.color)
apply_order!(obstacle::Obstacle) = nothing

Ball(position::AbstractVector{T}) where T<:Number = Obstacle(position=position)

"""
    WorldState

Data clump holding all information about the game.
"""
struct WorldState{R<:AbstractVector{<:Robot{<:Number}}, O<:AbstractVector{<:Obstacle{<:Number}}}
    myPlayers::R # indices of myPlayers
    opponents::R # indices of opponents
    obstacles::O
    # time::Number
    # max_time::Number
    # v_off::Number
    # player::AbstractVector{<:Robot{<:T}}
    # all::AbstractVector{<:Object{<:T}}
end
# WorldState(myPlayers, opponents, obstacles) = WorldState(myPlayers, opponents, obstacles, vcat(myPlayers, opponents), vcat(myPlayers, opponents, obstacles))
function WorldState(ws::WorldState)
    myPlayers = Robot.(ws.myPlayers)
    opponents = Robot.(ws.opponents)
    obstacles = Obstacle.(ws.obstacles)
    return WorldState(myPlayers, opponents, obstacles)
end

# struct WorldStateIter
#     ws::WorldState
#     start::Integer

#     iterate_obstacles::Bool
#     exclude::Object
# end


text:
# ---
# title: JuliaRL\_A2CGAE\_CartPole
# cover: assets/JuliaRL_A2CGAE_CartPole.png
# description: A2CGAE applied to CartPole
# date: 2021-05-22
# author: "[Sriram](https://github.com/sriram13m)"
# ---

#+ tangle=true
using ReinforcementLearning
using StableRNGs
using Flux
using Flux.Losses

function RL.Experiment(
    ::Val{:JuliaRL},
    ::Val{:A2CGAE},
    ::Val{:CartPole},
    ::Nothing;
    seed = 123,
)
    rng = StableRNG(seed)
    N_ENV = 16
    UPDATE_FREQ = 10
    env = MultiThreadEnv([
        CartPoleEnv(; T = Float32, rng = StableRNG(hash(seed + i))) for i in 1:N_ENV
    ])
    ns, na = length(state(env[1])), length(action_space(env[1]))
    RLBase.reset!(env, is_force = true)
    agent = Agent(
        policy = QBasedPolicy(
            learner = A2CGAELearner(
                approximator = ActorCritic(
                    actor = NeuralNetworkApproximator(
                        model = Chain(
                            Dense(ns, 256, relu; init = glorot_uniform(rng)),
                            Dense(256, na; init = glorot_uniform(rng)),
                        ),
                        optimizer = ADAM(1e-3),
                    ),
                    critic = NeuralNetworkApproximator(
                        model = Chain(
                            Dense(ns, 256, relu; init = glorot_uniform(rng)),
                            Dense(256, 1; init = glorot_uniform(rng)),
                        ),
                        optimizer = ADAM(1e-3),
                    ),
                ) |> gpu,
                γ = 0.99f0,
                λ = 0.97f0,
                actor_loss_weight = 1.0f0,
                critic_loss_weight = 0.5f0,
                entropy_loss_weight = 0.001f0,
                update_freq = UPDATE_FREQ,
            ),
            explorer = BatchExplorer(GumbelSoftmaxExplorer(;)),
        ),
        trajectory = CircularArraySARTTrajectory(;
            capacity = UPDATE_FREQ,
            state = Matrix{Float32} => (ns, N_ENV),
            action = Vector{Int} => (N_ENV,),
            reward = Vector{Float32} => (N_ENV,),
            terminal = Vector{Bool} => (N_ENV,),
        ),
    )
    stop_condition = StopAfterStep(50_000, is_show_progress=!haskey(ENV, "CI"))
    hook = TotalBatchRewardPerEpisode(N_ENV)
    Experiment(agent, env, stop_condition, hook, "# A2CGAE with CartPole")
end

#+ tangle=false
using Plots
using Statistics
pyplot() #hide
ex = E`JuliaRL_A2CGAE_CartPole`
run(ex)
n = minimum(map(length, ex.hook.rewards))
m = mean([@view(x[1:n]) for x in ex.hook.rewards])
s = std([@view(x[1:n]) for x in ex.hook.rewards])
plot(m,ribbon=s)
savefig("assets/JuliaRL_A2CGAE_CartPole.png") #hide

# ![](assets/JuliaRL_A2CGAE_CartPole.png)


text:
# ------------ Runtime ------------ #

(a::Type{Array{K, N}})(init, d::Diffed) where {K, N} = a(init, unwrap(d))
(f::Colon)(v::K, d::Diffed{K, NoChange}) where K = f(v, unwrap(d))
isless(v::K, d::Diffed{K, NoChange}) where K = isless(v, unwrap(d))
Base.Pair(v, d::Diffed) = Base.Pair(v, unwrap(d))
Base.Pair(d::Diffed, v) = Base.Pair(unwrap(d), v)

# ------------ Abstract ------------ #

@abstract DiffInterpreter convert(args...) = propagate(args...)
@abstract DiffInterpreter iterate(args...) = propagate(args...)
@abstract DiffInterpreter getfield(args...) = propagate(args...)
@abstract DiffInterpreter typeassert(args...) = propagate(args...)
@abstract DiffInterpreter sleep(args...) = propagate(args...)
@abstract DiffInterpreter collect(args...) = propagate(args...)
@abstract DiffInterpreter (::Colon)(args...) = propagate(args...)
@abstract DiffInterpreter Base.:(===)(args...) = propagate(args...)
@abstract DiffInterpreter (::Core.IntrinsicFunction)(args...) = propagate(args...)
@abstract DiffInterpreter (::Type{Pair})(args...) = propagate(args...)


text:
# Pointers with address space information

#
# Address spaces
#

abstract type AddressSpace end

module AS

using AMDGPUnative
import AMDGPUnative: AddressSpace

struct Generic  <: AddressSpace end
struct Global   <: AddressSpace end
struct Region   <: AddressSpace end
struct Local    <: AddressSpace end
struct Constant <: AddressSpace end
struct Private  <: AddressSpace end

end

#
# Device pointer
#

struct DevicePtr{T,A}
    ptr::Ptr{T}

    # inner constructors, fully parameterized
    DevicePtr{T,A}(ptr::Ptr{T}) where {T,A<:AddressSpace} = new(ptr)
end

# outer constructors, partially parameterized
DevicePtr{T}(ptr::Ptr{T}) where {T} = DevicePtr{T,AS.Generic}(ptr)

# outer constructors, non-parameterized
DevicePtr(ptr::Ptr{T})              where {T} = DevicePtr{T,AS.Generic}(ptr)

Base.show(io::IO, dp::DevicePtr{T,AS}) where {T,AS} =
    print(io, AS.name.name, " Device", pointer(dp))

## getters

Base.pointer(p::DevicePtr) = p.ptr

Base.eltype(::Type{<:DevicePtr{T}}) where {T} = T

addrspace(x::DevicePtr) = addrspace(typeof(x))
addrspace(::Type{DevicePtr{T,A}}) where {T,A} = A

## conversions

# between regular and device pointers
## simple conversions disallowed
Base.convert(::Type{Ptr{T}}, p::DevicePtr{T})        where {T} = throw(InexactError(:convert, Ptr{T}, p))
Base.convert(::Type{<:DevicePtr{T}}, p::Ptr{T})      where {T} = throw(InexactError(:convert, DevicePtr{T}, p))
## unsafe ones are allowed
Base.unsafe_convert(::Type{Ptr{T}}, p::DevicePtr{T}) where {T} = pointer(p)

# defer conversions to DevicePtr to unsafe_convert
Base.cconvert(::Type{<:DevicePtr}, x) = x

# between device pointers
Base.convert(::Type{<:DevicePtr}, p::DevicePtr)                         = throw(InexactError(:convert, DevicePtr, p))
Base.convert(::Type{DevicePtr{T,A}}, p::DevicePtr{T,A})   where {T,A}   = p
Base.unsafe_convert(::Type{DevicePtr{T,A}}, p::DevicePtr) where {T,A}   = DevicePtr{T,A}(reinterpret(Ptr{T}, pointer(p)))
## identical addrspaces
Base.convert(::Type{DevicePtr{T,A}}, p::DevicePtr{U,A}) where {T,U,A} = Base.unsafe_convert(DevicePtr{T,A}, p)
## convert to & from generic
Base.convert(::Type{DevicePtr{T,AS.Generic}}, p::DevicePtr)               where {T}     = Base.unsafe_convert(DevicePtr{T,AS.Generic}, p)
Base.convert(::Type{DevicePtr{T,A}}, p::DevicePtr{U,AS.Generic})          where {T,U,A} = Base.unsafe_convert(DevicePtr{T,A}, p)
Base.convert(::Type{DevicePtr{T,AS.Generic}}, p::DevicePtr{T,AS.Generic}) where {T}     = p  # avoid ambiguities
## unspecified, preserve source addrspace
Base.convert(::Type{DevicePtr{T}}, p::DevicePtr{U,A}) where {T,U,A} = Base.unsafe_convert(DevicePtr{T,A}, p)

## limited pointer arithmetic & comparison

Base.:(==)(a::DevicePtr, b::DevicePtr) = pointer(a) == pointer(b) && addrspace(a) == addrspace(b)

Base.isless(x::DevicePtr, y::DevicePtr) = Base.isless(pointer(x), pointer(y))
Base.:(-)(x::DevicePtr, y::DevicePtr)   = pointer(x) - pointer(y)

Base.:(+)(x::DevicePtr{T,A}, y::Integer) where {T,A} = DevicePtr{T,A}(pointer(x) + y)
Base.:(-)(x::DevicePtr{T,A}, y::Integer) where {T,A} = DevicePtr{T,A}(pointer(x) - y)
Base.:(+)(x::Integer, y::DevicePtr) = y + x

## memory operations

@static if Base.libllvm_version < v"7.0"
    # Old values (LLVM 6)
    Base.convert(::Type{Int}, ::Type{AS.Private})  = 0
    Base.convert(::Type{Int}, ::Type{AS.Global})   = 1
    Base.convert(::Type{Int}, ::Type{AS.Constant}) = 2
    Base.convert(::Type{Int}, ::Type{AS.Local})    = 3
    Base.convert(::Type{Int}, ::Type{AS.Generic})  = 4
    Base.convert(::Type{Int}, ::Type{AS.Region})   = 5
else
    # New values (LLVM 7+)
    Base.convert(::Type{Int}, ::Type{AS.Generic})  = 0
    Base.convert(::Type{Int}, ::Type{AS.Global})   = 1
    Base.convert(::Type{Int}, ::Type{AS.Region})   = 2
    Base.convert(::Type{Int}, ::Type{AS.Local})    = 3
    Base.convert(::Type{Int}, ::Type{AS.Constant}) = 4
    Base.convert(::Type{Int}, ::Type{AS.Private})  = 5
end

function tbaa_make_child(name::String, constant::Bool=false; ctx::LLVM.Context=JuliaContext())
    tbaa_root = MDNode([MDString("roctbaa", ctx)], ctx)
    tbaa_struct_type =
        MDNode([MDString("roctbaa_$name", ctx),
                tbaa_root,
                LLVM.ConstantInt(0, ctx)], ctx)
    tbaa_access_tag =
        MDNode([tbaa_struct_type,
                tbaa_struct_type,
                LLVM.ConstantInt(0, ctx),
                LLVM.ConstantInt(constant ? 1 : 0, ctx)], ctx)

    return tbaa_access_tag
end

tbaa_addrspace(as::Type{<:AddressSpace}) = tbaa_make_child(lowercase(String(as.name.name)))

@generated function Base.unsafe_load(p::DevicePtr{T,A}, i::Integer=1,
                                     ::Val{align}=Val(1)) where {T,A,align}
    eltyp = convert(LLVMType, T)

    T_int = convert(LLVMType, Int)
    T_ptr = convert(LLVMType, Ptr{T})

    T_actual_ptr = LLVM.PointerType(eltyp, convert(Int, A))

    # create a function
    param_types = [T_ptr, T_int]
    llvm_f, _ = create_function(eltyp, param_types)

    # generate IR
    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)

        ptr = inttoptr!(builder, parameters(llvm_f)[1], T_actual_ptr)

        ptr = gep!(builder, ptr, [parameters(llvm_f)[2]])
        ld = load!(builder, ptr)

        if A != AS.Generic
            metadata(ld)[LLVM.MD_tbaa] = tbaa_addrspace(A)
        end
        alignment!(ld, align)

        ret!(builder, ld)
    end

    call_function(llvm_f, T, Tuple{Ptr{T}, Int}, :((pointer(p), Int(i-one(i)))))
end

@generated function Base.unsafe_store!(p::DevicePtr{T,A}, x, i::Integer=1,
                                       ::Val{align}=Val(1)) where {T,A,align}
    eltyp = convert(LLVMType, T)

    T_int = convert(LLVMType, Int)
    T_ptr = convert(LLVMType, Ptr{T})

    T_actual_ptr = LLVM.PointerType(eltyp, convert(Int, A))

    # create a function
    param_types = [T_ptr, eltyp, T_int]
    llvm_f, _ = create_function(LLVM.VoidType(JuliaContext()), param_types)

    # generate IR
    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)

        ptr = inttoptr!(builder, parameters(llvm_f)[1], T_actual_ptr)

        ptr = gep!(builder, ptr, [parameters(llvm_f)[3]])
        val = parameters(llvm_f)[2]
        st = store!(builder, val, ptr)

        if A != AS.Generic
            metadata(st)[LLVM.MD_tbaa] = tbaa_addrspace(A)
        end
        alignment!(st, align)

        ret!(builder)
    end

    call_function(llvm_f, Cvoid, Tuple{Ptr{T}, T, Int}, :((pointer(p), convert(T,x), Int(i-one(i)))))
end


text:
using Ipopt
include("../include.jl")

nlp_raw2 = CUTEstModel("NET4")

my_pars = Class_parameters()
one_phase_solve(nlp_raw2, my_pars);
finalize(nlp_raw2)

ORG_STDOUT = STDOUT
file = open("net4.txt", "w")
redirect_stdout(file)

IpoptSolve(nlp_raw2);
my_pars = Class_parameters()
one_phase_solve(nlp_raw2, my_pars);

STDOUT = ORG_STDOUT

close(file)

finalize(nlp_raw2)


text:
# Test of Floyd-Warshall's algorithm

using Graphs
using Base.Test

eweights = [  0.  1.  5. Inf;
            Inf   0.  3.  8.;
            Inf   3.  0.  2.;
             6. Inf  Inf  0. ]

dists0 = [  0. 1.  4. 6.;
           11. 0.  3. 5.;
            8. 3.  0. 2.;
            6. 7. 10. 0.  ]

nexts0 = [1 2 2 2; 3 2 3 3; 4 2 3 4; 1 1 1 4]

dists = copy(eweights)
nexts = fill(-1, (4, 4))

floyd_warshall!(dists, nexts)

@test dists == dists0
@test nexts == nexts0

@test floyd_warshall(eweights) == dists0
@test_throws ArgumentError floyd_warshall!(rand(3,5))
@test_throws ArgumentError floyd_warshall!(rand(3,5),rand(1:7,5,5))
@test_throws ArgumentError floyd_warshall!(rand(5,5),rand(1:7,3,5))


text:
# This script illustrates backpropagation through a simulation, via Flux,
# DifferentialEquations and ForwardDiff. Flux's ADAM optimiser and training loop
# are then used to optimise parameters of the simulation.

using OrdinaryDiffEq, Plots

# ––––––––––––––––––––––––––––––––––––––––––––––––––––– #
#                      ODE setup                        #
# ––––––––––––––––––––––––––––––––––––––––––––––––––––– #

# The ODE
function lotka_volterra(du,u,p,t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = (α - β*y)x
  du[2] = dy = (δ*x - γ)y
end

const initial_pop = 1

# Solve the ODE with a given set of parameters, to see how the predator/prey
# populations behave over time.
function trajectory(predator, prey, tfinal = 10)
  params = [predator..., prey...]
  T = eltype(params)
  u0 = T.([initial_pop, initial_pop])
  tspan = (T(0), T(tfinal))
  prob = ODEProblem(lotka_volterra, u0, tspan, params)
  solve(prob, Tsit5(), dtmin = 1e-4)
end

# See an example solution
plot(trajectory([1.8, 1.5], [1.2, 3]), ylim=(0,6))

# For now, our loss is the deviation from the initial population;
# we are optimising for stability.
function stability(sol::ODESolution)
  sol.retcode != :Success && return zero(sol.u[1][1])
  series = sol.(range(0, stop = 10, length = 50))
  sum(x -> sum(x -> (x - initial_pop)^2, x), series)/length(series)
end

stability(predator, prey, tfinal = 100) =
  stability(trajectory(predator, prey, tfinal))

# Preview the loss
stability([1.8, 1.5], [1.2, 3])

# ––––––––––––––––––––––––––––––––––––––––––––––––––––– #
#                      Autodiff                         #
# ––––––––––––––––––––––––––––––––––––––––––––––––––––– #

using Flux
import Flux.Tracker: TrackedVector, @grad, track
using ForwardDiff: Dual, value, partials

# Provide a custom gradient for `stability`, which will use forward-mode AD
# internally.
stability(predator::TrackedVector, prey) = track(stability, predator, prey)

@grad function stability(predator, prey)
  s = stability(Tracker.data(predator) .+ [Dual(0,1,0),Dual(0,0,1)], prey)
  value(s), Δ -> (Δ .* partials(s), nothing)
end

# Now we can take gradients w.r.t. the parameters.

pred = [2.2, 1.0]
prey = [2, 3]
Tracker.derivative(p -> stability(p, prey), pred)

# ––––––––––––––––––––––––––––––––––––––––––––––––––––– #
#               Optimising Parameters                   #
# ––––––––––––––––––––––––––––––––––––––––––––––––––––– #

predator = param([2.2, 1.0])
prey = [2, 3]

data = Iterators.repeated((), 100)
opt = ADAM([predator], 0.1)
cb = () ->
  display(plot(trajectory(Flux.data(predator), prey), ylim=(0,6)))

# Display the ODE with the current parameter values.
cb()

# Running this in Juno will generate an animation of the ODE over time.
Flux.train!(() -> stability(predator, prey),
            data, opt, cb = cb)

max_stars_repo_name: NTimmons/model-zoo

max_issues_repo_name: NTimmons/model-zoo

max_forks_repo_name: NTimmons/model-zoo