from __future__ import print_function
from __future__ import division
import numpy
from scipy import optimize


class Objective(object):
    def __init__(self, n, alpha, sigma, S, alpha_fixed, sigma_fixed):
        self.n = n
        self.alpha = alpha
        self.sigma = sigma
        self.S = S
        self.alpha_fixed = alpha_fixed
        self.sigma_fixed = sigma_fixed

    def __call__(self, phi):
        n_v = len(self.S)
        A, Sigma_e = self.make_matrix(phi)
        Sigma = self.Sigma(A, Sigma_e)
        e = (self.S - Sigma) * numpy.tri(n_v)
        e.shape = (n_v * n_v,)
        return e

    def make_matrix(self, x):
        n = self.n
        it = iter(x)
        A = numpy.zeros([n, n])
        Sigma_e = numpy.zeros([n, n])
        for i, j in self.alpha:
            A[i, j] = next(it)
        for i, j, val in self.alpha_fixed:
            A[i, j] = val
        for i, j in self.sigma:
            Sigma_e[i, j] = Sigma_e[j, i] = next(it)
        for i, j, val in self.sigma_fixed:
            Sigma_e[i, j] = Sigma_e[j, i] = val
        return A, Sigma_e

    def Sigma(self, A, Sigma_e):
        n_v = len(self.S)
        I = numpy.identity(self.n)
        U = numpy.zeros((n_v, self.n))
        U[:, :n_v] = numpy.identity(n_v)
        T = numpy.linalg.inv(I - A)
        return numpy.dot(
            numpy.dot(
                numpy.dot(
                    numpy.dot(U, T),
                    Sigma_e),
                T.T),
            U.T)


def gfi(Sigma, S):
    n = len(Sigma)
    I = numpy.identity(n)
    SigmaInv = numpy.linalg.inv(Sigma)
    SigmaS = numpy.dot(SigmaInv, S)
    denom = numpy.trace(numpy.dot(SigmaS, SigmaS.T))
    numer = numpy.trace(numpy.dot(SigmaS - I, (SigmaS - I).T))
    return 1 - numer / denom


# n:   観測変数の数
# p:   推定する母数の数
# gfi: 適合度
def agfi(n, p, gfi):
    df = 0.5 * n * (n + 1) - p # 自由度
    denom = 2 * df
    numer = n * (n + 1) * (1 - gfi)
    return 1 - numer / denom


def sem(n, alpha, sigma, S, alpha_fixed=None, sigma_fixed=None):
    if alpha_fixed is None:
        alpha_fixed = []
    if sigma_fixed is None:
        sigma_fixed = []
    x0 = numpy.ones(len(alpha) + len(sigma)) / 10
    obj = Objective(n, alpha, sigma, S, alpha_fixed, sigma_fixed)
    if len(x0) > 0:
        sol = optimize.leastsq(obj, x0)[0]
        A, Sigma_e = obj.make_matrix(sol)
    else:
        A, Sigma_e = obj.make_matrix(x0)
    Sigma = obj.Sigma(A, Sigma_e)

    _gfi = gfi(Sigma, S)
    p = len(alpha) + len(sigma) # θの次数

    return A, Sigma_e, _gfi, agfi(n, p, _gfi)

@@ -1758,17 +1758,77 @@
 denom%0A%0A%0A
-#
+def agfi(n, p, gfi):%0A    %22%22%22%0A    AGFI%E3%82%92%E8%A8%88%E7%AE%97%E3%81%99%E3%82%8B%0A%0A    Args:%0A       
  n:   %E8%A6%B3%E6%B8%AC
@@ -1832,17 +1832,23 @@
  %E8%A6%B3%E6%B8%AC%E5%A4%89%E6%95%B0%E3%81%AE%E6%95%B0%0A
-#
+       
  p:   %E6%8E%A8%E5%AE%9A
@@ -1854,17 +1854,23 @@
 %E5%AE%9A%E3%81%99%E3%82%8B%E6%AF%8D%E6%95%B0%E3%81%AE%E6%95%B0%0A
-#
+       
  gfi: %E9%81%A9%E5%90%88
@@ -1871,36 +1871,64 @@
 fi: %E9%81%A9%E5%90%88%E5%BA%A6%0A
-def agfi(n, p, gfi):
+%0A    Returns:%0A        AGFI(%E8%87%AA%E7%94%B1%E5%BA%A6%E8%AA%BF%E6%95%B4%E6%B8%88%E3%81%BF%E9%81%A9%E5%90%88%E5%BA%A6%E6%8C%87%E6%A8%99)%0A    %22%22%22
 %0A    df 

Add docstring