Skip to content

Commit

Permalink
Added support for frequency bounds.
Browse files Browse the repository at this point in the history
  • Loading branch information
Jade Cheng committed Jul 27, 2017
1 parent 1b374cd commit 84361a0
Show file tree
Hide file tree
Showing 10 changed files with 108 additions and 37 deletions.
16 changes: 15 additions & 1 deletion src/cpax/jade.improver.hpp
View file
Expand Up @@ -48,13 +48,16 @@ namespace jade
const matrix_type & fb, ///< The 1-F matrix.
const matrix_type & qfa, ///< The Q*F matrix.
const matrix_type & qfb, ///< The Q*(1-F) matrix.
const matrix_type * fif) ///< The Fin-force matrix.
const matrix_type * fif, ///< The Fin-force matrix.
const bool frb) ///< Using frequency-bounds.
{
assert(verification_type::validate_gqf_sizes(g, q, fa));
assert(verification_type::validate_gqf_sizes(g, q, fb));
assert(verification_type::validate_q(q));
assert(verification_type::validate_f(fa));
assert(nullptr == fif || !frb);

const auto I = g.get_height();
const auto K = fa.get_height();
const auto J = fa.get_width();
assert(nullptr == fif || verification_type::
Expand All @@ -68,6 +71,9 @@ namespace jade
matrix_type derivative_vec (K, 1);
matrix_type hessian_mat (K, K);

const auto frb_delta = value_type(1.0) /
(value_type(2 * I) + value_type(1.0));

for (size_t j = 0; j < J; j++)
{
const auto f_column = fa.copy_column(j);
Expand All @@ -93,6 +99,14 @@ namespace jade
b_vec[k + K] = value_type(0);
}
}
else if (frb)
{
for (size_t k = 0; k < K; k++)
{
b_vec[k + 0] -= frb_delta;
b_vec[k + K] -= frb_delta;
}
}

const auto sqp_q = -hessian_mat;
const auto sqp_a = -a_mat;
Expand Down
11 changes: 8 additions & 3 deletions src/cpax/jade.main.cpp
View file
Expand Up @@ -40,8 +40,13 @@ OPTIONS
file identifying an assignment of components
for each individual and a range of Q values for
each component
--frequency-bounds,-frb indicates the algorithm applies bounds between
1 / (2n + 1) and 1 - (1 / (2n + 1)) for allele
frequencies, where n is the number of
individuals; without this flag, this bounds
are set to 0 and 1
--ksize,-k indicates the next argument is the number of
components; this value must be at least two
components; this value must be at least one
--max-iterations,-mi indicates the next argument is the maximum
number of iterations to execute the algorithm;
this value must be greater than or equal to
Expand Down Expand Up @@ -99,7 +104,7 @@ DESCRIPTION
[Notation]
K := Number of Components
This value must be greater than or equal to two.
This value must be greater than or equal to one.
I := Number of Individuals
This value must be greater than or equal to two.
Expand All @@ -123,7 +128,7 @@ DESCRIPTION
Q := [I x K] Proportion Matrix
This matrix consists of floating-point values ranging from 0 to 1. Each
row sums to 1.0, and each row contains more than one distinct value.
row sums to 1.0.
[Matrix File Format]
Expand Down
3 changes: 2 additions & 1 deletion src/cpax/jade.optimizer.hpp
View file
Expand Up @@ -54,6 +54,7 @@ namespace jade
const auto fg = settings.get_fg();
const auto fif = settings.get_fif();
const auto & g = settings.get_g();
const auto frb = opts.is_frb();

const stopwatch sw1;

Expand Down Expand Up @@ -87,7 +88,7 @@ namespace jade

if (!opts.is_fixed_f())
{
fa = improver_type::improve_f(g, q, fa, fb, qfa, qfb, fif);
fa = improver_type::improve_f(g, q, fa, fb, qfa, qfb, fif, frb);
_compute_fb(fa, fb);
matrix_type::gemm(q, fa, qfa);
matrix_type::gemm(q, fb, qfb);
Expand Down
23 changes: 23 additions & 0 deletions src/cpax/jade.options.hpp
View file
Expand Up @@ -57,6 +57,7 @@ namespace jade
, _qin (a.read<std::string>("--qin", "-qi"))
, _qout (a.read<std::string>("--qout", "-qo"))
, _seed (a.read("--seed", "-s", std::random_device()()))
, _frb (a.read_flag("--frequency-bounds", "-frb"))
, _fixed_f (a.read_flag("--fixed-f", "-ff"))
, _fixed_q (a.read_flag("--fixed-q", "-fq"))
, _quiet (a.read_flag("--quiet", "-q"))
Expand Down Expand Up @@ -103,6 +104,19 @@ namespace jade
throw error() << "invalid specification of --fixed-q option "
<< "and --force option";

if (is_frb())
{
if (is_fin_force_specified())
throw error()
<< "invalid specification of --fin-force "
<< "and --frequency-bounds options";

if (is_fixed_f() && is_fin_specified())
throw error()
<< "invalid specification of --fixed-f, --fin, "
<< "and --frequency-bounds options";
}

if (_num_threads == 0)
throw error() << "invalid number of threads specified for "
<< "--num-threads option: " << _num_threads;
Expand Down Expand Up @@ -222,6 +236,14 @@ namespace jade
return !std::isnan(_epsilon);
}

///
/// \return True if the frequency-bounds option is specified.
///
inline bool is_frb() const
{
return _frb;
}

///
/// \return True if the fin option is specified.
///
Expand Down Expand Up @@ -334,6 +356,7 @@ namespace jade
const seed_type _seed;

// options without arguments
const bool _frb;
const bool _fixed_f;
const bool _fixed_q;
const bool _quiet;
Expand Down
37 changes: 11 additions & 26 deletions src/lib/jade.verification.hpp
View file
Expand Up @@ -80,9 +80,9 @@ namespace jade
const auto K = f.get_height();
const auto J = f.get_width();

if (K < 2)
if (K < 1)
throw error() << "invalid F matrix size " << f.get_size_str()
<< " does not contain at least two components";
<< " does not contain at least one component";

if (J < 2)
throw error() << "invalid F matrix size " << f.get_size_str()
Expand Down Expand Up @@ -232,44 +232,29 @@ namespace jade
throw error() << "invalid Q matrix size " << q.get_size_str()
<< " does not contain at least two individuals";

if (K < 2)
if (K < 1)
throw error() << "invalid Q matrix size " << q.get_size_str()
<< " does not contain at least two components";
<< " does not contain at least one component";

for (size_t i = 0; i < I; i++)
{
auto sum = value_type(0);
auto min = q(i, 0);
auto max = q(i, 0);

for (size_t k = 0; k < K; k++)
{
const auto q_ik = q(i, k);

if (q_ik >= value_type(0) && q_ik <= value_type(1))
{
sum += q_ik;
min = std::min(min, q_ik);
max = std::max(max, q_ik);
continue;
}
if (q_ik < value_type(0) || q_ik > value_type(1))
throw error()
<< "invalid Q matrix cell [" << i+1 << "," << k+1
<< "] (" << q_ik << ") is not between 0 and 1";

throw error()
<< "invalid Q matrix cell [" << i+1 << "," << k+1
<< "] (" << q_ik << ") is not between 0 and 1";
sum += q_ik;
}

if (std::fabs(sum - value_type(1)) <= epsilon)
{
if (std::fabs(max - min) >= epsilon)
continue;

if (std::fabs(sum - value_type(1)) > epsilon)
throw error() << "invalid Q matrix row " << i+1
<< " has only one value (" << min << ")";
}

throw error() << "invalid Q matrix row " << i+1
<< " does not sum to 1 (" << sum << ")";
<< " does not sum to 1 (" << sum << ")";
}

return true;
Expand Down
2 changes: 1 addition & 1 deletion src/lib/jade.version.hpp
View file
Expand Up @@ -30,7 +30,7 @@ namespace jade
char_type const * const title, ///< The program name.
ostream_type & out) ///< The output stream.
{
static const values current { 0, 0, 6416, 41394, 2015, 2017 };
static const values current { 0, 0, 6416, 42877, 2015, 2017 };

out << "ohana/" << title << ' '
<< current.major << '.'
Expand Down
16 changes: 15 additions & 1 deletion src/qpas/jade.improver.hpp
View file
Expand Up @@ -48,13 +48,16 @@ namespace jade
const matrix_type & fb, ///< The 1-F matrix.
const matrix_type & qfa, ///< The Q*F matrix.
const matrix_type & qfb, ///< The Q*(1-F) matrix.
const matrix_type * fif) ///< The Fin-force matrix.
const matrix_type * fif, ///< The Fin-force matrix.
const bool frb) ///< Using frequency-bounds.
{
assert(verification_type::validate_gqf_sizes(g, q, fa));
assert(verification_type::validate_gqf_sizes(g, q, fb));
assert(verification_type::validate_q(q));
assert(verification_type::validate_f(fa));
assert(nullptr == fif || !frb);

const auto I = g.get_height();
const auto K = fa.get_height();
const auto J = fa.get_width();
assert(nullptr == fif || verification_type::
Expand All @@ -67,6 +70,9 @@ namespace jade
matrix_type derivative_vec (K, 1);
matrix_type hessian_mat (K, K);

const auto frb_delta = value_type(1.0) /
(value_type(2 * I) + value_type(1.0));

for (size_t j = 0; j < J; j++)
{
const auto f_column = fa.copy_column(j);
Expand All @@ -92,6 +98,14 @@ namespace jade
b_vec[k + K] = value_type(0);
}
}
else if (frb)
{
for (size_t k = 0; k < K; k++)
{
b_vec[k + 0] -= frb_delta;
b_vec[k + K] -= frb_delta;
}
}

std::vector<size_t> active_set { 0 };
matrix_type delta_vec (K, 1);
Expand Down
11 changes: 8 additions & 3 deletions src/qpas/jade.main.cpp
View file
Expand Up @@ -40,8 +40,13 @@ OPTIONS
file identifying an assignment of components
for each individual and a range of Q values for
each component
--frequency-bounds,-frb indicates the algorithm applies bounds between
1 / (2n + 1) and 1 - (1 / (2n + 1)) for allele
frequencies, where n is the number of
individuals; without this flag, this bounds
are set to 0 and 1
--ksize,-k indicates the next argument is the number of
components; this value must be at least two
components; this value must be at least one
--max-iterations,-mi indicates the next argument is the maximum
number of iterations to execute the algorithm;
this value must be greater than or equal to
Expand Down Expand Up @@ -99,7 +104,7 @@ DESCRIPTION
[Notation]
K := Number of Components
This value must be greater than or equal to two.
This value must be greater than or equal to one.
I := Number of Individuals
This value must be greater than or equal to two.
Expand All @@ -123,7 +128,7 @@ DESCRIPTION
Q := [I x K] Proportion Matrix
This matrix consists of floating-point values ranging from 0 to 1. Each
row sums to 1.0, and each row contains more than one distinct value.
row sums to 1.0.
[Matrix File Format]
Expand Down
3 changes: 2 additions & 1 deletion src/qpas/jade.optimizer.hpp
View file
Expand Up @@ -54,6 +54,7 @@ namespace jade
const auto fg = settings.get_fg();
const auto fif = settings.get_fif();
const auto & g = settings.get_g();
const auto frb = opts.is_frb();

const stopwatch sw1;

Expand Down Expand Up @@ -87,7 +88,7 @@ namespace jade

if (!opts.is_fixed_f())
{
fa = improver_type::improve_f(g, q, fa, fb, qfa, qfb, fif);
fa = improver_type::improve_f(g, q, fa, fb, qfa, qfb, fif, frb);
_compute_fb(fa, fb);
matrix_type::gemm(q, fa, qfa);
matrix_type::gemm(q, fb, qfb);
Expand Down
23 changes: 23 additions & 0 deletions src/qpas/jade.options.hpp
View file
Expand Up @@ -57,6 +57,7 @@ namespace jade
, _qin (a.read<std::string>("--qin", "-qi"))
, _qout (a.read<std::string>("--qout", "-qo"))
, _seed (a.read("--seed", "-s", std::random_device()()))
, _frb (a.read_flag("--frequency-bounds", "-frb"))
, _fixed_f (a.read_flag("--fixed-f", "-ff"))
, _fixed_q (a.read_flag("--fixed-q", "-fq"))
, _quiet (a.read_flag("--quiet", "-q"))
Expand Down Expand Up @@ -103,6 +104,19 @@ namespace jade
throw error() << "invalid specification of --fixed-q option "
<< "and --force option";

if (is_frb())
{
if (is_fin_force_specified())
throw error()
<< "invalid specification of --fin-force "
<< "and --frequency-bounds options";

if (is_fixed_f() && is_fin_specified())
throw error()
<< "invalid specification of --fixed-f, --fin, "
<< "and --frequency-bounds options";
}

if (_num_threads == 0)
throw error() << "invalid number of threads specified for "
<< "--num-threads option: " << _num_threads;
Expand Down Expand Up @@ -222,6 +236,14 @@ namespace jade
return !std::isnan(_epsilon);
}

///
/// \return True if the frequency-bounds option is specified.
///
inline bool is_frb() const
{
return _frb;
}

///
/// \return True if the fin option is specified.
///
Expand Down Expand Up @@ -334,6 +356,7 @@ namespace jade
const seed_type _seed;

// options without arguments
const bool _frb;
const bool _fixed_f;
const bool _fixed_q;
const bool _quiet;
Expand Down

0 comments on commit 84361a0

Please sign in to comment.