DynProg.h

Go to the documentation of this file.

/*
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or
 * (at your option) any later version.
 *
 * Written (W) 1999-2009 Gunnar Raetsch
 * Written (W) 1999-2009 Soeren Sonnenburg
 * Written (W) 2008-2009 Jonas Behr
 * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society
 */

#ifndef __CDYNPROG_H__
#define __CDYNPROG_H__

#include <shogun/mathematics/Math.h>
#include <shogun/lib/common.h>
#include <shogun/base/SGObject.h>
#include <shogun/io/SGIO.h>
#include <shogun/lib/config.h>
#include <shogun/structure/PlifMatrix.h>
#include <shogun/structure/PlifBase.h>
#include <shogun/structure/Plif.h>
#include <shogun/structure/IntronList.h>
#include <shogun/structure/SegmentLoss.h>
#include <shogun/features/StringFeatures.h>
#include <shogun/features/SparseFeatures.h>
#include <shogun/distributions/Distribution.h>
#include <shogun/lib/DynamicArray.h>
#include <shogun/lib/DynamicObjectArray.h>
#include <shogun/lib/Time.h>

#include <stdio.h>
#include <limits.h>

namespace shogun
{
    template <class T> class CSparseFeatures;
    class CIntronList;
    class CPlifMatrix;
    class CSegmentLoss;

    template <class T> class CDynamicArray;

//#define DYNPROG_TIMING

#ifdef USE_BIGSTATES
typedef uint16_t T_STATES ;
#else
typedef uint8_t T_STATES ;
#endif
typedef T_STATES* P_STATES ;

#ifndef DOXYGEN_SHOULD_SKIP_THIS

struct segment_loss_struct
{
    int32_t maxlookback;
    int32_t seqlen;
    int32_t *segments_changed;
    float64_t *num_segment_id;
    int32_t *length_segment_id ;
};
#endif

class CDynProg : public CSGObject
{
public:
    CDynProg(int32_t p_num_svms=8);
    virtual ~CDynProg();

    // model related functions
    void set_num_states(int32_t N);

    int32_t get_num_states();

    int32_t get_num_svms();

    void init_content_svm_value_array(const int32_t p_num_svms);

    void init_tiling_data(int32_t* probe_pos, float64_t* intensities, const int32_t num_probes);

    void precompute_tiling_plifs(CPlif** PEN, const int32_t* tiling_plif_ids, const int32_t num_tiling_plifs);

    void resize_lin_feat(int32_t num_new_feat);
    void set_p_vector(SGVector<float64_t> p);

    void set_q_vector(SGVector<float64_t> q);

    void set_a(SGMatrix<float64_t> a);

    void set_a_id(SGMatrix<int32_t> a);

    void set_a_trans_matrix(SGMatrix<float64_t> a_trans);

    void init_mod_words_array(SGMatrix<int32_t> p_mod_words_array);

    bool check_svm_arrays();

    void set_observation_matrix(SGNDArray<float64_t> seq);

    int32_t get_num_positions();

    void set_content_type_array(SGMatrix<float64_t> seg_path);

    void set_pos(SGVector<int32_t> pos);

    void set_orf_info(SGMatrix<int32_t> orf_info);

    void set_gene_string(SGVector<char> genestr);


    void set_dict_weights(SGMatrix<float64_t> dictionary_weights);

    void best_path_set_segment_loss(SGMatrix<float64_t> segment_loss);

    void best_path_set_segment_ids_mask(int32_t* segment_ids, float64_t* segment_mask, int32_t m);

    void set_sparse_features(CSparseFeatures<float64_t>* seq_sparse1, CSparseFeatures<float64_t>* seq_sparse2);

    void set_plif_matrices(CPlifMatrix* pm);

    // best_path result retrieval functions
    SGVector<float64_t> get_scores();

    SGMatrix<int32_t> get_states();

    SGMatrix<int32_t> get_positions();


    void compute_nbest_paths(int32_t max_num_signals,
                         bool use_orf, int16_t nbest, bool with_loss, bool with_multiple_sequences);


    void best_path_trans_deriv(
            int32_t* my_state_seq, int32_t *my_pos_seq,
            int32_t my_seq_len, const float64_t *seq_array, int32_t max_num_signals);

    // additional best_path_trans_deriv functions
    void set_my_state_seq(int32_t* my_state_seq);

    void set_my_pos_seq(int32_t* my_pos_seq);

    void get_path_scores(float64_t** my_scores, int32_t* seq_len);

    void get_path_losses(float64_t** my_losses, int32_t* seq_len);


    inline T_STATES get_N() const
    {
        return m_N ;
    }

    inline void set_q(T_STATES offset, float64_t value)
    {
        m_end_state_distribution_q[offset]=value;
    }

    inline void set_p(T_STATES offset, float64_t value)
    {
        m_initial_state_distribution_p[offset]=value;
    }

    inline void set_a(T_STATES line_, T_STATES column, float64_t value)
    {
      m_transition_matrix_a.element(line_,column)=value; // look also best_path!
    }

    inline float64_t get_q(T_STATES offset) const
    {
        return m_end_state_distribution_q[offset];
    }

    inline float64_t get_q_deriv(T_STATES offset) const
    {
        return m_end_state_distribution_q_deriv[offset];
    }

    inline float64_t get_p(T_STATES offset) const
    {
        return m_initial_state_distribution_p[offset];
    }

    inline float64_t get_p_deriv(T_STATES offset) const
    {
        return m_initial_state_distribution_p_deriv[offset];
    }

    void precompute_content_values();

    inline float64_t* get_lin_feat(int32_t & dim1, int32_t & dim2)
    {
        m_lin_feat.get_array_size(dim1, dim2);
        return m_lin_feat.get_array();
    }
    inline void set_lin_feat(float64_t* p_lin_feat, int32_t p_num_svms, int32_t p_seq_len)
    {
      m_lin_feat.set_array(p_lin_feat, p_num_svms, p_seq_len, true, true);
    }
    void create_word_string();

    void precompute_stop_codons();

    inline float64_t get_a(T_STATES line_, T_STATES column) const
    {
      return m_transition_matrix_a.element(line_, column); // look also best_path()!
    }

    inline float64_t get_a_deriv(T_STATES line_, T_STATES column) const
    {
      return m_transition_matrix_a_deriv.element(line_, column); // look also best_path()!
    }

    void set_intron_list(CIntronList* intron_list, int32_t num_plifs);

    CSegmentLoss* get_segment_loss_object()
    {
        return m_seg_loss_obj;
    }

    void long_transition_settings(bool use_long_transitions, int32_t threshold, int32_t max_len)
    {
        m_long_transitions = use_long_transitions;
        m_long_transition_threshold = threshold;
        SG_DEBUG("ignoring max_len\n") ;
        //m_long_transition_max = max_len;
    }

protected:

    /* helper functions */

    void lookup_content_svm_values(const int32_t from_state,
        const int32_t to_state, const int32_t from_pos, const int32_t to_pos,
        float64_t* svm_values, int32_t frame);

    inline void lookup_tiling_plif_values(const int32_t from_state,
        const int32_t to_state, const int32_t len, float64_t* svm_values);

    inline int32_t find_frame(const int32_t from_state);

    inline int32_t raw_intensities_interval_query(
        const int32_t from_pos, const int32_t to_pos, float64_t* intensities, int32_t type);

#ifndef DOXYGEN_SHOULD_SKIP_THIS

    struct svm_values_struct
    {
        int32_t maxlookback;
        int32_t seqlen;

        int32_t* start_pos;
        float64_t ** svm_values_unnormalized;
        float64_t * svm_values;
        bool *** word_used;
        int32_t **num_unique_words;
    };
#endif // DOXYGEN_SHOULD_SKIP_THIS

    bool extend_orf(int32_t orf_from, int32_t orf_to, int32_t start, int32_t &last_pos, int32_t to);

    virtual const char* get_name() const { return "DynProg"; }

private:

    T_STATES trans_list_len;
    T_STATES **trans_list_forward;
    T_STATES *trans_list_forward_cnt;
    float64_t **trans_list_forward_val;
    int32_t **trans_list_forward_id;
    bool mem_initialized;

#ifdef DYNPROG_TIMING
    CTime MyTime;
    CTime MyTime2;
    CTime MyTime3;

    float64_t segment_init_time;
    float64_t segment_pos_time;
    float64_t segment_clean_time;
    float64_t segment_extend_time;
    float64_t orf_time;
    float64_t content_time;
    float64_t content_penalty_time;
    float64_t content_svm_values_time ;
    float64_t content_plifs_time ;
    float64_t svm_init_time;
    float64_t svm_pos_time;
    float64_t inner_loop_time;
    float64_t inner_loop_max_time ;
    float64_t svm_clean_time;
    float64_t long_transition_time ;
#endif


protected:

    int32_t m_N;

    CDynamicArray<int32_t> m_transition_matrix_a_id; // 2d
    CDynamicArray<float64_t> m_transition_matrix_a; // 2d
    CDynamicArray<float64_t> m_transition_matrix_a_deriv; // 2d

    CDynamicArray<float64_t> m_initial_state_distribution_p;
    CDynamicArray<float64_t> m_initial_state_distribution_p_deriv;

    CDynamicArray<float64_t> m_end_state_distribution_q;
    CDynamicArray<float64_t> m_end_state_distribution_q_deriv;


    int32_t m_num_degrees;
    int32_t m_num_svms;

    CDynamicArray<int32_t> m_word_degree;
    CDynamicArray<int32_t> m_cum_num_words;
    int32_t * m_cum_num_words_array;
    CDynamicArray<int32_t> m_num_words;
    int32_t* m_num_words_array;
    CDynamicArray<int32_t> m_mod_words; // 2d
    int32_t* m_mod_words_array;
    CDynamicArray<bool> m_sign_words;
    bool* m_sign_words_array;
    CDynamicArray<int32_t> m_string_words;
    int32_t* m_string_words_array;

//  CDynamicArray<int32_t> m_svm_pos_start;
    CDynamicArray<int32_t> m_num_unique_words;
    bool m_svm_arrays_clean;
    int32_t m_max_a_id;

    // input arguments
    CDynamicArray<float64_t> m_observation_matrix; //3d
    CDynamicArray<int32_t> m_pos;
    int32_t m_seq_len;
    CDynamicArray<int32_t> m_orf_info; // 2d
    CDynamicArray<float64_t> m_segment_sum_weights; // 2d
    CDynamicObjectArray m_plif_list; // CPlifBase*
    CDynamicObjectArray m_PEN; // 2d, CPlifBase*
    CDynamicObjectArray m_PEN_state_signals; // 2d, CPlifBase*
    CDynamicArray<char> m_genestr;
    uint16_t*** m_wordstr;
    CDynamicArray<float64_t> m_dict_weights; // 2d
    CDynamicArray<float64_t> m_segment_loss; // 3d
    CDynamicArray<int32_t> m_segment_ids;   
    CDynamicArray<float64_t> m_segment_mask;    
    CDynamicArray<int32_t> m_my_state_seq;
    CDynamicArray<int32_t> m_my_pos_seq;
    CDynamicArray<float64_t> m_my_scores;
    CDynamicArray<float64_t> m_my_losses;

    CSegmentLoss* m_seg_loss_obj;

    // output arguments
    CDynamicArray<float64_t> m_scores;
    CDynamicArray<int32_t> m_states; // 2d
    CDynamicArray<int32_t> m_positions; // 2d

    CSparseFeatures<float64_t>* m_seq_sparse1;
    CSparseFeatures<float64_t>* m_seq_sparse2;
    CPlifMatrix* m_plif_matrices;

    CDynamicArray<bool> m_genestr_stop;

    CIntronList* m_intron_list;

    int32_t m_num_intron_plifs;

    CDynamicArray<float64_t> m_lin_feat; // 2d

    float64_t *m_raw_intensities;
    int32_t* m_probe_pos;
    int32_t* m_num_probes_cum;
    int32_t* m_num_lin_feat_plifs_cum;
    int32_t m_num_raw_data;

    bool m_long_transitions ;
    int32_t m_long_transition_threshold  ;
    //int32_t m_long_transition_max ;

    static int32_t word_degree_default[4];

    static int32_t cum_num_words_default[5];

    static int32_t frame_plifs[3];

    static int32_t num_words_default[4];

    static int32_t mod_words_default[32];

    static bool sign_words_default[16];

    static int32_t string_words_default[16];
};
}
#endif