using System;
using System.Collections.Generic;
using System.Drawing;
using System.Linq;

namespace Cg1Model
{
    public abstract class Shape
    {
        public Vector[] Points { get; protected set; }
       
        public abstract bool IsVisible { get; }
       

        protected static Vector[] SelectPointsToArray(IEnumerable<Vector> points)
        {
            return points.Select(p => p.Clone()).ToArray();
        }

        protected Shape(Vector[] points)
        {
            Points = points;
        }

        public void ProcessAllPoints(Action<Vector> action)
        {
            foreach (var point in Points)
            {
                action(point);
            }
        }
       
        protected virtual void OnProcessCenter(Action<Vector> action)
        {

        }

        public void ProcessAllPointsInCenter(Action<Vector> action)
        {

            var center = Center;
            foreach (var point in Points)
            {
                point.Move(-center);
                action(point);
                point.Move(center);
            }
            OnProcessCenter(action);
        }

        public void RotateX(double angle)
        {
            ProcessAllPoints(point => point.RotateX(angle));
        }

        public void RotateY(double angle)
        {
            ProcessAllPoints(point => point.RotateY(angle));
        }

        public void RotateZ(double angle)
        {
            ProcessAllPoints(point => point.RotateZ(angle));
        }

        public void Translate(Vector vector)
        {
            ProcessAllPoints(point => point.Move(vector));
        }

        public void ReflectX()
        {
            ProcessAllPoints(point => point.ReflectX());
        }

        public void ReflectY()
        {
            ProcessAllPoints(point => point.ReflectY());
        }

        public void ReflectZ()
        {
            ProcessAllPoints(point => point.ReflectZ());
        }

        public void RotateCenterX(double angle)
        {
            ProcessAllPointsInCenter(vector => vector.RotateX(angle));
        }

        public void RotateCenterY(double angle)
        {
            ProcessAllPointsInCenter(vector => vector.RotateY(angle));
        }

        public void RotateCenterZ(double angle)
        {
            ProcessAllPointsInCenter(vector => vector.RotateZ(angle));
        }


        public void ReflectCenterX()
        {
            ProcessAllPointsInCenter(vector => vector.ReflectX());
        }

        public void ReflectCenterY()
        {
            ProcessAllPointsInCenter(vector => vector.ReflectY());
        }

        public void ReflectCenterZ()
        {
            ProcessAllPointsInCenter(vector => vector.ReflectZ());
        }
       

        public T Clone<T>() where T : Shape
        {
            var clone = (T) MemberwiseClone();
            clone.Points = clone.Points.SelectToArray(p => p.Clone());
            return clone;
        }

        public Vector Center
        {
            get
            {
                return Points.Aggregate((sum, v) => sum + v)/Points.Length;
            }
        }


        public abstract void Draw(Graphics g);

        public void Move(double dx, double dy, double dz)
        {
            var delta = new Vector(dx, dy, dz);
            ProcessAllPoints(point => point.Move(delta));
        }
    }
}

open System

type Solver(f:float[]->float,n,A,B,epsilon) = class 
    let DichotomySolve (f:float->float)=
        (* single parameter optimization for coordinate descent method *)
        if B <= A then
            failwith "Error"
        else
        let middle a b = (a+b)/2.
        let l = epsilon
        let epsilon = l / 10.;
        let rec func (a,b) =
            let middle = middle a b
            let lambda, mu = middle - epsilon, middle + epsilon
            if b - a < l then
                middle
            else
                let flambda, fmu = f lambda, f mu
                if f lambda < f mu then func(a, mu) else func(lambda, b)
        func (A,B)
    let ( .* ) (a:float) (vector:float[]) =
        vector |> Array.map (fun x->x * a)
    let ( .+ ) (a:float) (vector:float[]) =
        vector |> Array.map (fun x->x + a)
    let ( |- ) (X:float[]) (Y:float[]) =
        Array.map2(fun x y -> x - y) X Y
    let ( |+ ) (X:float[]) (Y:float[]) =
        Array.map2(fun x y -> x + y) X Y
    let ( |* ) (X:float[]) (Y:float[]) =
        Array.map2(fun x y -> x * y) X Y |> Array.sum         
    let Norm vector:float =
        vector |* vector |> sqrt
    let ( ./) (vector:float[]) (a:float) =
        (1./a) .* vector     
    member this.Solve =
        let d = Array.init n <| fun i -> Array.init n <| fun index -> if index = i then 1. else 0.
        let newLambda (y:float[]) j lambda = f <| y.[j] .+ lambda  .* d.[j]
        let newD (newLambda:float[]) =             
             let a = Array.init n <| fun j ->
                 if (newLambda.[j] = 0.) then d.[j] else
                    let summand i = newLambda.[i] .* d.[i]
                    let summands = Array.init n (fun i -> summand(i))
                    let res = Seq.skip <| j <| Seq.ofArray summands
                    Seq.fold (fun acc elm -> acc |+ elm) <| Array.zeroCreate n <| res
             let b = Array.init n <| fun j ->
                 if (j = 0) then a.[j] else
                     let add = Array.init <| j  <| fun i -> (a.[j-1] |* d.[i]) .* d.[i]
                     let summand =  Array.fold (fun acc elm -> acc |+ elm)  (Array.zeroCreate n)  add
                     a.[j] |- summand
             let result = b |> Array.map (fun b -> b ./ Norm b )
             let proof = result.[0] |* result.[1]
             d
        let rec solve y d =
            let x = y
            let newLambdas = Array.zeroCreate n (* Array.init n <| fun j -> DichotomySolve <| newLambda y j *)
            let mutable my = y
            for i in [0..n-1] do
                let newLambda = DichotomySolve <| newLambda my i
                my.[i] <- my.[i] + newLambda
            let newX = my
            let delta = newX |- x;
            let norm = delta |> Norm

            if norm < epsilon then newX else solve <| newX <| newD newLambdas
        let starty = Array.zeroCreate n
        starty.[1] <- 3.;
        solve <| starty <|  d
end

[<EntryPoint>]
let main argv =
    let f (x : float[]) = (x.[0] - 2.)**4. + (x.[0] - 2. * x.[1])**2.//(x.[0] - 1.)**2. + (x.[1] - 2.)**2.// + (x.[2] + 2.)**2.
    let a = new Solver(f, 2, -10.,10.,0.01)   
    printfn "%A" a.Solve
    Console.ReadKey() |> ignore
    0 // exit code