"""C08 构造教学模型。标准库、确定性输出；所有最优性范围见各条输入。

执行后写入脚本相邻 JSON。未调用行情、求解器或外部资料。
"""
from itertools import product
from math import sqrt
from pathlib import Path
import json


def dot(a, b):
    return sum(x*y for x, y in zip(a, b))


def quad(w, s):
    return sum(w[i]*s[i][j]*w[j] for i in range(len(w)) for j in range(len(w)))


def cov(a, b):
    if len(a) < 2 or len(a) != len(b):
        raise ValueError('协方差需要至少两个对齐观测')
    ma, mb = sum(a)/len(a), sum(b)/len(b)
    return sum((x-ma)*(y-mb) for x, y in zip(a, b))/(len(a)-1)


def simplex(n, units):
    for k in product(range(units+1), repeat=n):
        if sum(k) == units:
            yield [v/units for v in k]


def c078():
    inp = {'assets':['A','B','C','D'], 'box':[0,70], 'step':10,
           'constraints':[['总资金',[1,1,1,1],'=',100],['行业甲',[1,1,0,0],'>=',60],['行业乙',[0,0,1,1],'>=',50]],
           'relaxed_industry_b_min':40, 'unit':'净资产百分点'}
    grid = list(product(range(inp['box'][0],inp['box'][1]+1,inp['step']), repeat=4))
    def check(w, mask, bmin=None):
        vals=[]
        for index,(_,coefficients,op,bound) in enumerate(inp['constraints']):
            if index==2 and bmin is not None:
                bound=bmin
            value=dot(w,coefficients)
            vals.append(value==bound if op=='=' else value>=bound)
        return all(vals[index] for index in range(3) if mask&(1<<index))
    subsets=[]
    for mask in range(8):
        witnesses=[list(w) for w in grid if check(w,mask)]
        subsets.append({'mask':mask,'kept':[inp['constraints'][i][0] for i in range(3) if mask&(1<<i)],
                        'feasible_count':len(witnesses),'witness':witnesses[0] if witnesses else None})
    witness=next(list(w) for w in grid if check(w,7,inp['relaxed_industry_b_min']))
    return {'input':inp,'output':{'matrix':[[name,*a,op,b] for name,a,op,b in inp['constraints']],
        'bounds':[['单资产下限合计',4*inp['box'][0]],['单资产上限合计',4*inp['box'][1]],['两行业下限和',sum(r[3] for r in inp['constraints'][1:])],['可用总资金',inp['constraints'][0][3]]],
        'subsets':subsets,'minimal_conflict_relative_to_boxes':all(x['witness'] is not None for x in subsets[:-1]) and subsets[-1]['witness'] is None,
        'relaxed_witness':witness,'relaxed_totals':[sum(witness),sum(witness[:2]),sum(witness[2:])]}}


def c079():
    inp={'nav':100,'names':['A','B','C','D'],'positions':[60,40,-55,-45], 'cash':100,
         'returns':[.05,-.02,.04,-.03],'cash_rate':.001,'short_fee':.002,'trades':[-5,5,5,-5],'trade_fee':.001,
         'contrasts':{'half_short':{'positions':[60,40,-27.5,-22.5],'cash':50},'borrow_cash':{'positions':[120,80,-27.5,-22.5],'cash':-50}},
         'unit':'万元；单期收益率小数'}
    h=inp['positions']; pnl=[v*r for v,r in zip(h,inp['returns'])]; marked=[v+p for v,p in zip(h,pnl)]
    interest=inp['cash']*inp['cash_rate']; borrow=sum(max(-v,0) for v in h)*inp['short_fee']
    fee=sum(abs(v) for v in inp['trades'])*inp['trade_fee']; c1=inp['cash']+interest-borrow
    c2=c1-sum(inp['trades'])-fee; final=[x+y for x,y in zip(marked,inp['trades'])]; nav=sum(final)+c2
    def exposure(v,c):
        n=sum(v)+c
        return {'long':sum(max(x,0) for x in v),'short':sum(max(-x,0) for x in v),'gross':sum(abs(x) for x in v),'net':sum(v),
                'nav':n,'gross_weight':sum(abs(x) for x in v)/n,'net_weight':sum(v)/n,'cash_weight':c/n}
    return {'input':inp,'output':{'positions':[[name,h[i],h[i]/100,inp['returns'][i],pnl[i],marked[i],inp['trades'][i],final[i]] for i,name in enumerate(inp['names'])],
      'cash_ledger':[['期初现金',inp['cash']],['现金利息',interest],['借券持仓费',-borrow],['交易现金流',-sum(inp['trades'])],['交易费用',-fee],['期末现金',c2]],
      'stages':{'initial':exposure(h,inp['cash']),'marked_after_holding_cost':exposure(marked,c1),'after_trade':exposure(final,c2)},
      'nav_bridge':[100,sum(pnl),interest,-borrow,-fee,nav], 'self_finance_residual':sum(inp['trades'])+(c2-c1)+fee,
      'contrast_half_short':exposure(inp['contrasts']['half_short']['positions'],inp['contrasts']['half_short']['cash']),
      'contrast_borrow_cash':exposure(inp['contrasts']['borrow_cash']['positions'],inp['contrasts']['borrow_cash']['cash'])}}


def c080():
    inp={'beta':[1.6,1.2,1.6,.4],'fixed_long':[.25,.25],'beta_perturbation':[.2,0,-.1,0], 'unit':'净资产权重小数'}
    a,b=inp['fixed_long']; beta=inp['beta']; long_beta=dot([a,b],beta[:2]); total_short=-(a+b)
    c=(-long_beta-beta[3]*total_short)/(beta[2]-beta[3]); d=total_short-c
    v=-long_beta/(beta[2]+beta[3]); solutions=[('仅金额中性',[a,b,total_short/2,total_short/2]),('金额与beta同时中性',[a,b,c,d]),('仅beta中性',[a,b,v,v])]
    pert=[x+y for x,y in zip(beta,inp['beta_perturbation'])]
    rows=[[n,*w,sum(w),dot(w,beta),sum(abs(x) for x in w),dot(w,pert),dot(w,inp['beta_perturbation'])] for n,w in solutions]
    return {'input':inp,'output':{'exposures':[[chr(65+i),beta[i],pert[i]] for i in range(4)],'solutions':rows,'equations':{'long_beta':long_beta,'short_sum':total_short,'short_beta_required':-long_beta}}}


def c081():
    inp={'current':[45,15,20,20],'target':[25,15,30,30],'frozen_index':0,'step':5,'group':[0,1],'caps':[40,50], 'unit':'净资产百分点，价格及总值固定'}
    rows=[]
    candidates=[[round(v*100,10) for v in w] for w in simplex(4,100//inp['step'])]
    fixed=inp['frozen_index']; group=inp['group']
    for frozen,cap in [(True,40),(True,50),(False,40)]:
        feasible=[w for w in candidates if (not frozen or w[fixed]==inp['current'][fixed]) and sum(w[k] for k in group)<=cap]
        best=min(feasible,key=lambda w:(sum((a-b)**2 for a,b in zip(w,inp['target'])),w)) if feasible else None
        rows.append({'frozen':frozen,'cap':cap,'remaining_group_room':cap-inp['current'][fixed] if frozen and fixed in group else cap,
                     'feasible_count':len(feasible),'witness':best,'target_residual':None if best is None else [a-b for a,b in zip(best,inp['target'])]})
    return {'input':inp,'output':{'scenarios':rows,'variable_matrix':[[chr(65+k),'固定数量折合权重' if k==inp['frozen_index'] else '自由变量',inp['current'][k],inp['target'][k]] for k in range(4)]}}


def c082():
    inp={'current':[55,25,10,10],'lower':[0,0,20,15],'upper':[35,30,30,30],'step':5,'one_side_caps':[15,20], 'unit':'净资产百分点，无现金变动、无费用'}
    sell=[max(x-u,0) for x,u in zip(inp['current'],inp['upper'])]; buy=[max(l-x,0) for x,l in zip(inp['current'],inp['lower'])]
    feasible=[[round(x*100) for x in w] for w in simplex(4,100//inp['step']) if all(l-1e-9<=x*100<=u+1e-9 for x,l,u in zip(w,inp['lower'],inp['upper']))]
    def l1(w): return sum(abs(a-b) for a,b in zip(w,inp['current']))
    w=min(feasible,key=lambda w:(l1(w),sum((x-25)**2 for x in w),w)); trades=[a-b for a,b in zip(w,inp['current'])]
    minimum=l1(w)/2
    return {'input':inp,'output':{'assets':[[chr(65+i),inp['current'][i],inp['lower'][i],inp['upper'][i],sell[i],buy[i],w[i],trades[i]] for i in range(4)],
      'lower_bound':max(sum(sell),sum(buy)), 'minimum_one_side':minimum,'minimum_l1':l1(w),'witness':w,
      'cap_tests':[[c,2*c,c>=minimum] for c in inp['one_side_caps']],'feasible_count_without_turnover':len(feasible)}}


def c083():
    inp={'mu':[.01,.03,.05,.07],'current':[.2,.5,.8],'curvature':.08,'buy_cost':.01,'sell_cost':.006,'bounds':[0,1],
         'threshold_mu':[.033,.034,.04,.05,.051],'threshold_current':.5,'unit':'单期净资产比例；解析连续解'}
    def solve(mu,h):
        k,cb,cs=inp['curvature'],inp['buy_cost'],inp['sell_cost']; lo=k*h-cs; hi=k*h+cb
        x=h if lo<=mu<=hi else ((mu+cs)/k if mu<lo else (mu-cb)/k)
        x=max(inp['bounds'][0],min(inp['bounds'][1],x)); t=x-h; cost=cb*max(t,0)+cs*max(-t,0); risk=.5*k*x*x
        return [mu,h,lo,hi,x,t,cost,risk,mu*x-risk-cost]
    return {'input':inp,'output':{'grid':[solve(mu,h) for h in inp['current'] for mu in inp['mu']],
      'threshold_checks':[solve(mu,inp['threshold_current']) for mu in inp['threshold_mu']],
      'zero_cost_contrast':[[mu,max(inp['bounds'][0],min(inp['bounds'][1],mu/inp['curvature']))] for mu in inp['mu']]}}


def c084():
    inp={'linear_return':.1,'quadratic_penalty':.06,'industry_cap':.4,'gammas':[0,.02,.05,.052,.1], 'unit':'单期净资产比例；x为行业甲，1-x为其他，连续解析解'}
    a,b,cap=inp['linear_return'],inp['quadratic_penalty'],inp['industry_cap']; base=lambda x:a*x-b*x*x
    hard=min(cap,max(0,a/(2*b))); rows=[]
    for g in inp['gammas']:
        candidates=[0,cap,1,max(0,min(cap,a/(2*b))),max(cap,min(1,(a-g)/(2*b)))]; x=max(candidates,key=lambda x:base(x)-g*max(x-cap,0))
        violation=max(0,x-cap); rows.append([g,x,1-x,violation,base(x),g*violation,base(x)-g*violation,base(cap)])
    return {'input':inp,'output':{'soft_solutions':rows,'hard':[hard,1-hard,0,base(hard)],'threshold_gamma':a-2*b*cap,'objective_candidates':['区间端点','上限折点','折点两侧驻点']}}


def c085():
    inp={'names':['A','B','C'],'panel':[['p1',-1,-1,None],['p2',1,1,None],['p3',-1,None,-1],['p4',1,None,1],['p5',None,-1,1],['p6',None,1,-1],['p7',-.5,-.5,-.5],['p8',.5,.5,.5]],
         'witness':[-1,1,1], 'unit':'日收益百分点；样本协方差分母n-1'}
    m=[[0.0]*3 for _ in range(3)]; pairs=[]
    for i in range(3):
        for j in range(i,3):
            rows=[r for r in inp['panel'] if r[i+1] is not None and r[j+1] is not None]
            x=[r[i+1] for r in rows]; y=[r[j+1] for r in rows]; v=cov(x,y); m[i][j]=m[j][i]=v
            corr=v/sqrt(cov(x,x)*cov(y,y)); pairs.append([inp['names'][i]+inp['names'][j],[r[0] for r in rows],len(rows),sum(x)/len(x),sum(y)/len(y),v,corr])
    common=[r for r in inp['panel'] if all(x is not None for x in r[1:])]; unified=[[cov([r[i+1] for r in common],[r[j+1] for r in common]) for j in range(3)] for i in range(3)]
    w=inp['witness']; terms=[[inp['names'][i]+inp['names'][j],w[i]*m[i][j]*w[j]*(1 if i==j else 2)] for i in range(3) for j in range(i,3)]
    return {'input':inp,'output':{'pairs':pairs,'pairwise_matrix':m,'negative_terms':terms,'negative_quadratic':quad(w,m),'common_ids':[r[0] for r in common], 'common_matrix':unified,'common_quadratic':quad(w,unified)}}


def c086():
    inp={'sample':[[-2,-1.9],[-1,-.9],[1,1.1],[2,1.9]],'perturbed_last_b':2.3,'lambda_grid':[0,.25,.5,1], 'unit':'日收益百分点，协方差百分点平方；样本分母n-1；两资产非负全投资'}
    def matrix(rows): return [[cov([r[i] for r in rows],[r[j] for r in rows]) for j in range(2)] for i in range(2)]
    s=matrix(inp['sample']); pert=[r[:] for r in inp['sample']]; pert[-1][1]=inp['perturbed_last_b']; sp=matrix(pert)
    def stats(m,lam):
        mu=(m[0][0]+m[1][1])/2; a=(1-lam)*m[0][0]+lam*mu; b=(1-lam)*m[1][1]+lam*mu; c=(1-lam)*m[0][1]
        den=a+b-2*c; x=.5 if abs(den)<1e-14 else max(0,min(1,(b-c)/den)); eig=[(a+b-sqrt((a-b)**2+4*c*c))/2,(a+b+sqrt((a-b)**2+4*c*c))/2]
        return {'matrix':[[a,c],[c,b]],'eigenvalues':eig,'weights':[x,1-x],'variance':quad([x,1-x],[[a,c],[c,b]])}
    rows=[]; details=[]
    for lam in inp['lambda_grid']:
        a,b=stats(s,lam),stats(sp,lam); drift=sum(abs(x-y) for x,y in zip(a['weights'],b['weights'])); details.append({'lambda':lam,'original':a,'perturbed':b,'l1_drift':drift})
        rows.append([lam,*a['eigenvalues'],*a['weights'],*b['weights'],drift,a['variance']])
    mu=(s[0][0]+s[1][1])/2; mup=(sp[0][0]+sp[1][1])/2
    return {'input':inp,'output':{'sample_matrix':s,'perturbed_matrix':sp,'diagonal_target':[[mu,0],[0,mu]],'perturbed_target':[[mup,0],[0,mup]],'solutions':rows,'details':details,'perturbed_sample':pert}}


def c087():
    inp={'centers':[.08,.075,.06,.04],'errors':[.04,.01,.01,.005],'diagonal_variance':[.04,.04,.0225,.01], 'gamma':.5,'units':4,'max_weight':.75,'perturbed_center_a':.07,'unit':'单期收益及权重小数；四分之一网格'}
    weights=[w for w in simplex(4,inp['units']) if max(w)<=inp['max_weight']]
    def evaluate(w,centers,errors):
        nominal=dot(w,centers); discount=dot([abs(x) for x in w],errors); risk=inp['gamma']*dot([x*x for x in w],inp['diagonal_variance'])
        return [*w,nominal,discount,nominal-discount,risk,nominal-risk,nominal-discount-risk]
    rows=[evaluate(w,inp['centers'],inp['errors']) for w in weights]
    centers=inp['centers'][:]; centers[0]=inp['perturbed_center_a']; pr=[evaluate(w,centers,inp['errors']) for w in weights]
    doubled=[evaluate(w,inp['centers'],[2*x for x in inp['errors']]) for w in weights]
    return {'input':inp,'output':{'all_feasible_grid':rows,'count':len(rows),'nominal_best':max(rows,key=lambda r:r[8]),'robust_best':max(rows,key=lambda r:r[9]),
       'perturbed_nominal_best':max(pr,key=lambda r:r[8]),'perturbed_robust_best':max(pr,key=lambda r:r[9]),'double_error_robust_best':max(doubled,key=lambda r:r[9])}}


def c088():
    inp={'inventory':[80,60,40],'prices':[10,20,25],'cash':100,'capacities':[[10,5,4],[15,5,4],[15,10,4],[20,10,4],[20,10,4]],'target':2000,'contrast_targets':[1800,2300], 'unit':'股、元；固定价格、整股、零费、立即可用的教学假设'}
    def upper(start):
        room=[sum(day[i] for day in inp['capacities'][start:]) for i in range(3)]; sold=[min(x,c) for x,c in zip(inp['inventory'],room)]
        return [start+1,*room,*sold,inp['cash']+dot(sold,inp['prices'])]
    bounds=[upper(s) for s in range(5)]; inv=inp['inventory'][:]; cash=inp['cash']; days=[]
    for d,caps in enumerate(inp['capacities']):
        sold=[]
        for i in range(3):
            need=max(0,inp['target']-cash); q=min(inv[i],caps[i],int(need//inp['prices'][i])); sold.append(q); inv[i]-=q; cash+=q*inp['prices'][i]
        days.append([d+1,*caps,*sold,*inv,cash])
    checks=[]
    for target in [inp['target'],*inp['contrast_targets']]:
        starts=[r[0] for r in bounds if r[-1]>=target]; checks.append([target,bounds[0][-1],max(starts) if starts else None,max(0,target-bounds[0][-1])])
    return {'input':inp,'output':{'daily_path':days,'start_day_upper_bounds':bounds,'deadline_tests':checks,'final_inventory':inv,'final_cash':cash,
            'conservation':dot(inv,inp['prices'])+cash,'initial_nav':dot(inp['inventory'],inp['prices'])+inp['cash']}}


def build():
    funcs=[c078,c079,c080,c081,c082,c083,c084,c085,c086,c087,c088]
    return {f'DJ4-{78+i:03d}':f() for i,f in enumerate(funcs)}


if __name__ == '__main__':
    result=build()
    Path(__file__).with_suffix('.json').write_text(json.dumps(result,ensure_ascii=False,indent=2,allow_nan=False)+'\n',encoding='utf-8',newline='\n')