Skip to content

Commit

Permalink
Finalized the t-walk implementation.
Browse files Browse the repository at this point in the history
  • Loading branch information
@robbert-harms
robbert-harms committed Oct 24, 2018
1 parent 1780454 commit f4a0425
Show file tree
Hide file tree
Showing 2 changed files with 126 additions and 71 deletions.
59 changes: 57 additions & 2 deletions mot/optimize/__init__.py
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
from mot.lib.cl_function import SimpleCLFunction
from mot.configuration import CLRuntimeInfo
from mot.lib.kernel_data import Array, Zeros, LocalMemory
from mot.lib.kernel_data import Array
from mot.library_functions import Powell, Subplex, NMSimplex, LevenbergMarquardt
from mot.optimize.base import OptimizeResults

Expand All @@ -13,7 +13,9 @@

def minimize(func, x0, data=None, method=None, nmr_observations=None, cl_runtime_info=None, options=None,
jacobian_func=None):
"""Minimization of scalar function of one or more variables.
"""Minimization of one or more variables.
For an easy wrapper of function maximization, see :func:`maximize`.
Args:
func (mot.lib.cl_function.CLFunction): A CL function with the signature:
Expand Down Expand Up @@ -93,6 +95,59 @@ def minimize(func, x0, data=None, method=None, nmr_observations=None, cl_runtime
raise ValueError('Could not find the specified method "{}".'.format(method))


def maximize(func, x0, nmr_observations, **kwargs):
"""Maximization of a function.
This wraps the objective function to take the negative of the computed values and passes it then on to one
of the minimization routines.
Args:
func (mot.lib.cl_function.CLFunction): A CL function with the signature:
.. code-block:: c
double <func_name>(local const mot_float_type* const x,
void* data,
local mot_float_type* objective_list);
The objective list needs to be filled when the provided pointer is not null. It should contain
the function values for each observation. This list is used by non-linear least-squares routines,
and will be squared by the least-square optimizer. This is only used by the ``Levenberg-Marquardt`` routine.
x0 (ndarray): Initial guess. Array of real elements of size (n, p), for 'n' problems and 'p'
independent variables.
nmr_observations (int): the number of observations returned by the optimization function.
**kwargs: see :func:`minimize`.
"""
wrapped_func = SimpleCLFunction.from_string('''
double _negate_''' + func.get_cl_function_name() + '''(
local mot_float_type* x,
void* data,
local mot_float_type* objective_list){
double return_val = ''' + func.get_cl_function_name() + '''(x, data, objective_list);
if(objective_list){
const uint nmr_observations = ''' + str(nmr_observations) + ''';
uint local_id = get_local_id(0);
uint workgroup_size = get_local_size(0);
uint observation_ind;
for(uint i = 0; i < (nmr_observations + workgroup_size - 1) / workgroup_size; i++){
observation_ind = i * workgroup_size + local_id;
if(observation_ind < nmr_observations){
objective_list[observation_ind] *= -1;
}
}
}
return -return_val;
}
''', dependencies=[func])
kwargs['nmr_observations'] = nmr_observations
return minimize(wrapped_func, x0, **kwargs)


def get_minimizer_options(method):
"""Return a dictionary with the default options for the given minimization method.
Expand Down
138 changes: 69 additions & 69 deletions mot/sample/t_walk.py
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,3 @@
import warnings

import numpy as np

from mot.lib.cl_function import SimpleCLFunction
Expand Down Expand Up @@ -69,8 +67,6 @@ def __init__(self, ll_func, log_prior_func, x0, x1, finalize_proposal_func=None,
"""
super().__init__(ll_func, log_prior_func, x0, **kwargs)

warnings.warn('The t-walk sampler is released as testing version only. It may not yet function 100% properly.')

self._subset_size = subset_size
self._param_choose_prob = min(self._nmr_params, self._subset_size) / (float(self._nmr_params))
self._walk_scale = walk_scale
Expand All @@ -96,6 +92,8 @@ def _get_mcmc_method_kernel_data(self):
'x1_position': Array(self._x1, 'mot_float_type', mode='rw', ensure_zero_copy=True),
'x1_log_likelihood': Array(self._x1_log_likelihood, 'mot_float_type', mode='rw', ensure_zero_copy=True),
'x1_log_prior': Array(self._x1_log_prior, 'mot_float_type', mode='rw', ensure_zero_copy=True),
'scratch_mft': LocalMemory('mot_float_type', self._nmr_params + 2),
'scratch_int': LocalMemory('int', self._nmr_params + 4),
}, '_twalk_data')

def _get_state_update_cl_func(self, nmr_samples, thinning, return_output):
Expand All @@ -105,23 +103,26 @@ def _get_state_update_cl_func(self, nmr_samples, thinning, return_output):
void* data,
ulong current_iteration,
void* rng_data,
char chain_ind,
int chain_ind,
global mot_float_type* current_position,
global mot_float_type* current_log_likelihood,
global mot_float_type* current_log_prior){
_twalk_data* twalk_data = (_twalk_data*)method_data;
bool is_first_work_item = get_local_id(0) == 0;
local char kernel_ind;
local bool params_selector[''' + str(self._nmr_params) + '''];
local uint nmr_params_selected;
local int* kernel_ind = twalk_data->scratch_int + 1;
local int* nmr_params_selected = twalk_data->scratch_int + 2;
local int* proposal_accepted = twalk_data->scratch_int + 3;
local int* params_selector = twalk_data->scratch_int + 4;
local mot_float_type proposal[''' + str(self._nmr_params) + '''];
local mot_float_type proposal_ll;
local mot_float_type proposal_lprior;
local bool proposal_accepted;
proposal_accepted = false;
local mot_float_type* proposal_ll = twalk_data->scratch_mft;
local mot_float_type* proposal_lprior = twalk_data->scratch_mft + 1;
local mot_float_type* proposal = twalk_data->scratch_mft + 2;
*proposal_accepted = false;
*proposal_ll = 0;
*proposal_lprior = 0;
global mot_float_type* main_chain;
global mot_float_type* helper_chain;
Expand All @@ -132,22 +133,22 @@ def _get_state_update_cl_func(self, nmr_samples, thinning, return_output):
float r = frand(rng_data);
if(r < ''' + str(self._move_probabilities[0]) + '''){
kernel_ind = 0;
*kernel_ind = 0;
}
else if(r < ''' + str(sum(self._move_probabilities[:2])) + '''){
kernel_ind = 1;
*kernel_ind = 1;
}
else if(r < ''' + str(sum(self._move_probabilities[:3])) + '''){
kernel_ind = 2;
*kernel_ind = 2;
}
else{
kernel_ind = 3;
*kernel_ind = 3;
}
nmr_params_selected = 0;
*nmr_params_selected = 0;
for(uint i = 0; i < ''' + str(self._nmr_params) + '''; i++){
params_selector[i] = frand(rng_data) < ''' + str(self._param_choose_prob) + ''';
nmr_params_selected++;
(*nmr_params_selected)++;
}
}
barrier(CLK_LOCAL_MEM_FENCE);
Expand All @@ -165,41 +166,41 @@ def _get_state_update_cl_func(self, nmr_samples, thinning, return_output):
main_lprior = current_log_prior;
}
if(kernel_ind == 0){
if(*kernel_ind == 0){
_twalk_walk_move(main_chain, helper_chain, main_ll, main_lprior,
proposal, &proposal_ll, &proposal_lprior, &proposal_accepted,
params_selector, nmr_params_selected, data, rng_data);
proposal, proposal_ll, proposal_lprior, proposal_accepted,
params_selector, *nmr_params_selected, data, rng_data);
}
else if(kernel_ind == 1){
else if(*kernel_ind == 1){
_twalk_traverse_move(main_chain, helper_chain, main_ll, main_lprior,
proposal, &proposal_ll, &proposal_lprior, &proposal_accepted,
params_selector, nmr_params_selected, data, rng_data);
proposal, proposal_ll, proposal_lprior, proposal_accepted,
params_selector, *nmr_params_selected, data, rng_data);
}
else if(kernel_ind == 2){
else if(*kernel_ind == 2){
_twalk_hop_move(main_chain, helper_chain, main_ll, main_lprior,
proposal, &proposal_ll, &proposal_lprior, &proposal_accepted,
params_selector, nmr_params_selected, data, rng_data);
proposal, proposal_ll, proposal_lprior, proposal_accepted,
params_selector, *nmr_params_selected, data, rng_data);
}
else{
_twalk_blow_move(main_chain, helper_chain, main_ll, main_lprior,
proposal, &proposal_ll, &proposal_lprior, &proposal_accepted,
params_selector, nmr_params_selected, data, rng_data);
proposal, proposal_ll, proposal_lprior, proposal_accepted,
params_selector, *nmr_params_selected, data, rng_data);
}
if(is_first_work_item && proposal_accepted){
if(is_first_work_item && *proposal_accepted){
if(chain_ind == 0){
for(uint k = 0; k < ''' + str(self._nmr_params) + '''; k++){
twalk_data->x1_position[k] = proposal[k];
}
*twalk_data->x1_log_likelihood = proposal_ll;
*twalk_data->x1_log_prior = proposal_lprior;
*twalk_data->x1_log_likelihood = *proposal_ll;
*twalk_data->x1_log_prior = *proposal_lprior;
}
else{
for(uint k = 0; k < ''' + str(self._nmr_params) + '''; k++){
current_position[k] = proposal[k];
}
*current_log_likelihood = proposal_ll;
*current_log_prior = proposal_lprior;
*current_log_likelihood = *proposal_ll;
*current_log_prior = *proposal_lprior;
}
}
barrier(CLK_LOCAL_MEM_FENCE);
Expand All @@ -214,21 +215,20 @@ def _get_state_update_cl_func(self, nmr_samples, thinning, return_output):
global mot_float_type* current_log_likelihood,
global mot_float_type* current_log_prior){
local char chain_ind;
if(get_local_id(0) == 0){
chain_ind = (char)(frand(rng_data) > 0.5);
}
barrier(CLK_LOCAL_MEM_FENCE);
_twalk_data* twalk_data = (_twalk_data*)method_data;
do{
_twalk_advance_chain(method_data, data, current_iteration, rng_data, chain_ind,
current_position, current_log_likelihood, current_log_prior);
local int* chain_ind = twalk_data->scratch_int;
*chain_ind = 0;
while(*chain_ind != 1){
if(get_local_id(0) == 0){
chain_ind = (char)(frand(rng_data) > 0.5);
*chain_ind = frand(rng_data) > 0.5;
}
barrier(CLK_LOCAL_MEM_FENCE);
} while(chain_ind != 1);
_twalk_advance_chain(method_data, data, current_iteration, rng_data, *chain_ind,
current_position, current_log_likelihood, current_log_prior);
}
}
''', dependencies=[self._finalize_proposal_func, self._get_walk_move(), self._get_traverse_move(),
Expand All @@ -241,7 +241,7 @@ def _get_walk_move(self):
global mot_float_type* main_chain,
global mot_float_type* helper_chain,
local mot_float_type* proposal,
local bool* params_selector,
local int* params_selector,
void* rng_data){
float u, a;
Expand All @@ -268,9 +268,9 @@ def _get_walk_move(self):
local mot_float_type* proposal,
local mot_float_type* proposal_ll,
local mot_float_type* proposal_lprior,
local bool* proposal_accepted,
local bool* params_selector,
uint nmr_params_selected,
local int* proposal_accepted,
local int* params_selector,
int nmr_params_selected,
void* data,
void* rng_data){
Expand Down Expand Up @@ -312,7 +312,7 @@ def _get_traverse_move(self):
global mot_float_type* main_chain,
global mot_float_type* helper_chain,
local mot_float_type* proposal,
local bool* params_selector,
local int* params_selector,
mot_float_type beta,
void* rng_data){
Expand All @@ -334,9 +334,9 @@ def _get_traverse_move(self):
local mot_float_type* proposal,
local mot_float_type* proposal_ll,
local mot_float_type* proposal_lprior,
local bool* proposal_accepted,
local bool* params_selector,
uint nmr_params_selected,
local int* proposal_accepted,
local int* params_selector,
int nmr_params_selected,
void* data,
void* rng_data){
Expand Down Expand Up @@ -377,7 +377,7 @@ def _get_hop_move(self):
global mot_float_type* main_chain,
global mot_float_type* helper_chain,
local mot_float_type* proposal,
local bool* params_selector,
local int* params_selector,
void* rng_data){
mot_float_type sigma = 0;
Expand All @@ -400,8 +400,8 @@ def _get_hop_move(self):
global mot_float_type* proposal,
local mot_float_type* main_chain,
global mot_float_type* helper_chain,
local bool* params_selector,
uint nmr_params_selected){
local int* params_selector,
int nmr_params_selected){
mot_float_type sigma = 0;
double sum = 0;
Expand All @@ -425,8 +425,8 @@ def _get_hop_move(self):
local mot_float_type* proposal,
global mot_float_type* main_chain,
global mot_float_type* helper_chain,
local bool* params_selector,
uint nmr_params_selected){
local int* params_selector,
int nmr_params_selected){
mot_float_type sigma = 0;
double sum = 0;
Expand All @@ -453,9 +453,9 @@ def _get_hop_move(self):
local mot_float_type* proposal,
local mot_float_type* proposal_ll,
local mot_float_type* proposal_lprior,
local bool* proposal_accepted,
local bool* params_selector,
uint nmr_params_selected,
local int* proposal_accepted,
local int* params_selector,
int nmr_params_selected,
void* data,
void* rng_data){
Expand Down Expand Up @@ -493,7 +493,7 @@ def _get_blow_move(self):
global mot_float_type* main_chain,
global mot_float_type* helper_chain,
local mot_float_type* proposal,
local bool* params_selector,
local int* params_selector,
void* rng_data){
mot_float_type sigma = 0;
Expand All @@ -515,8 +515,8 @@ def _get_blow_move(self):
global mot_float_type* proposal,
local mot_float_type* main_chain,
global mot_float_type* helper_chain,
local bool* params_selector,
uint nmr_params_selected){
local int* params_selector,
int nmr_params_selected){
mot_float_type sigma = 0;
double sum = 0;
Expand All @@ -537,8 +537,8 @@ def _get_blow_move(self):
local mot_float_type* proposal,
global mot_float_type* main_chain,
global mot_float_type* helper_chain,
local bool* params_selector,
uint nmr_params_selected){
local int* params_selector,
int nmr_params_selected){
mot_float_type sigma = 0;
double sum = 0;
Expand All @@ -563,9 +563,9 @@ def _get_blow_move(self):
local mot_float_type* proposal,
local mot_float_type* proposal_ll,
local mot_float_type* proposal_lprior,
local bool* proposal_accepted,
local bool* params_selector,
uint nmr_params_selected,
local int* proposal_accepted,
local int* params_selector,
int nmr_params_selected,
void* data,
void* rng_data){
Expand Down

0 comments on commit f4a0425

Please sign in to comment.