#!/bin/env python3
# coding: utf-8
import sys
import json

# perturb the vertices of a polytope

filename = sys.argv[1]

with open(filename,'r') as infile:
    original = json.load(infile)
    
## space for your own config options
vertexcount = 20

## pick a name
#name = "perturbed " + original['name']
#filename = "perturbed-" + filename
name = "perturbed"
filename = "perturbed.json"

## this is where you define your polytope by vertex set
vertices = []

import random
for vtex in original['vertices']:
    x = vtex[0] + random.gauss(0,0.1)
    y = vtex[1] + random.gauss(0,0.1)
    z = vtex[2] + random.gauss(0,0.1)
    vertices.append([x,y,z])

## last, choose whether to put the polytope in a nice position.
rescale = False

## that was it :)
## the rest should hopefully work out of the box

# filter extraneous points
from scipy.spatial import ConvexHull
hull = ConvexHull(vertices)
extremepoints = hull.simplices.flatten()
vertices = [vertices[index] for index in range(len(vertices)) if index in extremepoints]

polytope = {}
polytope['name'] = name
hull = ConvexHull(vertices)
unorientedfaces = hull.simplices.tolist()
polytope['edges']=[]
for i in range(vertexcount):
    for j in range(i+1,vertexcount):
        seen = False
        for f in unorientedfaces:
            if i in f and j in f:
                seen = True
        if seen:
            polytope['edges'].append([i,j])

centroid = [0,0,0]
for v in vertices:
    centroid[0] = centroid[0] + v[0]/vertexcount
    centroid[1] = centroid[1] + v[1]/vertexcount
    centroid[2] = centroid[2] + v[2]/vertexcount

mins = [1e9, 1e9, 1e9]
maxs = [-1e9, -1e9, -1e9]
for i in range(3):
    for v in vertices:
        if mins[i] > v[i]:
            mins[i] = v[i]
        if maxs[i] < v[i]:
            maxs[i] = v[i]
basicrescaling = []
for v in vertices:
    newpt = [0,0,0]
    for i in range(3):
        newpt[i] = (v[i]-centroid[i])/(maxs[i]-mins[i])
    basicrescaling.append(newpt)

if rescale:
    import numpy as np
    u,s,vh = np.linalg.svd(basicrescaling,full_matrices=False)

    polytope['vertices'] = ((u)@vh).tolist()
else:
    polytope['vertices'] = vertices

# polyhedron viewer expects all faces to have their vertices
# listed in the same order (clockwise/counterclockwise).
# making that happen right here.
import numpy
orientedfaces = []
for f in unorientedfaces:
    Vone = basicrescaling[f[0]]
    Vtwo = basicrescaling[f[1]]
    Vthree = basicrescaling[f[2]]
    
    det = numpy.linalg.det(numpy.array([Vone,Vtwo,Vthree]))
    if det == 0:
        print("this shouldn't happen :(")
    else:
        if det > 0:
            orientedfaces.append(f)
        else:
            orientedfaces.append([f[0],f[2],f[1]])

polytope['faces'] = orientedfaces

with open(filename, 'w') as outfile:
            json.dump(polytope, outfile)
