1 //
2 //  permutation.c
3 //  Symmetry
4 //
5 //  Created by Marcus Johansson on 01/02/15.
6 //  Copyright (c) 2015 Marcus Johansson.
7 //
8 //  Distributed under the MIT License ( See LICENSE file or copy at http://opensource.org/licenses/MIT )
9 //
10 
11 #include <stdlib.h>
12 #include <string.h>
13 
14 #include "msym.h"
15 #include "permutation.h"
16 #include "linalg.h"
17 
18 #include "debug.h"
19 
20 msym_error_t setPermutationCycles(msym_permutation_t *perm);
21 
freePermutationData(msym_permutation_t * perm)22 void freePermutationData(msym_permutation_t *perm){
23     if(perm != NULL){
24         free(perm->c);
25         free(perm->p);
26     }
27 }
28 
findPermutation(msym_symmetry_operation_t * sop,int l,double (* v[l])[3],msym_thresholds_t * t,msym_permutation_t * perm)29 msym_error_t findPermutation(msym_symmetry_operation_t *sop, int l, double (*v[l])[3], msym_thresholds_t *t, msym_permutation_t *perm){
30     msym_error_t ret = MSYM_SUCCESS;
31     double m[3][3];
32     symmetryOperationMatrix(sop, m);
33 
34     perm->p = malloc(sizeof(int[l]));
35     memset(perm->p, -1, sizeof(int[l])); //TODO: 2s complement
36     perm->p_length = l;
37 
38     for(int i = 0; i < l;i++){
39         int j;
40         double r[3];
41         mvmul(*v[i], m, r);
42         for(j = 0;j < l;j++){
43             if(vequal(r, *v[j],t->permutation)){
44                 perm->p[i] = j;
45                 break;
46             }
47         }
48         if(j == l) {
49             char buf[16];
50             symmetryOperationName(sop, sizeof(buf), buf);
51             msymSetErrorDetails("Unable to determine permutation for symmetry operation %s",buf);
52             ret = MSYM_PERMUTATION_ERROR;
53             goto err;
54         }
55     }
56     if(MSYM_SUCCESS != (ret = setPermutationCycles(perm))) goto err;
57 
58     return ret;
59 
60 err:
61     free(perm->p);
62     return ret;
63 }
64 
65 
66 typedef struct _perm_subgroup {
67     int sopsl;
68     int *sops;
69     int subgroup[2];
70 } perm_subgroup_t;
71 
72 
73 
findPermutationSubgroups(int l,msym_permutation_t perm[l],int sgmax,msym_symmetry_operation_t * sops,int * subgroupl,msym_subgroup_t ** subgroup)74 msym_error_t findPermutationSubgroups(int l, msym_permutation_t perm[l], int sgmax, msym_symmetry_operation_t *sops, int *subgroupl, msym_subgroup_t **subgroup){
75     msym_error_t ret = MSYM_SUCCESS;
76     perm_subgroup_t *group = calloc(l, sizeof(perm_subgroup_t));
77 
78     int *isops = malloc(sizeof(int[l]));
79     int *msops = malloc(sizeof(int[l]));
80     int gl = 0, glm = 0;
81 
82     for(int i = 0;i < l;i++){
83         if((sops[i].power == 1 && (sops[i].type == PROPER_ROTATION || sops[i].type == IMPROPER_ROTATION)) || sops[i].type == INVERSION || sops[i].type == REFLECTION){
84             msym_permutation_cycle_t* c = perm[i].c;
85             glm = gl;
86             memset(msops, 0, sizeof(int[l]));
87             group[gl].sopsl = c->l;
88             group[gl].sops = realloc(group[gl].sops, sizeof(int[c->l]));
89             memset(group[gl].sops,0,sizeof(int[c->l]));
90             //group[gl].sops = calloc(c->l, sizeof(int));
91             group[gl].subgroup[0] = group[gl].subgroup[1] = -1;
92             for(int next = c->s, j = 0;j < c->l;j++){
93                 msops[next] = 1;
94                 group[gl].sops[j] = next;
95                 next = perm[i].p[next];
96             }
97 
98             int n = 0;
99             for(int k = 0;k < l && n < l;k++){
100                 if(msops[k]){
101                     group[gl].sops[n] = k;
102                     n++;
103                 }
104             }
105             gl += n < l;
106         }
107     }
108 
109     if(glm == gl){
110         free(group[gl].sops);
111     }
112 
113     for(int i = 0;i < gl && gl < sgmax;i++){
114         for(int j = i+1;j < gl && gl < sgmax;j++){
115             int minl = group[i].sopsl < group[j].sopsl ? group[i].sopsl : group[j].sopsl;
116             if(0 == memcmp(group[i].sops,group[j].sops,sizeof(int[minl]))) continue;
117 
118             int n = 0;
119             memset(isops, 0, sizeof(int[l]));
120             memset(msops, 0, sizeof(int[l]));
121 
122             for(int k = 0;k < group[i].sopsl;k++){
123                 int s = group[i].sops[k];
124                 msops[s] = 1;
125                 isops[k] = s;
126             }
127             n = group[i].sopsl;
128             for(int k = 0;k < group[j].sopsl;k++){
129                 int s = group[j].sops[k];
130                 if(msops[s] == 0){
131                     msops[s] = 1;
132                     isops[n] = s;
133                     n++;
134                 }
135             }
136             for(int p = 0;p < n && n < l;p++){
137                 for(int q = 0;q < n && n < l;q++){
138                     int next = perm[isops[p]].p[isops[q]];
139                     if(msops[next] == 0){
140                         msops[next] = 1;
141                         isops[n] = next;
142                         n++;
143                     }
144                 }
145             }
146 
147             if(n < l && n > 1) {
148                 n = 0;
149                 memset(isops, 0, sizeof(int[l]));
150 
151                 for(int k = 0;k < l;k++){
152                     if(msops[k]){
153                         isops[n] = k;
154                         n++;
155                     }
156                 }
157                 int f;
158                 for(f = 0;f < gl;f++){
159                     if(group[f].sopsl == n && 0 == memcmp(group[f].sops, isops, sizeof(int[n]))){
160                         break;
161                     }
162                 }
163                 if(f == gl){
164                     group = realloc(group, sizeof(perm_subgroup_t[gl+1]));
165                     group[gl].sopsl = n;
166                     group[gl].sops = malloc(sizeof(int[n]));
167                     memcpy(group[gl].sops, isops, sizeof(int[n]));
168                     group[gl].subgroup[0] = i;
169                     group[gl].subgroup[1] = j;
170                     gl++;
171                 }
172             }
173         }
174     }
175 
176     msym_subgroup_t *mgroup = calloc(gl, sizeof(msym_subgroup_t));
177     for(int i = 0;i < gl;i++){
178         mgroup[i].sops = calloc(group[i].sopsl, sizeof(msym_symmetry_operation_t *));
179         mgroup[i].order = group[i].sopsl;
180         mgroup[i].generators[0] = group[i].subgroup[0] < 0 ? NULL : &mgroup[group[i].subgroup[0]];
181         mgroup[i].generators[1] = group[i].subgroup[1] < 0 ? NULL : &mgroup[group[i].subgroup[1]];
182 
183         for(int j = 0;j < group[i].sopsl;j++){
184             mgroup[i].sops[j] = &sops[group[i].sops[j]];
185         }
186     }
187 
188     *subgroup = mgroup;
189     *subgroupl = gl;
190 
191 //err:
192     for(int i = 0;i < gl;i++){
193         free(group[i].sops);
194     }
195     free(group);
196     free(isops);
197     free(msops);
198     return ret;
199 }
200 
findSymmetryOperationPermutations(int l,msym_symmetry_operation_t sops[l],msym_thresholds_t * t,msym_permutation_t ** rperm)201 msym_error_t findSymmetryOperationPermutations(int l, msym_symmetry_operation_t sops[l], msym_thresholds_t *t, msym_permutation_t **rperm){
202 
203     msym_error_t ret = MSYM_SUCCESS;
204     //Don't block allocate this, it's a pain to keep track of the pointers
205     msym_permutation_t *permutations = malloc(sizeof(msym_permutation_t[l]));
206 
207     for(int i = 0; i < l;i++){
208         permutations[i].p = malloc(sizeof(int[l]));
209         memset(permutations[i].p, -1, sizeof(int[l]));
210         permutations[i].p_length = l;
211     }
212 
213     double (*msops)[3][3] = malloc(sizeof(double[l][3][3]));
214 
215     for(int i = 0; i < l;i++){
216         symmetryOperationMatrix(&sops[i], msops[i]);
217     }
218 
219     for(int i = 0; i < l;i++){
220         if((sops[i].type == PROPER_ROTATION && sops[i].order == 0) || sops[i].type == IDENTITY){
221             for(int j = 0;j < l;j++) permutations[i].p[j] = j;
222         } else {
223             for(int j = 0; j < l;j++){
224                 int k;
225                 double rsop[3][3];
226                 mmmul(msops[i], msops[j], rsop);
227                 for(k = 0;k < l;k++){
228                     if(mequal(rsop,msops[k],t->permutation)){
229                         permutations[i].p[j] = k;
230                         break;
231                     }
232                 }
233                 if(k == l){
234                     char buf1[16];
235                     char buf2[16];
236                     symmetryOperationName(&sops[i], sizeof(buf1), buf1);
237                     symmetryOperationName(&sops[j], sizeof(buf2), buf2);
238                     msymSetErrorDetails("Unable to determine permutation for symmetry operation %s and %s",buf1, buf2);
239                     ret = MSYM_PERMUTATION_ERROR;
240                     goto err;
241                 }
242             }
243         }
244     }
245 
246     for(int i = 0; i < l;i++){
247         if(MSYM_SUCCESS != (ret = setPermutationCycles(&permutations[i]))) goto err;
248     }
249 
250     free(msops);
251     *rperm = permutations;
252 
253     return ret;
254 
255 err:
256     free(msops);
257     for(int i = 0; i < l;i++){
258         free(permutations[i].p);
259     }
260     free(permutations);
261     *rperm = NULL;
262     return ret;
263 
264 }
265 
setPermutationCycles(msym_permutation_t * perm)266 msym_error_t setPermutationCycles(msym_permutation_t *perm){
267     msym_error_t ret = MSYM_SUCCESS;
268     int l = perm->p_length;
269     int *icycle = malloc(sizeof(int[l]));
270     int *pcycle = malloc(sizeof(int[l]));
271     int *lcycle = malloc(sizeof(int[l]));
272 
273     int cl = 0;
274     memset(icycle, -1,sizeof(int[l])); //TODO: 2s complement
275     memset(lcycle,  0,sizeof(int[l]));
276 
277     perm->c = NULL;
278     perm->c_length = 0;
279 
280     for(int i = 0; i < l;i++){
281         if(icycle[i] >= 0) continue;
282         lcycle[cl] = 1;
283         pcycle[cl] = i;
284         icycle[i] = cl;
285         for(int next = perm->p[i], loop = 0; next != i;next = perm->p[next]){
286             if(loop++ > l) {
287                 msymSetErrorDetails("Encountered loop when determining permutation cycle");
288                 ret = MSYM_PERMUTATION_ERROR;
289                 goto err;
290             }
291             icycle[next] = cl;
292             lcycle[cl]++;
293         };
294         cl++;
295     }
296     perm->c_length = cl;
297     perm->c = malloc(sizeof(msym_permutation_cycle_t[cl]));
298     for(int c = 0; c < cl;c++){
299         perm->c[c].l = lcycle[c];
300         perm->c[c].s = pcycle[c];
301     }
302 
303 
304 err:
305     free(icycle);
306     free(pcycle);
307     free(lcycle);
308     return ret;
309 }
310 
311 
312 //We need these as doubles later, so might as well, even though we could represent these with ALOT less memory
permutationMatrix(msym_permutation_t * perm,double m[perm->p_length][perm->p_length])313 void permutationMatrix(msym_permutation_t *perm, double m[perm->p_length][perm->p_length]){
314     memset(m, 0, sizeof(double[perm->p_length][perm->p_length]));
315     for(int i = 0;i < perm->p_length;i++){
316         //m[i][perm->p[i]] = 1.0;
317         m[perm->p[i]][i] = 1.0;
318     }
319 }
320 
321