F# performance exploration code

fsharp_perf_test.fs

// small perf tests for doing floating point format conversions
// full credit to http://www.palladiumconsulting.com/2014/09/little-performance-explorations-f/ as the original source for most of this

open System.IO
open System.Runtime.InteropServices
open System.Diagnostics

type IbmFloat32 = 
    struct
        val Bits: uint32
    end
    new (b:uint32) = { Bits = b } 

[<StructLayout(LayoutKind.Explicit)>]
type IntFloat =
    struct
        [<FieldOffset(0)>] val mutable f: float32
        [<FieldOffset(0)>] val mutable i: uint32
    end

// test: remove inlining here
let inline ToFloat32 (i:uint32) =
    let mutable x = IntFloat()
    x.i <- i;
    x.f 

// test: remove inlining here
let inline ieeeOfPieces fr exp sgn = 
    uint32(fr >>> 9) ||| uint32(exp <<< 23) ||| sgn 
    |> ToFloat32

// test: remove inlining here
let inline ToIeee (ibm:IbmFloat32) : float32 =
    let mutable fr = ibm.Bits

    let sgn = fr &&& 0x80000000u // save sign
    fr <- fr <<< 1 // shift sign out

    let exp = int(fr >>> 25) // save exponent
    fr <- fr <<< 7 // shift exponent out
 
    if fr = 0u then // short-circuit, return properly signed zero
        ieeeOfPieces 0u 0 sgn
    else
        // adjust exponent from base 16 offset 64 radix point
        // before first digit to base 2 offset 127 radix point
        // after first digit 
        // (exp - 64) * 4 + 127 - 1 
        let mutable exp = (exp - 64) * 4 + 127 - 1

        // (re)normalize, 3 times max for normalized input */
        while fr < 0x80000000u do exp <- exp - 1; fr <- fr <<< 1
    
        if exp <= 0 then // underflow
            let fr = 
                if exp < -24 then
                    0u // complete underflow, return properly signed zero
                else
                    fr >>> -exp // partial underflow, return denormalized number
            ieeeOfPieces fr 0 sgn
        elif exp >= 255 then // overflow - return infinity
            ieeeOfPieces 0u 255 sgn
        else // just a plain old number - remove the assumed high bit
            ieeeOfPieces (fr <<< 1) exp sgn

let b v shift = (int v) <<< (8 * shift)

let readBigEndianInt32 (a:byte[]) i = 
    b a.[i] 3 ||| 
    b a.[i+1] 2 ||| 
    b a.[i+2] 1 ||| 
    b a.[i+3] 0
    |> uint32
 
let time f = Stopwatch.StartNew() |> (fun sw -> (f(), sw.Elapsed.TotalMilliseconds))
 
let bench () =
    let fileHeaderLength = 3200 + 400 // SEG-Y text and binary headers
    let nTraces = 64860 //500 // number of traces in the file
    let headerLength = 240 // each SEGY trace has a 240 byte header at the beginning, which we ignore
    let nSamples = 1501 // Each trace has 1501 4 byte samples, each an IBM floating point number
    let traceLength = headerLength + (nSamples * 4) // = 6244 length in bytes of each trace, including its header
 
    let buffer = File.ReadAllBytes @"C:\users\latkin\documents\fssandbox\filt_mig.sgy"

    ////////////////////
    // original approach
    ////////////////////
    let mapi' f (a:'T[]) = 
        for i in 0..a.Length-1 do a.[i] <- f i a.[i]
        a
    let convertOrig() =
        seq { for n in 0 .. nTraces-1 -> fileHeaderLength + n * traceLength + headerLength }
        |> Seq.fold
            (fun (samples:float32[]) traceOffset -> 
                samples |> mapi' (fun i s -> s + ToIeee(IbmFloat32(readBigEndianInt32 buffer (traceOffset + i*4)))))
            (Array.zeroCreate nSamples)
        |> ignore

    ////////////////////////////
    // fully imperative approach
    ////////////////////////////
    let convertImperative() =
        let samples : float32[] = Array.zeroCreate nSamples
        for n in 0 .. nTraces-1 do
            let traceOffset = fileHeaderLength + n * traceLength + headerLength        
            for i in 0 .. nSamples-1 do
                samples.[i] <- samples.[i] + ToIeee(IbmFloat32(readBigEndianInt32 buffer (traceOffset + i*4)))
    
    ///////////////////////////
    // inline the mapi function
    ///////////////////////////
    let inline mapi'' f (a:'T[]) = 
        for i in 0..a.Length-1 do a.[i] <- f i a.[i]
        a
    let convertInlineMapi() =
        seq { for n in 0 .. nTraces-1 -> fileHeaderLength + n * traceLength + headerLength }
        |> Seq.fold
            (fun (samples:float32[]) traceOffset -> 
                samples |> mapi'' (fun i s -> s + ToIeee(IbmFloat32(readBigEndianInt32 buffer (traceOffset + i*4)))))
            (Array.zeroCreate nSamples)
        |> ignore

    /////////////////////////
    // remove mapi altogether
    /////////////////////////
    let convertNoMapi() =      
        seq { for n in 0 .. nTraces-1 -> fileHeaderLength + n * traceLength + headerLength }
        |> Seq.fold
            (fun (samples:float32[]) traceOffset -> 
                for i in 0..samples.Length-1 do
                    samples.[i] <- samples.[i] + ToIeee(IbmFloat32(readBigEndianInt32 buffer (traceOffset + i*4)))
                samples)
            (Array.zeroCreate nSamples)
        |> ignore

    printfn "convertOrig        %f ms" ((time>>snd) convertOrig)
    printfn "convertImperative  %f ms" ((time>>snd) convertImperative)
    printfn "convertInlineMapi  %f ms" ((time>>snd) convertInlineMapi)
    printfn "convertNoMapi      %f ms" ((time>>snd) convertNoMapi)

    System.Console.ReadKey() |> ignore
    0

bench () |> ignore