//
// Graphs for Random Walking
// Nicholas Gorski (nicholasgorski@ngorski.com)
//

// Options
numberPaths1D = 5 // number of paths to be drawn (1D)
pathLengthMin1D = 50 // minimum length of a path to be drawn (1D)
pathLengthMax1D = 100 // maximum length of a path to be drawn (1D)
stepBasedLine1D = 0 // draw each step differently (1D), either on (1) or off (0)

numberPaths2D = 1 // number of paths to be drawn (2D)
pathLengthMin2D = 50 // minimum length of a path to be drawn (2D)
pathLengthMax2D = 1000 // maximum length of a path to be drawn (2D)
stepBasedLine2D = 1 // draw each step differently (2D), either on (1) or off (0)
shakeLine2D = 1 // offset duplicate points (2D), either on (1) or off (0)

numberPaths3D = 1 // number of paths to be drawn (3D)
pathLengthMin3D = 300 // minimum length of a path to be drawn (3D)
pathLengthMax3D = 1000 // maximum length of a path to be drawn (3D)
stepBasedLine3D = 1 // draw each step differently (3D), either on (1) or off (0)
shakeLine3D = 1 // offset duplicate points (3D), either on (1) or off (0)

requireTermination3D = 0 // require a 3D walk reach origin, either on (1) or off (0) 

// warning: enabling the previous flag may cause the process to take
// enourmous amounts of time, depending on length bounds; patience is required!



// Edit anything under this at your own risk!
////////////////////////////////////////////////////////////////////////////////////////////////////

// setup:
f = gcf() // get rid of any initial graphs, will be created manually
delete(f)

lines(0) // allow uninterrupted disp

function cm = customcolormap(n)

    cm = bonecolormap(max(3, n))
    cmColors = jetcolormap(max(3, n))

    for i = 1:n
        for j = 1:3            
            cm(i, j) = cmColors(i, j) * (1.0 - cm(i, j)) * 0.9
        end
    end

endfunction

// random seed
rand('seed', getdate('s'))

//////////////////////////////////////////////////////////////////////
// One-dimensional Walk
//////////////////////////////////////////////////////////////////////

// figure setup
fig = scf()
fig.auto_resize = "on"
fig.anti_aliasing = "16x"

if stepBasedLine1D == 1 then
    colorMap = []
else
    fig.color_map = [customcolormap(numberPaths1D); 0, 0, 0]
end

// axis setup
axis = gca()
xtitle("Random Walk 1D", "Step", "Position")

drawlater() // disable drawing while generating image, for speed
furthestWalk = 0 // for drawing the origin line later

drawnPaths = 0
while drawnPaths <> numberPaths1D

    // start at origin
    position = 0
    positions = [position]

    while size(positions, 'r') - 1 == 0 | position <> 0

        // add new value
        r = rand(1, 'uniform')
        if r <= 0.5 then
            position = position + 1
        else
            position = position - 1
        end
        
        positions = [positions; position]

        // check current length
        if size(positions, 'r') - 1 > pathLengthMax1D then
            break // too long already, stop walking
        end

    end

    // check final path length
    pathLength = size(positions, 'r') - 1
    
    if pathLength > pathLengthMax1D then
        continue // final path too long
    end

    if pathLength < pathLengthMin1D then
        continue // final path too short
    end

    if pathLength > furthestWalk then
        furthestWalk = size(positions, 'r') - 1 // longest walk yet
    end

    // draw path
    if stepBasedLine1D == 1 then

        styleBase = size(colorMap, "r")

        colorMap = [colorMap; customcolormap(pathLength + 1)]
        fig.color_map = colorMap

        for i = 1:pathLength
            plot2d([i - 1, i], ..
                    [positions(i), positions(i + 1)], ..
                    style = styleBase + i)
        end

    else

        plot2d(0:pathLength, positions(:), style = drawnPaths + 1)

    end

    drawnPaths = drawnPaths + 1

    disp("Path drawn (#, Steps):")
    disp([drawnPaths, pathLength])

end

if stepBasedLine1D == 1 then
    fig.color_map = [fig.color_map; 0, 0, 0]
end

// origin line
plot2d([0, furthestWalk], [0, 0], style = size(fig.color_map, "r"))

drawnow()
disp("Completed 1D walk.")

//////////////////////////////////////////////////////////////////////
// Two-dimensional Walk
//////////////////////////////////////////////////////////////////////

// figure setup
fig = scf()
fig.auto_resize = "on"
fig.anti_aliasing = "16x"
if stepBasedLine2D == 1 then
    colorMap = []
else
    fig.color_map = [customcolormap(numberPaths2D); 0, 0, 0]
end

// axis setup
axis = gca()
xtitle("Random Walk 2D", "Position X", "Position Y")

drawlater() // disable drawing while generating image, for speed

drawnPaths = 0
while drawnPaths <> numberPaths2D

    // start at origin
    position = [0, 0]
    realPositions = [position]
    positions = [position]

    while size(positions, 'r') - 1 == 0 | ..
            position(1) <> 0 | position(2) <> 0

        r = 4 * rand(1, 'uniform')
        if r <= 1 then
            position(1) = position(1) + 1
        elseif r <= 2 then
            position(1) = position(1) - 1
        elseif r <= 3 then
            position(2) = position(2) + 1
        else
            position(2) = position(2) - 1
        end
        
        realPositions = [realPositions; position]
        positions = [positions; position]

        // check current length
        if size(positions, 'r') - 1 > pathLengthMax2D then
            break // too long already, stop walking
        end

    end
    
    // check final path length
    pathLength = size(positions, 'r') - 1
    
    if pathLength > pathLengthMax2D then
        continue // final path too long
    end

    if pathLength < pathLengthMin2D then
        continue // final path too short
    end

    if shakeLine2D == 1 then

        // shake draw points slightly to avoid total overlap
        duplicateOccurrences = []
        for i = 1:(pathLength + 1)

            // get point
            position = realPositions(i, [1, 2])

            if position(1) == 0 & position(2) == 0 then
                continue // don't mess with origin
            end

            // find all other occurrences
            pointOccurrences = []
            for j = i:(pathLength + 1)
                if realPositions(j, [1, 2]) == position then
                    pointOccurrences = [pointOccurrences, j]
                end
            end

            if length(pointOccurrences) <> 1 then
                duplicateOccurrences = [duplicateOccurrences, ..
                                        pointOccurrences]
            end

        end

        offsetAngles = [%pi / 4, 3 * %pi / 4, .. // offset by 45 degrees
                        5 * %pi / 4, 7 * %pi / 4]
        for i = 1:length(duplicateOccurrences)

            offsetAngleBase = offsetAngles(round(1 + 3 * rand(1, 'uniform')))
            offsetAngleMin = offsetAngleBase - 0.9 * %pi / 4
            offsetAngleMax = offsetAngleBase + 0.9 * %pi / 4
            offsetAngle = offsetAngleMin + ..
                            rand(1, 'uniform') * (offsetAngleMax - offsetAngleMin)

            offsetRadius = sqrt(2) * 0.111605 * rand(1, 'normal')

            offset = offsetRadius * [cos(offsetAngle), sin(offsetAngle)]

            index = duplicateOccurrences(i)
            positions(index, [1, 2]) = realPositions(index, [1, 2]) + offset

        end

    end

    // draw path
    if stepBasedLine2D == 1 then

        styleBase = size(colorMap, "r")

        colorMap = [colorMap; customcolormap(pathLength)]
        fig.color_map = colorMap

        for i = 1:pathLength
            plot2d([positions(i, 1), positions(i + 1, 1)], ..
                    [positions(i, 2), positions(i + 1, 2)], ..
                    style = styleBase + i)
        end

    else

        plot2d(positions(:, 1), positions(:, 2), style = drawnPaths + 1)

    end

    drawnPaths = drawnPaths + 1

    disp("Path drawn (#, Steps):")
    disp([drawnPaths, pathLength])

end


// grid lines
plot([0, 0], [0, 0], "blackpentagram")

drawnow()
disp("Completed 2D walk.")

//////////////////////////////////////////////////////////////////////
// Three-dimensional Walk
//////////////////////////////////////////////////////////////////////

// figure setup
fig = scf()
fig.auto_resize = "on"
fig.anti_aliasing = "16x"
if stepBasedLine3D == 1 then
    colorMap = [0, 0, 0]
else
    fig.color_map = [0, 0, 0; customcolormap(numberPaths3D)]
end

// axis setup
axis = gca()
axis.hidden_axis_color = 1
xtitle("Random Walk 3D", "Position X", "Position Y", "Position Z")

drawlater() // disable drawing while generating image, for speed

// view and grid lines
viewTheta = 20
viewPhi = 70

majorVec = [sqrt(2) / 4, -sqrt(2) / 4]
minorVec = [-majorVec(1), -majorVec(2)]

param3d(majorVec, minorVec, minorVec, viewTheta, viewPhi)
param3d(minorVec, majorVec, minorVec, viewTheta, viewPhi)
param3d(minorVec, minorVec, majorVec, viewTheta, viewPhi)

drawnPaths = 0
while drawnPaths <> numberPaths3D

    // start at origin
    position = [0, 0, 0]
    realPositions = [position]
    positions = [position]

    while size(positions, 'r') - 1 == 0 | ..
            position(1) <> 0 | position(2) <> 0 | position(3) <> 0

        r = 6 * rand(1, 'uniform')
        if r <= 1 then
            position(1) = position(1) + 1
        elseif r <= 2 then
            position(1) = position(1) - 1
        elseif r <= 3 then
            position(2) = position(2) + 1
        elseif r <= 4 then
            position(2) = position(2) - 1
        elseif r <= 5 then
            position(3) = position(3) + 1
        else
            position(3) = position(3) - 1
        end

        realPositions = [realPositions; position]
        positions = [positions; position]

        // check current length
        if size(positions, 'r') - 1 > pathLengthMax3D then
            break // too long already, stop walking
        end

    end

    // check final path length
    pathLength = size(positions, 'r') - 1

    // (allow paths that don't terminate, not guaranteed in 3D)
    if requireTermination3D == 1 & pathLength > pathLengthMax3D then
        continue // final path too long
    end

    if pathLength < pathLengthMin3D then
        continue // final path too short
    end

    if shakeLine3D == 1 then

        // shake draw points slightly to avoid total overlap
        duplicateOccurrences = []
        for i = 1:(pathLength + 1)

            // get point
            position = realPositions(i, [1, 2, 3])

            if position(1) == 0 & position(2) == 0 & position(3) == 0 then
                continue // don't mess with origin
            end

            // find all other occurrences
            pointOccurrences = []
            for j = i:(pathLength + 1)
                if realPositions(j, [1, 2, 3]) == position then
                    pointOccurrences = [pointOccurrences, j]
                end
            end

            if length(pointOccurrences) <> 1 then
                duplicateOccurrences = [duplicateOccurrences, ..
                                        pointOccurrences]
            end

        end

        offsetAngles = [%pi / 4, 3 * %pi / 4, .. // offset 45 degrees
                        5 * %pi / 4, 7 * %pi / 4]
        for i = 1:length(duplicateOccurrences)

            offsetAngle1Base = offsetAngles(round(1 + 3 * rand(1, 'uniform')))
            offsetAngle1Min = offsetAngle1Base - 0.9 * %pi / 4
            offsetAngle1Max = offsetAngle1Base + 0.9 * %pi / 4
            offsetAngle1 = offsetAngle1Min + ..
                            rand(1, 'uniform') * (offsetAngle1Max - offsetAngle1Min)

            offsetAngle2Base = offsetAngles(round(1 + 3 * rand(1, 'uniform')))
            offsetAngle2Min = offsetAngle2Base - 0.9 * %pi / 4
            offsetAngle2Max = offsetAngle2Base + 0.9 * %pi / 4
            offsetAngle2 = offsetAngle2Min + ..
                            rand(1, 'uniform') * (offsetAngle2Max - offsetAngle2Min)

            offsetRadius = sqrt(2) * 0.111605 * rand(1, 'normal')

            cosOffsetAngle1 = cos(offsetAngle1)
            sinOffsetAngle1 = sin(offsetAngle1)
            cosOffsetAngle2 = cos(offsetAngle2)
            sinOffsetAngle2 = sin(offsetAngle2)

            offset = [cosOffsetAngle1 * cosOffsetAngle2 * offsetRadius, ..
                        sinOffsetAngle2 * offsetRadius, ..
                        sinOffsetAngle1 * cosOffsetAngle2 * offsetRadius]

            index = duplicateOccurrences(i)
            positions(index, [1, 2, 3]) = realPositions(index, [1, 2, 3]) + offset

        end

    end

    // draw path
    if stepBasedLine3D == 1 then

        styleBase = size(colorMap, "r")

        colorMap = [colorMap; customcolormap(pathLength)]
        fig.color_map = colorMap

        for i = 1:pathLength
            param3d([positions(i, 1), positions(i + 1, 1)], ..
                        [positions(i, 2), positions(i + 1, 2)], ..
                        [positions(i, 3), positions(i + 1, 3)], ..
                        viewTheta, viewPhi)

            axis.children(1).foreground = styleBase + i
        end

    else

        param3d(positions(:, 1), positions(:, 2), positions(:, 3), ..
                    viewTheta, viewPhi)

        axis.children(1).foreground = drawnPaths + 2

    end

    drawnPaths = drawnPaths + 1

    disp("Path drawn (#, Steps, Reached Origin):")
    disp([drawnPaths, pathLength, ..
            realPositions($, 1) == 0 & ..
                realPositions($, 2) == 0 & realPositions($, 3) == 0])

end

drawnow()
disp("Completed 3D walk.")
