Translating C function with a loop to Haskell

Question

I am trying to reason how to convert an imperative style program into a functional one like Haskell.

The function is:

void calcPerim(point polycake[], int v, int y, double *perim1, double *perim2){
    int next = 0;
    int index = 0;
    point points[2];
    *perim1 = 0.0;
    *perim2 = 0.0;

    for(int i = 0; i < v; i++)
    {
        next = (i + 1) % v;
        if(polycake[i].y < y && polycake[next].y < y)
            (*perim1) += distance(polycake[i], polycake[next]);
        else if(polycake[i].y > y && polycake[next].y > y)
            (*perim2) += distance(polycake[i], polycake[next]);
        else
        {
            points[index] = intersectPoint(polycake[i], polycake[next], y);
            if(polycake[i].y < y)
            {
                (*perim1) += distance(polycake[i], points[index]);
                (*perim2) += distance(polycake[next],points[index]);
            }
            else
            {
                (*perim2) += distance(polycake[i], points[index]);
                (*perim1) += distance(polycake[next],points[index]);
            }
            index++;
        }
    }

    (*perim1) += distance(points[0], points[1]);
    (*perim2) += distance(points[0], points[1]);
}

I am finding it difficult to understand how I can turn this into a functional approach when it is updating two variables at the same time in some cases. Would it make sense when translating this into recursion to pass in a tuple (perim1, perim2)?

Answer 1

You can give this a try, it is a direct translation of your C algorithm to Haskell

data Point = Point {x :: Float, y :: Float}


calcPerim :: [Point] -> Int -> Int -> (Float, Float)
calcPerim ls v some_y = 
    let (x:xs)       = take v ls
        r            = zip (x:xs) (xs ++ [x])
        (u, c, o, _) = foldl someFunction (0, 0, [], fromIntegral some_y :: Float) r
        points_0     = last o
        points_1     = o !! ((length o) - 2)
        answer       = (u + (distance points_0 points_1), c + (distance points_0 points_1))
    in  answer


someFunction :: (Float, Float, [Point], Float) -> (Point, Point) -> (Float, Float, [Point], Float)
someFunction (perim_1, perim_2, points, some_y) (i, nxt)
    | y i < some_y && y nxt < some_y    = (perim_1 + (distance i nxt), perim_2, points, some_y)
    | y i > some_y && y nxt > some_y    = (perim_1, perim_2 + (distance i nxt), points, some_y)
    | y i < some_y                      = (perim_1 + (distance i temp_pt), perim_2 + (distance nxt temp_pt), temp_pt:points, some_y)
    | otherwise                         = (perim_1 + (distance nxt temp_pt), perim_2 + (distance i temp_pt), temp_pt:points, some_y)
    where temp_pt = intersection i nxt some_y


distance :: Point -> Point -> Float
distance p q = undefined


intersection :: Point -> Point -> Float -> Point
intersection p q f = undefined

I didn't run it. Not sure if I used the right fold.

Answer 2

It's rarely a good idea to first write a C version and then try to translate it to Haskell.

Instead, consider what you're trying to do, rather than how you're trying to do it.

It appears that given a series of point representing a polygon and a horizontal line at height y, you want to divide it into two polygons at line y and return the perimeter of both. The algorithm assumes the polygon is convex on the vertical axis:

You're doing this by:

1. Dividing the segments into those entirely over and entirely under y
2. Segments that cross y are split into two parts, the one above and the one below y, indicated by red dots.
3. Adding the intersection line between the two split points (cyan) to both polygons.

We can just implement that logic directly, rather than trying to emulate the iterative approach. Here's an example:

type Length = Double
type Point = (Double, Double)
type Segment = (Point, Point)

-- Check whether a segment is over, under or on the line given by y
segmentCompare :: Double -> Segment -> Ordering
segmentCompare y (p,q) =                                                                          
    case () of                                                                                        
        _ | all (`isUnder` y) [p,q] -> LT
        _ | all (`isOver` y)  [p,q] -> GT
        _ -> EQ

-- Partition a list into (lt, eq, gt) based on f
partition3 :: (Segment -> Ordering) -> [Segment] -> ([Segment], [Segment], [Segment])
partition3 f = p' ([], [], [])
  where
    p' (lt, eq, gt) (x:xs) =
        case f x of
            LT -> p' (x:lt, eq, gt) xs
            EQ -> p' (lt, x:eq, gt) xs
            GT -> p' (lt, eq, x:gt) xs
    p' result [] = result

-- Split a crossing segment into an under part and over part, and return middle
divvy :: Double -> Segment -> (Segment, Segment, Point)
divvy y (start, end) =
    if start `isUnder` y
    then ((start, middle), (middle, end), middle)
    else ((middle, end), (start, middle), middle)
  where
    middle = intersectPoint y (start, end)

-- Split a polygon in two, or Nothing if it's not convex enough
splitPolygon :: Double -> [Point] -> Maybe ([Segment], [Segment])
splitPolygon y list = do
    let (under, crossing, over) = partition3 (segmentCompare y) pairs
    case crossing of
        -- No lines cross. Simple.
        [] -> return (under, over)
        -- Two segments cross. Divide them up.
        [(p1,p2),(q1,q2)] ->
            let (u1, o1, mid1) = divvy y (p1,p2)
                (u2, o2, mid2) = divvy y (q1, q2)
                split = (mid1, mid2) :: Segment
            in return (split:u1:u2:under, split:o1:o2:over)
        -- More segments cross. Algorithm doesn't work.
        rest -> fail "Can't split polygons concave at y"
  where
    pairs = zip list (drop 1 $ cycle list) :: [Segment]

-- Your original function that sums the perimeter of both polygons
calcPerim :: Double -> [Point] -> Maybe (Length, Length)
calcPerim y list = do
    (under, over) <- (splitPolygon y list :: Maybe ([Segment], [Segment]))
    return (sumSegments under, sumSegments over)

-- Self explanatory helpers
distance :: Segment -> Length
distance ((ax, ay), (bx, by)) = sqrt $ (bx-ax)^2 + (by-ay)^2

intersectPoint :: Double -> Segment -> Point
intersectPoint y ((px, py), (qx, qy)) =
    let slope     = (qx-px)/(qy-py)
        intercept = qy - slope*qx
        x = (y - intercept)/slope
    in
        if slope /= 0
        then (x,y)
        else (px, y)

sumSegments :: [Segment] -> Length
sumSegments = sum . map distance

isUnder :: Point -> Double -> Bool
isUnder (_, py) y = py < y
isOver (_, py) y = py > y

Answer 3

It might be a good idea to not translate it straight to Haskell but rather first to C++, which already allows to you structure it in a much more functional way.

First thing, as Cirdec commented, this function doesn't really take perim1 as arguments – those are “output arguments” as Fortran people would say, i.e. they're really results. Also, the v parameter seems to be basically just length of the input array. So in C++ you can reduce it to:

std::pair<double, double> calcPerim(std::vector <point> polycake, int y){
   double perim1 = 0, perim2 = 0;
   ...
   return std::make_pair(perim1, perim2);
}

Now, you have this mutating for loop. In a functional language, the general approach would be to replace that with recursion. For this, you need to make all mutable-state variables function parameters. That includes i, index, points and the perim accumulators (so they're back, in a way... but now as input arguments). You don't need next (which is anyways re-computed from scratch in each iteration).

std::pair<double, double> calcPerim_rec
                   ( std::vector<point> polycake, int y
                   , int i, int index, std::array<point,2> points
                   , double perim1Acc, double perim2Acc ){
   ...
}

...to be used by

std::pair<double, double> calcPerim(std::vector<point> polycake, int y){
   return calcPerim_rec(polycake, y, 0, 0, {}, 0, 0);
}

The recursive function looks very similar to your original loop body; you just need to phrase the end condition:

std::pair<double, double> calcPerim_rec
                   ( std::vector<point> polycake, int y
                   , int i, int index, std::array<point,2> points
                   , double perim1Acc, double perim2Acc ){
    if (i < polycake.length()) {
        int next = (i + 1) % polycake.length();
        if(polycake[i].y < y && polycake[next].y < y)
            perim1Acc += distance(polycake[i], polycake[next]);
        else if(polycake[i].y > y && polycake[next].y > y)
            perim2Acc += distance(polycake[i], polycake[next]);
        else
        {
            points[index] = intersectPoint(polycake[i], polycake[next], y);
            if(polycake[i].y < y)
            {
                perim1Acc += distance(polycake[i], points[index]);
                perim2Acc += distance(polycake[next],points[index]);
            }
            else
            {
                perim2Acc += distance(polycake[i], points[index]);
                perim1Acc += distance(polycake[next],points[index]);
            }
            ++index;
        }
        ++i;
        return calcPerim_rec
                 ( polycake, y, i, index, points, perim1Acc, perim2Acc );
    } else {
       perim1Acc += distance(points[0], points[1]);
       perim2Acc += distance(points[0], points[1]);
       return std::make_pair(perim1Acc, perim2Acc);
    }
}

There's still quite a bit of mutability going on, but we've already encapsulated it to happen all on local variables of the recursion function call, instead of variables lying around during the loop execution. And each of these variables is only updated once, followed by the recursive call, so you can just skip the mutation and simply pass a value plus update to the recursive call:

std::pair<double, double> calcPerim_rec
                   ( std::vector<point> polycake, int y
                   , int i, int index, std::array<point,2> points
                   , double perim1Acc, double perim2Acc ){
    if (i < polycake.length()) {
        int next = (i + 1) % polycake.length();
        if(polycake[i].y < y && polycake[next].y < y)
            return calcPerim_rec
             ( polycake, y, i+1, index, points
             , perim1Acc + distance(polycake[i], polycake[next])
             , perim2Acc
             );
        else if(polycake[i].y > y && polycake[next].y > y)
            return calcPerim_rec
             ( polycake, y, i+1, index, points
             , perim1Acc
             , perim2Acc + distance(polycake[i], polycake[next])
             );
        else
        {
            points[index] = intersectPoint(polycake[i], polycake[next], y);
            if(polycake[i].y < y)
            {
                return calcPerim_rec
                 ( polycake, y, i+1, index+1
                 , points
                 , perim1Acc + distance(polycake[i], points[index])
                 , perim2Acc + distance(polycake[next],points[index])
                 );
            }
            else
            {
                return calcPerim_rec
                 ( polycake, y, i+1, index+1
                 , points
                 , perim1Acc + distance(polycake[i], points[index])
                 , perim2Acc + distance(polycake[next],points[index])
                 );
            }
        }
    } else {
       return std::make_pair( perim1Acc + distance(points[0], points[1])
                            , perim2Acc + distance(points[0], points[1]) );
    }
}

Well, quite a bit of awkward passing-on of parameters, and we still have a mutation of points – but essentially, the code can now be translated to Haskell.

import Data.Vector (Vector, (!), length) as V

calcPerim_rec :: Vector Point -> Int -> Int -> Int -> Int -> [Point] -> (Double, Double) -> (Double, Double)
calcPerim_rec polycake y i index points (perim1Acc, perim2Acc)
 = if i < V.length polycake
    then
     let next = (i + 1) `mod` V.length polycake
     in if yCoord (polycake!i) < y && yCoord (polycake!next) < y
         then calcPerim_rec polycake v y (i+1) index points
                (perim1Acc + distance (polycake!i) (polycake!next)
                perim2Acc
         else
          if yCoord (polycake!i) > y && yCoord (polycake!next) > y)
           then calcPerim_rec polycake v y (i+1) index points
                     perim1Acc
                     (perim2Acc + distance (polycake!i) (polycake!next))
           else
            let points' = points ++ [intersectPoint (polycake!i) (polycake!next) y]
            in if yCoord (polycake!i) < y
                then calcPerim_rec polycake v y (i+1) (index+1)
                       points'
                       (perim1Acc + distance (polycake!i) (points!!index))
                       (perim2Acc + distance (polycake!next) (points!!index))
                else calcPerim_rec polycake v y (i+1) (index+1)
                       points'
                       (perim1Acc + distance (polycake!i) points!!index))
                       (perim2Acc + distance (polycake!next) points!!index))
    else ( perim1Acc + distance (points!!0) (points!!1)
         , perim2Acc + distance (points!!0) (points!!1) )

There's a lot here that could be stylistically improved, but it should in essence work.

A good first thing to actually make it idiomatic is to try and get rid of indices. Indices are strongly eschewed in Haskell, and can often be avoided when you properly work with lists instead of arrays.