15.4.980  src/lib/src/random-sample.pkg

## random-sample.pkg
## Author: Matthias Blume (blume@tti-c.org)

# Compiled by:
#     src/lib/std/standard.lib

#   Randomized linear-time selection from an unordered sample.


###            "Did you ever observe to whom the accidents happen?
###             Chance favors only the prepared mind."
###
###                                   -- Louis Pasteur


stipulate
    package f8b =  eight_byte_float;                                    # eight_byte_float      is from   src/lib/std/eight-byte-float.pkg
herein

    package   random_sample
    : (weak)
    api {

        # WARNING: Each of the functions exported from this module
        # modifies its argument rw_vector by (partially) sorting it.

        # Select the i-th order statistic:
        #
        random_selection  :   (Rw_Vector( Float), Int) -> Float;
        random_selection' :   (rw_vector_slice::Slice( Float ), Int) -> Float;

        # Calculate the median:
        #    if N is odd, then this is the (floor (N/2))th order statistic
        #    otherwise it is the average of (N/2-1)th and (N/2)th
        #
        median:   Rw_Vector( Float ) -> Float;
        median' : rw_vector_slice::Slice( Float ) -> Float;

    }
    {
        infix my 90  @@@ ;    my (@@@)        = unsafe::rw_vector::get;
        infix my 40  <-  ;    fun (a, i) <- x = unsafe::rw_vector::set (a, i, x);

        # Initialize random number generator:
        #
        rand = random::make_random_number_generator (123, 73256);


        # Select i-th order statistic from unsorted rw_vector with
        # starting point p and ending point r (inclusive):
        #
        fun random_selection_0 (a: Rw_Vector( Float ), p, r, i)
            =
            {   fun x + y = unt::to_int_x (unt::(+) (unt::from_int x, unt::from_int y));
                fun x - y = unt::to_int_x (unt::(-) (unt::from_int x, unt::from_int y));

                # Random partition: 
                #
                fun rp (p, r)
                    =
                    {   fun sw (i, j) = { t=a@@@i;  (a, i)<-a@@@j; (a, j)<-t; };

                        q = random::range (p, r) rand;
                        qv = a@@@q;

                        if   (q != p)
                             (a, q)<-a@@@p;
                             (a, p)<-qv;
                        fi;

                        fun up i
                            =
                            if (i>r or qv < a@@@i ) i; else up (i+1);fi;

                        fun dn i
                            =
                            if (i>=p and qv < a@@@i ) dn (i - 1); else i;fi;

                        fun lp (i, j)
                            =
                            {   my (i, j) = (up i, dn j);

                                if   (i > j)
                                     q' = i - 1;
                                     sw (p, q');
                                     (q', qv);
                                else
                                     sw (i, j);
                                     lp (i+1, j - 1);
                                fi;
                            };

                     lp (p+1, r); };

                # Random selection:
                #
                fun rs (p, r)
                    =
                    if  (p==r)
                         a@@@r;
                    else
                         my (q, qv) = rp (p, r);

                         if   (i==q)
                              qv;
                         else
                              if (i < q   )   rs (p, q - 1);
                                         else   rs (q+1, r);   fi;
                         fi;
                    fi;
             rs (p, r);
        };

        fun random_selection (a, i)
            =
            random_selection_0 (a, 0, rw_vector::length a - 1, i);

        fun random_selection' (s, i)
            =
            {   (rw_vector_slice::burst_slice s)
                    ->
                    (a, p, l);

                random_selection_0 (a, p, p+l - 1, p+i);
            };

        fun median0 (a, p, len)
            =
            {   mid = p + len / 2;
                r = p + len - 1;
                m0 = random_selection_0 (a, p, r, mid);

                fun l (i, m)
                    =
                    if (i>=mid)   m;
                    else l       (i+1, f8b::max (a@@@i, m));
                    fi;

                if (len % 2 == 1)   m0;
                else                (l (p+1, a@@@p) + m0) / 2.0;
                fi;
            };

        fun median  a =  median0 (a, 0, rw_vector::length a);
        fun median' s =  median0 (rw_vector_slice::burst_slice s);

    };
end;