Changeset 1324 for applications/robust/robustlib.h
- Timestamp:
- 04/07/11 18:33:51 (13 years ago)
- Files:
-
- 1 modified
Legend:
- Unmodified
- Added
- Removed
-
applications/robust/robustlib.h
r1320 r1324 31 31 class emlig; 32 32 33 33 34 /* 34 35 class t_simplex … … 61 62 this->value = value; 62 63 multiplicity = 1; 64 } 65 }; 66 67 class simplex 68 { 69 public: 70 71 set<vertex*> vertices; 72 73 double probability; 74 75 vector<list<pair<double,double>>> gamma_parameters; 76 77 simplex(set<vertex*> vertices) 78 { 79 this->vertices.insert(vertices.begin(),vertices.end()); 80 probability = 0; 81 } 82 83 simplex(vertex* vertex) 84 { 85 this->vertices.insert(vertex); 86 probability = 0; 63 87 } 64 88 }; … … 164 188 165 189 /// List of triangulation polyhedrons of the polyhedron given by their relative vertices. 166 map<double,set<vertex*>> triangulation;190 set<simplex*> triangulation; 167 191 168 192 /// A list of relative addresses serving for Hasse diagram construction. … … 275 299 vertices.insert(this); 276 300 277 set<vertex*> vert_simplex; 278 279 vert_simplex.insert(this); 280 281 triangulation.insert(pair<double,set<vertex*>>(0,vert_simplex)); 301 simplex* vert_simplex = new simplex(vertices); 302 303 triangulation.insert(vert_simplex); 282 304 } 283 305 … … 330 352 } 331 353 332 double integrate_simplex(s et<vertex*>simplex, char c);354 double integrate_simplex(simplex* simplex, char c); 333 355 334 356 }; 357 335 358 336 359 … … 774 797 cur_par_toprow->probability = 0.0; 775 798 776 map<double,set<vertex*>> new_triangulation;777 778 for( map<double,set<vertex*>>::iterator t_ref = current_parent->triangulation.begin();t_ref!=current_parent->triangulation.end();t_ref++)799 //set<simplex*> new_triangulation; 800 801 for(set<simplex*>::iterator s_ref = current_parent->triangulation.begin();s_ref!=current_parent->triangulation.end();s_ref++) 779 802 { 780 double cur_prob = cur_par_toprow->integrate_simplex((* t_ref).second,'C');803 double cur_prob = cur_par_toprow->integrate_simplex((*s_ref),'C'); 781 804 782 805 cur_par_toprow->probability += cur_prob; 783 806 784 new_triangulation.insert(pair<double,set<vertex*>>(cur_prob,(*t_ref).second));807 //new_triangulation.insert(pair<double,set<vertex*>>(cur_prob,(*t_ref).second)); 785 808 } 786 809 787 current_parent->triangulation.clear();788 current_parent->triangulation.insert(new_triangulation.begin(),new_triangulation.end());810 //current_parent->triangulation.clear(); 811 //current_parent->triangulation.insert(new_triangulation.begin(),new_triangulation.end()); 789 812 } 790 813 } … … 879 902 cur_par_toprow->probability = 0.0; 880 903 881 map<double,set<vertex*>> new_triangulation;904 //map<double,set<vertex*>> new_triangulation; 882 905 883 for( map<double,set<vertex*>>::iterator t_ref = current_parent->triangulation.begin();t_ref!=current_parent->triangulation.end();t_ref++)906 for(set<simplex*>::iterator s_ref = current_parent->triangulation.begin();s_ref!=current_parent->triangulation.end();s_ref++) 884 907 { 885 double cur_prob = cur_par_toprow->integrate_simplex((* t_ref).second,'C');908 double cur_prob = cur_par_toprow->integrate_simplex((*s_ref),'C'); 886 909 887 910 cur_par_toprow->probability += cur_prob; 888 911 889 new_triangulation.insert(pair<double,set<vertex*>>(cur_prob,(*t_ref).second));912 //new_triangulation.insert(pair<double,set<vertex*>>(cur_prob,(*t_ref).second)); 890 913 } 891 914 892 current_parent->triangulation.clear();893 current_parent->triangulation.insert(new_triangulation.begin(),new_triangulation.end());915 //current_parent->triangulation.clear(); 916 //current_parent->triangulation.insert(new_triangulation.begin(),new_triangulation.end()); 894 917 } 895 918 … … 960 983 for(int i = 0;i<statistic.size();i++) 961 984 { 962 int zero = 0;963 int one = 0;964 int two = 0;985 //int zero = 0; 986 //int one = 0; 987 //int two = 0; 965 988 966 989 for(polyhedron* horiz_ref = statistic.rows[i];horiz_ref!=statistic.get_end();horiz_ref=horiz_ref->next_poly) 967 990 { 991 992 993 if(i==statistic.size()-1) 994 { 995 cout << ((toprow*)horiz_ref)->condition_sum << " " << ((toprow*)horiz_ref)->probability << endl; 996 cout << "Order:" << ((toprow*)horiz_ref)->condition_order << endl; 997 } 968 998 969 999 /* 970 if(i==statistic.size()-1)971 {972 //cout << ((toprow*)horiz_ref)->condition_sum << " " << ((toprow*)horiz_ref)->probability << endl;973 cout << "Order:" << ((toprow*)horiz_ref)->condition_order << endl;974 }975 1000 if(i==0) 976 1001 { … … 979 1004 */ 980 1005 1006 /* 981 1007 char* string = "Checkpoint"; 982 1008 1009 983 1010 if((*horiz_ref).parentconditions.size()==0) 984 1011 { … … 993 1020 two++; 994 1021 } 1022 */ 995 1023 996 1024 } … … 1070 1098 void add_and_remove_condition(vec toadd, vec toremove) 1071 1099 { 1072 step_me(0);1100 //step_me(0); 1073 1101 1074 1102 likelihood_value = numeric_limits<double>::max(); … … 1712 1740 mat sample_mat(int n) 1713 1741 { 1714 mat sample_mat;1715 map<double,toprow*> ordered_toprows;1742 1743 1716 1744 1717 1745 /// \TODO tady je to spatne, tady nesmi byt conditions.size(), viz RARX.bayes() 1718 1746 if(conditions.size()-2-number_of_parameters>=0) 1719 1747 { 1720 double sum = 0; 1721 1722 for(polyhedron* top_ref = statistic.rows[number_of_parameters-1];top_ref!=statistic.end_poly;top_ref=top_ref->next_poly) 1748 mat sample_mat; 1749 map<double,toprow*> ordered_toprows; 1750 double sum_a = 0; 1751 1752 1753 for(polyhedron* top_ref = statistic.rows[number_of_parameters];top_ref!=statistic.end_poly;top_ref=top_ref->next_poly) 1723 1754 { 1724 1755 toprow* current_top = (toprow*)top_ref; 1725 1756 1726 sum+=current_top->probability; 1727 ordered_toprows.insert(pair<double,toprow*>(sum,current_top)); 1728 } 1757 sum_a+=current_top->probability; 1758 ordered_toprows.insert(pair<double,toprow*>(sum_a,current_top)); 1759 } 1760 1761 while(sample_mat.cols()<n) 1762 { 1763 double rnumber = randu()*sum_a; 1764 1765 // cout << "RND:" << rnumber << endl; 1766 1767 // This could be more efficient (log n), but map::upper_bound() doesn't let me dereference returned iterator 1768 toprow* sampled_toprow; 1769 for(map<double,toprow*>::iterator top_ref = ordered_toprows.begin();top_ref!=ordered_toprows.end();top_ref++) 1770 { 1771 // cout << "CDF:"<< (*top_ref).first << endl; 1772 1773 if((*top_ref).first >= rnumber) 1774 { 1775 sampled_toprow = (*top_ref).second; 1776 break; 1777 } 1778 } 1779 1780 // cout << &sampled_toprow << ";"; 1781 1782 rnumber = randu(); 1783 1784 set<simplex*>::iterator s_ref; 1785 double sum_b = 0; 1786 1787 for(s_ref = sampled_toprow->triangulation.begin();s_ref!=sampled_toprow->triangulation.end();s_ref++) 1788 { 1789 sum_b += (*s_ref)->probability; 1790 1791 if(sum_b/sampled_toprow->probability >= rnumber) 1792 break; 1793 } 1794 1795 // cout << &(*tri_ref) << endl; 1796 1797 int dimension = (*s_ref)->vertices.size()-1; 1798 1799 mat jacobian(dimension,dimension); 1800 vec gradient = sampled_toprow->condition_sum.right(dimension); 1801 1802 vertex* base_vert = *(*s_ref)->vertices.begin(); 1803 1804 int row_count = 0; 1805 1806 for(set<vertex*>::iterator vert_ref = ++(*s_ref)->vertices.begin();vert_ref!=(*s_ref)->vertices.end();vert_ref++) 1807 { 1808 jacobian.set_row(row_count,(*vert_ref)->get_coordinates()-base_vert->get_coordinates()); 1809 1810 row_count++; 1811 } 1812 1813 // cout << gradient << endl; 1814 // cout << jacobian << endl; 1815 1816 gradient = jacobian*gradient; 1817 1818 //cout << gradient << endl; 1819 1820 mat rotation_matrix = eye(dimension); 1821 1822 double squared_length = 1; 1823 1824 for(int i = 1;i<dimension;i++) 1825 { 1826 double t = gradient[i]/sqrt(squared_length); 1827 1828 double sin_theta = t/sqrt(1+pow(t,2)); 1829 double cos_theta = 1/sqrt(1+pow(t,2)); 1830 1831 mat partial_rotation = eye(dimension); 1832 1833 partial_rotation.set(0,0,cos_theta); 1834 partial_rotation.set(i,i,cos_theta); 1835 1836 partial_rotation.set(0,i,sin_theta); 1837 partial_rotation.set(i,0,-sin_theta); 1838 1839 rotation_matrix = rotation_matrix*partial_rotation; 1840 1841 squared_length += pow(gradient[i],2); 1842 } 1843 1844 // cout << rotation_matrix << endl; 1845 // cout << gradient << endl; 1846 1847 vec minima = numeric_limits<double>::max()*ones(number_of_parameters); 1848 vec maxima = -numeric_limits<double>::max()*ones(number_of_parameters); 1849 1850 vertex* min_grad; 1851 vertex* max_grad; 1852 1853 for(set<vertex*>::iterator vert_ref = (*s_ref)->vertices.begin();vert_ref!=(*s_ref)->vertices.end();vert_ref++) 1854 { 1855 vec new_coords = rotation_matrix*jacobian*(*vert_ref)->get_coordinates(); 1856 1857 // cout << new_coords << endl; 1858 1859 for(int j = 0;j<new_coords.size();j++) 1860 { 1861 if(new_coords[j]>maxima[j]) 1862 { 1863 maxima[j] = new_coords[j]; 1864 1865 if(j==0) 1866 { 1867 max_grad = (*vert_ref); 1868 } 1869 } 1870 1871 if(new_coords[j]<minima[j]) 1872 { 1873 minima[j] = new_coords[j]; 1874 1875 if(j==0) 1876 { 1877 min_grad = (*vert_ref); 1878 } 1879 } 1880 } 1881 } 1882 1883 vec sample_coordinates; 1884 1885 for(int j = 0;j<number_of_parameters;j++) 1886 { 1887 rnumber = randu(); 1888 1889 double coordinate; 1890 1891 if(j==0) 1892 { 1893 double pos_value = sampled_toprow->condition_sum.right(sampled_toprow->condition_sum.size()-1)*max_grad->get_coordinates()-sampled_toprow->condition_sum[0]; 1894 1895 1896 //coordinate = log(rnumber*(exp(sampled_toprow->condition_sum.size()-1)*max_grad->get_coordinates()-sampled_toprow->condition_sum[0])-exp(sampled_toprow->condition_sum.right(sampled_toprow->condition_sum.size()-1)*min_grad->get_coordinates()-sampled_toprow->condition_sum[0]))+ 1897 } 1898 } 1899 1900 // cout << sampled_toprow->condition_sum.right(sampled_toprow->condition_sum.size()-1)*min_grad->get_coordinates()-sampled_toprow->condition_sum[0] << endl; 1901 // cout << sampled_toprow->condition_sum.right(sampled_toprow->condition_sum.size()-1)*max_grad->get_coordinates()-sampled_toprow->condition_sum[0] << endl; 1902 1903 } 1904 1905 return sample_mat; 1729 1906 } 1730 1907 else 1731 1908 { 1732 1909 throw new exception("You are trying to sample from density that is not determined (parameters can't be integrated out)!"); 1733 } 1734 1735 while(sample_mat.cols()<n) 1736 { 1737 double rnumber = randu(); 1738 1739 toprow* sampled_toprow = ordered_toprows.upper_bound(rnumber)->second; 1740 1741 rnumber = randu(); 1742 1743 map<double,set<vertex*>>::iterator tri_ref = sampled_toprow->triangulation.begin() 1744 double sum = 0; 1745 1746 for(tri_ref = sampled_toprow->triangulation.begin();tri_ref!=sampled_toprow->triangulation.end();tri_ref++) 1747 { 1748 sum += (*tri_ref).first; 1749 1750 if(sum/sampled_toprow->probability >= rnumber) 1751 break; 1752 } 1753 1754 1755 } 1756 1757 return sample_mat; 1910 1911 return 0; 1912 } 1913 1914 1758 1915 } 1759 1916 … … 1920 2077 // This method guarantees that each polyhedron is already triangulated, therefore its triangulation 1921 2078 // is only one set of vertices and it is the set of all its vertices. 1922 set<vertex*> t_simplex1; 1923 set<vertex*> t_simplex2; 1924 1925 t_simplex1.insert(current_copy1->vertices.begin(),current_copy1->vertices.end()); 1926 t_simplex2.insert(current_copy2->vertices.begin(),current_copy2->vertices.end()); 2079 simplex* t_simplex1 = new simplex(current_copy1->vertices); 2080 simplex* t_simplex2 = new simplex(current_copy2->vertices); 1927 2081 1928 current_copy1->triangulation.insert( pair<double,set<vertex*>>(0,t_simplex1));1929 current_copy2->triangulation.insert( pair<double,set<vertex*>>(0,t_simplex2));2082 current_copy1->triangulation.insert(t_simplex1); 2083 current_copy2->triangulation.insert(t_simplex2); 1930 2084 1931 2085 // Now we have copied the polyhedron and we have to copy all of its relations. Because we are copying